Context Navigation

surf_inlandsis_mod.F90 @ 5018

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

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

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

Download in other formats:

Original Format