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

Index: ../trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp        (revision 15015)
+++ ../trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp        (revision 15016)
@@ -11,29 +11,16 @@
 #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,12 +37,11 @@
 };
 
 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,6 +51,9 @@
 
 void GiaDeflectionCorex( IssmDouble* pwi, IssmDouble* pdwidt, GiaDeflectionCoreArgs* arguments){
 
+       /*intermediary: */
+       int i;
+
        /*output: */
        IssmDouble wi=0;
        IssmDouble dwidt=0;
@@ -86,6 +75,15 @@
        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,6 +125,10 @@
        /*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,20 +142,27 @@
        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,8 +171,9 @@
        /*}}}*/
 
        /*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,4 +189,12 @@
                *pdwidt=dwidt;
        }
 
+       /*Free ressources: */
+       xDelete<IssmDouble>(blockt_time);
+       xDelete<IssmDouble>(blockt_bi);
+       xDelete<IssmDouble>(blockt_dmi);
+       xDelete<IssmDouble>(blocky_zhload);
+       xDelete<IssmDouble>(blocko_hload);
+
 }
+      
Index: ../trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f
===================================================================
--- ../trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f      (revision 15015)
+++ ../trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f      (revision 15016)
@@ -1,10 +1,8 @@
-      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,9 +10,7 @@
 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,7 +22,6 @@
       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)