Context Navigation

surf_inlandsis_mod.F90 @ 3814

Last change on this file since 3814 was 3792, checked in by evignon, 4 years ago
Ajout de INLANDSIS, nouvelle interface entre LMDZ et la neige de SISVAT Etienne, 04/01/2021
File size: 60.9 KB

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

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

Download in other formats:

Original Format