Changeset 19172 for issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

Timestamp:

03/04/15 11:09:39 (10 years ago)

Author:

lemorzad

Message:

new method to interpolate present day temp and prec from delta O18

File:

• issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp ¶

@@ -713 +713 @@
         this->inputs->AddInput(NewPrecipitationInput);
         
+        this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
+        this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
+
+        /*clean-up*/
+        delete gauss;
+}
+/*}}}*/
+void       Penta::Delta18opdParameterization(void){/*{{{*/
+        /*Are we on the base? If not, return*/
+        if(!IsOnBase()) return;
+
+        int        i;
+        IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
+        IssmDouble TemperaturesPresentday[NUMVERTICES][12];
+        IssmDouble PrecipitationsPresentday[NUMVERTICES][12];
+        IssmDouble tmp[NUMVERTICES];
+        IssmDouble Delta18oTime;
+        IssmDouble dpermil;
+        IssmDouble time,yts;
+        this->parameters->FindParam(&time,TimeEnum);
+        this->parameters->FindParam(&yts,ConstantsYtsEnum);
+
+        /*Get some pdd parameters*/
+        dpermil=matpar-> GetMaterialParameter(SurfaceforcingsDpermilEnum);
+
+        /*Recover present day temperature and precipitation*/
+        Input*     input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input);
+        Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input2);
+        GaussPenta* gauss=new GaussPenta();
+        for(int month=0;month<12;month++) {
+                for(int iv=0;iv<NUMVERTICES;iv++) {
+                        gauss->GaussVertex(iv);
+                        input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
+                        input2->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
+                }
+        }
+
+        /*Recover interpolation parameters at time t*/
+        this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time);
+
+        /*Compute the temperature and precipitation*/
+        for(int iv=0;iv<NUMVERTICES;iv++){
+          ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,
+                                        &PrecipitationsPresentday[iv][0], &TemperaturesPresentday[iv][0], 
+                                        &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
+        }
+
+        /*Update inputs*/ 
+        TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
+        TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
+        for (int imonth=0;imonth<12;imonth++) {
+                for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
+                PentaInput* newmonthinput1 = new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
+                NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
+
+                for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
+                PentaInput* newmonthinput2 = new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
+                NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
+        }
+        NewTemperatureInput->Configure(this->parameters);
+        NewPrecipitationInput->Configure(this->parameters);
+
+        this->inputs->AddInput(NewTemperatureInput);
+        this->inputs->AddInput(NewPrecipitationInput);
+
         this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
         this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);