Context Navigation

inlandsis.F @ 3792

Last change on this file since 3792 was 3792, checked in by evignon, 3 years ago
Ajout de INLANDSIS, nouvelle interface entre LMDZ et la neige de SISVAT Etienne, 04/01/2021
File size: 54.6 KB

Line
1	subroutine INLANDSIS(SnoMod,BloMod,jjtime)
2
3	USE dimphy
4
5	!--------------------------------------------------------------------------+
6	! INLANDSIS module \|
7	! Simplified SISVAT module, containing ice and snow processes for \|
8	! ice-covered surfaces \|
9	! version MARv3, november 2020 \|
10	! SubRoutine INLANDSIS contains the fortran 77 code of the \|
11	! Soil/Ice Snow Vegetation Atmosphere Transfer Scheme \|
12	! \|
13	!--------------------------------------------------------------------------+
14	! PARAMETERS: klonv: Total Number of columns = \|
15	! ^^^^^^^^^^ = Total Number of continental grid boxes \|
16	! X Number of Mosaic Cell per grid box \|
17	! \|
18	! INPUT: daHost : Date Host Model \|
19	! ^^^^^ \|
20	! \|
21	! INPUT: LSmask : 1: Land MASK \|
22	! ^^^^^ 0: Sea MASK \|
23	! isotSV = 0,...,12: Soil Type \|
24	! 0: Water, Liquid (Sea, Lake) \|
25	! 12: Water, Solid (Ice) \|
26	! \|
27	! INPUT: coszSV : Cosine of the Sun Zenithal Distance [-] \|
28	! ^^^^^ sol_SV : Surface Downward Solar Radiation [W/m2] \|
29	! IRd_SV : Surface Downward Longwave Radiation [W/m2] \|
30	! drr_SV : Rain Intensity [kg/m2/s] \|
31	! dsn_SV : Snow Intensity [mm w.e./s] \|
32	! dsnbSV : Snow Intensity, Drift Fraction [-] \|
33	! dbs_SV : Drift Amount [mm w.e.] \|
34	! za__SV : Surface Boundary Layer (SBL) Height [m] \|
35	! VV__SV :(SBL Top) Wind Velocity [m/s] \|
36	! TaT_SV : SBL Top Temperature [K] \|
37	! rhT_SV : SBL Top Air Density [kg/m3] \|
38	! QaT_SV : SBL Top Specific Humidity [kg/kg] \|
39	! qsnoSV : SBL Mean Snow Content [kg/kg] \|
40	! alb0SV : Soil Basic Albedo [-] \|
41	! slopSV : Surface Slope [-] \|
42	! dt__SV : Time Step [s] \|
43	! \|
44	! INPUT / isnoSV = total Nb of Ice/Snow Layers \|
45	! OUTPUT: ispiSV = 0,...,nsno: Uppermost Superimposed Ice Layer \|
46	! ^^^^^^ iiceSV = total Nb of Ice Layers \|
47	! istoSV = 0,...,5 : Snow History (see istdSV data) \|
48	! \|
49	! INPUT / alb_SV : Surface Albedo [-] \|
50	! OUTPUT: emi_SV : Surface Emissivity [-] \|
51	! ^^^^^^ IRs_SV : Soil IR Flux (negative) [W/m2] \|
52	! LMO_SV : Monin-Obukhov Scale [m] \|
53	! us__SV : Friction Velocity [m/s] \|
54	! uts_SV : Temperature Turbulent Scale [m/s] \|
55	! uqs_SV : Specific Humidity Velocity [m/s] \|
56	! uss_SV : Blowing Snow Turbulent Scale [m/s] \|
57	! usthSV : Blowing Snow Erosion Threshold [m/s] \|
58	! Z0m_SV : Momentum Roughness Length [m] \|
59	! Z0mmSV : Momentum Roughness Length (time mean) [m] \|
60	! Z0mnSV : Momentum Roughness Length (instantaneous)[m] \|
61	! Z0SaSV : Sastrugi Roughness Length [m] \|
62	! Z0e_SV : Erosion Snow Roughness Length [m] \|
63	! Z0emSV : Erosion Snow Roughness Length (time mean) [m] \|
64	! Z0enSV : Erosion Snow Roughness Length (instantaneous)[m] \|
65	! Z0roSV : Subgrid Topo Roughness Length [m] \|
66	! Z0h_SV : Heat Roughness Length [m] \|
67	! TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)\|
68	! & Snow Temperatures (layers 1,2,...,nsno) [K] \|
69	! ro__SV : Soil/Snow Volumic Mass [kg/m3] \|
70	! eta_SV : Soil/Snow Water Content [m3/m3] \|
71	! G1snSV : snow dendricity/sphericity \|
72	! G2snSV : snow sphericity/grain size \|
73	! dzsnSV : Snow Layer Thickness [m] \|
74	! agsnSV : Snow Age [day] \|
75	! BufsSV : Snow Buffer Layer [kg/m2] .OR. [mm] \|
76	! BrosSV : Snow Buffer Layer Density [kg/m3] \|
77	! BG1sSV : Snow Buffer Layer Dendricity / Sphericity [-] \|
78	! BG2sSV : Snow Buffer Layer Sphericity / Size [-] [0.1 mm] \|
79	! rusnSV : Surficial Water [kg/m2] .OR. [mm] \|
80	! \|
81	! OUTPUT: no__SV : OUTPUT file Unit Number [-] \|
82	! ^^^^^^ i___SV : OUTPUT point i Coordinate [-] \|
83	! j___SV : OUTPUT point j Coordinate [-] \|
84	! n___SV : OUTPUT point n Coordinate [-] \|
85	! lwriSV : OUTPUT point vec Index [-] \|
86	! \|
87	! OUTPUT: IRu_SV : Upward IR Flux (+, upw., effective) [K] \|
88	! ^^^^^^ hSalSV : Saltating Layer Height [m] \|
89	! qSalSV : Saltating Snow Concentration [kg/kg] \|
90	! RnofSV : RunOFF Intensity [kg/m2/s] \|
91	! \|
92	! Internal Variables: \|
93	! ^^^^^^^^^^^^^^^^^^ \|
94	! NLaysv = New Snow Layer Switch [-] \|
95	! albisv : Snow/Ice/Water/Soil Integrated Albedo [-] \|
96	! SoSosv : Absorbed Solar Radiation by Surfac.(Normaliz)[-] \|
97	! TBr_sv : Brightness Temperature [K] \|
98	! IRupsv : Upward IR Flux (-, upw.) [W/m2] \|
99	! ram_sv : Aerodynamic Resistance for Momentum [s/m] \|
100	! rah_sv : Aerodynamic Resistance for Heat [s/m] \|
101	! Evp_sv : Evaporation [kg/m2] \|
102	! EvT_sv : Evapotranspiration [kg/m2] \|
103	! HSs_sv : Surface Sensible Heat Flux + => absorb.[W/m2] \|
104	! HLs_sv : Surface Latent Heat Flux + => absorb.[W/m2] \|
105	! Lx_H2O : Latent Heat of Vaporization/Sublimation [J/kg] \|
106	! Tsrfsv : Surface Temperature [K] \|
107	! sEX_sv : Verticaly Integrated Extinction Coefficient [-] \|
108	! LSdzsv : Vertical Discretization Factor [-] \|
109	! = 1. Soil \|
110	! = 1000. Ocean \|
111	! z_snsv : Snow Pack Thickness [m] \|
112	! zzsnsv : Snow Pack Thickness [m] \|
113	! albssv : Soil Albedo [-] \|
114	! Eso_sv : Soil+Snow Emissivity [-] \|
115	! Khydsv : Soil Hydraulic Conductivity [m/s] \|
116	! \|
117	! ETSo_0 : Snow/Soil Energy Power, before Forcing [W/m2] \|
118	! ETSo_1 : Snow/Soil Energy Power, after Forcing [W/m2] \|
119	! ETSo_d : Snow/Soil Energy Power Forcing [W/m2] \|
120	! EqSn_0 : Snow Energy, before Phase Change [J/m2] \|
121	! EqSn_1 : Snow Energy, after Phase Change [J/m2] \|
122	! EqSn_d : Snow Energy, net Forcing [J/m2] \|
123	! Enrsvd : SVAT Energy Power Forcing [W/m2] \|
124	! Enrbal : SVAT Energy Balance [W/m2] \|
125	! Wats_0 : Soil Water, before Forcing [mm] \|
126	! Wats_1 : Soil Water, after Forcing [mm] \|
127	! Wats_d : Soil Water Forcing [mm] \|
128	! SIWm_0 : Snow initial Mass [mm w.e.] \|
129	! SIWm_1 : Snow final Mass [mm w.e.] \|
130	! SIWa_i : Snow Atmos. initial Forcing [mm w.e.] \|
131	! SIWa_f : Snow Atmos. final Forcing(noConsumed)[mm w.e.] \|
132	! SIWe_i : SnowErosion initial Forcing [mm w.e.] \|
133	! SIWe_f : SnowErosion final Forcing(noConsumed)[mm w.e.] \|
134	! SIsubl : Snow sublimed/deposed Mass [mm w.e.] \|
135	! SImelt : Snow Melted Mass [mm w.e.] \|
136	! SIrnof : Surficial Water + Run OFF Change [mm w.e.] \|
137	! SIvAcr : Sea-Ice vertical Acretion [mm w.e.] \|
138	! Watsvd : SVAT Water Forcing [mm] \|
139	! Watbal : SVAT Water Balance [W/m2] \|
140	! \|
141	! vk2 : Square of Von Karman Constant [-] \|
142	! sqrCm0 : Factor of Neutral Drag Coeffic.Momentum [s/m] \|
143	! sqrCh0 : Factor of Neutral Drag Coeffic.Heat [s/m] \|
144	! EmiSol : Soil Emissivity [-] \|
145	! EmiSno : Snow Emissivity [-] \|
146	! EmiWat : Water Emissivity [-] \|
147	! Z0mLnd : Land Roughness Length [m] \|
148	! sqrrZ0 : ut/u \|
149	! f_eff : Marticorena & B. 1995 JGR (20) \|
150	! A_Fact : Fundamental * Roughness \|
151	! Z0mBSn : BSnow Roughness Length [m] \|
152	! Z0mBS0 : Mimimum BSnow Roughness Length (blown* ) [m] \|
153	! Z0m_Sn : Snow Roughness Length (surface) [m] \|
154	! Z0m_S0 : Mimimum Snow Roughness Length [m] \|
155	! Z0m_S1 : Maximum Snow Roughness Length [m] \|
156	! Z0_GIM : Minimum GIMEX Roughness Length [m] \|
157	! Z0_ICE : Sea Ice ISW Roughness Length [m] \|
158	! \|
159	! \|
160	!--------------------------------------------------------------------------+
161
162
163
164	! Global Variables
165	! ================
166
167
168	USE VARphy
169	USE VAR_SV
170	USE VARdSV
171	USE VAR0SV
172	USE VARxSV
173	USE VARySV
174	USE VARtSV
175	USE surface_data, only: iflag_tsurf_inlandsis
176
177
178	IMPLICIT NONE
179
180	logical SnoMod
181	logical BloMod
182	integer jjtime
183
184
185	! Internal Variables
186	! ==================
187
188	! Non Local
189	! ---------
190
191	real TBr_sv(klonv) ! Brightness Temperature
192	real IRdwsv(klonv) ! DOWNward IR Flux
193	real IRupsv(klonv) ! UPward IR Flux
194	real d_Bufs,Bufs_N ! Buffer Snow Layer Increment
195	real Buf_ro,Bros_N ! Buffer Snow Layer Density
196	real BufPro ! Buffer Snow Layer Density
197	real Buf_G1,BG1__N ! Buffer Snow Layer Dendr/Sphe[-]
198	real Buf_G2,BG2__N ! Buffer Snow Layer Spher/Size[-]
199	real Bdzssv(klonv) ! Buffer Snow Layer Thickness
200	real z_snsv(klonv) ! Snow-Ice, current Thickness
201
202
203
204	! Local
205	! -----
206
207	integer iwr
208	integer ikl ,isn ,isl ,ist !
209	integer ist__s,ist__w ! Soil/Water Body Identifier
210	integer growth ! Seasonal Mask
211	integer LISmsk ! Land+Ice / Open Sea Mask
212	integer LSnMsk ! Snow-Ice / No Snow-Ice Mask
213	integer IceMsk,IcIndx(klonv) ! Ice / No Ice Mask
214	integer SnoMsk ! Snow / No Snow Mask
215
216	real roSMin,roSMax,roSn_1,roSn_2,roSn_3 ! Fallen Snow Density (PAHAUT)
217	real Dendr1,Dendr2,Dendr3 ! Fallen Snow Dendric.(GIRAUD)
218	real Spher1,Spher2,Spher3,Spher4 ! Fallen Snow Spheric.(GIRAUD)
219	real Polair ! Polar Snow Switch
220	real PorSno,Por_BS,Salt_f,PorRef !
221	c #sw real PorVol,rWater !
222	c #sw real rusNEW,rdzNEW,etaNEW !
223	real ro_new !
224	real TaPole ! Maximum Polar Temperature
225	real T__Min ! Minimum realistic Temperature
226	real EmiSol ! Emissivity of Soil
227	real EmiSno ! Emissivity of Snow
228	real EmiWat ! Emissivity of a Water Area
229	real vk2 ! Square of Von Karman Constant
230	real u2star !(u)*2
231	real Z0mLnd ! Land Roughness Length
232	c #ZN real sqrrZ0 ! ut/u
233	real f_eff ! Marticorena & B. 1995 JGR (20)
234	real A_Fact ! Fundamental * Roughness
235	real Z0m_nu ! Smooth R Snow Roughness Length
236	real Z0mBSn ! BSnow Roughness Length
237	real Z0mBS0 ! Mimimum BSnow Roughness Length
238	real Z0m_S0 ! Mimimum Snow Roughness Length
239	real Z0m_S1 ! Maximum Snow Roughness Length
240	c #SZ real Z0Sa_N ! Regime Snow Roughness Length
241	c #SZ real Z0SaSi ! 1.IF Rgm Snow Roughness Length
242	c #GL real Z0_GIM ! Mimimum GIMEX Roughness Length
243	real Z0_ICE ! Ice ISW Roughness Length
244	real Z0m_Sn,Z0m_90 ! Snow Surface Roughness Length
245	real SnoWat ! Snow Layer Switch
246	c #RN real rstar,alors !
247	c #RN real rstar0,rstar1,rstar2 !
248	real SameOK ! 1. => Same Type of Grains
249	real G1same ! Averaged G1, same Grains
250	real G2same ! Averaged G2, same Grains
251	real typ__1 ! 1. => Lay1 Type: Dendritic
252	real zroNEW ! dz X ro, if fresh Snow
253	real G1_NEW ! G1, if fresh Snow
254	real G2_NEW ! G2, if fresh Snow
255	real zroOLD ! dz X ro, if old Snow
256	real G1_OLD ! G1, if old Snow
257	real G2_OLD ! G2, if old Snow
258	real SizNEW ! Size, if fresh Snow
259	real SphNEW ! Spheric.,if fresh Snow
260	real SizOLD ! Size, if old Snow
261	real SphOLD ! Spheric.,if old Snow
262	real Siz_av ! Averaged Grain Size
263	real Sph_av ! Averaged Grain Spher.
264	real Den_av ! Averaged Grain Dendr.
265	real DendOK ! 1. => Average is Dendr.
266	real G1diff ! Averaged G1, diff. Grains
267	real G2diff ! Averaged G2, diff. Grains
268	real G1 ! Averaged G1
269	real G2 ! Averaged G2
270	real param ! Polynomial fit z0=f(T)
271	real Z0_obs ! Fit Z0_obs=f(T) (m)
272	real tamin ! min T of linear fit (K)
273	real tamax ! max T of linear fit (K)
274	real coefa,coefb,coefc,coefd ! Coefs for z0=f(T)
275	real ta1,ta2,ta3 ! Air temperature thresholds
276	real z01,z02,z03 ! z0 thresholds
277	real tt_c,vv_c ! Critical param.
278	real tt_tmp,vv_tmp,vv_virt ! Temporary variables
279	logical density_kotlyakov ! .true. if Kotlyakov 1961
280	real e_prad,e1pRad,A_Rad0,absg_V,absgnI,exdRad ! variables for SoSosv calculations
281	real zm1, zm2, coefslope ! variables for surface temperature extrapolation
282
283
284
285	! Internal DATA
286	! =============
287
288	data T__Min / 200.00/ ! Minimum realistic Temperature
289	data TaPole / 263.15/ ! Maximum Polar Temperature
290	data roSMin / 30. / ! Minimum Snow Density
291	data roSMax / 400. / ! Max Fresh Snow Density
292	data tt_c / -2.0 / ! Critical Temp. (degC)
293	data vv_c / 14.3 / ! Critical Wind speed (m/s)
294	data roSn_1 / 109. / ! Fall.Sno.Density, Indep. Param.
295	data roSn_2 / 6. / ! Fall.Sno.Density, Temper.Param.
296	data roSn_3 / 26. / ! Fall.Sno.Density, Wind Param.
297	data Dendr1 / 17.12/ ! Fall.Sno.Dendric.,Wind 1/Param.
298	data Dendr2 / 128. / ! Fall.Sno.Dendric.,Wind 2/Param.
299	data Dendr3 / -20. / ! Fall.Sno.Dendric.,Indep. Param.
300	data Spher1 / 7.87/ ! Fall.Sno.Spheric.,Wind 1/Param.
301	data Spher2 / 38. / ! Fall.Sno.Spheric.,Wind 2/Param.
302	data Spher3 / 50. / ! Fall.Sno.Spheric.,Wind 3/Param.
303	data Spher4 / 90. / ! Fall.Sno.Spheric.,Indep. Param.
304	data EmiSol / 0.99999999/ ! 0.94Emissivity of Soil
305	data EmiWat / 0.99999999/ ! Emissivity of a Water Area
306	data EmiSno / 0.99999999/ ! Emissivity of Snow
307
308	! DATA Emissivities ! Pielke, 1984, pp. 383,409
309
310	data Z0mBS0 / 0.5e-6/ ! MINimum Snow Roughness Length
311	! for Momentum if Blowing Snow
312	! Gallee et al. 2001 BLM 99 (19)
313	data Z0m_S0/ 0.00005/ ! MINimum Snow Roughness Length
314	! MegaDunes included
315	data Z0m_S1/ 0.030 / ! MAXimum Snow Roughness Length
316	! (Sastrugis)
317	c #GL data Z0_GIM/ 0.0013/ ! Ice Min Z0 = 0.0013 m (Broeke)
318	! ! Old Ice Z0 = 0.0500 m (Bruce)
319	! ! 0.0500 m (Smeets)
320	! ! 0.1200 m (Broeke)
321	data Z0_ICE/ 0.0010/ ! Sea-Ice Z0 = 0.0010 m (Andreas)
322	! ! (Ice Station Weddel -- ISW)
323	vk2 = vonKrm * vonKrm ! Square of Von Karman Constant
324
325
326	! BEGIN.main.
327	! ===========================
328
329
330
331
332	! "Soil" Humidity of Water Bodies
333	! ===============================
334
335	DO ikl=1,knonv
336
337	ist = isotSV(ikl) ! Soil Type
338	ist__s = min(ist, 1) ! 1 => Soil
339	ist__w = 1 - ist__s ! 1 => Water Body
340	DO isl=-nsol,0
341	eta_SV(ikl,isl) = eta_SV(ikl,isl) * ist__s ! Soil
342	. + etadSV(ist) * ist__w ! Water Body
343	END DO
344
345
346	! Vertical Discretization Factor
347	! ==============================
348
349	LSdzsv(ikl) = ist__s ! Soil
350	. + OcndSV * ist__w ! Water Body
351	END DO
352
353
354	! Blowing Particles Threshold Friction velocity
355	! =============================================
356
357	c #AE usthSV(ikl) = 1.0e+2
358	! END DO
359	!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
360
361
362
363
364	! Contribution of Snow to the Surface Snow Pack
365	! =============================================
366
367	IF (SnoMod) THEN
368
369
370
371	C +--Blowing Snow
372	C + ------------
373
374	IF (BloMod) then
375	if (klonv.eq.1) then
376	if(isnoSV(1).ge.2 .and.
377	. TsisSV(1,max(1,isnoSV(1)))<273. .and.
378	. ro__SV(1,max(1,isnoSV(1)))<500. .and.
379	. eta_SV(1,max(1,isnoSV(1)))<epsi) then
380	C + **********
381	call SISVAT_BSn
382	endif
383	else
384	call SISVAT_BSn
385	C + **********
386	endif
387	ENDIF
388
389
390
391	C +--Buffer Layer
392	C + ------------
393
394	DO ikl=1,knonv
395	c BufsSV(ikl) [mm w.e.] i.e, i.e., [kg/m2]
396	d_Bufs = max(dsn_SV(ikl) *dt__SV,0.) !
397	dsn_SV(ikl) = 0. !
398	Bufs_N = BufsSV(ikl) +d_Bufs !
399
400
401	C +--Snow Density
402	C + ^^^^^^^^^^^^
403	Polair = zero
404	c #NP Polair = max(zero, !
405	c #NP. sign(unun,TaPole !
406	c #NP. -TaT_SV(ikl))) !
407	Polair = max(zero, !
408	. sign(unun,TaPole !
409	. -TaT_SV(ikl))) !
410	Buf_ro = max( rosMin, ! Fallen Snow Density
411	. roSn_1+roSn_2* (TaT_SV(ikl)-TfSnow) ! [kg/m3]
412	. +roSn_3*sqrt( VV__SV(ikl))) ! Pahaut (CEN), Etienne: use wind speed at first model level instead of 10m wind
413	c #NP BufPro = max( rosMin, ! Fallen Snow Density
414	c #NP. 104. *sqrt( max( VV10SV(ikl)-6.0,0.0))) ! Kotlyakov (1961)
415
416	density_kotlyakov = .true.
417	c #AC density_kotlyakov = .false. !C.Agosta snow densisty as if BS is on b
418	c #BS density_kotlyakov = .false. !C.Amory BS 2018
419	C + ... Fallen Snow Density, Adapted for Antarctica
420	if (density_kotlyakov) then
421	tt_tmp = TaT_SV(ikl)-TfSnow
422	!vv_tmp = VV10SV(ikl)
423	vv_tmp=VV__SV(ikl) ! Etienne: use wind speed at first model level instead of 10m wind
424	C + ... [ A compromise between
425	C + ... Kotlyakov (1961) and Lenaerts (2012, JGR, Part1) ]
426	if (tt_tmp.ge.-10) then
427	BufPro = max( rosMin,
428	. 104. *sqrt( max( vv_tmp-6.0,0.0))) ! Kotlyakov (1961)
429	else
430	vv_virt = (tt_cvv_tmp+vv_c(tt_tmp+10))
431	. /(tt_c+tt_tmp+10)
432	BufPro = 104. *sqrt( max( vv_virt-6.0,0.0))
433	endif
434	else
435	C + ... [ density derived from observations of the first 50cm of
436	C + ... snow - cf. Rajashree Datta - and multiplied by 0.8 ]
437	C + ... C. Agosta, 2016-09
438	cc #SD BufPro = 149.2 + 6.84VV10SV(ikl) + 0.48Tsrfsv(ikl)
439	cc #SD BufPro = 125 + 14VV10SV(ikl) + 0.6Tsrfsv(ikl) !MAJ CK and CAm
440	! BufPro = 200 + 21 * VV10SV(ikl)!CK 29/07/19
441	BufPro = 200 + 21 * VV__SV(ikl)!Etienne: use wind speed at first model level instead of 10m wind
442	endif
443
444	Bros_N = (1. - Polair) * Buf_ro ! Temperate Snow
445	. + Polair * BufPro ! Polar Snow
446
447	Bros_N = max( 20.,max(rosMin, Bros_N))
448	Bros_N = min(400.,min(rosMax-1,Bros_N)) ! for dz_min in SISVAT_zSn
449
450
451	! Density of deposited blown snow
452	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
453
454	c #BS Bros_N = frsno
455	c #BS ro_new = ro__SV(ikl,max(1,isnoSV(ikl)))
456	c #BS ro_new = max(Bros_N,min(roBdSV,ro_new))
457	c #BS Fac = 1-((ro__SV(ikl,max(1,isnoSV(ikl)))
458	c #BS. -roBdSV)/(500.-roBdSV))
459	c #BS Fac = max(0.,min(1.,Fac))
460	c #BS dsnbSV(ikl) = Fac*dsnbSV(ikl)
461	c #BS Bros_N = Bros_N * (1.0-dsnbSV(ikl))
462	c #BS. + ro_new * dsnbSV(ikl)
463
464
465
466	! Time averaged Density of deposited blown Snow
467	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
468
469	BrosSV(ikl) =(Bros_N * d_Bufs !
470	. +BrosSV(ikl)* BufsSV(ikl))!
471	. / max(epsi,Bufs_N) !
472
473
474	C +-- S.Falling Snow Properties (computed as in SISVAT_zAg)
475	C + ^^^^^^^^^^^^^^^^^^^^^^^
476	Buf_G1 = max(-G1_dSV, ! Temperate Snow
477	. min(Dendr1*VV__SV(ikl)-Dendr2, ! Dendricity
478	. Dendr3 )) !
479	Buf_G2 = min( Spher4, ! Temperate Snow
480	. max(Spher1*VV__SV(ikl)+Spher2, ! Sphericity
481	. Spher3 )) !
482	Buf_G1 = (1. - Polair) * Buf_G1 ! Temperate Snow
483	. + Polair * G1_dSV ! Polar Snow
484	Buf_G2 = (1. - Polair) * Buf_G2 ! Temperate Snow
485	. + Polair * ADSdSV ! Polar Snow
486	G1 = Buf_G1 ! NO Blown Snow
487	G2 = Buf_G2 ! NO Blown Snow
488
489
490	! S.1. Meme Type de Neige / same Grain Type
491	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
492	c #BS SameOK = max(zero,
493	c #BS. sign(unun, Buf_G1 *G1_dSV
494	c #BS. - eps_21 ))
495	c #BS G1same = ((1.0-dsnbSV(ikl))Buf_G1+dsnbSV(ikl) G1_dSV)
496	c #BS G2same = ((1.0-dsnbSV(ikl))Buf_G2+dsnbSV(ikl) ADSdSV)
497	! Blowing Snow Properties: G1_dSV, ADSdSV
498
499	! S.2. Types differents / differents Types
500	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
501	c #BS typ__1 = max(zero,sign(unun,epsi-Buf_G1)) ! =1.=> Dendritic
502	c #BS zroNEW = typ__1 *(1.0-dsnbSV(ikl)) ! fract.Dendr.Lay.
503	c #BS. + (1.-typ__1) * dsnbSV(ikl) !
504	c #BS G1_NEW = typ__1 *Buf_G1 ! G1 of Dendr.Lay.
505	c #BS. + (1.-typ__1) *G1_dSV !
506	c #BS G2_NEW = typ__1 *Buf_G2 ! G2 of Dendr.Lay.
507	c #BS. + (1.-typ__1) *ADSdSV !
508	c #BS zroOLD = (1.-typ__1) *(1.0-dsnbSV(ikl)) ! fract.Spher.Lay.
509	c #BS. + typ__1 * dsnbSV(ikl) !
510	c #BS G1_OLD = (1.-typ__1) *Buf_G1 ! G1 of Spher.Lay.
511	c #BS. + typ__1 *G1_dSV !
512	c #BS G2_OLD = (1.-typ__1) *Buf_G2 ! G2 of Spher.Lay.
513	c #BS. + typ__1 *ADSdSV !
514	c #BS SizNEW = -G1_NEW *DDcdSV/G1_dSV ! Size Dendr.Lay.
515	c #BS. +(1.+G1_NEW /G1_dSV) !
516	c #BS. (G2_NEW DScdSV/G1_dSV !
517	c #BS. +(1.-G2_NEW /G1_dSV)*DFcdSV) !
518	c #BS SphNEW = G2_NEW /G1_dSV ! Spher.Dendr.Lay.
519	c #BS SizOLD = G2_OLD ! Size Spher.Lay.
520	c #BS SphOLD = G1_OLD /G1_dSV ! Spher.Spher.Lay.
521	c #BS Siz_av = (zroNEWSizNEW+zroOLDSizOLD) ! Averaged Size
522	c #BS Sph_av = min( zroNEWSphNEW+zroOLDSphOLD !
523	c #BS. , unun) ! Averaged Sphericity
524	c #BS Den_av = min((Siz_av -( Sph_av *DScdSV !
525	c #BS. +(1.-Sph_av)*DFcdSV)) !
526	c #BS. / (DDcdSV -( Sph_av *DScdSV !
527	c #BS. +(1.-Sph_av)*DFcdSV)) !
528	c #BS. , unun) !
529	c #BS DendOK = max(zero, !
530	c #BS. sign(unun, Sph_av *DScdSV ! Small Grains
531	c #BS. +(1.-Sph_av)*DFcdSV ! Faceted Grains
532	c #BS. - Siz_av )) !
533	C +... REMARQUE: le type moyen (dendritique ou non) depend
534	C + ^^^^^^^^ de la comparaison avec le diametre optique
535	C + d'une neige recente de dendricite nulle
536	C +... REMARK: the mean type (dendritic or not) depends
537	C + ^^^^^^ on the comparaison with the optical diameter
538	C + of a recent snow having zero dendricity
539
540	c #BS G1diff =( -DendOK *Den_av
541	c #BS. +(1.-DendOK)Sph_av) G1_dSV
542	c #BS G2diff = DendOK Sph_av G1_dSV
543	c #BS. +(1.-DendOK)*Siz_av
544	c #BS G1 = SameOK *G1same
545	c #BS. +(1.-SameOK)*G1diff
546	c #BS G2 = SameOK *G2same
547	c #BS. +(1.-SameOK)*G2diff
548
549
550
551	! S.1. Meme Type de Neige / same Grain Type
552	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
553	SameOK = max(zero,
554	. sign(unun, Buf_G1 *BG1sSV(ikl)
555	. - eps_21 ))
556	G1same = (d_BufsBuf_G1+BufsSV(ikl)BG1sSV(ikl))
557	. /max(epsi,Bufs_N)
558	G2same = (d_BufsBuf_G2+BufsSV(ikl)BG2sSV(ikl))
559	. /max(epsi,Bufs_N)
560
561	! S.2. Types differents / differents Types
562	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
563
564	typ__1 = max(zero,sign(unun,epsi-Buf_G1)) ! =1.=> Dendritic
565	zroNEW =( typ__1 *d_Bufs ! fract.Dendr.Lay.
566	. + (1.-typ__1) *BufsSV(ikl)) !
567	. /max(epsi,Bufs_N) !
568	G1_NEW = typ__1 *Buf_G1 ! G1 of Dendr.Lay.
569	. + (1.-typ__1) *BG1sSV(ikl) !
570	G2_NEW = typ__1 *Buf_G2 ! G2 of Dendr.Lay.
571	. + (1.-typ__1) *BG2sSV(ikl) !
572	zroOLD =((1.-typ__1) *d_Bufs ! fract.Spher.Lay.
573	. + typ__1 *BufsSV(ikl)) !
574	. /max(epsi,Bufs_N) !
575	G1_OLD = (1.-typ__1) *Buf_G1 ! G1 of Spher.Lay.
576	. + typ__1 *BG1sSV(ikl) !
577	G2_OLD = (1.-typ__1) *Buf_G2 ! G2 of Spher.Lay.
578	. + typ__1 *BG2sSV(ikl) !
579	SizNEW = -G1_NEW *DDcdSV/G1_dSV ! Size Dendr.Lay.
580	. +(1.+G1_NEW /G1_dSV) !
581	. (G2_NEW DScdSV/G1_dSV !
582	. +(1.-G2_NEW /G1_dSV)*DFcdSV) !
583	SphNEW = G2_NEW /G1_dSV ! Spher.Dendr.Lay.
584	SizOLD = G2_OLD ! Size Spher.Lay.
585	SphOLD = G1_OLD /G1_dSV ! Spher.Spher.Lay.
586	Siz_av = ( zroNEW SizNEW+zroOLDSizOLD) ! Averaged Size
587	Sph_av = min( zroNEW SphNEW+zroOLDSphOLD !
588	. , unun ) ! Averaged Sphericity
589	Den_av = min((Siz_av - ( Sph_av *DScdSV !
590	. +(1.-Sph_av)*DFcdSV)) !
591	. / (DDcdSV - ( Sph_av *DScdSV !
592	. +(1.-Sph_av)*DFcdSV)) !
593	. , unun )!
594	DendOK = max(zero, !
595	. sign(unun, Sph_av *DScdSV ! Small Grains
596	. +(1.-Sph_av)*DFcdSV ! Faceted Grains
597	. - Siz_av )) !
598	C +... REMARQUE: le type moyen (dendritique ou non) depend
599	C + ^^^^^^^^ de la comparaison avec le diametre optique
600	C + d'une neige recente de dendricite nulle
601	C +... REMARK: the mean type (dendritic or not) depends
602	C + ^^^^^^ on the comparaison with the optical diameter
603	C + of a recent snow having zero dendricity
604
605	G1diff =( -DendOK *Den_av
606	. +(1.-DendOK)Sph_av) G1_dSV
607	G2diff = DendOK Sph_av G1_dSV
608	. +(1.-DendOK)*Siz_av
609	G1 = SameOK *G1same
610	. +(1.-SameOK)*G1diff
611	G2 = SameOK *G2same
612	. +(1.-SameOK)*G2diff
613
614	BG1sSV(ikl) = G1 !
615	. * Bufs_N/max(epsi,Bufs_N) !
616	BG2sSV(ikl) = G2 !
617	. * Bufs_N/max(epsi,Bufs_N) !
618
619
620	C +--Update of Buffer Layer Content & Decision about creating a new snow layer
621	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
622	BufsSV(ikl) = Bufs_N ! [mm w.e.]
623	NLaysv(ikl) = min(unun, !
624	. max(zero, ! Allows to create
625	. sign(unun,BufsSV(ikl) ! a new snow Layer
626	. -SMndSV )) ! if Buffer > SMndSV
627	. max(zero, ! Except if Erosion
628	. sign(unun,0.50 ! dominates
629	. -dsnbSV(ikl))) !
630	. +max(zero, ! Allows to create
631	. sign(unun,BufsSV(ikl) ! a new snow Layer
632	. -SMndSV3.00))) ! is Buffer > SMndSV3
633	Bdzssv(ikl) = 1.e-3BufsSV(ikl)ro_Wat ! [mm w.e.] -> [m w.e.]
634	. /max(epsi,BrosSV(ikl))!& [m w.e.] -> [m]
635
636
637
638	END DO
639
640
641
642	! Snow Pack Discretization
643	! ========================
644
645	! **********
646	call SISVAT_zSn
647	! **********
648
649	! **********
650	! #ve call SISVAT_wEq('_zSn ',0)
651	! **********
652
653
654
655	! Add a new Snow Layer
656	! ====================
657
658	DO ikl=1,knonv
659
660	isnoSV(ikl) = isnoSV(ikl) +NLaysv(ikl)
661	isn = isnoSV(ikl)
662	dzsnSV(ikl,isn) = dzsnSV(ikl,isn) * (1-NLaysv(ikl))
663	. + Bdzssv(ikl) * NLaysv(ikl)
664	TsisSV(ikl,isn) = TsisSV(ikl,isn) * (1-NLaysv(ikl))
665	. + min(TaT_SV(ikl),Tf_Sno) *NLaysv(ikl)
666
667	ro__SV(ikl,isn) = ro__SV(ikl,isn) * (1-NLaysv(ikl))
668	. + Brossv(ikl) * NLaysv(ikl)
669	eta_SV(ikl,isn) = eta_SV(ikl,isn) * (1-NLaysv(ikl)) ! + 0.
670	agsnSV(ikl,isn) = agsnSV(ikl,isn) * (1-NLaysv(ikl)) ! + 0.
671	G1snSV(ikl,isn) = G1snSV(ikl,isn) * (1-NLaysv(ikl))
672	. + BG1ssv(ikl) * NLaysv(ikl)
673	G2snSV(ikl,isn) = G2snSV(ikl,isn) * (1-NLaysv(ikl))
674	. + BG2ssv(ikl) * NLaysv(ikl)
675	istoSV(ikl,isn) = istoSV(ikl,isn) * (1-NLaysv(ikl))
676	. + max(zer0,sign(un_1,TaT_SV(ikl)
677	. -Tf_Sno-eps_21)) * istdSV(2)
678	. * NLaysv(ikl)
679	BufsSV(ikl) = BufsSV(ikl) * (1-NLaysv(ikl))
680	NLaysv(ikl) = 0
681
682
683	END DO
684
685
686	! Snow Pack Thickness
687	! -------------------
688
689	DO ikl=1,knonv
690	z_snsv(ikl) = 0.0
691	END DO
692	DO isn=1,nsno
693	DO ikl=1,knonv
694	z_snsv(ikl) = z_snsv(ikl) + dzsnSV(ikl,isn)
695	zzsnsv(ikl,isn) = z_snsv(ikl)
696	END DO
697	END DO
698
699
700
701	END IF
702
703
704
705	! Soil Albedo: Soil Humidity Correction
706	! ==========================================
707
708	! REFERENCE: McCumber and Pielke (1981), Pielke (1984)
709	! ^^^^^^^^^
710	DO ikl=1,knonv
711	albssv(ikl) =
712	. alb0SV(ikl) *(1.0-min(half,eta_SV( ikl,0)
713	. /etadSV(isotSV(ikl))))
714	! REMARK: Albedo of Water Surfaces (isotSV=0):
715	! ^^^^^^ alb0SV := 2 X effective value, while
716	! eta_SV := etadSV
717	END DO
718
719
720	! Snow Pack Optical Properties
721	! ============================
722
723	IF (SnoMod) THEN
724
725	! ******
726	call SnOptP(jjtime)
727	! ******
728
729	ELSE
730	DO ikl=1,knonv
731	sEX_sv(ikl,1) = 1.0
732	sEX_sv(ikl,0) = 0.0
733	albisv(ikl) = albssv(ikl)
734	END DO
735	END IF
736
737
738
739	! Soil optical properties
740	! =============================
741	!Etienne: as in inlandis we do not call vgopt, we need to define
742	!the albedo alb_SV and to calculate the
743	!absorbed Solar Radiation by Surfac (Normaliz)[-] SoSosv
744
745
746	DO ikl=1,klonv
747
748	e_pRad = 2.5 * coszSV(ikl) ! exponential argument,
749	! V/nIR radiation partitioning,
750	! DR97, 2, eqn (2.53) & (2.54)
751	e1pRad = 1.-exp(-e_pRad) ! exponential, V/nIR Rad. Part.
752	exdRad= 1.
753
754	! Visible Part of the Solar Radiation Spectrum (V, 0.4--0.7mi.m)
755	A_Rad0 = 0.25 + 0.697 * e1pRad ! Absorbed Vis. Radiation
756	absg_V = (1.-albisv(ikl))(A_Rad0exdRad) !
757
758	! Near-IR Part of the Solar Radiation Spectrum (nIR, 0.7--2.8mi.m)
759
760	A_Rad0 = 0.80 + 0.185 * e1pRad ! Absorbed nIR. Radiation
761	absgnI = (1.-albisv(ikl))(A_Rad0exdRad) !
762
763	SoSosv(ikl) = (absg_V+absgnI)*0.5d0
764
765	alb_SV(ikl) = albisv(ikl)
766
767	END DO
768
769	! **********
770	! #ve call SISVAT_wEq('SnOptP',0)
771	! **********
772
773
774	! Surface Emissivity (Etienne: simplified calculation for landice)
775	! =============================================================
776	!
777	DO ikl=1,knonv
778	LSnMsk = min( 1,isnoSV(ikl))
779	Eso_sv(ikl)= EmiSol(1-LSnMsk)+EmiSnoLSnMsk ! Sol+Sno Emissivity
780	emi_SV(ikl)= EmiSol(1-LSnMsk) + EmiSnoLSnMsk
781	END DO
782
783
784
785
786	! Upward IR (INPUT, from previous time step)
787	! ===================================================================
788
789	DO ikl=1,knonv
790	! #e1 Enrsvd(ikl) = - IRs_SV(ikl)
791	IRupsv(ikl) = IRs_SV(ikl)
792	END DO
793
794
795	! Turbulence
796	! ==========
797
798	! Latent Heat of Vaporization/Sublimation
799	! ---------------------------------------
800
801	DO ikl=1,knonv
802	SnoWat = min(isnoSV(ikl),0)
803	Lx_H2O(ikl) =
804	. (1.-SnoWat) * LhvH2O
805	. + SnoWat (LhsH2O (1.-eta_SV(ikl,isnoSV(ikl)))
806	. +LhvH2O * eta_SV(ikl,isnoSV(ikl)) )
807	END DO
808
809
810
811
812	! Aerodynamic Resistance
813	! ^^^^^^^^^^^^^^^^^^^^^^
814
815
816	DO ikl=1,knonv
817	ram_sv(ikl) = 1./(cdM_SV(ikl)*max(VV__SV(ikl),eps6))
818	rah_sv(ikl) = 1./(cdH_SV(ikl)*max(VV__SV(ikl),eps6))
819	END DO
820
821
822
823	! Soil Energy Balance
824	! =====================
825
826
827	if (iflag_tsurf_inlandsis .eq. 0) then
828
829	call SISVAT_TSo
830
831	else
832
833	call SISVAT_TS2
834
835	end if
836
837
838	! **********
839	! #ve call SISVAT_wEq('_TSo ',0)
840	! **********
841
842
843
844	! Soil Water Potential
845	! ------------------------
846
847	DO isl=-nsol,0
848	DO ikl=1,knonv
849	ist = isotSV(ikl) ! Soil Type
850	psi_sv(ikl,isl) = psidSV(ist) ! DR97, Eqn.(3.34)
851	. *(etadSV(ist) /max(eps6,eta_SV(ikl,isl))) !
852	. **bCHdSV(ist) !
853
854
855	! Soil Hydraulic Conductivity
856	! ---------------------------
857
858	Khydsv(ikl,isl) = s2__SV(ist) ! DR97, Eqn.(3.35)
859	. (eta_SV(ikl,isl)(2.bCHdSV(ist)+3.)) !
860	END DO
861	END DO
862
863
864	! Melting / Refreezing in the Snow Pack
865	! =====================================
866
867	IF (SnoMod) THEN
868
869	! **********
870	call SISVAT_qSn
871	! **********
872
873	! **********
874	! #ve call SISVAT_wEq('_qSn ',0)
875	! **********
876
877
878	! Snow Pack Thickness
879	! -------------------
880
881	DO ikl=1,knonv
882	z_snsv(ikl) = 0.0
883	END DO
884	DO isn=1,nsno
885	DO ikl=1,knonv
886	z_snsv(ikl) = z_snsv(ikl) + dzsnSV(ikl,isn)
887	zzsnsv(ikl,isn) = z_snsv(ikl)
888	END DO
889	END DO
890
891
892	! Energy in Excess is added to the first Soil Layer
893	! -------------------------------------------------
894	DO ikl=1,knonv
895	z_snsv(ikl) = max(zer0,
896	. sign(un_1,eps6-z_snsv(ikl)))
897	TsisSV(ikl,0) = TsisSV(ikl,0) + EExcsv(ikl)
898	. /(rocsSV(isotSV(ikl))
899	. +rcwdSV*eta_SV(ikl,0))
900	EExcsv(ikl) = 0.
901	END DO
902
903
904	END IF
905
906
907	! Soil Water Balance
908	! =====================
909
910	! **********
911	call SISVAT_qSo
912	! #m0. (Wats_0,Wats_1,Wats_d)
913	! **********
914
915
916	! Surface Fluxes
917	! =====================
918
919	DO ikl=1,knonv
920	IRdwsv(ikl)=IRd_SV(ikl)*Eso_sv(ikl) ! Downward IR
921	! IRdwsv(ikl)=tau_sv(ikl) IRd_SV(ikl)Eso_sv(ikl) ! Downward IR
922	! . +(1.0-tau_sv(ikl))IRd_SV(ikl)Evg_sv(ikl) ! ! Etienne, remove vegetation component
923	IRupsv(ikl) = IRupsv(ikl) ! Upward IR
924	IRu_SV(ikl) = -IRupsv(ikl) ! Upward IR
925	. +IRd_SV(ikl) ! (effective)
926	. -IRdwsv(ikl) ! (positive)
927
928	TBr_sv(ikl) =sqrt(sqrt(IRu_SV(ikl)/StefBo)) ! Brightness
929	! ! Temperature
930	uts_SV(ikl) = (HSv_sv(ikl) +HSs_sv(ikl)) ! uT
931	. /(rhT_SV(ikl) *cp) !
932	uqs_SV(ikl) = (HLv_sv(ikl) +HLs_sv(ikl)) ! uq
933	. /(rhT_SV(ikl) *LhvH2O) !
934	LMO_SV(ikl) = TaT_SV(ikl)(us__SV(ikl)*3)
935	. /gravit/uts_SV(ikl)/vonKrm ! MO length
936
937
938	! Surface Temperature
939	! ^^^^^^^^^^^^^^^^^^^^
940	! Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
941
942	! Etienne: extrapolation from the two uppermost levels:
943
944	if (isnoSV(ikl) >=2) then
945	zm1=-dzsnSV(ikl,isnoSV(ikl))/2.
946	zm2=-(dzsnSV(ikl,isnoSV(ikl)) + dzsnSV(ikl,isnoSV(ikl)-1)/2.)
947	else if (isnoSV(ikl) .EQ. 1) then
948	zm1=-dzsnSV(ikl,isnoSV(ikl))/2.
949	zm2=-(dzsnSV(ikl,isnoSV(ikl))+dz_dSV(0)/2.)
950	else
951	zm1=-dz_dSV(0)/2.
952	zm2=-(dz_dSV(0)+dz_dSV(-1)/2.)
953
954	end if
955
956	coefslope=(TsisSV(ikl,isnoSV(ikl))-TsisSV(ikl,isnoSV(ikl)-1))
957	. /(zm1-zm2)
958	Tsrfsv(ikl)=TsisSV(ikl,isnoSV(ikl))+coefslope*(0. - zm1)
959
960
961	END DO
962
963
964	! Snow Pack Properties (sphericity, dendricity, size)
965	! ===================================================
966
967	IF (SnoMod) THEN
968
969	! **********
970	call SISVAT_GSn
971	! **********
972
973	! **********
974	! #ve call SISVAT_wEq('_GSn ',0)
975	! **********
976
977
978
979	END IF
980
981
982	! Roughness Length for next time step
983	!====================================
984
985	! Note that in INLANDSIS, we treat only ice covered surfaces so calculation
986	! of z0 is much simpler (no subgrid fraction of ocean or land)
987	! old calculations are commented below
988
989
990	C +--Roughness Length for Momentum
991	C + -----------------------------
992
993	C +--Land+Sea-Ice / Ice-free Sea Mask
994	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
995	DO ikl=1,klonv
996	IcIndx(ikl) = 0
997	ENDDO
998	DO isn=1,nsno
999	DO ikl=1,klonv
1000	IcIndx(ikl) = max(IcIndx(ikl),
1001	. isn*max(0,
1002	. sign(1,
1003	. int(ro__SV(ikl,isn)-900.))))
1004	ENDDO
1005	ENDDO
1006
1007	DO ikl=1,klonv
1008	LISmsk = 1. ! in inlandsis, land only
1009	IceMsk = max(0,sign(1 ,IcIndx(ikl)-1) )
1010	SnoMsk = max(min(isnoSV(ikl)-iiceSV(ikl),1),0)
1011
1012
1013
1014	Z0mLnd =max( Z0_ICE , 5.e-5 ) ! Min set := Z0 on *
1015
1016	C +--Z0 Smooth Regime over Snow (Andreas 1995, CRREL Report 95-16, p. 8)
1017	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1018	Z0m_nu = 5.e-5 ! z0s~(10-d)*exp(-vonkar/sqrt(1.1e-03))
1019
1020	C +--Z0 Saltat.Regime over Snow (Gallee et al., 2001, BLM 99 (19) p.11)
1021	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1022	u2star = us__SV(ikl) *us__SV(ikl)
1023	Z0mBSn = u2star *0.536e-3 - 61.8e-6
1024	Z0mBSn = max(Z0mBS0 ,Z0mBSn)
1025
1026	C +--Z0 Smooth + Saltat. Regime
1027	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1028	Z0enSV(ikl) = Z0m_nu
1029	. + Z0mBSn
1030
1031	C +--Rough Snow Surface Roughness Length (Typical Value)
1032	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1033	c #tz Z0m_Sn = 0.250e-3 ! Andreas 1995, CRREL Report 95-16, fig.1&p.2
1034	! z0r~(10-d)*exp(-vonkar/sqrt(1.5e-03))-5.e-5
1035	Z0m_Sn = 2.000e-3 ! Calibration of MAR
1036	c #TZ Z0m_Sn = 1.000e-3 ! Exemple Tuning in RACMO
1037	c #TZ Z0m_Sn = 0.500e-3 ! Exemple Tuning in MAR
1038
1039	C +--Rough Snow Surface Roughness Length (Variable Sastrugi Height)
1040	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1041	A_Fact = 1.0000 ! Andreas et al., 2004, p.4
1042	! ams.confex.com/ams/pdfpapers/68601.pdf
1043
1044	! Parameterization of z0 dependance on Temperature (C. Amory, 2017)
1045	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1046	! Z0=f(T) deduced from observations, Adelie Land, dec2012-dec2013
1047	coefa = 0.1658 !0.1862 !Ant
1048	coefb = -50.3869 !-55.7718 !Ant
1049	ta1 = 253.15 !255. Ant
1050	ta2 = 273.15
1051	ta3 = 273.15+3
1052	z01 = exp(coefa*ta1 + coefb) !~0.2 ! ~0.25 mm
1053	z02 = exp(coefa*ta2 + coefb) !~6 !~7 mm
1054	z03 = z01
1055	coefc = log(z03/z02)/(ta3-ta2)
1056	coefd = log(z03)-coefc*ta3
1057	if (TaT_SV(ikl) .lt. ta1) then
1058	Z0_obs = z01
1059	else if (TaT_SV(ikl).ge.ta1 .and. TaT_SV(ikl).lt.ta2) then
1060	Z0_obs = exp(coefa*TaT_SV(ikl) + coefb)
1061	else if (TaT_SV(ikl).ge.ta2 .and. TaT_SV(ikl).lt.ta3) then
1062	! if st > 0, melting induce smooth surface
1063	Z0_obs = exp(coefc*TaT_SV(ikl) + coefd)
1064	else
1065	Z0_obs = z03
1066	endif
1067
1068
1069	! pour le moment, on choisit une valeur fixe
1070	Z0_obs = 1.000e-3
1071
1072	cCA Snow roughness lenght deduced from observations
1073	cCA (parametrization if no Blowing Snow)
1074	cCA ----------------------------------- C. Agosta 09-2016 -----
1075	cCA Substract Z0enSV(ikl) because re-added later in Z0mnSV(ikl)
1076	Z0m_Sn = Z0_obs - Z0enSV(ikl)
1077	cCA -----------------------------------------------------------
1078
1079	param = Z0_obs/1. ! param(s) \| 1.(m/s)=TUNING
1080
1081	c #SZ Z0Sa_N = (us__SV(ikl) -0.2)*param ! 0.0001=TUNING
1082	c #SZ. * max(zero,sign(unun,TfSnow-eps9
1083	c #SZ. -TsisSV(ikl , isnoSV(ikl))))
1084	!!#SZ Z0SaSi = max(zero,sign(unun,Z0Sa_N ))! 1 if erosion
1085	c #SZ Z0SaSi = max(zero,sign(unun,zero -eps9 -uss_SV(ikl)))!
1086	c #SZ Z0Sa_N = max(zero, Z0Sa_N)
1087	c #SZ Z0SaSV(ikl) =
1088	c #SZ. max(Z0SaSV(ikl) ,Z0SaSV(ikl)
1089	c #SZ. + Z0SaSi(Z0Sa_N-Z0SaSV(ikl))exp(-dt__SV/43200.))
1090	c #SZ. - min(dz0_SV(ikl) , Z0SaSV(ikl))
1091
1092	c #SZ A_Fact = Z0SaSV(ikl) * 5.0/0.15 ! A=5 if h~10cm
1093	C +... CAUTION: The influence of the sastrugi direction is not yet included
1094
1095	c #SZ Z0m_Sn = Z0SaSV(ikl) !
1096	c #SZ. - Z0m_nu !
1097
1098	C +--Z0 Saltat.Regime over Snow (Shao & Lin, 1999, BLM 91 (46) p.222)
1099	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1100	c #ZN sqrrZ0 = usthSV(ikl)/max( us__SV(ikl),0.001)
1101	c #ZN sqrrZ0 = min( sqrrZ0 ,0.999)
1102	c #ZN Z0mBSn = 0.55 0.55 exp(-sqrrZ0 *sqrrZ0)
1103	c #ZN. us__SV(ikl) us__SV(ikl)grvinv0.5
1104
1105	C +--Z0 Smooth + Saltat. Regime (Shao & Lin, 1999, BLM 91 (46) p.222)
1106	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1107	c #ZN Z0enSV(ikl) = (Z0m_nu ** sqrrZ0 )
1108	c #ZN. * (Z0mBSn **(1.-sqrrZ0))
1109	c #ZN Z0enSV(ikl) = max(Z0enSV(ikl), Z0m_nu)
1110
1111	C +--Z0 Smooth Regime over Snow (Andreas etAl., 2004
1112	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ams.confex.com/ams/pdfpapers/68601.pdf)
1113	c #ZA Z0m_nu = 0.135*akmol / max(us__SV(ikl) , epsi)
1114
1115	C +--Z0 Saltat.Regime over Snow (Andreas etAl., 2004
1116	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ams.confex.com/ams/pdfpapers/68601.pdf)
1117	c #ZA Z0mBSn = 0.035u2star grvinv
1118
1119	C +--Z0 Smooth + Saltat. Regime (Andreas etAl., 2004
1120	! ( used by Erosion) ams.confex.com/ams/pdfpapers/68601.pdf)
1121	! ^^^^^^^^^^^^^^^^^^^^^^^^^^
1122	c #ZA Z0enSV(ikl) = Z0m_nu
1123	c #ZA. + Z0mBSn
1124
1125	C +--Z0 Rough Regime over Snow (Andreas etAl., 2004
1126	C + (.NOT. used by Erosion) ams.confex.com/ams/pdfpapers/68601.pdf)
1127	! ^^^^^^^^^^^^^^^^^^^^^^^^^^
1128	!!#ZA u2star = (us__SV(ikl) -0.1800) / 0.1
1129	!!#ZA Z0m_Sn =A_FactZ0mBSn exp(-u2star*u2star)
1130	c #ZA Z0m_90 =(10.-0.025*VVs_SV(ikl)/5.)
1131	c #ZA. exp(-0.4/sqrt(.00275+.00001max(0.,VVs_SV(ikl)-5.)))
1132	c #ZA Z0m_Sn = DDs_SV(ikl)* Z0m_90 / 45.
1133	c #ZA. - DDs_SV(ikl)DDs_SV(ikl) Z0m_90 /(90.*90.)
1134
1135	C +--Z0 (Erosion) over Snow (instantaneous or time average)
1136	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1137	Z0e_SV(ikl) = Z0enSV(ikl)
1138	Z0e_SV(ikl) = Z0emSV(ikl)
1139
1140	C +--Momentum Roughness Length
1141	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ! Contribution of
1142	Z0mnSV(ikl) = Z0mLnd ! land Form
1143	. + (Z0m_Sn ! Sastrugi Form
1144	. + Z0enSV(ikl)) *SnoMsk ! Snow Erosion
1145
1146
1147	C +--GIS Roughness Length
1148	C + ^^^^^^^^^^^^^^^^^^^^^
1149	c #GL Z0mnSV(ikl) =
1150	c #GL. (1-LSmask(ikl)) * Z0mnSV(ikl)
1151	c #GL. + LSmask(ikl) * max(Z0mnSV(ikl),max(Z0_GIM,
1152	c #GL. Z0_GIM+
1153	c #GL. (0.0032-Z0_GIM)*(ro__SV(ikl,isnoSV(ikl))-600.) !
1154	c #GL. /(920.00 -600.))) !
1155
1156	C +--Mom. Roughness Length, Instantaneous OR Box Moving Average in Time
1157	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1158	Z0m_SV(ikl) = Z0mnSV(ikl) ! Z0mnSV instant.
1159	! Z0m_SV(ikl) = Z0mmSV(ikl) ! Z0mnSV Average
1160
1161	C +--Corrected Threshold Friction Velocity before Erosion ! Marticorena and
1162	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! Bergametti 1995
1163	! not used anymore since Marticorena and Bergametti disabled !CK 18/07/2018
1164	cc #BS Z0e_SV(ikl) = min(Z0m_SV(ikl),Z0e_SV(ikl)) !
1165	cc #MB f_eff= log(0.35(0.1 /Z0e_SV(ikl))*0.8) ! JGR 100
1166	cc #MB f_eff=1.-(log( Z0m_SV(ikl)/Z0e_SV(ikl) ))! (20) p. 16420
1167	cc #MB. /(max( f_eff ,epsi ))! p.16426 2nd ?
1168	cc #MB f_eff= max( f_eff ,epsi )! CONTROL
1169	cc #MB f_eff=1.0 -(1.0 - f_eff) /5.00 ! TUNING
1170	cc #MB f_eff= min( f_eff ,1.00 )!
1171	cc #MB usthSV(ikl) = usthSV(ikl)/f_eff !
1172
1173
1174
1175	C +--Roughness Length for Scalars
1176	C + ----------------------------
1177
1178	Z0hnSV(ikl) = Z0mnSV(ikl)/ 7.4
1179	c #SH Z0hnSV(ikl) = Z0mnSV(ikl)/100.0
1180	C + Z0h = Z0m /100.0 over the Sahel
1181	C + (Taylor & Clark, QJRMS 127,p864)
1182
1183	c #RN rstar = Z0mnSV(ikl) * us__SV(ikl) / akmol
1184	c #RN rstar = max(epsi,min(rstar,thous))
1185	c #RN alors = log(rstar)
1186	c #RN rstar0 = 1.250e0 * max(zero,sign(unun,0.135e0 - rstar))
1187	c #RN. +(1. - max(zero,sign(unun,0.135e0 - rstar)))
1188	c #RN. (0.149e0 max(zero,sign(unun,2.500e0 - rstar))
1189	c #RN. + 0.317e0
1190	c #RN. *(1. - max(zero,sign(unun,2.500e0 - rstar))))
1191	c #RN rstar1 = 0. * max(zero,sign(unun,0.135e0 - rstar))
1192	c #RN. +(1. - max(zero,sign(unun,0.135e0 - rstar)))
1193	c #RN. (-0.55e0 max(zero,sign(unun,2.500e0 - rstar))
1194	c #RN. - 0.565
1195	c #RN. *(1. - max(zero,sign(unun,2.500e0 - rstar))))
1196	c #RN rstar2 = 0. * max(zero,sign(unun,0.135e0 - rstar))
1197	c #RN. +(1. - max(zero,sign(unun,0.135e0 - rstar)))
1198	c #RN. (0. max(zero,sign(unun,2.500e0 - rstar))
1199	c #RN. - 0.183
1200	c #RN. *(unun - max(zero,sign(unun,2.500e0 - rstar))))
1201
1202	cXF #RN does not work over bare ice
1203	cXF MAR is then too warm and not enough melt
1204
1205	c #RN if(ro__SV(ikl,isnoSV(ikl))>50
1206	c #RN. .and.ro__SV(ikl,isnoSV(ikl))<roSdSV)then
1207
1208	c #RN Z0hnSV(ikl) = max(zero
1209	c #RN. , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi))
1210	c #RN. * exp(rstar0+rstar1alors+rstar2alors*alors)
1211	c #RN. * 0.001e0 + Z0hnSV(ikl) * ( 1. - max(zero
1212	c #RN. , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi)))
1213
1214	c #RN endif
1215
1216	Z0h_SV(ikl) = Z0hnSV(ikl)
1217	! Z0h_SV(ikl) = Z0hmSV(ikl)
1218
1219
1220	c #MT Z0m_SV(ikl) = max(2.0e-6 ,Z0m_SV(ikl)) ! Min Z0_m (Garrat Scheme)
1221	! Z0m_SV(ikl) = min(Z0m_SV(ikl),za__SV(ikl)*0.3333)
1222
1223
1224	END DO
1225
1226
1227	return
1228	end
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239

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

Download in other formats:

Original Format