Context Navigation

surf_inlandsis_mod.f90 @ 5353

Last change on this file since 5353 was 5296, checked in by abarral, 6 weeks ago
Turn compbl.h into a module Move calcul_REGDYN.h to obsolete Create phys_constants_mod.f90
File size: 62.2 KB

Line
1	MODULE surf_inlandsis_mod
2
3	IMPLICIT NONE; PRIVATE
4	PUBLIC surf_inlandsis, get_soil_levels, SISVAT_ini, sisvatetat0, sisvatredem
5
6	CONTAINS
7
8
9	SUBROUTINE surf_inlandsis(knon, rlon, rlat, ikl2i, itime, dtime, debut, lafin, &
10	rmu0, swdown, lwdown, albedo_old, pexner, ps, p1lay, &
11	precip_rain, precip_snow, &
12	zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, &
13	rugos, snow_cont_air, alb_soil, alt, slope, cloudf, &
14	radsol, qsol, tsoil, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, &
15	AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
16	runoff_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat, dflux_s,dflux_l, &
17	tsurf_new, alb1, alb2, alb3, alb6, emis_new, z0m, z0h, qsurf)
18
19	! \| \|
20	! \| SubRoutine surf_inlandsis: Interfacing Lmdz AND Sisvat's Ice and Snow\|
21	! \| (INLANDSIS) \|
22	! \| SISVAT (Soil/Ice Snow Vegetation Atmosphere Transfer Scheme) \|
23	! \| surface scheme of the Modele Atmospherique Regional (MAR) \|
24	! \| Author: Heinz Juergen Punge, LSCE June 2009 \|
25	! \| based on the MAR-SISVAT interface by Hubert Gallee \|
26	! \| Updated by Etienne Vignon, Cecile Agosta \|
27	! \| \|
28	! +------------------------------------------------------------------------+
29	! \|
30	! \| In the current setup, SISVAT is used only to model the land ice \|
31	! \| part of the surface; hence it is called with the compressed variables\|
32	! \| from pbl_surface, and only by the surf_landice routine. \|
33	! \| \|
34	! \| In this interface it is assumed that the partitioning of the soil, \|
35	! \| and hence the number of grid points is constant during a simulation, \|
36	! \| hence eg. snow properties remain stored in the global SISVAT \|
37	! \| variables between the calls and don't need to be handed over as \|
38	! \| arguments. When the partitioning is supposed to change, make sure to \|
39	! \| update the variables. \|
40	! \| \|
41	! \| INPUT (via MODULES VARxSV, VARySV, VARtSV ...) \|
42	! \| ^^^^^ xxxxSV: SISVAT/LMDZ interfacing variables \|
43	! \| \|
44	! +------------------------------------------------------------------------+
45
46	USE dimphy
47	USE VAR_SV
48	USE VARdSV
49	USE VARxSV
50	USE VARySV
51	USE VARtSV
52	USE VARphy
53	USE surface_data, only : iflag_tsurf_inlandsis, SnoMod, BloMod, ok_outfor
54	USE dimsoil_mod_h, ONLY: nsoilmx, nsnowmx, nsismx
55
56	IMPLICIT NONE
57
58	! +--Global Variables
59	! + ================
60	! Input Variables for SISVAT
61	INTEGER, INTENT(IN) :: knon
62	INTEGER, INTENT(IN) :: itime
63	REAL, INTENT(IN) :: dtime
64	LOGICAL, INTENT(IN) :: debut ! true if first step
65	LOGICAL, INTENT(IN) :: lafin ! true if last step
66
67	INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i ! Index Decompression
68	REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat
69	REAL, DIMENSION(klon), INTENT(IN) :: rmu0 ! cos sol. zenith angle
70	REAL, DIMENSION(klon), INTENT(IN) :: swdown !
71	REAL, DIMENSION(klon), INTENT(IN) :: lwdown !
72	REAL, DIMENSION(klon), INTENT(IN) :: albedo_old
73	REAL, DIMENSION(klon), INTENT(IN) :: pexner ! Exner potential
74	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
75	REAL, DIMENSION(klon), INTENT(IN) :: zsl_height, wind_velo
76	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum, ps, p1lay
77	REAL, DIMENSION(klon), INTENT(IN) :: dens_air, tsurf
78	REAL, DIMENSION(klon), INTENT(IN) :: rugos
79	REAL, DIMENSION(klon), INTENT(IN) :: snow_cont_air
80	REAL, DIMENSION(klon), INTENT(IN) :: alb_soil, slope
81	REAL, DIMENSION(klon), INTENT(IN) :: alt ! surface elevation
82	REAL, DIMENSION(klon), INTENT(IN) :: cloudf
83	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ
84	REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ
85	REAL, DIMENSION(klon), INTENT(IN) :: cdragm, cdragh
86	REAL, DIMENSION(klon), INTENT(IN) :: ustar ! friction velocity
87
88	! Variables exchanged between LMDZ and SISVAT
89	REAL, DIMENSION(klon), INTENT(IN) :: radsol ! Surface absorbed rad.
90	REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! Tot snow mass [kg/m2]
91	REAL, DIMENSION(klon), INTENT(INOUT) :: zfra ! snwo surface fraction [0-1]
92	REAL, DIMENSION(klon, nsoilmx), INTENT(OUT) :: tsoil ! Soil Temperature
93	REAL, DIMENSION(klon), INTENT(OUT) :: qsol ! Soil Water Content
94	REAL, DIMENSION(klon), INTENT(INOUT) :: z0m ! Momentum Roughn Lgt
95	REAL, DIMENSION(klon), INTENT(INOUT) :: z0h ! Momentum Roughn Lgt
96
97	! Output Variables for LMDZ
98	REAL, DIMENSION(klon), INTENT(OUT) :: alb1 ! Albedo SW
99	REAL, DIMENSION(klon), INTENT(OUT) :: alb2, alb3 ! Albedo NIR and LW
100	REAL, DIMENSION(klon,6), INTENT(OUT) :: alb6 ! 6 band Albedo
101	REAL, DIMENSION(klon), INTENT(OUT) :: emis_new ! Surface Emissivity
102	REAL, DIMENSION(klon), INTENT(OUT) :: runoff_lic ! Runoff
103	REAL, DIMENSION(klon), INTENT(OUT) :: ffonte ! enthalpy flux due to surface melting
104	REAL, DIMENSION(klon), INTENT(OUT) :: fqfonte ! water flux due to surface melting
105	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s ! d/dT sens. ht flux
106	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_l ! d/dT latent ht flux
107	REAL, DIMENSION(klon), INTENT(OUT) :: fluxsens ! Sensible ht flux
108	REAL, DIMENSION(klon), INTENT(OUT) :: fluxlat ! Latent heat flux
109	REAL, DIMENSION(klon), INTENT(OUT) :: evap ! Evaporation
110	REAL, DIMENSION(klon), INTENT(OUT) :: erod ! Erosion of surface snow (flux)
111	REAL, DIMENSION(klon), INTENT(OUT) :: agesno ! Snow age (top layer)
112	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! Surface Temperature
113	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf ! Surface Humidity
114
115	! Specific INLANDIS outputs
116	REAL, DIMENSION(klon), INTENT(OUT) :: qsnow ! Total H2O snow[kg/m2]
117	REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt ! Snow height (m)
118	REAL, DIMENSION(klon), INTENT(OUT) :: to_ice ! Snow passed to ice
119	REAL, DIMENSION(klon), INTENT(OUT) :: sissnow ! Snow in model (kg/m2)
120
121	! +--Internal Variables
122	! + ===================
123
124	CHARACTER(len = 20) :: fn_outfor ! Name for output file
125	CHARACTER (len = 80) :: abort_message
126	CHARACTER (len = 20) :: modname = 'surf_inlandsis_mod'
127
128	INTEGER :: i, ig, ikl, isl, isn, nt
129	INTEGER :: gp_outfor, un_outfor
130	REAL, PARAMETER :: f1 = 0.5
131	REAL, PARAMETER :: sn_upp = 10000., sn_low = 500.
132	REAL, PARAMETER :: sn_add = 400., sn_div = 2.
133	! snow mass upper,lower limit,
134	! added mass/division lowest layer
135	REAL, PARAMETER :: c1_zuo = 12.960e+4, c2_zuo = 2.160e+6
136	REAL, PARAMETER :: c3_zuo = 1.400e+2, czemin = 1.e-3
137	! Parameters for drainage
138	! c1_zuo/ 2.796e+4/,c2_zuo/2.160e+6/,c3_zuo/1.400e+2/ ! Tuning
139	! +... Run Off Parameters
140	! + 864001.5 day ...25 days (Modif. ETH Camp: 86400*0.3day)
141	! + (Zuo and Oerlemans 1996, J.Glacio. 42, 305--317)
142
143	REAL, DIMENSION(klon) :: eps0SL ! surface Emissivity
144	REAL :: zsigma, Ua_min, Us_min, lati
145	REAL, PARAMETER :: cdmax=0.05
146	REAL :: lambda ! Par. soil discret.
147	REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2 ! Soil layer thicknesses
148	!$OMP THREADPRIVATE(dz1,dz2)
149	LOGICAL, SAVE :: firstcall
150	!$OMP THREADPRIVATE(firstcall)
151
152	INTEGER :: iso
153	LOGICAL :: file_exists
154	CHARACTER(len = 20) :: fichnom
155	LOGICAL :: is_init_domec
156	! CA initialization
157	! dz_profil_15 : 1 m in 15 layers [m]
158	real, parameter :: dz_profil_15(15) = (/0.005, 0.01, 0.015, 0.02, 0.03, 0.04, 0.05, &
159	0.06, 0.07, 0.08, 0.09, 0.1, 0.12, 0.14, 0.17/)
160	! mean_temp : mean annual surface temperature [K]
161	real, dimension(klon) :: mean_temp
162	! mean_dens : mean surface density [kg/m3]
163	real, dimension(klon) :: mean_dens
164	! lat_scale : temperature lapse rate against latitude [K degree-1]
165	real :: lat_scale
166	! sh_scale : temperature lapse rate against altitude [K km-1]
167	real :: sh_scale
168	! variables for density profile
169	! E0, E1 : exponent
170	real :: E0, E1
171	! depth at which 550 kg m-3 is reached [m]
172	real :: z550
173	! depths of snow layers
174	real :: depth, snow_depth, distup
175	! number of initial snow layers
176	integer :: nb_snow_layer
177	! For density calc.
178	real :: alpha0, alpha1, ln_smb
179	! theoritical densities [kg m-3]
180	real :: rho0, rho1, rho1_550
181	! constants for density profile
182	! C0, C1 : constant, 0.07 for z <= 550 kg m-3
183	real, parameter :: C0 = 0.07
184	real, parameter :: C1 = 0.03
185	! rho_i : ice density [kg m-3]
186	real, parameter :: rho_ice = 917.
187	! E_c : activation energy [J mol-1]
188	real, parameter :: E_c = 60000.
189	! E_g : activation energy [J mol-1]
190	real, parameter :: E_g = 42400.
191	! R : gas constant [J mol-1 K-1]
192	real, parameter :: R = 8.3144621
193
194
195
196
197
198	! + PROGRAM START
199	! + -----------------------------------------
200
201	zsigma = 1000.
202	dt__SV = dtime
203
204	IF (debut) THEN
205	firstcall = .TRUE.
206	INI_SV = .false.
207	ELSE
208	firstcall = .false.
209	INI_SV = .true.
210	END IF
211
212	IF (ok_outfor) THEN
213	un_outfor = 51 ! unit number for point output file
214	gp_outfor = 1 ! grid point number for point output 1 for 1D, 273 for zoom-nudg DC
215	fn_outfor = 'outfor_SV.dat'
216	END IF
217
218	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
219	! + INITIALISATION: BEGIN +++
220	! + -----------------------------------------
221	IF (firstcall) THEN
222
223	! +--Array size
224	! + -----------------------
225
226	klonv = klon
227	knonv = knon
228	write(, ) 'ikl, lon and lat in INLANDSIS'
229
230	DO ikl = 1, knon
231	i=ikl2i(ikl)
232	write(, ) 'ikl=', ikl, 'rlon=', rlon(i), 'rlat=', rlat(i)
233	END DO
234
235	! +--Variables initizialisation
236	! + ---------------------------
237
238	CALL INIT_VARtSV
239	CALL INIT_VARxSV
240	CALL INIT_VARySV
241
242
243
244	! +--Surface Fall Line Slope
245	! + -----------------------
246	IF (SnoMod) THEN
247	DO ikl = 1, knon
248	slopSV(ikl) = slope(ikl)
249	SWf_SV(ikl) = & ! Normalized Decay of the
250	exp(-dt__SV & ! Surficial Water Content
251	/ (c1_zuo & !(Zuo and Oerlemans 1996,
252	+ c2_zuo * exp(-c3_zuo * abs(slopSV(ikl))))) ! J.Glacio. 42, 305--317)
253	END DO
254	END IF
255
256
257
258	! +--Soil layer thickness . Compute soil discretization (as for LMDZ)
259	! + ----------------------------------------------------------------
260	! write(*,'(/a)') 'Start SISVAT init: soil discretization ', nsoilmx
261	CALL get_soil_levels(dz1, dz2, lambda)
262
263	lambSV = lambda
264	dz1_SV(1:knon, 1:) = 0.
265	dz2_SV(1:knon, 1:) = 0.
266
267	DO isl = -nsol, 0
268	dz_dSV(isl) = 0.5e-3 * dz2(1 - isl) ! Soil layer thickness
269	DO ikl = 1, knon
270	dz1_SV(ikl, isl) = dz1(1 - isl) !1.e-3*
271	dz2_SV(ikl, isl) = dz2(1 - isl) !1.e-3*
272	END DO
273	END DO
274
275
276	! Set variables
277	! =============
278	DO ikl = 1, knon
279	! LSmask : Land/Sea Mask
280	LSmask(ikl) = 1
281	! isotSV : Soil Type -> 12 = ice
282	isotSV(ikl) = 12
283	! iWaFSV : Soil Drainage (1,0)=(y,n)
284	iWaFSV(ikl) = 1
285	! eps0SL : Surface Emissivity
286	eps0SL(ikl) = 1.
287	! alb0SV : Soil Albedo
288	alb0SV(ikl) = alb_soil(ikl)
289	! Tsf_SV : Surface Temperature, must be bellow freezing
290	Tsf_SV(ikl) = min(temp_air(ikl), TfSnow)
291	END DO
292
293	! +--Initialization of soil and snow variables in case startsis is not read
294	! + ----------------------------------------------------------------------
295
296	is_init_domec=.FALSE.
297
298
299	IF (is_init_domec) THEN
300	! Coarse initilization inspired from vertcical profiles at Dome C,
301	! Antarctic Plateaui (10m of snow, 19 levels)
302
303	DO ikl = 1,knon
304	! + Soil
305	DO isl = -nsol,0
306	TsisSV(ikl,isl) = min(tsoil(ikl,1+nsol),TfSnow-0.2) !temp_air(ikl)
307	!tsoil(ikl,1-isl) Soil Temperature
308	!TsisSV(ikl,isl) = min(temp_air(ikl),TfSnow-0.2)
309	eta_SV(ikl,isl) = epsi !etasoil(ikl,1-isl) Soil Water[m3/m3]
310	ro__SV(ikl,isl) = rhoIce !rosoil(ikl,1-isl) volumic mass
311	END DO
312
313
314	! Snow
315	isnoSV(ikl) = 19
316	istoSV(ikl, 1:isnoSV(ikl)) = 100
317	ro__SV(ikl, 1:isnoSV(ikl)) = 350.
318	eta_SV(ikl, 1:isnoSV(ikl)) = epsi
319	TsisSV(ikl, 1:isnoSV(ikl)) = min(tsoil(ikl, 1), TfSnow - 0.2)
320	G1snSV(ikl, 1:isnoSV(ikl)) = 99.
321	G2snSV(ikl, 1:isnoSV(ikl)) = 2.
322	agsnSV(ikl, 1:isnoSV(ikl)) = 50.
323	dzsnSV(ikl, 19) = 0.015
324	dzsnSV(ikl, 18) = 0.015
325	dzsnSV(ikl, 17) = 0.020
326	dzsnSV(ikl, 16) = 0.030
327	dzsnSV(ikl, 15) = 0.040
328	dzsnSV(ikl, 14) = 0.060
329	dzsnSV(ikl, 13) = 0.080
330	dzsnSV(ikl, 12) = 0.110
331	dzsnSV(ikl, 11) = 0.150
332	dzsnSV(ikl, 10) = 0.200
333	dzsnSV(ikl, 9) = 0.300
334	dzsnSV(ikl, 8) = 0.420
335	dzsnSV(ikl, 7) = 0.780
336	dzsnSV(ikl, 6) = 1.020
337	dzsnSV(ikl, 5) = 0.980
338	dzsnSV(ikl, 4) = 1.020
339	dzsnSV(ikl, 3) = 3.970
340	dzsnSV(ikl, 2) = 1.020
341	dzsnSV(ikl, 1) = 1.020
342
343	END DO
344	ELSE
345
346	! Initilialisation with climatological temperature and density
347	! profiles as in MAR. Methodology developed by Cecile Agosta
348
349	! initialize with 0., for unused snow layers
350	dzsnSV = 0.
351	G1snSV = 0.
352	G2snSV = 0.
353	istoSV = 0
354	TsisSV = 0.
355
356
357	! initialize mean variables (unrealistic)
358	mean_temp = TfSnow
359	mean_dens = 300.
360	! loop on grid cells
361	DO ikl = 1, knon
362	lati=rlat(ikl2i(ikl))
363	! approximations for mean_temp and mean_dens
364	! from Feulner et al., 2013 (DOI: 10.1175/JCLI-D-12-00636.1)
365	! Fig. 3 and 5 : the lapse rate vs. latitude at high latitude is about 0.55 °C °lat-1
366	! with a moist-adiabatic lapse rate of 5 °C km-1 everywhere except for Antarctica,
367	! for Antarctica, a dry-adiabatic lapse rate of 9.8 °C km-1 is assumed.
368	if (lati > 60.) then
369	! CA todo : add longitude bounds
370	! Greenland mean temperature : function of altitude and latitude
371	! for altitudes 0. to 1000. m, lat_scale varies from 0.9 to 0.75 °C °lat-1
372	lat_scale = (0.75 - 0.9) / 1000. * alt(ikl) + 0.9
373	lat_scale = max(min(lat_scale, 0.9), 0.75)
374	! sh_scale equals the environmental lapse rate : 6.5 °C km-1
375	sh_scale = 6.5
376	mean_temp(ikl) = TfSnow + 1.5 - sh_scale * alt(ikl) / 1000. - lat_scale * (lati - 60.)
377	! surface density: Fausto et al. 2018, https://doi.org/10.3389/feart.2018.00051
378	mean_dens(ikl) = 315.
379	else if (lati < -60.) then
380	! Antarctica mean temperature : function of altitude and latitude
381	! for altitudes 0. to 500. m, lat_scale varies from 1.3 to 0.6 °C °lat-1
382	lat_scale = (0.6 - 1.3) / 500. * alt(ikl) + 1.3
383	lat_scale = max(min(lat_scale, 1.3), 0.6)
384	! for altitudes 0. to 500. m, sh_scale varies from 6.5 to 9.8 °C km-1
385	sh_scale = (9.8 - 6.5) / 500. * alt(ikl) + 6.5
386	sh_scale = max(min(sh_scale, 9.8), 6.5)
387	mean_temp(ikl) = TfSnow - 7. - sh_scale * alt(ikl) / 1000. + lat_scale * (lati + 60.)
388	! Antarctica surface density : function of mean annual temperature
389	! surface density of 350. kg m-3 at Dome C and 450. kg m-3 at Prud'homme (Agosta et al. 2013)
390	! 350 kg m-3 is a typical value for the Antarctic plateau around 3200 m.
391	! Weinhart et al 2020 https://doi.org/10.5194/tc-14-3663-2020 and Sugiyama et al. 2011 oi: 10.3189/2012JoG11J201
392	! 320 kg m-3 is reached at Dome A, 4100 m a.s.l.
393	! Dome C : st_ant_param(3233, -75.1) = -47.7
394	! Dumont d'Urville : st_ant_param(0, -66.66) = -15.7
395	mean_dens(ikl) = (450. - 320.) / (-15.7 + 47.7) * (mean_temp(ikl) - TfSnow + 15.7) + 450.
396	mean_dens(ikl) = min(450., max(320., mean_dens(ikl)))
397	else
398
399	! write(, ) 'Attention: temperature initialization is only defined for Greenland and Antarctica'
400
401	mean_dens(ikl) =350.
402	mean_temp(ikl) = min(tsoil(ikl,1),TfSnow-0.2)
403
404	!abort_message='temperature initialization is only defined for Greenland and Antarctica'
405	!CALL abort_physic(modname,abort_message,1)
406
407	end if
408
409	! mean_temp is defined for ice ground only
410	mean_temp(ikl) = min(mean_temp(ikl), TfSnow - 0.2)
411
412	! Soil layers
413	! ===========
414	DO isl = -nsol, 0
415	! TsisSV : Temperature [K]
416	TsisSV(ikl, isl) = mean_temp(ikl)
417	! eta_SV : Soil Water [m3/m3]
418	eta_SV(ikl, isl) = epsi
419	! ro__SV : Volumic Mass [kg/m3]
420	ro__SV(ikl, isl) = rhoIce
421	END DO
422
423	! Snow layers
424	! ===========
425	! snow_depth : initial snow depth
426	snow_depth = 20.
427	! nb_snow_layer : initial nb of snow layers
428	nb_snow_layer = 15
429	! isnoSV : total nb of snow layers
430	isnoSV(ikl) = nb_snow_layer
431	! depth : depth of each layer
432	depth = snow_depth
433	do isl = 1, nb_snow_layer
434	! dzsnSV : snow layer thickness
435	dzsnSV(ikl, isl) = max(0.01, snow_depth * dz_profil_15(nb_snow_layer - isl + 1))
436	! G1snSV : dendricity (<0) or sphericity (>0) : 99. = sperical
437	G1snSV(ikl, isl) = 99.
438	! G2snSV : Sphericity (>0) or Size [1/10 mm] : 2. = small grain size
439	G2snSV(ikl, isl) = 3.
440	agsnSV(ikl, isl) = 0.
441	istoSV(ikl, isl) = 0
442	! eta_SV : Liquid Water Content [m3/m3]
443	eta_SV(ikl, isl) = 0.
444	! distance to surface
445	depth = depth - dzsnSV(ikl,isl) / 2.
446	distup = min(1., max(0., depth / snow_depth))
447	! TsisSV : Temperature [K], square interpolation between Tsf_SV (surface) and mean_temp (bottom)
448	TsisSV(ikl, isl) = Tsf_SV(ikl) * (1. - distup*2) + mean_temp(ikl) distup**2
449	! firn density : densification formulas from :
450	! Ligtenberg et al 2011 eq. (6) (www.the-cryosphere.net/5/809/2011/)
451	! equivalent to Arthern et al. 2010 eq. (4) "Nabarro-Herring" (doi:10.1029/2009JF001306)
452	! Integration of the steady state equation
453	! ln_smb approximated as a function of temperature
454	ln_smb = max((mean_temp(ikl) - TfSnow) * 5. / 60. + 8., 3.)
455	! alpha0, alpha1 : correction coefficient as a function of ln_SMB from Ligtenberg 2011, adjusted for alpha1
456	alpha0 = max(1.435 - 0.151 * ln_smb, 0.25)
457	alpha1 = max(2.0111 - 0.2051 * ln_smb, 0.25)
458	E0 = C0 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha0
459	E1 = C1 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha1
460	z550 = log((rho_ice/mean_dens(ikl) - 1.)/(rho_ice/550. - 1.)) / E0
461	rho0 = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice
462	rho1 = exp(E1 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E1 * depth)) * rho_ice
463	if (depth <= z550) then
464	ro__SV(ikl, isl) = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice
465	else
466	ro__SV(ikl, isl) = exp(E1 * (depth - z550)) / (rho_ice / 550. - 1 + exp(E1 * (depth - z550))) * rho_ice
467	end if
468	depth = depth - dzsnSV(ikl,isl) / 2.
469
470	end do
471
472	END DO
473
474	END IF
475
476
477	! + Numerics paramaters, SISVAT_ini
478	! + ----------------------
479	CALL SISVAT_ini(knon)
480
481
482	! +--Read restart file
483	! + =================================================
484
485	INQUIRE(FILE = "startsis.nc", EXIST = file_exists)
486	IF (file_exists) THEN
487	CALL sisvatetat0("startsis.nc", ikl2i)
488	END IF
489
490
491
492	! +--Output ascii file
493	! + =================================================
494
495	! open output file
496	IF (ok_outfor) THEN
497	open(unit = un_outfor, status = 'replace', file = fn_outfor)
498	ikl = gp_outfor ! index sur la grille land ice
499	write(un_outfor, *) fn_outfor, ikl, dt__SV, rlon(ikl2i(ikl)), rlat(ikl2i(ikl))
500	write(un_outfor, *) 'nsnow - albedo - z0m - z0h , dz [m,30], temp [K,41], rho [kg/m3,41], eta [kg/kg,41] &
501	& G1 [-,30], G2 [-,30], agesnow [d,30], history [-,30], DOP [m,30]'
502	END IF
503
504	END IF ! firstcall
505	! +
506	! + +++ INITIALISATION: END +++
507	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
508
509
510
511	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
512	! + READ FORCINGS
513	! + ------------------------
514
515	! + Update Forcings for SISVAT given by the LMDZ model.
516	! +
517	DO ikl = 1, knon
518
519	! +--Atmospheric Forcing (INPUT)
520	! + ^^^^^^^^^^^^^^^^^^^ ^^^^^
521	za__SV(ikl) = zsl_height(ikl) ! surface layer height (fisr model level) [m]
522	Ua_min = 0.2 * sqrt(za__SV(ikl)) !
523	VV__SV(ikl) = max(Ua_min, wind_velo(ikl)) ! Wind velocity [m/s]
524	TaT_SV(ikl) = temp_air(ikl) ! BL top Temperature [K]
525	ExnrSV(ikl) = pexner(ikl) ! Exner potential
526	rhT_SV(ikl) = dens_air(ikl) ! Air density
527	QaT_SV(ikl) = spechum(ikl) ! Specific humidity
528	ps__SV(ikl) = ps(ikl) ! surface pressure [Pa]
529	p1l_SV(ikl) = p1lay(ikl) ! lowest atm. layer press[Pa]
530
531	! +--Surface properties
532	! + ^^^^^^^^^^^^^^^^^^
533
534	Z0m_SV(ikl) = z0m(ikl) ! Moment.Roughn.L.
535	Z0h_SV(ikl) = z0h(ikl) ! Moment.Roughn.L.
536
537	! +--Energy Fluxes (INPUT)
538	! + ^^^^^^^^^^^^^ ^^^^^
539	coszSV(ikl) = max(czemin, rmu0(ikl)) ! cos(zenith.Dist.)
540	sol_SV(ikl) = swdown(ikl) ! downward Solar
541	IRd_SV(ikl) = lwdown(ikl) ! downward IR
542	rsolSV(ikl) = radsol(ikl) ! surface absorbed rad.
543
544	! +--Water Fluxes (INPUT)
545	! + ^^^^^^^^^^^^^ ^^^^^
546	drr_SV(ikl) = precip_rain(ikl) ! Rain fall rate [kg/m2/s]
547	dsn_SV(ikl) = precip_snow(ikl) ! Snow fall rate [kg/m2/s]
548
549	! #BS dbs_SV(ikl) = blowSN(i,j,n)
550	! dbs_SV = Maximum potential erosion amount [kg/m2]
551	! => Upper bound for eroded snow mass
552	! uss_SV(ikl) = SLussl(i,j,n) ! uqs (only for Tv in sisvatesbl.f)
553	! #BS if(dsn_SV(ikl)>eps12.and.erprev(i,j,n).gt.eps9) then
554	! #BS dsnbSV(ikl) =1.0-min(qsHY(i,j,kB) !BS neglib. at kb ~100 magl)
555	! #BS. /max(qshy(i,j,mz),eps9),unun)
556	! #BS dsnbSV(ikl) = max(dsnbSV(ikl),erprev(i,j,n)/dsn_SV(ikl))
557	! #BS dsnbSV(ikl) = max(0.,min(1.,dsnbSV(ikl)))
558	! #BS else
559	! #BS dsnbSV(ikl) = 0.
560	! #BS endif
561	! dsnbSV is the drift fraction of deposited snow updated in sisvat.f
562	! will be used for characterizing the Buffer Layer
563	! (see update of Bros_N, G1same, G2same, zroOLD, zroNEW)
564	! #BS if(n==1) qbs_HY(i,j) = dsnbSV(ikl)
565	qsnoSV(ikl) = snow_cont_air(ikl)
566
567
568
569	! +--Soil/BL (INPUT)
570	! + ^^^^^^^ ^^^^^
571	alb0SV(ikl) = alb_soil(ikl) ! Soil background Albedo
572	AcoHSV(ikl) = AcoefH(ikl)
573	BcoHSV(ikl) = BcoefH(ikl)
574	AcoQSV(ikl) = AcoefQ(ikl)
575	BcoQSV(ikl) = BcoefQ(ikl)
576	cdH_SV(ikl) = min(cdragh(ikl),cdmax)
577	cdM_SV(ikl) = min(cdragm(ikl),cdmax)
578	rcdmSV(ikl) = sqrt(cdM_SV(ikl))
579	Us_min = 0.01
580	us__SV(ikl) = max(Us_min, ustar(ikl))
581	ram_sv(ikl) = 1. / (cdM_SV(ikl) * max(VV__SV(ikl), eps6))
582	rah_sv(ikl) = 1. / (cdH_SV(ikl) * max(VV__SV(ikl), eps6))
583
584	! +--Energy Fluxes (INPUT/OUTPUT)
585	! + ^^^^^^^^^^^^^ ^^^^^^^^^^^^
586	!IF (.not.firstcall) THEN
587	Tsrfsv(ikl) = tsurf(ikl) !hj 12 03 2010
588	cld_SV(ikl) = cloudf(ikl) ! Cloudiness
589	!END IF
590
591	END DO
592
593	!
594	! + +++ READ FORCINGS: END +++
595	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
596
597
598	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
599	! +--SISVAT EXECUTION
600	! + ----------------
601
602	call INLANDSIS(SnoMod, BloMod, 1)
603
604
605
606	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
607	! + RETURN RESULTS
608	! + --------------
609	! + Return (compressed) SISVAT variables to LMDZ
610	! +
611	DO ikl = 1, knon ! use only 1:knon (actual ice sheet..)
612	dflux_s(ikl) = dSdTSV(ikl) ! Sens.H.Flux T-Der.
613	dflux_l(ikl) = dLdTSV(ikl) ! Latn.H.Flux T-Der.
614	fluxsens(ikl) = HSs_sv(ikl) ! HS
615	fluxlat(ikl) = HLs_sv(ikl) ! HL
616	evap(ikl) = -1*HLs_sv(ikl) / LHvH2O ! Evaporation
617	erod(ikl) = 0.
618
619	IF (BloMod) THEN
620	! + Blowing snow
621
622	! SLussl(i,j,n)= 0.
623	! #BS SLussl(i,j,n)= !Effective erosion
624	! #BS. (- dbs_ER(ikl))/(dtrhT_SV(ikl)) !~uqs* from previous time step
625	! #BS blowSN(i,j,n)= dt*uss_SV(ikl) !New max. pot. Erosion [kg/m2]
626	! #BS. *rhT_SV(ikl) !(further bounded in sisvat_bsn.f)
627	! #BS erprev(i,j,n) = dbs_Er(ikl)/dt__SV
628	erod(ikl) = dbs_Er(ikl) / dt__SV
629	ENDIF
630
631	! + Check snow thickness, substract if too thick, add if too thin
632
633	sissnow(ikl) = 0. !()
634	DO isn = 1, isnoSV(ikl)
635	sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn)
636	END DO
637
638	IF (sissnow(ikl) .LE. sn_low) THEN !add snow
639	IF (isnoSV(ikl).GE.1) THEN
640	dzsnSV(ikl, 1) = dzsnSV(ikl, 1) + sn_add / max(ro__SV(ikl, 1), epsi)
641	toicSV(ikl) = toicSV(ikl) - sn_add
642	ELSE
643	write(, ) 'Attention, bare ice... point ', ikl
644	isnoSV(ikl) = 1
645	istoSV(ikl, 1) = 0
646	ro__SV(ikl, 1) = 350.
647	dzsnSV(ikl, 1) = sn_add / max(ro__SV(ikl, 1), epsi) ! 1.
648	eta_SV(ikl, 1) = epsi
649	TsisSV(ikl, 1) = min(TsisSV(ikl, 0), TfSnow - 0.2)
650	G1snSV(ikl, 1) = 0.
651	G2snSV(ikl, 1) = 0.3
652	agsnSV(ikl, 1) = 10.
653	toicSV(ikl) = toicSV(ikl) - sn_add
654	END IF
655	END IF
656
657	IF (sissnow(ikl) .ge. sn_upp) THEN !thinnen snow layer below
658	dzsnSV(ikl, 1) = dzsnSV(ikl, 1) / sn_div
659	toicSV(ikl) = toicSV(ikl) + dzsnSV(ikl, 1) * ro__SV(ikl, 1) / sn_div
660	END IF
661
662	sissnow(ikl) = 0.
663	qsnow(ikl) = 0.
664	snow(ikl) = 0.
665	snowhgt(ikl) = 0.
666
667	DO isn = 1, isnoSV(ikl)
668	sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn)
669	snowhgt(ikl) = snowhgt(ikl) + dzsnSV(ikl, isn)
670	! Etienne: check calc qsnow
671	qsnow(ikl) = qsnow(ikl) + rhoWat * eta_SV(ikl, isn) * dzsnSV(ikl, isn)
672	END DO
673
674	zfra(ikl) = max(min(isnoSV(ikl) - iiceSV(ikl), 1), 0)
675	! Etienne: comment following line
676	! snow(ikl) = sissnow(ikl)+toicSV(ikl)
677	snow(ikl) = sissnow(ikl)
678
679	to_ice(ikl) = toicSV(ikl)
680	runoff_lic(ikl) = RnofSV(ikl) ! RunOFF: intensity (flux due to melting + liquid precip)
681	fqfonte(ikl)= max(0., (wem_SV(ikl)-wer_SV(ikl))/dtime) ! net melting = melting - refreezing
682	ffonte(ikl)=fqfonte(ikl)*Lf_H2O
683
684	qsol(ikl) = 0.
685	DO isl = -nsol, 0
686	tsoil(ikl, 1 - isl) = TsisSV(ikl, isl) ! Soil Temperature
687	! Etienne: check calc qsol
688	qsol(ikl) = qsol(ikl) &
689	+ eta_SV(ikl, isl) * dz_dSV(isl)
690	END DO
691	agesno(ikl) = agsnSV(ikl, isnoSV(ikl)) ! [day]
692
693	alb1(ikl) = alb1sv(ikl) ! Albedo VIS
694	! alb2(ikl) = ((So1dSV - f1) * alb1sv(ikl) &
695	! & + So2dSV * alb2sv(ikl) + So3dSV * alb3sv(ikl)) / f1
696	alb2(ikl)=alb2sv(ikl)
697	! Albedo NIR
698	alb3(ikl) = alb3sv(ikl) ! Albedo FIR
699	! 6 band Albedo
700	alb6(ikl,:)=alb6sv(ikl,:)
701
702	tsurf_new(ikl) = Tsrfsv(ikl)
703
704	qsurf(ikl) = QsT_SV(ikl)
705	emis_new(ikl) = eps0SL(ikl)
706	z0m(ikl) = Z0m_SV(ikl)
707	z0h(ikl) = Z0h_SV(ikl)
708
709
710	END DO
711
712	IF (ok_outfor) THEN
713	ikl= gp_outfor
714	write(un_outfor, *) '+++++++++++', rlon(ikl2i(ikl)), rlat(ikl2i(ikl)),alt(ikl),'+++++++++++'
715	write(un_outfor, *) isnoSV(ikl), alb_SV(ikl), Z0m_SV(ikl), Z0h_SV(ikl),HSs_sv(ikl),HLs_sv(ikl),alb1(ikl),alb2(ikl)
716	write(un_outfor, *) dzsnSV(ikl, :)
717	write(un_outfor, *) TsisSV(ikl, :)
718	write(un_outfor, *) ro__SV(ikl, :)
719	write(un_outfor, *) eta_SV(ikl, :)
720	write(un_outfor, *) G1snSV(ikl, :)
721	write(un_outfor, *) G2snSV(ikl, :)
722	write(un_outfor, *) agsnSV(ikl, :)
723	write(un_outfor, *) istoSV(ikl, :)
724	write(un_outfor, *) DOPsnSV(ikl, :)
725	ENDIF
726
727
728
729	! + -----------------------------
730	! + END --- RETURN RESULTS
731	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
732	IF (lafin) THEN
733	fichnom = "restartsis.nc"
734	CALL sisvatredem("restartsis.nc", ikl2i, rlon, rlat)
735
736	IF (ok_outfor) THEN
737	close(unit = un_outfor)
738	END IF
739	END IF
740
741	END SUBROUTINE surf_inlandsis
742
743
744	!=======================================================================
745
746	SUBROUTINE get_soil_levels(dz1, dz2, lambda)
747	! ======================================================================
748	! Routine to compute the vertical discretization of the soil in analogy
749	! to LMDZ. In LMDZ it is done in soil.F, which is not used in the case
750	! of SISVAT, therefore it's needed here.
751	!
752	USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root
753	USE mod_phys_lmdz_para
754	USE VAR_SV
755	USE dimsoil_mod_h, ONLY: nsoilmx, nsnowmx, nsismx
756
757	REAL, DIMENSION(nsoilmx), INTENT(OUT) :: dz2, dz1
758	REAL, INTENT(OUT) :: lambda
759
760
761	!-----------------------------------------------------------------------
762	! Depthts:
763	! --------
764	REAL fz, rk, fz1, rk1, rk2
765	REAL min_period, dalph_soil
766	INTEGER ierr, jk
767
768	fz(rk) = fz1 * (dalph_soil**rk - 1.) / (dalph_soil - 1.)
769
770	! write(,)'Start soil level computation'
771	!-----------------------------------------------------------------------
772	! Calculation of some constants
773	! NB! These constants do not depend on the sub-surfaces
774	!-----------------------------------------------------------------------
775	!-----------------------------------------------------------------------
776	! ground levels
777	! grnd=z/l where l is the skin depth of the diurnal cycle:
778	!-----------------------------------------------------------------------
779
780	min_period = 1800. ! en secondes
781	dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
782	! !$OMP MASTER
783	! IF (is_mpi_root) THEN
784	! OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
785	! IF (ierr == 0) THEN ! Read file only if it exists
786	! READ(99,*) min_period
787	! READ(99,*) dalph_soil
788	! PRINT*,'Discretization for the soil model'
789	! PRINT*,'First level e-folding depth',min_period, &
790	! ' dalph',dalph_soil
791	! CLOSE(99)
792	! END IF
793	! ENDIF
794	! !$OMP END MASTER
795	! CALL bcast(min_period)
796	! CALL bcast(dalph_soil)
797
798	! la premiere couche represente un dixieme de cycle diurne
799	fz1 = SQRT(min_period / 3.14)
800
801	DO jk = 1, nsoilmx
802	rk1 = jk
803	rk2 = jk - 1
804	dz2(jk) = fz(rk1) - fz(rk2)
805	ENDDO
806	DO jk = 1, nsoilmx - 1
807	rk1 = jk + .5
808	rk2 = jk - .5
809	dz1(jk) = 1. / (fz(rk1) - fz(rk2))
810	ENDDO
811	lambda = fz(.5) * dz1(1)
812	DO jk = 1, nsoilmx
813	rk = jk
814	rk1 = jk + .5
815	rk2 = jk - .5
816	ENDDO
817
818	END SUBROUTINE get_soil_levels
819
820
821	!===========================================================================
822
823	SUBROUTINE SISVAT_ini(knon)
824
825	!C +------------------------------------------------------------------------+
826	!C \| MAR SISVAT_ini Jd 11-10-2007 MAR \|
827	!C \| SubRoutine SISVAT_ini generates non time dependant SISVAT parameters \|
828	!C +------------------------------------------------------------------------+
829	!C \| PARAMETERS: klonv: Total Number of columns = \|
830	!C \| ^^^^^^^^^^ = Total Number of continental grid boxes \|
831	!C \| X Number of Mosaic Cell per grid box \|
832	!C \| \|
833	!C \| INPUT: dt__SV : Time Step [s] \|
834	!C \| ^^^^^ dz_dSV : Layer Thickness [m] \|
835	!C \| \|
836	!C \| OUTPUT: [-] \|
837	!C \| ^^^^^^ rocsSV : Soil Contrib. to (ro c)_s exclud.Water [J/kg/K] \|
838	!C \| etamSV : Soil Minimum Humidity [m3/m3] \|
839	!C \| (based on a prescribed Soil Relative Humidity) \|
840	!C \| s1__SV : Factor of eta**( b+2) in Hydraul.Diffusiv. \|
841	!C \| s2__SV : Factor of eta**( b+2) in Hydraul.Conduct. \|
842	!C \| aKdtSV : KHyd: Piecewise Linear Profile: a * dt [m] \|
843	!C \| bKdtSV : KHyd: Piecewise Linear Profile: b * dt [m/s] \|
844	!C \| dzsnSV(0): Soil first Layer Thickness [m] \|
845	!C \| dzmiSV : Distance between two contiguous levels [m] \|
846	!C \| dz78SV : 7/8 (Layer Thickness) [m] \|
847	!C \| dz34SV : 3/4 (Layer Thickness) [m] \|
848	!C \| dz_8SV : 1/8 (Layer Thickness) [m] \|
849	!C \| dzAvSV : 1/8 dz_(i-1) + 3/4 dz_(i) + 1/8 dz_(i+1) [m] \|
850	!C \| dtz_SV : dt/dz [s/m] \|
851	!C \| OcndSV : Swab Ocean / Soil Ratio [-] \|
852	!C \| Implic : Implicit Parameter (0.5: Crank-Nicholson) \|
853	!C \| Explic : Explicit Parameter = 1.0 - Implic \|
854	!C \| \|
855	!C \| # OPTIONS: #ER: Richards Equation is not smoothed \|
856	!C \| # ^^^^^^^ #kd: De Ridder Discretization \|
857	!C \| # #SH: Hapex-Sahel Values !
858	!C \| \|
859	!C +------------------------------------------------------------------------+
860	!
861	!
862
863	!C +--Global Variables
864	!C + ================
865
866	USE dimphy
867	USE VARphy
868	USE VAR_SV
869	USE VARdSV
870	USE VAR0SV
871	USE VARxSV
872	USE VARtSV
873	USE VARxSV
874	USE VARySV
875	IMPLICIT NONE
876
877
878
879	!C +--Arguments
880	!C + ==================
881	INTEGER, INTENT(IN) :: knon
882
883	!C +--Internal Variables
884	!C + ==================
885
886	INTEGER :: ivt, ist, ikl, isl, isn, ikh
887	INTEGER :: misl_2, nisl_2
888	REAL :: d__eta, eta__1, eta__2, Khyd_1, Khyd_2
889	REAL, PARAMETER :: RHsMin = 0.001 ! Min.Soil Relative Humidity
890	REAL :: PsiMax ! Max.Soil Water Potential
891	REAL :: a_Khyd, b_Khyd ! Water conductivity
892
893
894	!c #WR REAL :: Khyd_x,Khyd_y
895
896
897
898	!C +--Non Time Dependant SISVAT parameters
899	!C + ====================================
900
901	!C +--Soil Discretization
902	!C + -------------------
903
904	!C +--Numerical Scheme Parameters
905	!C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^
906	Implic = 0.75 ! 0.5 <==> Crank-Nicholson
907	Explic = 1.00 - Implic !
908
909	!C +--Soil/Snow Layers Indices
910	!C + ^^^^^^^^^^^^^^^^^^^^^^^^
911	DO isl = -nsol, 0
912	islpSV(isl) = isl + 1
913	islpSV(isl) = min(islpSV(isl), 0)
914	islmSV(isl) = isl - 1
915	islmSV(isl) = max(-nsol, islmSV(isl))
916	END DO
917
918	DO isn = 1, nsno
919	isnpSV(isn) = isn + 1
920	isnpSV(isn) = min(isnpSV(isn), nsno)
921	END DO
922
923	!C +--Soil Layers Thicknesses
924	!C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
925	! Not used here as LMDZ method is applied, see SUBROUTINE get_soil_levels!
926	!c #kd IF (nsol.gt.4) THEN
927	!c #kd DO isl=-5,-nsol,-1
928	!c #kd dz_dSV(isl)= 1.
929	!c #kd END DO
930	!c #kd END IF
931	!
932	! IF (nsol.ne.4) THEN
933	! DO isl= 0,-nsol,-1
934	! misl_2 = -mod(isl,2)
935	! nisl_2 = -isl/2
936	! dz_dSV(isl)=(((1-misl_2) * 0.001
937	! . + misl_2 * 0.003) * 10*(nisl_2)) 4.
938	!C +... dz_dSV(0) = Hapex-Sahel Calibration: 4 mm
939	!
940	!c +SH dz_dSV(isl)=(((1-misl_2) * 0.001
941	!c +SH. + misl_2 * 0.003) * 10*(nisl_2)) 1.
942	!
943	!c #05 dz_dSV(isl)=(((1-misl_2) * 0.001
944	!c #05. + misl_2 * 0.008) * 10*(nisl_2)) 0.5
945	! END DO
946	! dz_dSV(0) = 0.001
947	! dz_dSV(-1) = dz_dSV(-1) - dz_dSV(0) + 0.004
948	! END IF
949
950	zz_dSV = 0.
951	DO isl = -nsol, 0
952	dzmiSV(isl) = 0.500 * (dz_dSV(isl) + dz_dSV(islmSV(isl)))
953	dziiSV(isl) = 0.500 * dz_dSV(isl) / dzmiSV(isl)
954	dzi_SV(isl) = 0.500 * dz_dSV(islmSV(isl)) / dzmiSV(isl)
955	dtz_SV(isl) = dt__SV / dz_dSV(isl)
956	dtz_SV2(isl) = 1. / dz_dSV(isl)
957	dz78SV(isl) = 0.875 * dz_dSV(isl)
958	dz34SV(isl) = 0.750 * dz_dSV(isl)
959	dz_8SV(isl) = 0.125 * dz_dSV(isl)
960	dzAvSV(isl) = 0.125 * dz_dSV(islmSV(isl)) &
961	& + 0.750 * dz_dSV(isl) &
962	& + 0.125 * dz_dSV(islpSV(isl))
963	zz_dSV = zz_dSV + dz_dSV(isl)
964	END DO
965	DO ikl = 1, knon !v
966	dzsnSV(ikl, 0) = dz_dSV(0)
967	END DO
968
969	!C +--Conversion to a 50 m Swab Ocean Discretization
970	!C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
971	OcndSV = 0.
972	DO isl = -nsol, 0
973	OcndSV = OcndSV + dz_dSV(isl)
974	END DO
975	OcndSV = 50. / OcndSV
976
977
978	!C +--Secondary Soil Parameters
979	!C + -------------------------------
980
981	DO ist = 0, nsot
982	rocsSV(ist) = (1.0 - etadSV(ist)) * 1.2E+6 ! Soil Contrib. to (ro c)_s
983	s1__SV(ist) = bCHdSV(ist) & ! Factor of (eta)**(b+2)
984	& * psidSV(ist) * Ks_dSV(ist) & ! in DR97, Eqn.(3.36)
985	& / (etadSV(ist)**(bCHdSV(ist) + 3.)) !
986	s2__SV(ist) = Ks_dSV(ist) & ! Factor of (eta)**(2b+3)
987	& / (etadSV(ist)*(2. bCHdSV(ist) + 3.)) ! in DR97, Eqn.(3.35)
988
989	!C +--Soil Minimum Humidity (from a prescribed minimum relative Humidity)
990	!C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
991	Psimax = -(log(RHsMin)) / 7.2E-5 ! DR97, Eqn 3.15 Inversion
992	etamSV(ist) = etadSV(ist) &
993	& * (PsiMax / psidSV(ist))**(-min(10., 1. / bCHdSV(ist)))
994	END DO
995	etamSV(12) = 0.
996
997	!C +--Piecewise Hydraulic Conductivity Profiles
998	!C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
999	DO ist = 0, nsot
1000
1001	d__eta = etadSV(ist) / nkhy
1002	eta__1 = 0.
1003	eta__2 = d__eta
1004	DO ikh = 0, nkhy
1005	Khyd_1 = s2__SV(ist) & ! DR97, Eqn.(3.35)
1006	& * (eta__1 *(2. bCHdSV(ist) + 3.)) !
1007	Khyd_2 = s2__SV(ist) &!
1008	& * (eta__2 *(2. bCHdSV(ist) + 3.)) !
1009
1010	a_Khyd = (Khyd_2 - Khyd_1) / d__eta !
1011	b_Khyd = Khyd_1 - a_Khyd * eta__1 !
1012	!c #WR Khyd_x = a_Khyd*eta__1 +b_Khyd !
1013	!c #WR Khyd_y = a_Khyd*eta__2 +b_Khyd !
1014	aKdtSV(ist, ikh) = a_Khyd * dt__SV !
1015	bKdtSV(ist, ikh) = b_Khyd * dt__SV !
1016
1017	eta__1 = eta__1 + d__eta
1018	eta__2 = eta__2 + d__eta
1019	END DO
1020	END DO
1021
1022	return
1023
1024	END SUBROUTINE SISVAT_ini
1025
1026
1027	!***************************************************************************
1028
1029	SUBROUTINE sisvatetat0 (fichnom, ikl2i)
1030	USE clesphys_mod_h
1031	USE dimphy
1032	USE mod_grid_phy_lmdz
1033	USE mod_phys_lmdz_para
1034	USE iostart
1035	USE VAR_SV
1036	USE VARdSV
1037	USE VARxSV
1038	USE VARtSV
1039	USE indice_sol_mod
1040	USE dimsoil_mod_h, ONLY: nsoilmx, nsnowmx, nsismx
1041	USE clesphys_mod_h
1042	USE compbl_mod_h
1043	IMPLICIT none
1044	!======================================================================
1045	! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
1046	! Objet: Lecture du fichier de conditions initiales pour SISVAT
1047	!======================================================================
1048
1049	CHARACTER(LEN = *) :: fichnom
1050
1051	INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
1052	REAL, DIMENSION(klon) :: rlon
1053	REAL, DIMENSION(klon) :: rlat
1054
1055	! les variables globales ecrites dans le fichier restart
1056	REAL, DIMENSION(klon) :: isno
1057	REAL, DIMENSION(klon) :: ispi
1058	REAL, DIMENSION(klon) :: iice
1059	REAL, DIMENSION(klon) :: rusn
1060	REAL, DIMENSION(klon, nsno) :: isto
1061
1062	REAL, DIMENSION(klon, nsismx) :: Tsis
1063	REAL, DIMENSION(klon, nsismx) :: eta
1064	REAL, DIMENSION(klon, nsismx) :: ro
1065
1066	REAL, DIMENSION(klon, nsno) :: dzsn
1067	REAL, DIMENSION(klon, nsno) :: G1sn
1068	REAL, DIMENSION(klon, nsno) :: G2sn
1069	REAL, DIMENSION(klon, nsno) :: agsn
1070
1071	REAL, DIMENSION(klon) :: toic
1072
1073	INTEGER :: isl, ikl, i, isn, errT, erreta, errro, errdz, snopts
1074	CHARACTER (len = 2) :: str2
1075	LOGICAL :: found
1076
1077	errT = 0
1078	errro = 0
1079	erreta = 0
1080	errdz = 0
1081	snopts = 0
1082	! Ouvrir le fichier contenant l'etat initial:
1083
1084	CALL open_startphy(fichnom)
1085
1086	! Lecture des latitudes, longitudes (coordonnees):
1087
1088	CALL get_field("latitude", rlat, found)
1089	CALL get_field("longitude", rlon, found)
1090
1091	CALL get_field("n_snows", isno, found)
1092	IF (.NOT. found) THEN
1093	PRINT*, 'phyetat0: Le champ <n_snows> est absent'
1094	PRINT *, 'fichier startsisvat non compatible avec sisvatetat0'
1095	ENDIF
1096
1097	CALL get_field("n_ice_top", ispi, found)
1098	CALL get_field("n_ice", iice, found)
1099	CALL get_field("surf_water", rusn, found)
1100
1101
1102	CALL get_field("to_ice", toic, found)
1103	IF (.NOT. found) THEN
1104	PRINT*, 'phyetat0: Le champ <to_ice> est absent'
1105	toic(:) = 0.
1106	ENDIF
1107
1108	DO isn = 1, nsno
1109	IF (isn.LE.99) THEN
1110	WRITE(str2, '(i2.2)') isn
1111	CALL get_field("AGESNOW" // str2, &
1112	agsn(:, isn), found)
1113	ELSE
1114	PRINT*, "Trop de couches"
1115	CALL abort
1116	ENDIF
1117	ENDDO
1118	DO isn = 1, nsno
1119	IF (isn.LE.99) THEN
1120	WRITE(str2, '(i2.2)') isn
1121	CALL get_field("DZSNOW" // str2, &
1122	dzsn(:, isn), found)
1123	ELSE
1124	PRINT*, "Trop de couches"
1125	CALL abort
1126	ENDIF
1127	ENDDO
1128	DO isn = 1, nsno
1129	IF (isn.LE.99) THEN
1130	WRITE(str2, '(i2.2)') isn
1131	CALL get_field("G2SNOW" // str2, &
1132	G2sn(:, isn), found)
1133	ELSE
1134	PRINT*, "Trop de couches"
1135	CALL abort
1136	ENDIF
1137	ENDDO
1138	DO isn = 1, nsno
1139	IF (isn.LE.99) THEN
1140	WRITE(str2, '(i2.2)') isn
1141	CALL get_field("G1SNOW" // str2, &
1142	G1sn(:, isn), found)
1143	ELSE
1144	PRINT*, "Trop de couches"
1145	CALL abort
1146	ENDIF
1147	ENDDO
1148	DO isn = 1, nsismx
1149	IF (isn.LE.99) THEN
1150	WRITE(str2, '(i2.2)') isn
1151	CALL get_field("ETA" // str2, &
1152	eta(:, isn), found)
1153	ELSE
1154	PRINT*, "Trop de couches"
1155	CALL abort
1156	ENDIF
1157	ENDDO
1158	DO isn = 1, nsismx
1159	IF (isn.LE.99) THEN
1160	WRITE(str2, '(i2.2)') isn
1161	CALL get_field("RO" // str2, &
1162	ro(:, isn), found)
1163	ELSE
1164	PRINT*, "Trop de couches"
1165	CALL abort
1166	ENDIF
1167	ENDDO
1168	DO isn = 1, nsismx
1169	IF (isn.LE.99) THEN
1170	WRITE(str2, '(i2.2)') isn
1171	CALL get_field("TSS" // str2, &
1172	Tsis(:, isn), found)
1173	ELSE
1174	PRINT*, "Trop de couches"
1175	CALL abort
1176	ENDIF
1177	ENDDO
1178	DO isn = 1, nsno
1179	IF (isn.LE.99) THEN
1180	WRITE(str2, '(i2.2)') isn
1181	CALL get_field("HISTORY" // str2, &
1182	isto(:, isn), found)
1183	ELSE
1184	PRINT*, "Trop de couches"
1185	CALL abort
1186	ENDIF
1187	ENDDO
1188	write(, )'Read ', fichnom, ' finished!!'
1189
1190	!*********************************************************************************
1191	! Compress restart file variables for SISVAT
1192
1193	DO ikl = 1, klon
1194	i = ikl2i(ikl)
1195	IF (i > 0) THEN
1196	isnoSV(ikl) = INT(isno(i)) ! Nb Snow/Ice Lay.
1197	ispiSV(ikl) = INT(ispi(i)) ! Nb Supr.Ice Lay.
1198	iiceSV(ikl) = INT(iice(i)) ! Nb Ice Lay.
1199
1200	DO isl = -nsol, 0
1201	ro__SV(ikl, isl) = ro(i, nsno + 1 - isl) !
1202	eta_SV(ikl, isl) = eta(i, nsno + 1 - isl) ! Soil Humidity
1203	!hjp 15/10/2010
1204	IF (eta_SV(ikl, isl) <= 1.e-6) THEN !hj check
1205	eta_SV(ikl, isl) = 1.e-6
1206	ENDIF
1207	TsisSV(ikl, isl) = Tsis(i, nsno + 1 - isl) ! Soil Temperature
1208	IF (TsisSV(ikl, isl) <= 1.) THEN !hj check
1209	! errT=errT+1
1210	TsisSV(ikl, isl) = 273.15 - 0.2 ! Etienne: negative temperature since soil is ice
1211	ENDIF
1212
1213	END DO
1214	write(, )'Copy histo', ikl
1215
1216	DO isn = 1, isnoSV(ikl) !nsno
1217	snopts = snopts + 1
1218	IF (isto(i, isn) > 10.) THEN !hj check
1219	write(, )'Irregular isto', ikl, i, isn, isto(i, isn)
1220	isto(i, isn) = 1.
1221	ENDIF
1222
1223	istoSV(ikl, isn) = INT(isto(i, isn)) ! Snow History
1224	ro__SV(ikl, isn) = ro(i, isn) ! [kg/m3]
1225	eta_SV(ikl, isn) = eta(i, isn) ! [m3/m3]
1226	TsisSV(ikl, isn) = Tsis(i, isn) ! [K]
1227
1228	IF (TsisSV(ikl, isn) <= 1.) THEN !hj check
1229	errT = errT + 1
1230	TsisSV(ikl, isn) = TsisSV(ikl, 0)
1231	ENDIF
1232	IF (TsisSV(ikl, isn) <= 1.) THEN !hj check
1233	TsisSV(ikl, isn) = 263.15
1234	ENDIF
1235	IF (eta_SV(ikl, isn) < 1.e-9) THEN !hj check
1236	eta_SV(ikl, isn) = 1.e-6
1237	erreta = erreta + 1
1238	ENDIF
1239	IF (ro__SV(ikl, isn) <= 10.) THEN !hj check
1240	ro__SV(ikl, isn) = 11.
1241	errro = errro + 1
1242	ENDIF
1243	write(, )ikl, i, isn, Tsis(i, isn), G1sn(i, isn)
1244	G1snSV(ikl, isn) = G1sn(i, isn) ! [-] [-]
1245	G2snSV(ikl, isn) = G2sn(i, isn) ! [-] [0.0001 m]
1246	dzsnSV(ikl, isn) = dzsn(i, isn) ! [m]
1247	agsnSV(ikl, isn) = agsn(i, isn) ! [day]
1248	END DO
1249	rusnSV(ikl) = rusn(i) ! Surficial Water
1250	toicSV(ikl) = toic(i) ! bilan snow to ice
1251	END IF
1252	END DO
1253
1254	END SUBROUTINE sisvatetat0
1255
1256
1257	!======================================================================
1258	SUBROUTINE sisvatredem (fichnom, ikl2i, rlon, rlat)
1259
1260
1261
1262	!======================================================================
1263	! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
1264	! Objet: Ecriture de l'etat de redemarrage pour SISVAT
1265	!======================================================================
1266	USE compbl_mod_h
1267	USE mod_grid_phy_lmdz
1268	USE mod_phys_lmdz_para
1269	USE iostart
1270	USE VAR_SV
1271	USE VARxSV
1272	USE VARySV !hj tmp 12 03 2010
1273	USE VARtSV
1274	USE indice_sol_mod
1275	USE dimphy
1276	USE dimsoil_mod_h, ONLY: nsoilmx, nsnowmx, nsismx
1277
1278	IMPLICIT none
1279
1280	!======================================================================
1281
1282	CHARACTER(LEN = *) :: fichnom
1283	INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
1284	REAL, DIMENSION(klon), INTENT(IN) :: rlon
1285	REAL, DIMENSION(klon), INTENT(IN) :: rlat
1286
1287	! les variables globales ecrites dans le fichier restart
1288	REAL, DIMENSION(klon) :: isno
1289	REAL, DIMENSION(klon) :: ispi
1290	REAL, DIMENSION(klon) :: iice
1291	REAL, DIMENSION(klon, nsnowmx) :: isto
1292
1293	REAL, DIMENSION(klon, nsismx) :: Tsis
1294	REAL, DIMENSION(klon, nsismx) :: eta
1295	REAL, DIMENSION(klon, nsnowmx) :: dzsn
1296	REAL, DIMENSION(klon, nsismx) :: ro
1297	REAL, DIMENSION(klon, nsnowmx) :: G1sn
1298	REAL, DIMENSION(klon, nsnowmx) :: G2sn
1299	REAL, DIMENSION(klon, nsnowmx) :: agsn
1300	REAL, DIMENSION(klon) :: IRs
1301	REAL, DIMENSION(klon) :: LMO
1302	REAL, DIMENSION(klon) :: rusn
1303	REAL, DIMENSION(klon) :: toic
1304	REAL, DIMENSION(klon) :: Bufs
1305	REAL, DIMENSION(klon) :: alb1, alb2, alb3
1306
1307	INTEGER isl, ikl, i, isn, ierr
1308	CHARACTER (len = 2) :: str2
1309	INTEGER :: pass
1310
1311	isno(:) = 0
1312	ispi(:) = 0
1313	iice(:) = 0
1314	IRs(:) = 0.
1315	LMO(:) = 0.
1316	eta(:, :) = 0.
1317	Tsis(:, :) = 0.
1318	isto(:, :) = 0
1319	ro(:, :) = 0.
1320	G1sn(:, :) = 0.
1321	G2sn(:, :) = 0.
1322	dzsn(:, :) = 0.
1323	agsn(:, :) = 0.
1324	rusn(:) = 0.
1325	toic(:) = 0.
1326	Bufs(:) = 0.
1327	alb1(:) = 0.
1328	alb2(:) = 0.
1329	alb3(:) = 0.
1330
1331	!***************************************************************************
1332	! Uncompress SISVAT output variables for storage
1333
1334	DO ikl = 1, klon
1335	i = ikl2i(ikl)
1336	IF (i > 0) THEN
1337	isno(i) = 1. * isnoSV(ikl) ! Nb Snow/Ice Lay.
1338	ispi(i) = 1. * ispiSV(ikl) ! Nb Supr.Ice Lay.
1339	iice(i) = 1. * iiceSV(ikl) ! Nb Ice Lay.
1340
1341	! IRs(i) = IRs_SV(ikl)
1342	! LMO(i) = LMO_SV(ikl)
1343
1344	DO isl = -nsol, 0 !
1345	eta(i, nsno + 1 - isl) = eta_SV(ikl, isl) ! Soil Humidity
1346	Tsis(i, nsno + 1 - isl) = TsisSV(ikl, isl) ! Soil Temperature
1347	ro(i, nsno + 1 - isl) = ro__SV(ikl, isl) ! [kg/m3]
1348	END DO
1349
1350	DO isn = 1, nsno
1351	isto(i, isn) = 1. * istoSV(ikl, isn) ! Snow History
1352	ro(i, isn) = ro__SV(ikl, isn) ! [kg/m3]
1353	eta(i, isn) = eta_SV(ikl, isn) ! [m3/m3]
1354	Tsis(i, isn) = TsisSV(ikl, isn) ! [K]
1355	G1sn(i, isn) = G1snSV(ikl, isn) ! [-] [-]
1356	G2sn(i, isn) = G2snSV(ikl, isn) ! [-] [0.0001 m]
1357	dzsn(i, isn) = dzsnSV(ikl, isn) ! [m]
1358	agsn(i, isn) = agsnSV(ikl, isn) ! [day]
1359	END DO
1360	rusn(i) = rusnSV(ikl) ! Surficial Water
1361	toic(i) = toicSV(ikl) ! to ice
1362	alb1(i) = alb1sv(ikl)
1363	alb2(i) = alb2sv(ikl)
1364	alb3(i) = alb3sv(ikl)
1365	! Bufs(i) = BufsSV(ikl)
1366	END IF
1367	END DO
1368
1369	CALL open_restartphy(fichnom)
1370
1371	DO pass = 1, 2
1372	CALL put_field(pass, "longitude", &
1373	"Longitudes de la grille physique", rlon)
1374	CALL put_field(pass, "latitude", "Latitudes de la grille physique", rlat)
1375
1376	CALL put_field(pass, "n_snows", "number of snow/ice layers", isno)
1377	CALL put_field(pass, "n_ice_top", "number of top ice layers", ispi)
1378	CALL put_field(pass, "n_ice", "number of ice layers", iice)
1379	CALL put_field(pass, "IR_soil", "Soil IR flux", IRs)
1380	CALL put_field(pass, "LMO", "Monin-Obukhov Scale", LMO)
1381	CALL put_field(pass, "surf_water", "Surficial water", rusn)
1382	CALL put_field(pass, "snow_buffer", "Snow buffer layer", Bufs)
1383	CALL put_field(pass, "alb_1", "albedo sw", alb1)
1384	CALL put_field(pass, "alb_2", "albedo nIR", alb2)
1385	CALL put_field(pass, "alb_3", "albedo fIR", alb3)
1386	CALL put_field(pass, "to_ice", "Snow passed to ice", toic)
1387
1388	DO isn = 1, nsno
1389	IF (isn.LE.99) THEN
1390	WRITE(str2, '(i2.2)') isn
1391	CALL put_field(pass, "AGESNOW" // str2, &
1392	"Age de la neige layer No." // str2, &
1393	agsn(:, isn))
1394	ELSE
1395	PRINT*, "Trop de couches"
1396	CALL abort
1397	ENDIF
1398	ENDDO
1399	DO isn = 1, nsno
1400	IF (isn.LE.99) THEN
1401	WRITE(str2, '(i2.2)') isn
1402	CALL put_field(pass, "DZSNOW" // str2, &
1403	"Snow/ice thickness layer No." // str2, &
1404	dzsn(:, isn))
1405	ELSE
1406	PRINT*, "Trop de couches"
1407	CALL abort
1408	ENDIF
1409	ENDDO
1410	DO isn = 1, nsno
1411	IF (isn.LE.99) THEN
1412	WRITE(str2, '(i2.2)') isn
1413	CALL put_field(pass, "G2SNOW" // str2, &
1414	"Snow Property 2, layer No." // str2, &
1415	G2sn(:, isn))
1416	ELSE
1417	PRINT*, "Trop de couches"
1418	CALL abort
1419	ENDIF
1420	ENDDO
1421	DO isn = 1, nsno
1422	IF (isn.LE.99) THEN
1423	WRITE(str2, '(i2.2)') isn
1424	CALL put_field(pass, "G1SNOW" // str2, &
1425	"Snow Property 1, layer No." // str2, &
1426	G1sn(:, isn))
1427	ELSE
1428	PRINT*, "Trop de couches"
1429	CALL abort
1430	ENDIF
1431	ENDDO
1432	DO isn = 1, nsismx
1433	IF (isn.LE.99) THEN
1434	WRITE(str2, '(i2.2)') isn
1435	CALL put_field(pass, "ETA" // str2, &
1436	"Soil/snow water content layer No." // str2, &
1437	eta(:, isn))
1438	ELSE
1439	PRINT*, "Trop de couches"
1440	CALL abort
1441	ENDIF
1442	ENDDO
1443	DO isn = 1, nsismx !nsno
1444	IF (isn.LE.99) THEN
1445	WRITE(str2, '(i2.2)') isn
1446	CALL put_field(pass, "RO" // str2, &
1447	"Snow density layer No." // str2, &
1448	ro(:, isn))
1449	ELSE
1450	PRINT*, "Trop de couches"
1451	CALL abort
1452	ENDIF
1453	ENDDO
1454	DO isn = 1, nsismx
1455	IF (isn.LE.99) THEN
1456	WRITE(str2, '(i2.2)') isn
1457	CALL put_field(pass, "TSS" // str2, &
1458	"Soil/snow temperature layer No." // str2, &
1459	Tsis(:, isn))
1460	ELSE
1461	PRINT*, "Trop de couches"
1462	CALL abort
1463	ENDIF
1464	ENDDO
1465	DO isn = 1, nsno
1466	IF (isn.LE.99) THEN
1467	WRITE(str2, '(i2.2)') isn
1468	CALL put_field(pass, "HISTORY" // str2, &
1469	"Snow history layer No." // str2, &
1470	isto(:, isn))
1471	ELSE
1472	PRINT*, "Trop de couches"
1473	CALL abort
1474	ENDIF
1475	ENDDO
1476
1477	CALL enddef_restartphy
1478	ENDDO
1479	CALL close_restartphy
1480
1481	END SUBROUTINE sisvatredem
1482
1483	END MODULE surf_inlandsis_mod

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

Download in other formats:

Original Format