source: issm/oecreview/Archive/25834-26739/ISSM-26100-26101.diff@ 26740

Last change on this file since 26740 was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739
File size: 6.9 KB
RevLine 
[26740]1Index: ../trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (revision 26100)
4+++ ../trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (revision 26101)
5@@ -43,113 +43,58 @@
6
7 void GiaDeflectionCorex( IssmDouble* pwi, IssmDouble* pdwidt, GiaDeflectionCoreArgs* arguments){
8
9- /*intermediary: */
10- int i;
11-
12- /*output: */
13- IssmDouble wi=0;
14- IssmDouble dwidt=0;
15-
16- /*inputs: {{{*/
17- /*constant: */
18- int idisk=1; // disk #
19- IssmDouble yts;
20-
21- /*coming from the model structure, runtime configurable:*/
22- /*gia solution parameters: */
23- int iedge=0;
24-
25- /*gia inputs: */
26- IssmDouble ri; //radial distance from center of disk to vertex i
27- IssmDouble re; //radius of disk
28- IssmDouble* hes; //loading history (in ice thickness)
29- IssmDouble* times; //loading history times
30- int numtimes; //loading history length
31- IssmDouble currenttime;
32- int Ntime; // number of times with load history
33- int Ntimm; // Ntime-1 : for slope/y-cept of load segments
34- int Ntimp; // Ntime+1 : for evaluation time
35- IssmDouble* blockt_time=NULL;
36- IssmDouble* blockt_bi=NULL;
37- IssmDouble* blockt_dmi=NULL;
38- IssmDouble* blocky_zhload=NULL;
39-
40- /*gia material parameters: */
41- IssmDouble lithosphere_shear_modulus;
42- IssmDouble lithosphere_density;
43- IssmDouble mantle_shear_modulus;
44- IssmDouble mantle_viscosity;
45- IssmDouble mantle_density;
46- IssmDouble lithosphere_thickness;
47-
48- /*ice properties: */
49- IssmDouble rho_ice;
50-
51- /*some debug info: */
52- int disk_id;
53-
54-/*}}}*/
55-
56 /*Recover material parameters and loading history: see GiaDeflectionCoreArgs for more details {{{*/
57- ri = arguments->ri;
58- re = arguments->re;
59- hes = arguments->hes;
60- times = arguments->times;
61- numtimes = arguments->numtimes;
62- currenttime = arguments->currenttime;
63- lithosphere_shear_modulus = arguments->lithosphere_shear_modulus;
64- lithosphere_density = arguments->lithosphere_density;
65- mantle_shear_modulus = arguments->mantle_shear_modulus;
66- mantle_viscosity = arguments->mantle_viscosity;
67- mantle_density = arguments->mantle_density;
68- lithosphere_thickness = arguments->lithosphere_thickness;
69- rho_ice = arguments->rho_ice;
70- disk_id = arguments->idisk;
71- iedge = arguments->iedge;
72- yts = arguments->yts;
73+ IssmDouble ri = arguments->ri; //radial distance from center of disk to vertex i
74+ IssmDouble re = arguments->re; //radius of disk
75+ IssmDouble *hes = arguments->hes; //loading history (in ice thickness)
76+ IssmDouble *times = arguments->times; //loading history times
77+ int numtimes = arguments->numtimes; //loading history length
78+ IssmDouble currenttime = arguments->currenttime;
79+ IssmDouble lithosphere_shear_modulus = arguments->lithosphere_shear_modulus;
80+ IssmDouble lithosphere_density = arguments->lithosphere_density;
81+ IssmDouble mantle_shear_modulus = arguments->mantle_shear_modulus;
82+ IssmDouble mantle_viscosity = arguments->mantle_viscosity;
83+ IssmDouble mantle_density = arguments->mantle_density;
84+ IssmDouble lithosphere_thickness = arguments->lithosphere_thickness;
85+ IssmDouble rho_ice = arguments->rho_ice;
86+ int disk_id = arguments->idisk;
87+ int iedge = arguments->iedge;
88+ IssmDouble yts = arguments->yts;
89
90 /*}}}*/
91
92 /*Modify inputs to match naruse code: */
93- Ntime=numtimes;
94- Ntimm=Ntime-1;
95- Ntimp=Ntime+1;
96+ int Ntime=numtimes; // number of times with load history
97+ int Ntimm=Ntime-1; // Ntime-1 : for slope/y-cept of load segments
98+ int Ntimp=Ntime+1; // Ntime+1 : for evaluation time
99
100 /*Prepare block inputs for fortran distme and what0 routines of the naruse code: {{{*/
101 /*Now, let's set pset from the data that we got in input to GiaDeflectionCorex: */
102- blockp_.pset[0]=lithosphere_thickness;
103- blockp_.pset[1]=mantle_viscosity;
104- blockp_.pset[2]=lithosphere_shear_modulus;
105- blockp_.pset[3]=mantle_shear_modulus;
106- blockp_.pset[4]=lithosphere_density;
107- blockp_.pset[5]=mantle_density;
108- blockp_.pset[6]=re;
109- blocko_.rhoi=rho_ice;
110+ blockp_.pset[0]=reCast<IssmPDouble>(lithosphere_thickness);
111+ blockp_.pset[1]=reCast<IssmPDouble>(mantle_viscosity);
112+ blockp_.pset[2]=reCast<IssmPDouble>(lithosphere_shear_modulus);
113+ blockp_.pset[3]=reCast<IssmPDouble>(mantle_shear_modulus);
114+ blockp_.pset[4]=reCast<IssmPDouble>(lithosphere_density);
115+ blockp_.pset[5]=reCast<IssmPDouble>(mantle_density);
116+ blockp_.pset[6]=reCast<IssmPDouble>(re);
117+ blocko_.rhoi=reCast<IssmPDouble>(rho_ice);
118+ blockrad_.distrad=reCast<IssmPDouble>(ri)/1000.0; // in km
119
120 /*loading history: */
121- blocky_zhload=xNew<IssmDouble>(Ntime);
122- for(i=0;i<Ntime;i++){
123- blocky_zhload[i]=hes[i];
124- }
125+ IssmPDouble* blocky_zhload=xNew<IssmPDouble>(Ntime);
126+ for(int i=0;i<Ntime;i++) blocky_zhload[i]=reCast<IssmPDouble>(hes[i]);
127
128 /*times in kyr: */
129- blockt_time=xNew<IssmDouble>(Ntimp);
130- for (i=0;i<Ntimp;i++){
131+ IssmPDouble* blockt_time=xNew<IssmPDouble>(Ntimp);
132+ for(int i=0;i<Ntimp;i++){
133 blockt_time[i]=times[i]/1000.0/yts;
134- if(i==numtimes-1)blockt_time[i]=times[numtimes-1]/1000.0/yts; // final loading time, same as evaluation time
135- if(i==numtimes)blockt_time[i]=times[numtimes-1]/1000.0/yts; // evaluation time
136+ if(i==numtimes-1) blockt_time[i]=reCast<IssmPDouble>(times[numtimes-1])/1000.0/yts; // final loading time, same as evaluation time
137+ if(i==numtimes) blockt_time[i]=reCast<IssmPDouble>(times[numtimes-1])/1000.0/yts; // evaluation time
138 }
139
140- /*bi: */
141- blockt_bi=xNew<IssmDouble>(Ntimm);
142+ IssmPDouble* blockt_bi=xNew<IssmPDouble>(Ntimm);
143+ IssmPDouble* blockt_dmi=xNew<IssmPDouble>(Ntimm);
144
145- /*dmi: */
146- blockt_dmi=xNew<IssmDouble>(Ntimm);
147-
148- /*radial distance of i-th element: */
149- blockrad_.distrad=ri/1000.0; // in km
150- /*}}}*/
151-
152 /*Call distme driver: */
153 distme_(&Ntime,&Ntimp,&Ntimm,blockt_time,blockt_bi,blockt_dmi,blocky_zhload);
154
155@@ -157,15 +102,13 @@
156 what0_(&iedge,&Ntimp,&Ntimm,blockt_time,blockt_bi,blockt_dmi);
157
158 /*output solution: */
159- wi = blocks_.aswokm_w;
160- dwidt = blocks_.aswokm_dwdt;
161- *pwi=wi;
162- *pdwidt=dwidt;
163+ *pwi=reCast<IssmDouble>(blocks_.aswokm_w);
164+ *pdwidt=reCast<IssmDouble>(blocks_.aswokm_dwdt);
165
166 /*Free ressources: */
167- xDelete<IssmDouble>(blockt_time);
168- xDelete<IssmDouble>(blockt_bi);
169- xDelete<IssmDouble>(blockt_dmi);
170- xDelete<IssmDouble>(blocky_zhload);
171+ xDelete<IssmPDouble>(blockt_time);
172+ xDelete<IssmPDouble>(blockt_bi);
173+ xDelete<IssmPDouble>(blockt_dmi);
174+ xDelete<IssmPDouble>(blocky_zhload);
175
176 }
Note: See TracBrowser for help on using the repository browser.