Index: ../trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp =================================================================== --- ../trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (revision 26100) +++ ../trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (revision 26101) @@ -43,113 +43,58 @@ void GiaDeflectionCorex( IssmDouble* pwi, IssmDouble* pdwidt, GiaDeflectionCoreArgs* arguments){ - /*intermediary: */ - int i; - - /*output: */ - IssmDouble wi=0; - IssmDouble dwidt=0; - - /*inputs: {{{*/ - /*constant: */ - int idisk=1; // disk # - IssmDouble yts; - - /*coming from the model structure, runtime configurable:*/ - /*gia solution parameters: */ - int iedge=0; - - /*gia inputs: */ - IssmDouble ri; //radial distance from center of disk to vertex i - IssmDouble re; //radius of disk - IssmDouble* hes; //loading history (in ice thickness) - IssmDouble* times; //loading history times - int numtimes; //loading history length - IssmDouble currenttime; - int Ntime; // number of times with load history - int Ntimm; // Ntime-1 : for slope/y-cept of load segments - int Ntimp; // Ntime+1 : for evaluation time - IssmDouble* blockt_time=NULL; - IssmDouble* blockt_bi=NULL; - IssmDouble* blockt_dmi=NULL; - IssmDouble* blocky_zhload=NULL; - - /*gia material parameters: */ - IssmDouble lithosphere_shear_modulus; - IssmDouble lithosphere_density; - IssmDouble mantle_shear_modulus; - IssmDouble mantle_viscosity; - IssmDouble mantle_density; - IssmDouble lithosphere_thickness; - - /*ice properties: */ - IssmDouble rho_ice; - - /*some debug info: */ - int disk_id; - -/*}}}*/ - /*Recover material parameters and loading history: see GiaDeflectionCoreArgs for more details {{{*/ - ri = arguments->ri; - re = arguments->re; - hes = arguments->hes; - times = arguments->times; - numtimes = arguments->numtimes; - currenttime = arguments->currenttime; - lithosphere_shear_modulus = arguments->lithosphere_shear_modulus; - lithosphere_density = arguments->lithosphere_density; - mantle_shear_modulus = arguments->mantle_shear_modulus; - mantle_viscosity = arguments->mantle_viscosity; - mantle_density = arguments->mantle_density; - lithosphere_thickness = arguments->lithosphere_thickness; - rho_ice = arguments->rho_ice; - disk_id = arguments->idisk; - iedge = arguments->iedge; - yts = arguments->yts; + IssmDouble ri = arguments->ri; //radial distance from center of disk to vertex i + IssmDouble re = arguments->re; //radius of disk + IssmDouble *hes = arguments->hes; //loading history (in ice thickness) + IssmDouble *times = arguments->times; //loading history times + int numtimes = arguments->numtimes; //loading history length + IssmDouble currenttime = arguments->currenttime; + IssmDouble lithosphere_shear_modulus = arguments->lithosphere_shear_modulus; + IssmDouble lithosphere_density = arguments->lithosphere_density; + IssmDouble mantle_shear_modulus = arguments->mantle_shear_modulus; + IssmDouble mantle_viscosity = arguments->mantle_viscosity; + IssmDouble mantle_density = arguments->mantle_density; + IssmDouble lithosphere_thickness = arguments->lithosphere_thickness; + IssmDouble rho_ice = arguments->rho_ice; + int disk_id = arguments->idisk; + int iedge = arguments->iedge; + IssmDouble yts = arguments->yts; /*}}}*/ /*Modify inputs to match naruse code: */ - Ntime=numtimes; - Ntimm=Ntime-1; - Ntimp=Ntime+1; + int Ntime=numtimes; // number of times with load history + int Ntimm=Ntime-1; // Ntime-1 : for slope/y-cept of load segments + int Ntimp=Ntime+1; // Ntime+1 : for evaluation time /*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: */ - blockp_.pset[0]=lithosphere_thickness; - blockp_.pset[1]=mantle_viscosity; - blockp_.pset[2]=lithosphere_shear_modulus; - blockp_.pset[3]=mantle_shear_modulus; - blockp_.pset[4]=lithosphere_density; - blockp_.pset[5]=mantle_density; - blockp_.pset[6]=re; - blocko_.rhoi=rho_ice; + blockp_.pset[0]=reCast(lithosphere_thickness); + blockp_.pset[1]=reCast(mantle_viscosity); + blockp_.pset[2]=reCast(lithosphere_shear_modulus); + blockp_.pset[3]=reCast(mantle_shear_modulus); + blockp_.pset[4]=reCast(lithosphere_density); + blockp_.pset[5]=reCast(mantle_density); + blockp_.pset[6]=reCast(re); + blocko_.rhoi=reCast(rho_ice); + blockrad_.distrad=reCast(ri)/1000.0; // in km /*loading history: */ - blocky_zhload=xNew(Ntime); - for(i=0;i(Ntime); + for(int i=0;i(hes[i]); /*times in kyr: */ - blockt_time=xNew(Ntimp); - for (i=0;i(Ntimp); + for(int i=0;i(times[numtimes-1])/1000.0/yts; // final loading time, same as evaluation time + if(i==numtimes) blockt_time[i]=reCast(times[numtimes-1])/1000.0/yts; // evaluation time } - /*bi: */ - blockt_bi=xNew(Ntimm); + IssmPDouble* blockt_bi=xNew(Ntimm); + IssmPDouble* blockt_dmi=xNew(Ntimm); - /*dmi: */ - blockt_dmi=xNew(Ntimm); - - /*radial distance of i-th element: */ - blockrad_.distrad=ri/1000.0; // in km - /*}}}*/ - /*Call distme driver: */ distme_(&Ntime,&Ntimp,&Ntimm,blockt_time,blockt_bi,blockt_dmi,blocky_zhload); @@ -157,15 +102,13 @@ what0_(&iedge,&Ntimp,&Ntimm,blockt_time,blockt_bi,blockt_dmi); /*output solution: */ - wi = blocks_.aswokm_w; - dwidt = blocks_.aswokm_dwdt; - *pwi=wi; - *pdwidt=dwidt; + *pwi=reCast(blocks_.aswokm_w); + *pdwidt=reCast(blocks_.aswokm_dwdt); /*Free ressources: */ - xDelete(blockt_time); - xDelete(blockt_bi); - xDelete(blockt_dmi); - xDelete(blocky_zhload); + xDelete(blockt_time); + xDelete(blockt_bi); + xDelete(blockt_dmi); + xDelete(blocky_zhload); }