Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sisvat_gsn.f90 @ 5452

Last change on this file since 5452 was 5246, checked in by abarral, 2 months ago
Convert fixed-form to free-form sources .F -> .{f,F}90 (WIP: some .F remain, will be handled in subsequent commits)
File size: 42.0 KB

Line
1
2	subroutine SISVAT_GSn
3
4	! +------------------------------------------------------------------------+
5	! \| MAR SISVAT_GSn 20-09-2003 MAR \|
6	! \| SubRoutine SISVAT_GSn simulates SNOW Metamorphism \|
7	! +------------------------------------------------------------------------+
8	! \| \|
9	! \| PARAMETERS: knonv: Total Number of columns = \|
10	! \| ^^^^^^^^^^ = Total Number of continental grid boxes \|
11	! \| X Number of Mosaic Cell per grid box \|
12	! \| \|
13	! \| INPUT / isnoSV = total Nb of Ice/Snow Layers \|
14	! \| OUTPUT: iiceSV = total Nb of Ice Layers \|
15	! \| ^^^^^^ istoSV = 0,...,5 : Snow History (see istdSV data) \|
16	! \| \|
17	! \| INPUT: TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)\|
18	! \| ^^^^^ & Snow Temperatures (layers 1,2,...,nsno) [K] \|
19	! \| ro__SV : Soil/Snow Volumic Mass [kg/m3] \|
20	! \| eta_SV : Soil/Snow Water Content [m3/m3] \|
21	! \| slopSV : Surface Slope [-] \|
22	! \| dzsnSV : Snow Layer Thickness [m] \|
23	! \| dt__SV2 : Time Step [s] \|
24	! \| \|
25	! \| INPUT / G1snSV : Dendricity (<0) or Sphericity (>0) of Snow Layer \|
26	! \| OUTPUT: G2snSV : Sphericity (>0) or Size of Snow Layer \|
27	! \| ^^^^^^ \|
28	! \| \|
29	! \| Formalisme adopte pour la Representation des Grains: \|
30	! \| Formalism for the Representation of Grains: \|
31	! \| ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ \|
32	! \| \|
33	! \| 1 - -1 Neige Fraiche \|
34	! \| / \ \| ------------- \|
35	! \| / \ \| Dendricite decrite par Dendricite \|
36	! \| / \ \| Dendricity et Sphericite \|
37	! \| / \ \| \|
38	! \| 2---------3 - 0 described by Dendricity \|
39	! \| and Sphericity \|
40	! \| \|---------\| \|
41	! \| 0 1 \|
42	! \| Sphericite \|
43	! \| Sphericity \|
44	! \| \|
45	! \| 4---------5 - \|
46	! \| \| \| \| \|
47	! \| \| \| \| Diametre (1/10eme de mm) (ou Taille) \|
48	! \| \| \| \| Diameter (1/10th of mm) (or Size ) \|
49	! \| \| \| \| \|
50	! \| \| \| \| Neige non dendritique \|
51	! \| 6---------7 - --------------------- \|
52	! \| decrite par Sphericite \|
53	! \| et Taille \|
54	! \| described by Sphericity \|
55	! \| and Size \|
56	! \| \|
57	! \| Les Variables du Modele: \|
58	! \| Model Variables: \|
59	! \| ^^^^^^^^^^^^^^^^^^^^^^^^ \|
60	! \| Cas Dendritique Cas non Dendritique \|
61	! \| \|
62	! \| G1snSV : Dendricite G1snSV : Sphericite \|
63	! \| G2snSV : Sphericite G2snSV : Taille (1/10e mm) \|
64	! \| Size \|
65	! \| \|
66	! \| Cas Dendritique/ Dendritic Case \|
67	! \| ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ \|
68	! \| Dendricite(Dendricity) G1snSV \|
69	! \| varie de -G1_dSV (-99 par defaut / etoile) a 0 \|
70	! \| division par -G1_dSV pour obtenir des valeurs entre 1 et 0 \|
71	! \| varies from -G1_dSV (default -99 / fresh snow) to 0 \|
72	! \| division by -G1_dSV to obtain values between 1 and 0 \|
73	! \| \|
74	! \| Sphericite(Sphericity) G2snSV \|
75	! \| varie de 0 (cas completement anguleux) \|
76	! \| a G1_dSV (99 par defaut, cas spherique) \|
77	! \| division par G1_dSV pour obtenir des valeurs entre 0 et 1 \|
78	! \| varies from 0 (full faceted) to G1_dSV \|
79	! \| \|
80	! \| Cas non Dendritique / non Dendritic Case \|
81	! \| ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ \|
82	! \| Sphericite(Sphericity) G1snSV \|
83	! \| varie de 0 (cas completement anguleux) \|
84	! \| a G1_dSV (99 par defaut, cas spherique) \|
85	! \| division par G1_dSV pour obtenir des valeurs entre 0 et 1 \|
86	! \| varies from 0 (full faceted) to G1_dSV \|
87	! \| \|
88	! \| Taille (Size) G2snSV \|
89	! \| superieure a ADSdSV (.4 mm) et ne fait que croitre \|
90	! \| greater than ADSdSV (.4 mm) always increases \|
91	! \| \|
92	! \| Exemples: Points caracteristiques des Figures ci-dessus \|
93	! \| ^^^^^^^^^ \|
94	! \| \|
95	! \| G1snSV G2snSV dendricite sphericite taille \|
96	! \| dendricity sphericity size \|
97	! \| ------------------------------------------------------------------ \|
98	! \| [1/10 mm] \|
99	! \| 1 -G1_dSV sph3SN 1 0.5 \|
100	! \| 2 0 0 0 0 \|
101	! \| 3 0 G1_dSV 0 1 \|
102	! \| 4 0 ADSdSV 0 4. \|
103	! \| 5 G1_dSV ADSdSV-vsphe1 1 3. \|
104	! \| 6 0 -- 0 -- \|
105	! \| 7 G1_dSV -- 1 -- \|
106	! \| \|
107	! \| par defaut: G1_dSV=99. \|
108	! \| sph3SN=50. \|
109	! \| ADSdSV= 4. \|
110	! \| vsphe1=1. \|
111	! \| \|
112	! \| Methode: \|
113	! \| ^^^^^^^^ \|
114	! \| 1. Evolution Types de Grains selon Lois de Brun et al. (1992): \|
115	! \| Grain metamorphism according to Brun et al. (1992): \|
116	! \| Plusieurs Cas sont a distiguer / the different Cases are: \|
117	! \| 1.1 Metamorphose Neige humide / wet Snow \|
118	! \| 1.2 Metamorphose Neige seche / dry Snow \|
119	! \| 1.2.1 Gradient faible / low Temperature Gradient \|
120	! \| 1.2.2 Gradient moyen / moderate Temperature Gradient \|
121	! \| 1.2.3 Gradient fort / high Temperature Gradient \|
122	! \| Dans chaque Cas on separe Neige Dendritique et non Dendritique \|
123	! \| le Passage Dendritique -> non Dendritique \|
124	! \| se fait lorsque G1snSV devient > 0 \|
125	! \| the Case of Dentritic or non Dendritic Snow is treated separately \|
126	! \| the Limit Dentritic -> non Dendritic is reached when G1snSV > 0 \|
127	! \| \|
128	! \| 2. Tassement: Loi de Viscosite adaptee selon le Type de Grains \|
129	! \| Snow Settling: Viscosity depends on the Grain Type \|
130	! \| \|
131	! \| 3. Update Variables historiques (cas non dendritique seulement) \|
132	! \| nhSNow defaut \|
133	! \| 0 Cas normal \|
134	! \| istdSV(1) 1 Grains anguleux / faceted cristal \|
135	! \| istdSV(2) 2 Grains ayant ete en presence d eau liquide \|
136	! \| mais n'ayant pas eu de caractere anguleux / \|
137	! \| liquid water and no faceted cristals before \|
138	! \| istdSV(3) 3 Grains ayant ete en presence d eau liquide \|
139	! \| ayant eu auparavant un caractere anguleux / \|
140	! \| liquid water and faceted cristals before \|
141	! \| \|
142	! \| REFER. : Brun et al. 1989, J. Glaciol 35 pp. 333--342 \|
143	! \| ^^^^^^^^ Brun et al. 1992, J. Glaciol 38 pp. 13-- 22 \|
144	! \| (CROCUS Model, adapted to MAR at CEN by H.Gallee) \|
145	! \| \|
146	! \| REFER. : Marbouty, D. 1980, J. Glaciol 26 pp. xxx--xxx \|
147	! \| ^^^^^^^^ (CROCUS Model, adapted to MAR at CEN by H.Gallee) \|
148	! \| (for angular shapes) \|
149	! \| \|
150	! \| Preprocessing Option: SISVAT IO (not always a standard preprocess.) \|
151	! \| ^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^ \|
152	! \| FILE \| CONTENT \|
153	! \| ~~~~~~~~~~~~~~~~~~~~~+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \|
154	! \| # SISVAT_GSn.vp \| #vp: OUTPUT/Verification: Snow Properties \|
155	! \| \| unit 47, SubRoutines SISVAT_zSn, _GSn \|
156	! \| # stdout \| #wp: OUTPUT/Verification: Snow Properties \|
157	! \| \| unit 6, SubRoutine SISVAT_GSn \|
158	! \| \|
159	! +------------------------------------------------------------------------+
160
161
162
163
164	! +--Global Variables
165	! + ================
166
167	use VARphy
168	use VAR_SV
169	use VARdSV
170	use VAR0SV
171	use VARxSV
172	use VARtSV
173
174
175	IMPLICIT NONE
176
177
178
179	! +--INPUT/OUTPUT
180	! + ------------
181
182
183	! +--OUTPUT
184	! + ------
185
186	integer :: dt__SV2
187
188
189	! +--Local Variables
190	! + ================
191
192	logical :: vector !
193	integer :: ikl !
194	integer :: isn ,isnp !
195	integer :: istoOK !
196	real :: G1_bak,G2_bak ! Old Values of G1, G2
197	real :: ro_dry(knonv, nsno) ! Dry Density [g/cm3]
198	real :: etaSno(knonv, nsno) ! Liquid Water Content [g/cm2]
199	real :: SnMass(knonv) ! Snow Mass [kg/m2]
200	real :: dTsndz ! Temperature Gradient
201	real :: sWater ! Water Content [%]
202	real :: exp1Wa !
203	real :: dDENDR ! Dendricity Increment
204	real :: DENDRn ! Normalized Dendricity
205	real :: SPHERn ! Normalized Sphericity
206	real :: Wet_OK ! Wet Metamorphism Switch
207	real :: OK__DE !
208	real :: OK__wd ! New G*, from wet Dendritic
209	real :: G1__wd ! New G1, from wet Dendritic
210	real :: G2__wd ! New G2, from wet Dendritic
211	real :: OKlowT !
212	real :: facVap !
213	real :: OK_ldd !
214	real :: G1_ldd !
215	real :: G2_ldd !
216	real :: DiamGx !
217	real :: DiamOK !
218	real :: No_Big !
219	real :: dSPHER !
220	real :: SPHER0 !
221	real :: SPHbig !
222	real :: G1_lds !
223	real :: OK_mdT !
224	real :: OKmidT !
225	real :: OKhigT !
226	real :: OK_mdd !
227	real :: G1_mdd !
228	real :: G2_mdd !
229	real :: G1_mds !
230	real :: OK_hdd !
231	real :: G1_hdd !
232	real :: G2_hdd !
233	real :: OK_hds !
234	real :: G1_hds !
235	real :: T1__OK,T2__OK !
236	real :: T3_xOK,T3__OK,T3_nOK !
237	real :: ro1_OK,ro2_OK !
238	real :: dT1_OK,dT2_OK,dT3xOK,dT3_OK !
239	real :: dT4xOK,dT4_OK,dT4nOK,AngSno !
240	real :: G2_hds,SphrOK,HISupd !
241	real :: H1a_OK,H1b_OK,H1__OK !
242	real :: H23aOK,H23bOK,H23_OK !
243	real :: H2__OK,H3__OK !
244	real :: H45_OK,H4__OK,H5__OK !
245	real :: ViscSn,OK_Liq,OK_Ang,OKxLiq !
246	real :: dSnMas,dzsnew,rosnew,rosmax,smb_old,smb_new
247	real :: zn_old,zn_new
248
249	real :: epsi5 ! Alpha ev67 single precision
250	real :: vdiam1 ! Small Grains Min.Diam.[.0001m]
251	real :: vdiam2 ! Spher.Variat.Max Diam. [mm]
252	real :: vdiam3 ! Min.Diam.\|Limit Spher. [mm]
253	real :: vdiam4 ! Min.Diam.\|Viscosity Change
254	real :: vsphe1 ! Max Sphericity
255	real :: vsphe2 ! Low T Metamorphism Coeff.
256	real :: vsphe3 ! Max.Sphericity (history=1)
257	real :: vsphe4 ! Min.Sphericity=>history=1
258	real :: vtang1,vtang2,vtang3,vtang4 ! Temperature Contribution
259	real :: vtang5,vtang6,vtang7,vtang8 !
260	real :: vtang9,vtanga,vtangb,vtangc !
261	real :: vrang1,vrang2 ! Density Contribution
262	real :: vgang1,vgang2,vgang3,vgang4 ! Grad(T) Contribution
263	real :: vgang5,vgang6,vgang7,vgang8 !
264	real :: vgang9,vganga,vgangb,vgangc !
265	real :: vgran6 ! Max.Sphericity for Settling
266	real :: vtelv1 ! Threshold \| history = 2, 3
267	real :: vvap1 ! Vapor Pressure Coefficient
268	real :: vvap2 ! Vapor Pressure Exponent
269	real :: vgrat1 ! Boundary weak/mid grad(T)
270	real :: vgrat2 ! Boundary mid/strong grad(T)
271	real :: vfi ! PHI, strong grad(T)
272	real :: vvisc1,vvisc2,vvisc3,vvisc4 ! Viscosity Coefficients
273	real :: vvisc5,vvisc6,vvisc7 ! id., wet Snow
274	real :: rovisc ! Wet Snow Density Influence
275	real :: vdz3 ! Maximum Layer Densification
276	real :: OK__ws ! New G2
277	real :: G1__ws ! New G1, from wet Spheric
278	real :: G2__ws ! New G2, from wet Spheric
279	real :: husi_0,husi_1,husi_2,husi_3 ! Constants for New G2
280	real :: vtail1,vtail2 ! Constants for New G2
281	real :: frac_j ! Time Step [Day]
282
283	real :: vdent1 ! Wet Snow Metamorphism
284	integer :: nvdent1 ! (Coefficients for
285	integer :: nvdent2 ! Dendricity)
286
287	! +--Snow Properties: IO
288	! + ~~~~~~~~~~~~~~~~~~~
289	! #vp real G_curr(18),Gcases(18)
290	! #vp common /GSnLOC/ Gcases
291	! #wp real D__MAX
292	! #wp common /GSnMAX/ D__MAX
293
294
295	! +--DATA
296	! + ====
297
298	data vector/.true./ ! Vectorization Switch
299	data vdent1/ 0.5e8/ ! Wet Snow Metamorphism
300	!XF tuned for Greenland (2.e8=old value)
301	data nvdent1/ 3 / ! (Coefficients for
302	data nvdent2/16 / ! Dendricity)
303
304	data husi_0 /20. / ! 10 * 2
305	data husi_1 / 0.23873 / ! (3/4) /pi
306	data husi_2 / 4.18880 / ! (4/3) *pi
307	data husi_3 / 0.33333 / ! 1/3
308	data vtail1 / 1.28e-08/ ! Wet Metamorphism
309	data vtail2 / 4.22e-10/ ! (NON Dendritic / Spheric)
310
311	data epsi5 / 1.0e-5 / !
312
313	data vdiam1 / 4.0 / ! Small Grains Min.Diameter
314
315	data vdiam2 / 0.5 / ! Spher.Variat.Max Diam.[mm]
316	data vdiam3 / 3.0 / ! Min.Diam.\|Limit Spher.[mm]
317	data vdiam4 / 2.0 / ! Min.Diam.\|Viscosity Change
318
319	data vsphe1 / 1.0 / ! Max Sphericity
320	data vsphe2 / 1.0e9 / ! Low T Metamorphism Coeff.
321	data vsphe3 / 0.5 / ! Max.Sphericity (history=1)
322	data vsphe4 / 0.1 / ! Min.Sphericity=>history=1
323
324	data vgran6 / 51. / ! Max.Sphericity for Settling
325	data vtelv1 / 5.e-1 / ! Threshold \| history = 2, 3
326
327	data vvap1 /-6.e3 / ! Vapor Pressure Coefficient
328	data vvap2 / 0.4 / ! Vapor Pressure Exponent
329
330	data vgrat1 /0.05 / ! Boundary weak/mid grad(T)
331	data vgrat2 /0.15 / ! Boundary mid/strong grad(T)
332	data vfi /0.09 / ! PHI, strong grad(T)
333
334	data vvisc1 / 0.70 / ! Viscosity Coefficients
335	data vvisc2 / 1.11e5 / !
336	data vvisc3 /23.00 / !
337	data vvisc4 / 0.10 / !
338	data vvisc5 / 1.00 / ! id., wet Snow
339	data vvisc6 / 2.00 / !
340	data vvisc7 /10.00 / !
341	data rovisc / 0.25 / ! Wet Snow Density Influence
342	data vdz3 / 0.30 / ! Maximum Layer Densification
343
344
345	! +--DATA (Coefficient Fonction fort Gradient Marbouty)
346	! + --------------------------------------------------
347
348	data vtang1 /40.0/ ! Temperature Contribution
349	data vtang2 / 6.0/ !
350	data vtang3 /22.0/ !
351	data vtang4 / 0.7/ !
352	data vtang5 / 0.3/ !
353	data vtang6 / 6.0/ !
354	data vtang7 / 1.0/ !
355	data vtang8 / 0.8/ !
356	data vtang9 /16.0/ !
357	data vtanga / 0.2/ !
358	data vtangb / 0.2/ !
359	data vtangc /18.0/ !
360
361	data vrang1 / 0.40/ ! Density Contribution
362	data vrang2 / 0.15/ !
363
364	data vgang1 / 0.70/ ! Grad(T) Contribution
365	data vgang2 / 0.25/ !
366	data vgang3 / 0.40/ !
367	data vgang4 / 0.50/ !
368	data vgang5 / 0.10/ !
369	data vgang6 / 0.15/ !
370	data vgang7 / 0.10/ !
371	data vgang8 / 0.55/ !
372	data vgang9 / 0.65/ !
373	data vganga / 0.20/ !
374	data vgangb / 0.85/ !
375	data vgangc / 0.15/ !
376
377	! #wp data D__MAX / 4.00/ !
378
379
380	! +-- 1. Metamorphoses dans les Strates
381	! + Metamorphism
382	! + ==============================
383
384	dt__SV2= dt__SV
385	frac_j = dt__SV2 / 86400. ! Time Step [Day]
386
387	zn4_SV = 0
388
389
390	! +-- 1.1 Initialisation: teneur en eau liquide et gradient de temperature
391	! + ------------------ liquid water content and temperature gradient
392
393	DO ikl=1,knonv
394	DO isn=1,isnoSV(ikl)
395
396	ro_dry(ikl,isn) = 1.e-3 *ro__SV(ikl,isn) & ! Dry Density
397	*(1. -eta_SV(ikl,isn)) ! [g/cm3]
398	etaSno(ikl,isn) = 1.e-1 *dzsnSV(ikl,isn) & ! Liquid Water
399	* ro__SV(ikl,isn) & ! Content [g/cm2]
400	* eta_SV(ikl,isn) !
401	END DO
402	END DO
403
404	!!$OMP PARALLEL DO default(firstprivate)
405	!!$OMP.shared (/xSISVAT_I/,/xSISVAT_R/,/SoR0SV/,/SoI0SV/,/Sn_dSV/)
406	DO ikl=1,knonv
407	DO isn=1,isnoSV(ikl)
408	isnp = min(isn+1,isnoSV(ikl))
409
410	dTsndz = abs( (TsisSV(ikl,isnp)-TsisSV(ikl,isn-1)) *2.e-2 &
411	/max(((dzsnSV(ikl,isnp)+dzsnSV(ikl,isn) ) &
412	*( isnp - isn) &
413	+(dzsnSV(ikl,isn )+dzsnSV(ikl,isn-1))),epsi))
414	! +... Factor 1.d-2 for Conversion K/m --> K/cm
415
416
417	! +-- 1.2 Metamorphose humide
418	! + Wet Snow Metamorphism
419	! + ---------------------
420
421	Wet_OK = max(zero,sign(unun,eta_SV(ikl,isn)-epsi))
422
423
424	! +-- Vitesse de diminution de la dendricite
425	! + Rate of the dendricity decrease
426	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
427	sWater=1.d-1ro__SV(ikl,isn)eta_SV(ikl,isn) &
428	/max(epsi,ro_dry(ikl,isn))
429	! +... sWater:Water Content [%]
430	! + 1.d-1= 1.d2(1->%) * 1.d-3(ro__SV*eta_SV:kg/m3->g/cm3)
431
432	exp1Wa= sWater**nvdent1
433	dDENDR=max(exp1Wa/nvdent2,vdent1*exp(vvap1/TfSnow))
434
435	! +-- 1.2.1 Cas dendritique/dendritic Case
436	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
437	OK__wd=max(zero, & !
438	sign(unun,-G1snSV(ikl,isn) & !
439	-epsi )) !
440
441	DENDRn=-G1snSV(ikl,isn)/G1_dSV ! Normalized Dendricity (+)
442	SPHERn= G2snSV(ikl,isn)/G1_dSV ! Normalized Sphericity
443	DENDRn= DENDRn -dDENDR *frac_j ! New Dendricity (+)
444	SPHERn= SPHERn +dDENDR *frac_j ! New Sphericity
445
446	OK__DE=max(zero, & ! IF 1.,
447	sign(unun, DENDRn & ! NO change
448	-epsi )) ! Dendr. -> Spheric
449
450	G1__wd=OK__DE * ( -DENDRn*G1_dSV) & ! Dendritic
451	+(1.-OK__DE)* min(G1_dSV,SPHERn*G1_dSV) ! Dendr. -> Spheric
452	G2__wd=OK__DE * min(G1_dSV,SPHERn*G1_dSV) & ! Spheric
453	+(1.-OK__DE)*(ADSdSV-min(SPHERn,vsphe1)) ! Spher. -> Size
454
455	! +-- 1.2.2 Cas non dendritique non completement spherique
456	! + Evolution de la Sphericite seulement.
457	! + Non dendritic and not completely spheric Case
458	! + Evolution of Sphericity only (not size)
459	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
460	OK__ws=max(zero, & !
461	sign(unun, G1_dSV & !
462	-epsi5 & !
463	-G1snSV(ikl,isn))) !
464
465	SPHERn= G1snSV(ikl,isn)/G1_dSV
466	SPHERn= SPHERn +dDENDR *frac_j
467	G1__ws= min(G1_dSV,SPHERn*G1_dSV)
468
469	! +-- 1.2.3 Cas non dendritique et spherique / non dendritic and spheric
470	! + Evolution de la Taille seulement / Evolution of Size only
471	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
472	G2__ws = husi_0 &
473	*( husi_1 &
474	(husi_2 ( G2snSV(ikl,isn)/husi_0)**3 &
475	+(vtail1 +vtail2 exp1Wa )dt__SV2)) &
476	** husi_3
477
478
479	! +-- 1.3 Metamorposes seches / Dry Metamorphism
480	! + --------------------------------------
481
482
483	! +-- 1.3.1 Calcul Metamorphose faible/low Gradient (0.00-0.05 deg/cm)
484	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
485	OKlowT=max(zero, & !
486	sign(unun, vgrat1 & !
487	-dTsndz )) !
488
489	facVap=exp(vvap1/TsisSV(ikl,isn))
490
491	! +-- 1.3.1.1 Cas dendritique / dendritic Case
492
493	OK_ldd=max(zero, & !
494	sign(unun,-G1snSV(ikl,isn) & !
495	-epsi )) !
496
497	DENDRn=-G1snSV(ikl,isn) /G1_dSV
498	SPHERn= G2snSV(ikl,isn) /G1_dSV
499	DENDRn= DENDRn-vdent1facVapfrac_j
500	SPHERn= SPHERn+vsphe2facVapfrac_j
501
502	OK__DE=max(zero, & ! IF 1.,
503	sign(unun, DENDRn & ! NO change
504	-epsi )) ! Dendr. -> Spheric
505
506	G1_ldd= OK__DE * ( -DENDRn*G1_dSV) & ! Dendritic
507	+(1.-OK__DE)* min(G1_dSV,SPHERn*G1_dSV) ! Dendr. -> Spheric
508	G2_ldd= OK__DE * min(G1_dSV,SPHERn*G1_dSV) & ! Spheric
509	+(1.-OK__DE)*(ADSdSV-min(SPHERn,vsphe1)) ! Spher. -> Size
510
511	! +-- 1.3.1.2 Cas non dendritique / non dendritic Case
512
513	SPHERn=G1snSV(ikl,isn)/G1_dSV
514	DiamGx=G2snSV(ikl,isn)*0.1
515
516	istoOK=min( abs(istoSV(ikl,isn)- &
517	istdSV(1) ),1) ! zero if istoSV = 1
518	DiamOK=max(zero, sign(unun,vdiam2-DiamGx))
519	No_Big= istoOK+DiamOK
520	No_Big=min(No_Big,unun)
521
522	dSPHER= vsphe2facVapfrac_j !
523	SPHER0= SPHERn+dSPHER ! small grains
524	SPHbig= SPHERn+dSPHER & ! big grains
525	*exp(min(zero,vdiam3-G2snSV(ikl,isn))) ! (history = 2 or 3)
526	SPHbig= min(vsphe3,SPHbig) ! limited sphericity
527	SPHERn= No_Big * SPHER0 &
528	+ (1.-No_Big)* SPHbig
529
530	G1_lds= min(G1_dSV,SPHERn*G1_dSV)
531
532	! +-- 1.3.2 Calcul Metamorphose Gradient Moyen/Moderate (0.05-0.15)
533	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
534	OK_mdT=max(zero, & !
535	sign(unun, vgrat2 & !
536	-dTsndz)) !
537	OKmidT= OK_mdT *(1.-OKlowT) !
538	OKhigT= (1. -OK_mdT) *(1.-OKlowT) !
539
540	facVap=vdent1*exp(vvap1/TsisSV(ikl,isn)) &
541	* (1.e2 dTsndz)*vvap2
542
543	! +-- 1.3.2.1 cas dendritique / dendritic case.
544
545	OK_mdd=max(zero, & !
546	sign(unun,-G1snSV(ikl,isn) & !
547	-epsi )) !
548
549	DENDRn=-G1snSV(ikl,isn)/G1_dSV
550	SPHERn= G2snSV(ikl,isn)/G1_dSV
551	DENDRn= DENDRn - facVap*frac_j
552	SPHERn= SPHERn - facVap*frac_j
553
554	OK__DE=max(zero, & ! IF 1.,
555	sign(unun, DENDRn & ! NO change
556	-epsi )) ! Dendr. -> Spheric
557
558	G1_mdd= OK__DE * ( -DENDRn*G1_dSV) & ! Dendritic
559	+(1.-OK__DE)* max(zero ,SPHERn*G1_dSV) ! Dendr. -> Spheric
560	G2_mdd= OK__DE * max(zero ,SPHERn*G1_dSV) & ! Spheric
561	+(1.-OK__DE)*(ADSdSV-max(SPHERn,zero )) ! Spher. -> Size
562
563	! +-- 1.3.2.2 Cas non dendritique / non dendritic Case
564
565	SPHERn=G1snSV(ikl,isn)/G1_dSV
566	SPHERn= SPHERn-facVap*frac_j
567	G1_mds=max(zero,SPHERn*G1_dSV)
568
569	! +-- 1.3.3 Calcul Metamorphose fort / high Gradient
570	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
571	facVap=vdent1*exp(vvap1/TsisSV(ikl,isn)) &
572	* (1.e2 dTsndz)*vvap2
573
574	! +-- 1.3.3.1 Cas dendritique / dendritic Case
575
576	OK_hdd=max(zero, & !
577	sign(unun,-G1snSV(ikl,isn) & !
578	-epsi )) !
579
580	DENDRn=-G1snSV(ikl,isn)/G1_dSV !
581	SPHERn= G2snSV(ikl,isn)/G1_dSV !
582	DENDRn= DENDRn - facVap*frac_j !
583	SPHERn= SPHERn - facVap*frac_j ! Non dendritic
584	! + ! and angular
585	OK__DE=max(zero, & ! IF 1.,
586	sign(unun, DENDRn & ! NO change
587	-epsi )) ! Dendr. -> Spheric
588
589	G1_hdd= OK__DE * ( -DENDRn*G1_dSV) & ! Dendritic
590	+(1.-OK__DE)* max(zero ,SPHERn*G1_dSV) ! Dendr. -> Spheric
591	G2_hdd= OK__DE * max(zero ,SPHERn*G1_dSV) & ! Spheric
592	+(1.-OK__DE)*(ADSdSV-max(SPHERn,zero )) ! Spher. -> Size
593
594	! +-- 1.3.3.2 Cas non dendritique non completement anguleux.
595	! + non dendritic and spericity gt. 0
596
597	OK_hds=max(zero, & !
598	sign(unun, G1snSV(ikl,isn) & !
599	-epsi )) !
600
601	SPHERn= G1snSV(ikl,isn)/G1_dSV
602	SPHERn= SPHERn - facVap*frac_j
603	G1_hds= max(zero,SPHERn*G1_dSV)
604
605	! +-- 1.3.3.3 Cas non dendritique et anguleux
606	! + dendritic and spericity = 0.
607
608	T1__OK = max(zero,sign(unun,TsisSV(ikl,isn)-TfSnow+vtang1))
609	T2__OK = max(zero,sign(unun,TsisSV(ikl,isn)-TfSnow+vtang2))
610	T3_xOK = max(zero,sign(unun,TsisSV(ikl,isn)-TfSnow+vtang3))
611	T3__OK = T3_xOK * (1. - T2__OK)
612	T3_nOK = (1. - T3_xOK) * (1. - T2__OK)
613	ro1_OK = max(zero,sign(unun,vrang1-ro_dry(ikl,isn)))
614	ro2_OK = max(zero,sign(unun,ro_dry(ikl,isn)-vrang2))
615	dT1_OK = max(zero,sign(unun,vgang1-dTsndz ))
616	dT2_OK = max(zero,sign(unun,vgang2-dTsndz ))
617	dT3xOK = max(zero,sign(unun,vgang3-dTsndz ))
618	dT3_OK = dT3xOK * (1. - dT2_OK)
619	dT4xOK = max(zero,sign(unun,vgang4-dTsndz ))
620	dT4_OK = dT4xOK * (1. - dT3_OK) &
621	* (1. - dT2_OK)
622	dT4nOK = (1. - dT4xOK) * (1. - dT3_OK) &
623	* (1. - dT2_OK)
624
625	! +-- Influence de la Temperature /Temperature Influence
626	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
627	AngSno = &
628	T1__OK & ! 11
629	(T2__OK(vtang4+vtang5*(TfSnow -TsisSV(ikl,isn)) & ! 12
630	/vtang6) & !
631	+T3__OK(vtang7-vtang8(TfSnow-vtang2-TsisSV(ikl,isn)) & ! 13
632	/vtang9) & !
633	+T3_nOK(vtanga-vtangb(TfSnow-vtang3-TsisSV(ikl,isn)) & ! 14
634	/vtangc)) & !
635
636	! +-- Influence de la Masse Volumique /Density Influence
637	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
638	* ro1_OK &
639	( ro2_OK(1. - (ro_dry(ikl,isn)-vrang2) & !
640	/(vrang1-vrang2)) & !
641	+1.-ro2_OK ) & !
642
643	! +-- Influence du Gradient de Temperature /Temperature Gradient Influence
644	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
645	( dT1_OK(dT2_OKvgang5(dTsndz-vgang6) & ! 15
646	/(vgang2-vgang6) & !
647	+dT3_OK*vgang7 & ! 16
648	+dT4_OK*vgang9 & ! 17
649	+dT4nOK*vgangb ) & ! 18
650	+1.-dT1_OK ) & !
651	+ ro1_OK &
652	* dT1_OK(dT3_OKvgang8*(dTsndz-vgang2) &
653	/(vgang3-vgang2) &
654	+dT4_OKvganga(dTsndz-vgang3) &
655	/(vgang4-vgang3) &
656	+dT4nOKvgangc(dTsndz-vgang4) &
657	/(vgang1-vgang4))
658
659	G2_hds = G2snSV(ikl,isn) + 1.d2 AngSnovfi *frac_j
660
661
662	! +--New Properties
663	! + --------------
664
665	G1_bak = G1snSV(ikl,isn)
666	G2_bak = G2snSV(ikl,isn)
667
668	G1snSV(ikl,isn) = Wet_OK * ( OK__wd *G1__wd & ! 1
669	+(1.-OK__wd)* OK__ws *G1__ws & ! 2
670	+(1.-OK__wd)(1.-OK__ws)G1_bak) & ! 3
671	+(1. - Wet_OK) & !
672	( OKlowT ( OK_ldd *G1_ldd & ! 4
673	+(1.-OK_ldd) *G1_lds) & ! 5
674	+ OKmidT ( OK_mdd G1_mdd & ! 6
675	+(1.-OK_mdd) *G1_mds) & ! 7
676	+ OKhigT ( OK_hdd G1_hdd & ! 8
677	+(1.-OK_hdd)* OK_hds *G1_hds & ! 9
678	+(1.-OK_hdd)(1.-OK_hds)G1_bak)) ! 10
679
680	!XF
681	if(G1snSV(ikl,isn)<0.1) &
682	G2_hds = G2snSV(ikl,isn) + 1.d1 AngSnovfi *frac_j
683	!XF
684
685
686	G2snSV(ikl,isn) = Wet_OK * ( OK__wd *G2__wd & ! 1
687	+(1.-OK__wd)* OK__ws *G2_bak & ! 2
688	+(1.-OK__wd)(1.-OK__ws)G2__ws) & ! 3
689	+(1. - Wet_OK) & !
690	( OKlowT ( OK_ldd *G2_ldd & ! 4
691	+(1.-OK_ldd) *G2_bak) & ! 5
692	+ OKmidT ( OK_mdd G2_mdd & ! 6
693	+(1.-OK_mdd) *G2_bak) & ! 7
694	+ OKhigT ( OK_hdd G2_hdd & ! 8
695	+(1.-OK_hdd)* OK_hds *G2_bak & ! 9
696	+(1.-OK_hdd)(1.-OK_hds)G2_hds)) ! 10
697
698	! +--Snow Properties: IO Set Up
699	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~
700	! #vp G_curr( 1) = Wet_OK * OK__wd
701	! #vp G_curr( 2) = Wet_OK (1.-OK__wd) OK__ws
702	! #vp G_curr( 3) = Wet_OK (1.-OK__wd)(1.-OK__ws)
703	! #vp G_curr( 4) = (1.-Wet_OK)* OKlowT * OK_ldd
704	! #vp G_curr( 5) = (1.-Wet_OK)* OKlowT *(1.-OK_ldd)
705	! #vp G_curr( 6) = (1.-Wet_OK)* OKmidT * OK_mdd
706	! #vp G_curr( 7) = (1.-Wet_OK)* OKmidT *(1.-OK_mdd)
707	! #vp G_curr( 8) = (1.-Wet_OK)* OKhigT * OK_hdd
708	! #vp G_curr( 9) = (1.-Wet_OK)* OKhigT (1.-OK_hdd) OK_hds
709	! #vp G_curr(10) = (1.-Wet_OK)* OKhigT (1.-OK_hdd)(1.-OK_hds)
710	! #vp G_curr(11) = T1__OK * G_curr(10)
711	! #vp G_curr(12) = T2__OK * G_curr(10)
712	! #vp G_curr(13) = T3__OK * G_curr(10)
713	! #vp G_curr(14) = T3_nOK * G_curr(10)
714	! #vp G_curr(15) = ro1_OK* dT1_OK * dT2_OK * G_curr(10)
715	! #vp G_curr(16) = ro1_OK* dT1_OK * dT3_OK * G_curr(10)
716	! #vp G_curr(17) = ro1_OK* dT1_OK * dT4_OK * G_curr(10)
717	! #vp G_curr(18) = ro1_OK* dT1_OK * dT4nOK * G_curr(10)
718
719	! #vp Gcases( 1) = max(Gcases( 1),G_curr( 1))
720	! #vp Gcases( 2) = max(Gcases( 2),G_curr( 2))
721	! #vp Gcases( 3) = max(Gcases( 3),G_curr( 3))
722	! #vp Gcases( 4) = max(Gcases( 4),G_curr( 4))
723	! #vp Gcases( 5) = max(Gcases( 5),G_curr( 5))
724	! #vp Gcases( 6) = max(Gcases( 6),G_curr( 6))
725	! #vp Gcases( 7) = max(Gcases( 7),G_curr( 7))
726	! #vp Gcases( 8) = max(Gcases( 8),G_curr( 8))
727	! #vp Gcases( 9) = max(Gcases( 9),G_curr( 9))
728	! #vp Gcases(10) = max(Gcases(10),G_curr(10))
729	! #vp Gcases(11) = max(Gcases(11),G_curr(11))
730	! #vp Gcases(12) = max(Gcases(12),G_curr(12))
731	! #vp Gcases(13) = max(Gcases(13),G_curr(13))
732	! #vp Gcases(14) = max(Gcases(14),G_curr(14))
733	! #vp Gcases(15) = max(Gcases(15),G_curr(15))
734	! #vp Gcases(16) = max(Gcases(16),G_curr(16))
735	! #vp Gcases(17) = max(Gcases(17),G_curr(17))
736	! #vp Gcases(18) = max(Gcases(18),G_curr(18))
737
738	! +--Snow Properties: IO
739	! + ~~~~~~~~~~~~~~~~~~~
740	! #vp IF (isn .le. isnoSV(ikl))
741	! #vp. write(47,471)isn ,isnoSV(ikl) ,
742	! #vp. TsisSV(ikl,isn),ro__SV(ikl,isn),eta_SV(ikl,isn),
743	! #vp. G1_bak ,G2_bak ,istoSV(ikl,isn),
744	! #vp. dTsndz,
745	! #vp. ( k ,k=1,18),
746	! #vp. (G_curr(k),k=1,18),
747	! #vp. (Gcases(k),k=1,18),
748	! #vp. Wet_OK,OK__wd,G1__wd,G2__wd,
749	! #vp. 1.-OK__wd,OK__ws,G1__ws,1.-OK__ws,G2__ws,
750	! #vp. 1.-Wet_OK,OKlowT,OK_ldd,G1_ldd, G2_ldd,
751	! #vp. 1.-OK_ldd,G1_lds,
752	! #vp. OKmidT,OK_mdd,G1_mdd, G1_mdd,
753	! #vp. 1.-OK_mdd,G1_mds,
754	! #vp. OKhigT,OK_hdd,G1_hdd, G2_hdd,
755	! #vp. 1.-OK_hdd,OK_hds, G1_hds,
756	! #vp. 1.-OK_hds,G2_hds,
757	! #vp. G1snSV(ikl,isn),
758	! #vp. G2snSV(ikl,isn)
759
760	END DO
761	END DO
762	!!$OMP END PARALLEL DO
763
764	! +-- 2. Mise a Jour Variables Historiques (Cas non dendritique)
765	! + Update of the historical Variables
766	! + =======================================================
767
768	IF (vector) THEN
769	!XF
770	DO ikl=1,knonv
771	DO isn=1,isnoSV(ikl)
772	SphrOK = max(zero,sign(unun, G1snSV(ikl,isn)))
773	H1a_OK = max(zero,sign(unun,vsphe4-G1snSV(ikl,isn)))
774	H1b_OK = 1 - min(1 , istoSV(ikl,isn))
775	H1__OK = H1a_OK*H1b_OK
776	H23aOK = max(zero,sign(unun,vsphe4-G1_dSV &
777	+G1snSV(ikl,isn)))
778	H23bOK = max(zero,sign(unun,etaSno(ikl,isn) &
779	/max(epsi,dzsnSV(ikl,isn)) &
780	-vtelv1 ))
781	H23_OK = H23aOK*H23bOK
782	H2__OK = 1 - min(1 , istoSV(ikl,isn))
783	H3__OK = 1 - min(1 , abs(istoSV(ikl,isn)-istdSV(1)))
784	H45_OK = max(zero,sign(unun,TfSnow-TsisSV(ikl,isn)+epsi))
785	H4__OK = 1 - min(1 , abs(istoSV(ikl,isn)-istdSV(2)))
786	H5__OK = 1 - min(1 , abs(istoSV(ikl,isn)-istdSV(3)))
787
788	HISupd = &
789	SphrOK(H1__OK istdSV(1) &
790	+(1.-H1__OK)* H23_OK (H2__OKistdSV(2) &
791	+H3__OK*istdSV(3)) &
792	+(1.-H1__OK)(1.-H23_OK) H45_OK(H4__OKistdSV(4) &
793	+H5__OK*istdSV(5)))
794	istoSV(ikl,isn) = HISupd + &
795	(1.-min(unun,HISupd)) *istoSV(ikl,isn)
796	END DO
797	END DO
798	ELSE
799
800
801	! +-- 2. Mise a Jour Variables Historiques (Cas non dendritique)
802	! + Update of the historical Variables
803	! + =======================================================
804
805	DO ikl=1,knonv
806	DO isn=iiceSV(ikl),isnoSV(ikl)
807	IF (G1snSV(ikl,isn).ge.0.) THEN
808	IF(G1snSV(ikl,isn).lt.vsphe4.and.istoSV(ikl,isn).eq.0) THEN
809	istoSV(ikl,isn)=istdSV(1)
810	ELSEIF(G1_dSV-G1snSV(ikl,isn) .lt.vsphe4.and. &
811	etaSno(ikl,isn)/dzsnSV(ikl,isn).gt.vtelv1) THEN
812	IF (istoSV(ikl,isn).eq.0) &
813	istoSV(ikl,isn)= istdSV(2)
814	IF (istoSV(ikl,isn).eq.istdSV(1)) &
815	istoSV(ikl,isn)= istdSV(3)
816	ELSEIF(TsisSV(ikl,isn).lt.TfSnow) THEN
817	IF (istoSV(ikl,isn).eq.istdSV(2)) &
818	istoSV(ikl,isn)= istdSV(4)
819	IF (istoSV(ikl,isn).eq.istdSV(3)) &
820	istoSV(ikl,isn)= istdSV(5)
821	END IF
822	END IF
823	END DO
824	END DO
825	END IF
826
827
828	! +-- 3. Tassement mecanique /mechanical Settlement
829	! + ==========================================
830
831	DO ikl=1,knonv
832	SnMass(ikl) = 0.
833	END DO
834	!XF
835	DO ikl=1,knonv
836
837	smb_old = 0.
838	zn_old = 0
839	DO isn = 1, isnoSV(ikl)
840	smb_old = smb_old + dzsnSV(ikl,isn) *ro__SV(ikl,isn)
841	zn_old = zn_old + dzsnSV(ikl,isn)
842	ENDDO
843
844	DO isn=isnoSV(ikl),1,-1
845	dSnMas = 100.dzsnSV(ikl,isn)ro_dry(ikl,isn)
846	SnMass(ikl)= SnMass(ikl)+0.5*dSnMas
847	ViscSn = vvisc1 *vvisc2 &
848	exp(vvisc3 ro_dry(ikl,isn) &
849	+vvisc4*abs(TfSnow-TsisSV(ikl,isn))) &
850	*ro_dry(ikl,isn)/rovisc
851
852	! +-- Changement de Viscosite si Teneur en Eau liquide
853	! + Change of the Viscosity if liquid Water Content
854	! + ------------------------------------------------
855
856	OK_Liq = max(zero,sign(unun,etaSno(ikl,isn)-epsi))
857	OK_Ang = max(zero,sign(unun,vgran6-G1snSV(ikl,isn))) &
858	*(1-min(1 , abs(istoSV(ikl,isn)-istdSV(1))))
859	! #wp IF (G1snSV(ikl,isn).gt.0..AND.G1snSV(ikl,isn).lt.vsphe4
860	! #wp. .AND.istoSV(ikl,isn).eq. 0)
861	! #wp. THEN
862	! #wp write(6,*) ikl,isn,' G1,G2,hist,OK_Ang ',
863	! #wp. G1snSV(ikl,isn), G2snSV(ikl,isn),istoSV(ikl,isn),OK_Ang
864	! #wp stop "Grains anguleux mal d?finis"
865	! #wp END IF
866	OKxLiq = max(zero,sign(unun,vtelv1-etaSno(ikl,isn) &
867	/max(epsi,dzsnSV(ikl,isn)))) &
868	* max(0 ,sign(1 ,istoSV(ikl,isn) &
869	-istdSV(1) ))
870	ViscSn = &
871	ViscSn( OK_Liq/(vvisc5+vvisc6etaSno(ikl,isn) &
872	/max(epsi,dzsnSV(ikl,isn))) &
873	+(1.-OK_Liq) ) &
874	( OK_Angexp(min(ADSdSV,G2snSV(ikl,isn)-vdiam4)) &
875	+(1.-OK_Ang) ) &
876	( OKxLiq vvisc7 &
877	+(1.-OKxLiq) )
878
879
880	! +-- Calcul nouvelle Epaisseur / new Thickness
881	! + -----------------------------------------
882
883	dzsnew = &
884	dzsnSV(ikl,isn) &
885	*max(vdz3, &
886	(unun-dt__SV2max(SnMass(ikl)cos(slopSV(ikl)),unun) &
887	/max(ViscSn ,epsi)))
888	rosnew = ro__SV(ikl,isn) *dzsnSV(ikl,isn) &
889	/max(1e-10,dzsnew)
890	rosmax = 1. /( (1. -eta_SV(ikl,isn)) /ro_Ice &
891	+ eta_SV(ikl,isn) /ro_Wat)
892	rosnew = min(rosnew ,rosmax)
893	dzsnew = dzsnSV(ikl,isn) *ro__SV(ikl,isn) &
894	/max(1e-10,rosnew)
895	ro__SV(ikl,isn)= rosnew
896	dzsnSV(ikl,isn)= dzsnew
897	ro_dry(ikl,isn)= ro__SV(ikl,isn)(1.-eta_SV(ikl,isn))1.e-3
898	! +... ro_dry: Dry Density (g/cm3)
899	! +
900	SnMass(ikl) = SnMass(ikl)+dSnMas*0.5
901	END DO
902
903	smb_new = 0.
904	DO isn = 1, isnoSV(ikl)
905	smb_new = smb_new + dzsnSV(ikl,isn) *ro__SV(ikl,isn)
906	ENDDO
907
908	isn=1
909	if (dzsnSV(ikl,isn)>0.and.ro__SV(ikl,isn)>0) then
910	dzsnSV(ikl,isn) = dzsnSV(ikl,isn) +0.9999*(smb_old-smb_new) &
911	/ ro__SV(ikl,isn)
912	endif
913
914	zn_new = 0
915	DO isn = 1, isnoSV(ikl)
916	zn_new = zn_new + dzsnSV(ikl,isn)
917	ENDDO
918	zn4_SV(ikl) = zn4_SV(ikl) + (zn_new - zn_old)
919
920	END DO
921
922
923
924	return
925	end subroutine sisvat_gsn

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format