Context Navigation

surf_inlandsis_mod.f90 @ 5413

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

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

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

Download in other formats:

Original Format