Context Navigation

surf_inlandsis_mod.F90 @ 5018

Last change on this file since 5018 was 3900, checked in by evignon, 3 years ago
Commission de la nouvelle interface LMDZ-SISVAT la prochaine commission consistera a supprimer l'ancien repertoire sisvat et a faire un peu de nettoyage.
File size: 62.3 KB

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

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

Download in other formats:

Original Format