source: issm/oecreview/Archive/14312-15392/ISSM-15015-15016.diff@ 15393

../trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp

@@ -11 +11 @@
 #include "../../toolkits/toolkits.h"
 #include "../InputUpdateFromConstantx/InputUpdateFromConstantx.h"
 
-#define Ntime 5 // number of times with load history 
-#define Ntimm 4 // Ntime-1 : for slope/y-cept of load segments 
-#define Ntimp 6 // Ntime+1 : for evaluation time  
 
 /*External blocks: {{{*/
 struct blockp{
         double pset[7];
 };
-
-struct blockt{
-        double time[Ntimp];
-        double bi[Ntimm];
-        double dmi[Ntimm];
-};
-
+        
 struct blocko{
         double rhoi;
-        double hload[Ntime];
 };
 
-struct blocky{
-        double zhload[Ntime];
-};
 
 struct blockn{
         int irate; 
@@ -50 +37 @@
 };
 
 extern "C" { 
-        int distme_( int* pidisk,int* piedge);
+        int distme_( int* pidisk,int* piedge, int* pNtime, int* pNtimp, int* pNtimm, int* pNafter,double* time,double* bi,double* dmi,double* zhload,double* hload);
+
         int what0_( int* pidisk,int* piedge);
         extern struct blockp blockp_;
-        extern struct blockt blockt_;
         extern struct blocko blocko_;
-        extern struct blocky blocky_;
         extern struct blockn blockn_;
         extern struct blockrad blockrad_;
         extern struct blocks blocks_;
@@ -65 +51 @@
 
 void GiaDeflectionCorex( IssmDouble* pwi, IssmDouble* pdwidt, GiaDeflectionCoreArgs* arguments){
 
+        /*intermediary: */
+        int i;
+
         /*output: */
         IssmDouble wi=0;
         IssmDouble dwidt=0;
@@ -86 +75 @@
         IssmDouble* times; //loading history times
         int numtimes; //loading history length
         IssmDouble currenttime;
+        int Ntime=5; // number of times with load history 
+        int Ntimm;  // Ntime-1 : for slope/y-cept of load segments 
+        int Ntimp; // Ntime+1 : for evaluation time  
+        int Nafter;
+        IssmDouble* blockt_time=NULL;
+        IssmDouble* blockt_bi=NULL;
+        IssmDouble* blockt_dmi=NULL;
+        IssmDouble* blocky_zhload=NULL;
+        IssmDouble* blocko_hload=NULL;
 
         /*gia material parameters: */
         IssmDouble lithosphere_shear_modulus;
@@ -127 +125 @@
         /*Modify inputs to match naruse code: */
         //from our model, irate comes in with values in [1,2], which maps into [0,1] in naruse: 
         irate=irate-1;
+        Nafter=1;
+        Ntime=numtimes;
+        Ntimm=Ntime-1;
+        Ntimp=Ntime+Nafter;
 
         /*Prepare block inputs for fortran distme and what0 routines of the naruse code: {{{*/
         /*Now, let's set pset from the data that we got in input to GiaDeflectionCorex: */
@@ -140 +142 @@
         blocko_.rhoi=rho_ice; 
 
         /*loading history: */
-        blocky_.zhload[0]=hes[0];
-        blocky_.zhload[1]=hes[1];
-        blocky_.zhload[2]=hes[2];
-        blocky_.zhload[3]=hes[3];
-        blocky_.zhload[4]=hes[numtimes-1];
+        blocky_zhload=xNew<IssmDouble>(numtimes);
+        for(i=0;i<numtimes;i++){
+        blocky_zhload[i]=hes[i];
+        }
 
         /*times: */
-        blockt_.time[0]=times[0]/1000.0/yts;    // in kyr
-        blockt_.time[1]=times[1]/1000.0/yts;    // in kyr
-        blockt_.time[2]=times[2]/1000.0/yts;    // in kyr
-        blockt_.time[3]=times[3]/1000.0/yts;    // in kyr
-        blockt_.time[4]=2500.0;                 // in kyr  
-        blockt_.time[5]=times[numtimes-1]/1000.0/yts;  // evaluation time, in kyr 
+        blockt_time=xNew<IssmDouble>(numtimes+1);
+        for (i=0;i<numtimes;i++){
+                blockt_time[i]=times[i]/1000.0/yts; //in kyr
+                if(i==numtimes-2)blockt_time[i]=2500.0;                 // in kyr  
+        }
+        
+        /*bi: */
+        blockt_bi=xNew<IssmDouble>(numtimes-1);
 
+        /*dmi: */
+        blockt_dmi=xNew<IssmDouble>(numtimes-1);
+        
+        /*hload:*/
+        blocko_hload=xNew<IssmDouble>(numtimes);
+
         /*irate: */
         blockn_.irate=irate; 
 
@@ -162 +171 @@
         /*}}}*/
 
         /*Call distme driver: */
-        distme_(&idisk,&iedge); 
+        distme_(&idisk,&iedge,&Ntime,&Ntimp,&Ntimm,&Nafter,blockt_time,blockt_bi,blockt_dmi,blocky_zhload,blocko_hload); 
 
+
         /*Call what0 driver: */
         what0_(&idisk,&iedge); 
 
@@ -179 +189 @@
                 *pdwidt=dwidt;
         }
 
+        /*Free ressources: */
+        xDelete<IssmDouble>(blockt_time);
+        xDelete<IssmDouble>(blockt_bi);
+        xDelete<IssmDouble>(blockt_dmi);
+        xDelete<IssmDouble>(blocky_zhload);
+        xDelete<IssmDouble>(blocko_hload);
+
 }
+      

../trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f

@@ -1 @@
-      subroutine distme(idisk,iedge)
+      subroutine distme(idisk,iedge,Ntime,Ntimp,Ntimm,Nafter,time,bi,
+     &dmi,zhload,hload)
       implicit double precision (a-h,o-y)
+      integer idisk,iedge,Ntime,Ntimp,Ntimm,Nafter
       parameter (N3G = 1)
-      parameter (Ntime = 5)
-      parameter (Ntimm = Ntime-1)
-      parameter (Nafter = 1)
-      parameter (Ntimp = Ntime + Nafter)
       double precision pset(7)
       double precision time(Ntimp),dmi(Ntimm),bi(Ntimm),dumbt(Ntimp)
       double precision hload(Ntime),qpat(Ntime),qt(Ntime)
@@ -12 +10 @@
 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
       common /blockp/ pset
       common /blockrad/ distrad 
-      common /blockt/ time,bi,dmi
-      common /blocky/ zhload
-      common /blocko/ rhoi,hload
+      common /blocko/ rhoi
       data g /9.832186d0/, yearco /3.15576d7/, eradm/6.371d6/
       data dpi /3.1415926535897932d0/, dzero/0.0d0/
 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
@@ -26 +22 @@
       dumbt(k) = time(k)
   776 continue
 c      write(6,*) time(1), time(2), time(3)
-c      write(6,*) zhload(1), zhload(2), distrad, rhoi
 c      write(6,*) pset(1), pset(2), pset(3), pset(4), pset(5), pset(6)
 c      write(6,*) pset(7)
 c      call dvecpr(time,Ntime,'::::: time @ distme.f :::::',79,0,0)