Context Navigation

surf_inlandsis_mod.F90 @ 5135

Last change on this file since 5135 was 5134, checked in by abarral, 16 months ago
Replace academic.h, alpale.h, comdissip.h, comdissipn.h, comdissnew.h by modules Remove unused clesph0.h
File size: 61.8 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	! + +++ READ FORCINGS: END +++
596	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
597
598
599	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
600	! +--SISVAT EXECUTION
601	! + ----------------
602
603	call INLANDSIS(SnoMod, BloMod, 1)
604
605
606
607	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
608	! + RETURN RESULTS
609	! + --------------
610	! + Return (compressed) SISVAT variables to LMDZ
611	! +
612	DO ikl = 1, knon ! use only 1:knon (actual ice sheet..)
613	dflux_s(ikl) = dSdTSV(ikl) ! Sens.H.Flux T-Der.
614	dflux_l(ikl) = dLdTSV(ikl) ! Latn.H.Flux T-Der.
615	fluxsens(ikl) = HSs_sv(ikl) ! HS
616	fluxlat(ikl) = HLs_sv(ikl) ! HL
617	evap(ikl) = -1*HLs_sv(ikl) / LHvH2O ! Evaporation
618	erod(ikl) = 0.
619
620	IF (BloMod) THEN
621	! + Blowing snow
622
623	! SLussl(i,j,n)= 0.
624	! #BS SLussl(i,j,n)= !Effective erosion
625	! #BS. (- dbs_ER(ikl))/(dtrhT_SV(ikl)) !~uqs* from previous time step
626	! #BS blowSN(i,j,n)= dt*uss_SV(ikl) !New max. pot. Erosion [kg/m2]
627	! #BS. *rhT_SV(ikl) !(further bounded in sisvat_bsn.f)
628	! #BS erprev(i,j,n) = dbs_Er(ikl)/dt__SV
629	erod(ikl) = dbs_Er(ikl) / dt__SV
630	ENDIF
631
632	! + Check snow thickness, substract if too thick, add if too thin
633
634	sissnow(ikl) = 0. !()
635	DO isn = 1, isnoSV(ikl)
636	sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn)
637	END DO
638
639	IF (sissnow(ikl) <= sn_low) THEN !add snow
640	IF (isnoSV(ikl)>=1) THEN
641	dzsnSV(ikl, 1) = dzsnSV(ikl, 1) + sn_add / max(ro__SV(ikl, 1), epsi)
642	toicSV(ikl) = toicSV(ikl) - sn_add
643	ELSE
644	WRITE(, ) 'Attention, bare ice... point ', ikl
645	isnoSV(ikl) = 1
646	istoSV(ikl, 1) = 0
647	ro__SV(ikl, 1) = 350.
648	dzsnSV(ikl, 1) = sn_add / max(ro__SV(ikl, 1), epsi) ! 1.
649	eta_SV(ikl, 1) = epsi
650	TsisSV(ikl, 1) = min(TsisSV(ikl, 0), TfSnow - 0.2)
651	G1snSV(ikl, 1) = 0.
652	G2snSV(ikl, 1) = 0.3
653	agsnSV(ikl, 1) = 10.
654	toicSV(ikl) = toicSV(ikl) - sn_add
655	END IF
656	END IF
657
658	IF (sissnow(ikl) >= sn_upp) THEN !thinnen snow layer below
659	dzsnSV(ikl, 1) = dzsnSV(ikl, 1) / sn_div
660	toicSV(ikl) = toicSV(ikl) + dzsnSV(ikl, 1) * ro__SV(ikl, 1) / sn_div
661	END IF
662
663	sissnow(ikl) = 0.
664	qsnow(ikl) = 0.
665	snow(ikl) = 0.
666	snowhgt(ikl) = 0.
667
668	DO isn = 1, isnoSV(ikl)
669	sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn)
670	snowhgt(ikl) = snowhgt(ikl) + dzsnSV(ikl, isn)
671	! Etienne: check calc qsnow
672	qsnow(ikl) = qsnow(ikl) + rhoWat * eta_SV(ikl, isn) * dzsnSV(ikl, isn)
673	END DO
674
675	zfra(ikl) = max(min(isnoSV(ikl) - iiceSV(ikl), 1), 0)
676	! Etienne: comment following line
677	! snow(ikl) = sissnow(ikl)+toicSV(ikl)
678	snow(ikl) = sissnow(ikl)
679
680	to_ice(ikl) = toicSV(ikl)
681	runoff_lic(ikl) = RnofSV(ikl) ! RunOFF: intensity (flux due to melting + liquid precip)
682	fqfonte(ikl)= max(0., (wem_SV(ikl)-wer_SV(ikl))/dtime) ! net melting = melting - refreezing
683	ffonte(ikl)=fqfonte(ikl)*Lf_H2O
684
685	qsol(ikl) = 0.
686	DO isl = -nsol, 0
687	tsoil(ikl, 1 - isl) = TsisSV(ikl, isl) ! Soil Temperature
688	! Etienne: check calc qsol
689	qsol(ikl) = qsol(ikl) &
690	+ eta_SV(ikl, isl) * dz_dSV(isl)
691	END DO
692	agesno(ikl) = agsnSV(ikl, isnoSV(ikl)) ! [day]
693
694	alb1(ikl) = alb1sv(ikl) ! Albedo VIS
695	! alb2(ikl) = ((So1dSV - f1) * alb1sv(ikl) &
696	! & + So2dSV * alb2sv(ikl) + So3dSV * alb3sv(ikl)) / f1
697	alb2(ikl)=alb2sv(ikl)
698	! Albedo NIR
699	alb3(ikl) = alb3sv(ikl) ! Albedo FIR
700	! 6 band Albedo
701	alb6(ikl,:)=alb6sv(ikl,:)
702
703	tsurf_new(ikl) = Tsrfsv(ikl)
704
705	qsurf(ikl) = QsT_SV(ikl)
706	emis_new(ikl) = eps0SL(ikl)
707	z0m(ikl) = Z0m_SV(ikl)
708	z0h(ikl) = Z0h_SV(ikl)
709
710
711	END DO
712
713	IF (ok_outfor) THEN
714	ikl= gp_outfor
715	WRITE(un_outfor, *) '+++++++++++', rlon(ikl2i(ikl)), rlat(ikl2i(ikl)),alt(ikl),'+++++++++++'
716	WRITE(un_outfor, *) isnoSV(ikl), alb_SV(ikl), Z0m_SV(ikl), Z0h_SV(ikl),HSs_sv(ikl),HLs_sv(ikl),alb1(ikl),alb2(ikl)
717	WRITE(un_outfor, *) dzsnSV(ikl, :)
718	WRITE(un_outfor, *) TsisSV(ikl, :)
719	WRITE(un_outfor, *) ro__SV(ikl, :)
720	WRITE(un_outfor, *) eta_SV(ikl, :)
721	WRITE(un_outfor, *) G1snSV(ikl, :)
722	WRITE(un_outfor, *) G2snSV(ikl, :)
723	WRITE(un_outfor, *) agsnSV(ikl, :)
724	WRITE(un_outfor, *) istoSV(ikl, :)
725	WRITE(un_outfor, *) DOPsnSV(ikl, :)
726	ENDIF
727
728
729
730	! + -----------------------------
731	! + END --- RETURN RESULTS
732	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
733	IF (lafin) THEN
734	fichnom = "restartsis.nc"
735	CALL sisvatredem("restartsis.nc", ikl2i, rlon, rlat)
736
737	IF (ok_outfor) THEN
738	close(unit = un_outfor)
739	END IF
740	END IF
741
742	END SUBROUTINE surf_inlandsis
743
744
745	!=======================================================================
746
747	SUBROUTINE get_soil_levels(dz1, dz2, lambda)
748	! ======================================================================
749	! Routine to compute the vertical discretization of the soil in analogy
750	! to LMDZ. In LMDZ it is done in soil.F, which is not used in the case
751	! of SISVAT, therefore it's needed here.
752
753	USE lmdz_phys_mpi_data, ONLY: is_mpi_root
754	USE lmdz_phys_para
755	USE VAR_SV
756
757
758	! INCLUDE "dimsoil.h"
759
760	REAL, DIMENSION(nsoilmx), INTENT(OUT) :: dz2, dz1
761	REAL, INTENT(OUT) :: lambda
762
763
764	!-----------------------------------------------------------------------
765	! Depthts:
766	! --------
767	REAL fz, rk, fz1, rk1, rk2
768	REAL min_period, dalph_soil
769	INTEGER ierr, jk
770
771	fz(rk) = fz1 * (dalph_soil**rk - 1.) / (dalph_soil - 1.)
772
773	! WRITE(,)'Start soil level computation'
774	!-----------------------------------------------------------------------
775	! Calculation of some constants
776	! NB! These constants do not depend on the sub-surfaces
777	!-----------------------------------------------------------------------
778	!-----------------------------------------------------------------------
779	! ground levels
780	! grnd=z/l where l is the skin depth of the diurnal cycle:
781	!-----------------------------------------------------------------------
782
783	min_period = 1800. ! en secondes
784	dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
785	!!! !$OMP MASTER
786	! IF (is_mpi_root) THEN
787	! OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
788	! IF (ierr == 0) THEN ! Read file only if it exists
789	! READ(99,*) min_period
790	! READ(99,*) dalph_soil
791	! PRINT*,'Discretization for the soil model'
792	! PRINT*,'First level e-folding depth',min_period, &
793	! ' dalph',dalph_soil
794	! CLOSE(99)
795	! END IF
796	! ENDIF
797	!!! !$OMP END MASTER
798	! CALL bcast(min_period)
799	! CALL bcast(dalph_soil)
800
801	! la premiere couche represente un dixieme de cycle diurne
802	fz1 = SQRT(min_period / 3.14)
803
804	DO jk = 1, nsoilmx
805	rk1 = jk
806	rk2 = jk - 1
807	dz2(jk) = fz(rk1) - fz(rk2)
808	ENDDO
809	DO jk = 1, nsoilmx - 1
810	rk1 = jk + .5
811	rk2 = jk - .5
812	dz1(jk) = 1. / (fz(rk1) - fz(rk2))
813	ENDDO
814	lambda = fz(.5) * dz1(1)
815	DO jk = 1, nsoilmx
816	rk = jk
817	rk1 = jk + .5
818	rk2 = jk - .5
819	ENDDO
820
821	END SUBROUTINE get_soil_levels
822
823
824	!===========================================================================
825
826	SUBROUTINE SISVAT_ini(knon)
827
828	!C +------------------------------------------------------------------------+
829	!C \| MAR SISVAT_ini Jd 11-10-2007 MAR \|
830	!C \| SubRoutine SISVAT_ini generates non time dependant SISVAT parameters \|
831	!C +------------------------------------------------------------------------+
832	!C \| PARAMETERS: klonv: Total Number of columns = \|
833	!C \| ^^^^^^^^^^ = Total Number of continental grid boxes \|
834	!C \| X Number of Mosaic Cell per grid box \|
835	!C \| \|
836	!C \| INPUT: dt__SV : Time Step [s] \|
837	!C \| ^^^^^ dz_dSV : Layer Thickness [m] \|
838	!C \| \|
839	!C \| OUTPUT: [-] \|
840	!C \| ^^^^^^ rocsSV : Soil Contrib. to (ro c)_s exclud.Water [J/kg/K] \|
841	!C \| etamSV : Soil Minimum Humidity [m3/m3] \|
842	!C \| (based on a prescribed Soil Relative Humidity) \|
843	!C \| s1__SV : Factor of eta**( b+2) in Hydraul.Diffusiv. \|
844	!C \| s2__SV : Factor of eta**( b+2) in Hydraul.Conduct. \|
845	!C \| aKdtSV : KHyd: Piecewise Linear Profile: a * dt [m] \|
846	!C \| bKdtSV : KHyd: Piecewise Linear Profile: b * dt [m/s] \|
847	!C \| dzsnSV(0): Soil first Layer Thickness [m] \|
848	!C \| dzmiSV : Distance between two contiguous levels [m] \|
849	!C \| dz78SV : 7/8 (Layer Thickness) [m] \|
850	!C \| dz34SV : 3/4 (Layer Thickness) [m] \|
851	!C \| dz_8SV : 1/8 (Layer Thickness) [m] \|
852	!C \| dzAvSV : 1/8 dz_(i-1) + 3/4 dz_(i) + 1/8 dz_(i+1) [m] \|
853	!C \| dtz_SV : dt/dz [s/m] \|
854	!C \| OcndSV : Swab Ocean / Soil Ratio [-] \|
855	!C \| Implic : Implicit Parameter (0.5: Crank-Nicholson) \|
856	!C \| Explic : Explicit Parameter = 1.0 - Implic \|
857	!C \| \|
858	!C \| # OPTIONS: #ER: Richards Equation is not smoothed \|
859	!C \| # ^^^^^^^ #kd: De Ridder Discretization \|
860	!C \| # #SH: Hapex-Sahel Values !
861	!C \| \|
862	!C +------------------------------------------------------------------------+
863
864
865	!C +--Global Variables
866	!C + ================
867
868	USE dimphy
869	USE VARphy
870	USE VAR_SV
871	USE VARdSV
872	USE VAR0SV
873	USE VARxSV
874	USE VARtSV
875	USE VARxSV
876	USE VARySV
877	IMPLICIT NONE
878
879
880
881	!C +--Arguments
882	!C + ==================
883	INTEGER, INTENT(IN) :: knon
884
885	!C +--Internal Variables
886	!C + ==================
887
888	INTEGER :: ivt, ist, ikl, isl, isn, ikh
889	INTEGER :: misl_2, nisl_2
890	REAL :: d__eta, eta__1, eta__2, Khyd_1, Khyd_2
891	REAL, PARAMETER :: RHsMin = 0.001 ! Min.Soil Relative Humidity
892	REAL :: PsiMax ! Max.Soil Water Potential
893	REAL :: a_Khyd, b_Khyd ! Water conductivity
894
895
896	!c #WR REAL :: Khyd_x,Khyd_y
897
898
899
900	!C +--Non Time Dependant SISVAT parameters
901	!C + ====================================
902
903	!C +--Soil Discretization
904	!C + -------------------
905
906	!C +--Numerical Scheme Parameters
907	!C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^
908	Implic = 0.75 ! 0.5 <==> Crank-Nicholson
909	Explic = 1.00 - Implic !
910
911	!C +--Soil/Snow Layers Indices
912	!C + ^^^^^^^^^^^^^^^^^^^^^^^^
913	DO isl = -nsol, 0
914	islpSV(isl) = isl + 1
915	islpSV(isl) = min(islpSV(isl), 0)
916	islmSV(isl) = isl - 1
917	islmSV(isl) = max(-nsol, islmSV(isl))
918	END DO
919
920	DO isn = 1, nsno
921	isnpSV(isn) = isn + 1
922	isnpSV(isn) = min(isnpSV(isn), nsno)
923	END DO
924
925	!C +--Soil Layers Thicknesses
926	!C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
927	! Not used here as LMDZ method is applied, see SUBROUTINE get_soil_levels!
928	!c #kd IF (nsol.gt.4) THEN
929	!c #kd DO isl=-5,-nsol,-1
930	!c #kd dz_dSV(isl)= 1.
931	!c #kd END DO
932	!c #kd END IF
933
934	! IF (nsol.NE.4) THEN
935	! DO isl= 0,-nsol,-1
936	! misl_2 = -mod(isl,2)
937	! nisl_2 = -isl/2
938	! dz_dSV(isl)=(((1-misl_2) * 0.001
939	! . + misl_2 * 0.003) * 10*(nisl_2)) 4.
940	!C +... dz_dSV(0) = Hapex-Sahel Calibration: 4 mm
941
942	!c +SH dz_dSV(isl)=(((1-misl_2) * 0.001
943	!c +SH. + misl_2 * 0.003) * 10*(nisl_2)) 1.
944
945	!c #05 dz_dSV(isl)=(((1-misl_2) * 0.001
946	!c #05. + misl_2 * 0.008) * 10*(nisl_2)) 0.5
947	! END DO
948	! dz_dSV(0) = 0.001
949	! dz_dSV(-1) = dz_dSV(-1) - dz_dSV(0) + 0.004
950	! END IF
951
952	zz_dSV = 0.
953	DO isl = -nsol, 0
954	dzmiSV(isl) = 0.500 * (dz_dSV(isl) + dz_dSV(islmSV(isl)))
955	dziiSV(isl) = 0.500 * dz_dSV(isl) / dzmiSV(isl)
956	dzi_SV(isl) = 0.500 * dz_dSV(islmSV(isl)) / dzmiSV(isl)
957	dtz_SV(isl) = dt__SV / dz_dSV(isl)
958	dtz_SV2(isl) = 1. / dz_dSV(isl)
959	dz78SV(isl) = 0.875 * dz_dSV(isl)
960	dz34SV(isl) = 0.750 * dz_dSV(isl)
961	dz_8SV(isl) = 0.125 * dz_dSV(isl)
962	dzAvSV(isl) = 0.125 * dz_dSV(islmSV(isl)) &
963	+ 0.750 * dz_dSV(isl) &
964	+ 0.125 * dz_dSV(islpSV(isl))
965	zz_dSV = zz_dSV + dz_dSV(isl)
966	END DO
967	DO ikl = 1, knon !v
968	dzsnSV(ikl, 0) = dz_dSV(0)
969	END DO
970
971	!C +--Conversion to a 50 m Swab Ocean Discretization
972	!C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
973	OcndSV = 0.
974	DO isl = -nsol, 0
975	OcndSV = OcndSV + dz_dSV(isl)
976	END DO
977	OcndSV = 50. / OcndSV
978
979
980	!C +--Secondary Soil Parameters
981	!C + -------------------------------
982
983	DO ist = 0, nsot
984	rocsSV(ist) = (1.0 - etadSV(ist)) * 1.2E+6 ! Soil Contrib. to (ro c)_s
985	s1__SV(ist) = bCHdSV(ist) & ! Factor of (eta)**(b+2)
986	* psidSV(ist) * Ks_dSV(ist) & ! in DR97, Eqn.(3.36)
987	/ (etadSV(ist)**(bCHdSV(ist) + 3.)) !
988	s2__SV(ist) = Ks_dSV(ist) & ! Factor of (eta)**(2b+3)
989	/ (etadSV(ist)*(2. bCHdSV(ist) + 3.)) ! in DR97, Eqn.(3.35)
990
991	!C +--Soil Minimum Humidity (from a prescribed minimum relative Humidity)
992	!C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
993	Psimax = -(log(RHsMin)) / 7.2E-5 ! DR97, Eqn 3.15 Inversion
994	etamSV(ist) = etadSV(ist) &
995	* (PsiMax / psidSV(ist))**(-min(10., 1. / bCHdSV(ist)))
996	END DO
997	etamSV(12) = 0.
998
999	!C +--Piecewise Hydraulic Conductivity Profiles
1000	!C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1001	DO ist = 0, nsot
1002
1003	d__eta = etadSV(ist) / nkhy
1004	eta__1 = 0.
1005	eta__2 = d__eta
1006	DO ikh = 0, nkhy
1007	Khyd_1 = s2__SV(ist) & ! DR97, Eqn.(3.35)
1008	* (eta__1 *(2. bCHdSV(ist) + 3.)) !
1009	Khyd_2 = s2__SV(ist) &!
1010	* (eta__2 *(2. bCHdSV(ist) + 3.)) !
1011
1012	a_Khyd = (Khyd_2 - Khyd_1) / d__eta !
1013	b_Khyd = Khyd_1 - a_Khyd * eta__1 !
1014	!c #WR Khyd_x = a_Khyd*eta__1 +b_Khyd !
1015	!c #WR Khyd_y = a_Khyd*eta__2 +b_Khyd !
1016	aKdtSV(ist, ikh) = a_Khyd * dt__SV !
1017	bKdtSV(ist, ikh) = b_Khyd * dt__SV !
1018
1019	eta__1 = eta__1 + d__eta
1020	eta__2 = eta__2 + d__eta
1021	END DO
1022	END DO
1023
1024
1025
1026	END SUBROUTINE SISVAT_ini
1027
1028
1029	!***************************************************************************
1030
1031	SUBROUTINE sisvatetat0(fichnom, ikl2i)
1032
1033	USE dimphy
1034	USE lmdz_grid_phy
1035	USE lmdz_phys_para
1036
1037	USE iostart
1038	USE VAR_SV
1039	USE VARdSV
1040	USE VARxSV
1041	USE VARtSV
1042	USE indice_sol_mod
1043
1044	IMPLICIT NONE
1045	!======================================================================
1046	! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
1047	! Objet: Lecture du fichier de conditions initiales pour SISVAT
1048	!======================================================================
1049	include "clesphys.h"
1050	include "compbl.h"
1051
1052	!======================================================================
1053	CHARACTER(LEN = *) :: fichnom
1054
1055	INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
1056	REAL, DIMENSION(klon) :: rlon
1057	REAL, DIMENSION(klon) :: rlat
1058
1059	! les variables globales ecrites dans le fichier restart
1060	REAL, DIMENSION(klon) :: isno
1061	REAL, DIMENSION(klon) :: ispi
1062	REAL, DIMENSION(klon) :: iice
1063	REAL, DIMENSION(klon) :: rusn
1064	REAL, DIMENSION(klon, nsno) :: isto
1065
1066	REAL, DIMENSION(klon, nsismx) :: Tsis
1067	REAL, DIMENSION(klon, nsismx) :: eta
1068	REAL, DIMENSION(klon, nsismx) :: ro
1069
1070	REAL, DIMENSION(klon, nsno) :: dzsn
1071	REAL, DIMENSION(klon, nsno) :: G1sn
1072	REAL, DIMENSION(klon, nsno) :: G2sn
1073	REAL, DIMENSION(klon, nsno) :: agsn
1074
1075	REAL, DIMENSION(klon) :: toic
1076
1077	INTEGER :: isl, ikl, i, isn, errT, erreta, errro, errdz, snopts
1078	CHARACTER (len = 2) :: str2
1079	LOGICAL :: found
1080
1081	errT = 0
1082	errro = 0
1083	erreta = 0
1084	errdz = 0
1085	snopts = 0
1086	! Ouvrir le fichier contenant l'etat initial:
1087
1088	CALL open_startphy(fichnom)
1089
1090	! Lecture des latitudes, longitudes (coordonnees):
1091
1092	CALL get_field("latitude", rlat, found)
1093	CALL get_field("longitude", rlon, found)
1094
1095	CALL get_field("n_snows", isno, found)
1096	IF (.NOT. found) THEN
1097	PRINT*, 'phyetat0: Le champ <n_snows> est absent'
1098	PRINT *, 'fichier startsisvat non compatible avec sisvatetat0'
1099	ENDIF
1100
1101	CALL get_field("n_ice_top", ispi, found)
1102	CALL get_field("n_ice", iice, found)
1103	CALL get_field("surf_water", rusn, found)
1104
1105
1106	CALL get_field("to_ice", toic, found)
1107	IF (.NOT. found) THEN
1108	PRINT*, 'phyetat0: Le champ <to_ice> est absent'
1109	toic(:) = 0.
1110	ENDIF
1111
1112	DO isn = 1, nsno
1113	IF (isn<=99) THEN
1114	WRITE(str2, '(i2.2)') isn
1115	CALL get_field("AGESNOW" // str2, &
1116	agsn(:, isn), found)
1117	ELSE
1118	PRINT*, "Trop de couches"
1119	CALL abort
1120	ENDIF
1121	ENDDO
1122	DO isn = 1, nsno
1123	IF (isn<=99) THEN
1124	WRITE(str2, '(i2.2)') isn
1125	CALL get_field("DZSNOW" // str2, &
1126	dzsn(:, isn), found)
1127	ELSE
1128	PRINT*, "Trop de couches"
1129	CALL abort
1130	ENDIF
1131	ENDDO
1132	DO isn = 1, nsno
1133	IF (isn<=99) THEN
1134	WRITE(str2, '(i2.2)') isn
1135	CALL get_field("G2SNOW" // str2, &
1136	G2sn(:, isn), found)
1137	ELSE
1138	PRINT*, "Trop de couches"
1139	CALL abort
1140	ENDIF
1141	ENDDO
1142	DO isn = 1, nsno
1143	IF (isn<=99) THEN
1144	WRITE(str2, '(i2.2)') isn
1145	CALL get_field("G1SNOW" // str2, &
1146	G1sn(:, isn), found)
1147	ELSE
1148	PRINT*, "Trop de couches"
1149	CALL abort
1150	ENDIF
1151	ENDDO
1152	DO isn = 1, nsismx
1153	IF (isn<=99) THEN
1154	WRITE(str2, '(i2.2)') isn
1155	CALL get_field("ETA" // str2, &
1156	eta(:, isn), found)
1157	ELSE
1158	PRINT*, "Trop de couches"
1159	CALL abort
1160	ENDIF
1161	ENDDO
1162	DO isn = 1, nsismx
1163	IF (isn<=99) THEN
1164	WRITE(str2, '(i2.2)') isn
1165	CALL get_field("RO" // str2, &
1166	ro(:, isn), found)
1167	ELSE
1168	PRINT*, "Trop de couches"
1169	CALL abort
1170	ENDIF
1171	ENDDO
1172	DO isn = 1, nsismx
1173	IF (isn<=99) THEN
1174	WRITE(str2, '(i2.2)') isn
1175	CALL get_field("TSS" // str2, &
1176	Tsis(:, isn), found)
1177	ELSE
1178	PRINT*, "Trop de couches"
1179	CALL abort
1180	ENDIF
1181	ENDDO
1182	DO isn = 1, nsno
1183	IF (isn<=99) THEN
1184	WRITE(str2, '(i2.2)') isn
1185	CALL get_field("HISTORY" // str2, &
1186	isto(:, isn), found)
1187	ELSE
1188	PRINT*, "Trop de couches"
1189	CALL abort
1190	ENDIF
1191	ENDDO
1192	WRITE(, )'Read ', fichnom, ' finished!!'
1193
1194	!*********************************************************************************
1195	! Compress restart file variables for SISVAT
1196
1197	DO ikl = 1, klon
1198	i = ikl2i(ikl)
1199	IF (i > 0) THEN
1200	isnoSV(ikl) = INT(isno(i)) ! Nb Snow/Ice Lay.
1201	ispiSV(ikl) = INT(ispi(i)) ! Nb Supr.Ice Lay.
1202	iiceSV(ikl) = INT(iice(i)) ! Nb Ice Lay.
1203
1204	DO isl = -nsol, 0
1205	ro__SV(ikl, isl) = ro(i, nsno + 1 - isl) !
1206	eta_SV(ikl, isl) = eta(i, nsno + 1 - isl) ! Soil Humidity
1207	!hjp 15/10/2010
1208	IF (eta_SV(ikl, isl) <= 1.e-6) THEN !hj check
1209	eta_SV(ikl, isl) = 1.e-6
1210	ENDIF
1211	TsisSV(ikl, isl) = Tsis(i, nsno + 1 - isl) ! Soil Temperature
1212	IF (TsisSV(ikl, isl) <= 1.) THEN !hj check
1213	! errT=errT+1
1214	TsisSV(ikl, isl) = 273.15 - 0.2 ! Etienne: negative temperature since soil is ice
1215	ENDIF
1216
1217	END DO
1218	WRITE(, )'Copy histo', ikl
1219
1220	DO isn = 1, isnoSV(ikl) !nsno
1221	snopts = snopts + 1
1222	IF (isto(i, isn) > 10.) THEN !hj check
1223	WRITE(, )'Irregular isto', ikl, i, isn, isto(i, isn)
1224	isto(i, isn) = 1.
1225	ENDIF
1226
1227	istoSV(ikl, isn) = INT(isto(i, isn)) ! Snow History
1228	ro__SV(ikl, isn) = ro(i, isn) ! [kg/m3]
1229	eta_SV(ikl, isn) = eta(i, isn) ! [m3/m3]
1230	TsisSV(ikl, isn) = Tsis(i, isn) ! [K]
1231
1232	IF (TsisSV(ikl, isn) <= 1.) THEN !hj check
1233	errT = errT + 1
1234	TsisSV(ikl, isn) = TsisSV(ikl, 0)
1235	ENDIF
1236	IF (TsisSV(ikl, isn) <= 1.) THEN !hj check
1237	TsisSV(ikl, isn) = 263.15
1238	ENDIF
1239	IF (eta_SV(ikl, isn) < 1.e-9) THEN !hj check
1240	eta_SV(ikl, isn) = 1.e-6
1241	erreta = erreta + 1
1242	ENDIF
1243	IF (ro__SV(ikl, isn) <= 10.) THEN !hj check
1244	ro__SV(ikl, isn) = 11.
1245	errro = errro + 1
1246	ENDIF
1247	WRITE(, )ikl, i, isn, Tsis(i, isn), G1sn(i, isn)
1248	G1snSV(ikl, isn) = G1sn(i, isn) ! [-] [-]
1249	G2snSV(ikl, isn) = G2sn(i, isn) ! [-] [0.0001 m]
1250	dzsnSV(ikl, isn) = dzsn(i, isn) ! [m]
1251	agsnSV(ikl, isn) = agsn(i, isn) ! [day]
1252	END DO
1253	rusnSV(ikl) = rusn(i) ! Surficial Water
1254	toicSV(ikl) = toic(i) ! bilan snow to ice
1255	END IF
1256	END DO
1257
1258	END SUBROUTINE sisvatetat0
1259
1260
1261	!======================================================================
1262	SUBROUTINE sisvatredem(fichnom, ikl2i, rlon, rlat)
1263
1264
1265
1266	!======================================================================
1267	! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
1268	! Objet: Ecriture de l'etat de redemarrage pour SISVAT
1269	!======================================================================
1270	USE lmdz_grid_phy
1271	USE lmdz_phys_para
1272	USE iostart
1273	USE VAR_SV
1274	USE VARxSV
1275	USE VARySV !hj tmp 12 03 2010
1276	USE VARtSV
1277	USE indice_sol_mod
1278	USE dimphy
1279
1280	IMPLICIT NONE
1281
1282	include "clesphys.h"
1283	include "compbl.h"
1284
1285	!======================================================================
1286
1287	CHARACTER(LEN = *) :: fichnom
1288	INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
1289	REAL, DIMENSION(klon), INTENT(IN) :: rlon
1290	REAL, DIMENSION(klon), INTENT(IN) :: rlat
1291
1292	! les variables globales ecrites dans le fichier restart
1293	REAL, DIMENSION(klon) :: isno
1294	REAL, DIMENSION(klon) :: ispi
1295	REAL, DIMENSION(klon) :: iice
1296	REAL, DIMENSION(klon, nsnowmx) :: isto
1297
1298	REAL, DIMENSION(klon, nsismx) :: Tsis
1299	REAL, DIMENSION(klon, nsismx) :: eta
1300	REAL, DIMENSION(klon, nsnowmx) :: dzsn
1301	REAL, DIMENSION(klon, nsismx) :: ro
1302	REAL, DIMENSION(klon, nsnowmx) :: G1sn
1303	REAL, DIMENSION(klon, nsnowmx) :: G2sn
1304	REAL, DIMENSION(klon, nsnowmx) :: agsn
1305	REAL, DIMENSION(klon) :: IRs
1306	REAL, DIMENSION(klon) :: LMO
1307	REAL, DIMENSION(klon) :: rusn
1308	REAL, DIMENSION(klon) :: toic
1309	REAL, DIMENSION(klon) :: Bufs
1310	REAL, DIMENSION(klon) :: alb1, alb2, alb3
1311
1312	INTEGER isl, ikl, i, isn, ierr
1313	CHARACTER (len = 2) :: str2
1314	INTEGER :: pass
1315
1316	isno(:) = 0
1317	ispi(:) = 0
1318	iice(:) = 0
1319	IRs(:) = 0.
1320	LMO(:) = 0.
1321	eta(:, :) = 0.
1322	Tsis(:, :) = 0.
1323	isto(:, :) = 0
1324	ro(:, :) = 0.
1325	G1sn(:, :) = 0.
1326	G2sn(:, :) = 0.
1327	dzsn(:, :) = 0.
1328	agsn(:, :) = 0.
1329	rusn(:) = 0.
1330	toic(:) = 0.
1331	Bufs(:) = 0.
1332	alb1(:) = 0.
1333	alb2(:) = 0.
1334	alb3(:) = 0.
1335
1336	!***************************************************************************
1337	! Uncompress SISVAT output variables for storage
1338
1339	DO ikl = 1, klon
1340	i = ikl2i(ikl)
1341	IF (i > 0) THEN
1342	isno(i) = 1. * isnoSV(ikl) ! Nb Snow/Ice Lay.
1343	ispi(i) = 1. * ispiSV(ikl) ! Nb Supr.Ice Lay.
1344	iice(i) = 1. * iiceSV(ikl) ! Nb Ice Lay.
1345
1346	! IRs(i) = IRs_SV(ikl)
1347	! LMO(i) = LMO_SV(ikl)
1348
1349	DO isl = -nsol, 0 !
1350	eta(i, nsno + 1 - isl) = eta_SV(ikl, isl) ! Soil Humidity
1351	Tsis(i, nsno + 1 - isl) = TsisSV(ikl, isl) ! Soil Temperature
1352	ro(i, nsno + 1 - isl) = ro__SV(ikl, isl) ! [kg/m3]
1353	END DO
1354
1355	DO isn = 1, nsno
1356	isto(i, isn) = 1. * istoSV(ikl, isn) ! Snow History
1357	ro(i, isn) = ro__SV(ikl, isn) ! [kg/m3]
1358	eta(i, isn) = eta_SV(ikl, isn) ! [m3/m3]
1359	Tsis(i, isn) = TsisSV(ikl, isn) ! [K]
1360	G1sn(i, isn) = G1snSV(ikl, isn) ! [-] [-]
1361	G2sn(i, isn) = G2snSV(ikl, isn) ! [-] [0.0001 m]
1362	dzsn(i, isn) = dzsnSV(ikl, isn) ! [m]
1363	agsn(i, isn) = agsnSV(ikl, isn) ! [day]
1364	END DO
1365	rusn(i) = rusnSV(ikl) ! Surficial Water
1366	toic(i) = toicSV(ikl) ! to ice
1367	alb1(i) = alb1sv(ikl)
1368	alb2(i) = alb2sv(ikl)
1369	alb3(i) = alb3sv(ikl)
1370	! Bufs(i) = BufsSV(ikl)
1371	END IF
1372	END DO
1373
1374	CALL open_restartphy(fichnom)
1375
1376	DO pass = 1, 2
1377	CALL put_field(pass, "longitude", &
1378	"Longitudes de la grille physique", rlon)
1379	CALL put_field(pass, "latitude", "Latitudes de la grille physique", rlat)
1380
1381	CALL put_field(pass, "n_snows", "number of snow/ice layers", isno)
1382	CALL put_field(pass, "n_ice_top", "number of top ice layers", ispi)
1383	CALL put_field(pass, "n_ice", "number of ice layers", iice)
1384	CALL put_field(pass, "IR_soil", "Soil IR flux", IRs)
1385	CALL put_field(pass, "LMO", "Monin-Obukhov Scale", LMO)
1386	CALL put_field(pass, "surf_water", "Surficial water", rusn)
1387	CALL put_field(pass, "snow_buffer", "Snow buffer layer", Bufs)
1388	CALL put_field(pass, "alb_1", "albedo sw", alb1)
1389	CALL put_field(pass, "alb_2", "albedo nIR", alb2)
1390	CALL put_field(pass, "alb_3", "albedo fIR", alb3)
1391	CALL put_field(pass, "to_ice", "Snow passed to ice", toic)
1392
1393	DO isn = 1, nsno
1394	IF (isn<=99) THEN
1395	WRITE(str2, '(i2.2)') isn
1396	CALL put_field(pass, "AGESNOW" // str2, &
1397	"Age de la neige layer No." // str2, &
1398	agsn(:, isn))
1399	ELSE
1400	PRINT*, "Trop de couches"
1401	CALL abort
1402	ENDIF
1403	ENDDO
1404	DO isn = 1, nsno
1405	IF (isn<=99) THEN
1406	WRITE(str2, '(i2.2)') isn
1407	CALL put_field(pass, "DZSNOW" // str2, &
1408	"Snow/ice thickness layer No." // str2, &
1409	dzsn(:, isn))
1410	ELSE
1411	PRINT*, "Trop de couches"
1412	CALL abort
1413	ENDIF
1414	ENDDO
1415	DO isn = 1, nsno
1416	IF (isn<=99) THEN
1417	WRITE(str2, '(i2.2)') isn
1418	CALL put_field(pass, "G2SNOW" // str2, &
1419	"Snow Property 2, layer No." // str2, &
1420	G2sn(:, isn))
1421	ELSE
1422	PRINT*, "Trop de couches"
1423	CALL abort
1424	ENDIF
1425	ENDDO
1426	DO isn = 1, nsno
1427	IF (isn<=99) THEN
1428	WRITE(str2, '(i2.2)') isn
1429	CALL put_field(pass, "G1SNOW" // str2, &
1430	"Snow Property 1, layer No." // str2, &
1431	G1sn(:, isn))
1432	ELSE
1433	PRINT*, "Trop de couches"
1434	CALL abort
1435	ENDIF
1436	ENDDO
1437	DO isn = 1, nsismx
1438	IF (isn<=99) THEN
1439	WRITE(str2, '(i2.2)') isn
1440	CALL put_field(pass, "ETA" // str2, &
1441	"Soil/snow water content layer No." // str2, &
1442	eta(:, isn))
1443	ELSE
1444	PRINT*, "Trop de couches"
1445	CALL abort
1446	ENDIF
1447	ENDDO
1448	DO isn = 1, nsismx !nsno
1449	IF (isn<=99) THEN
1450	WRITE(str2, '(i2.2)') isn
1451	CALL put_field(pass, "RO" // str2, &
1452	"Snow density layer No." // str2, &
1453	ro(:, isn))
1454	ELSE
1455	PRINT*, "Trop de couches"
1456	CALL abort
1457	ENDIF
1458	ENDDO
1459	DO isn = 1, nsismx
1460	IF (isn<=99) THEN
1461	WRITE(str2, '(i2.2)') isn
1462	CALL put_field(pass, "TSS" // str2, &
1463	"Soil/snow temperature layer No." // str2, &
1464	Tsis(:, isn))
1465	ELSE
1466	PRINT*, "Trop de couches"
1467	CALL abort
1468	ENDIF
1469	ENDDO
1470	DO isn = 1, nsno
1471	IF (isn<=99) THEN
1472	WRITE(str2, '(i2.2)') isn
1473	CALL put_field(pass, "HISTORY" // str2, &
1474	"Snow history layer No." // str2, &
1475	isto(:, isn))
1476	ELSE
1477	PRINT*, "Trop de couches"
1478	CALL abort
1479	ENDIF
1480	ENDDO
1481
1482	CALL enddef_restartphy
1483	ENDDO
1484	CALL close_restartphy
1485
1486	END SUBROUTINE sisvatredem
1487
1488	END MODULE surf_inlandsis_mod

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

Download in other formats:

Original Format