Context Navigation

surf_sisvat_mod.F90 @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

File size: 84.4 KB

Rev	Line
[1865]	1	MODULE surf_sisvat_mod
	2	USE dimphy
	3	IMPLICIT NONE
	4
	5	INTEGER, PARAMETER :: nsnowmx=35
	6	INTEGER, PARAMETER :: nsismx=46 ! = nsnowmx + nsoilmx
	7
	8	CONTAINS
	9
	10
	11
	12	SUBROUTINE surf_sisvat(knon,rlon,rlat, ikl2i, itime, dtime, debut, lafin, &
	13	rmu0, swdown, lwdown, pexner, ps, p1lay, &
	14	precip_rain, precip_snow, precip_snow_adv, snow_adv, &
	15	bl_height, wind_velo, temp_air, dens_air, spechum, tsurf, &
	16	rugos, snow_cont_air, alb_soil, slope, cloudf, &
	17	radsol, qsol, tsoil, snow, snowhgt, qsnow, to_ice, sissnow, agesno, &
	18	AcoefH, AcoefQ, BcoefH, BcoefQ, cdragh, &
	19	runoff_lic, evap, fluxsens, fluxlat, dflux_s, dflux_l, &
	20	tsurf_new, alb1, alb2, alb3, &
	21	emis_new, z0_new, qsurf)
	22
	23	! +------------------------------------------------------------------------+
	24	! \| \|
	25	! \| SubRoutine surf_sisvat: Interface between LMDZ and landice scheme \|
	26	! \| of the SISVAT (Soil/Ice Snow Vegetation Atmosphere Transfer Scheme)\|
	27	! \| \|
	28	! \| Author: Heinz Juergen Punge, LSCE June 2009 \|
	29	! \| based on the MAR-SISVAT interface by Hubert Gallee \|
	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	! \| ^^^^^ VegMod: SISVAT is set up when .T. \|
	46	! \| SnoMod: Snow Pack is set up when .T. \|
	47	! \| reaLBC: Update Bound.Condit.when .T. \|
	48	! \| \|
	49	! \| INPUT (via MODULES VARxSV, VARySV, VARtSV) \|
	50	! \| ^^^^^ xxxxSV: SISVAT/LMDZ interfacing variables \|
	51	! \| \|
	52	! \| Preprocessing Option: SISVAT PHYSICS \|
	53	! \| ^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^ \|
	54	! \| # #HY \|
	55	! \| # #SN: Snow Model \|
	56	! \| # #BS: Blowing Snow Parameterization \|
	57	! \| \|
	58	! \| # #DS: diffuse radiation differing from direct \|
	59	! \| (variable RADsod must still be included) \|
	60	! \| # #CP: SBL, Col de Porte \|
	61	! \| # #cp Solar Radiation, Col de Porte \|
	62	! \| # #AG: Snow Ageing, Col de Porte \|
	63	! \| \|
	64	! \| Preprocessing Option: SISVAT IO \|
	65	! \| ^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^ \|
	66	! \| FILE \| CONTENT \|
	67	! \| ~~~~~~~~~~~~~~~~~~~~~+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \|
	68	! \| # ANI.yyyymmdd.LAB.nc \| #NC: OUTPUT on NetCDF File (Stand Alone EXP.) \|
	69	! \| # SISVAT_iii_jjj_n \| #WV: OUTPUT on ASCII File (SISVAT Variables) \|
	70	! \| # SISVATtroubles \| #VF: OUTPUT on ASCII File (SISVAT Troubles) \|
	71	! \| \| \|
	72	! \| # \| #ES: OUTPUT/Verification: Energy Conservation \|
	73	! \| # \| #E2: OUTPUT/Verification: Energy Consrv.2e pt.\|
	74	! \| \| (no premature stop) \|
	75	! \| # \| #MW: OUTPUT/Verification: H2O Conservation \|
	76	! \| # \| #MS: OUTPUT/Verification: * Mass Conservation \|
	77	! \| # \| #MI: OUTPUT/Verification: SeaIce Conservation \|
	78	! \| \| \|
	79	! \| # SISVAT__zSn.OUT \| #as: OUTPUT/Verification: Snow Layers Agrega. \|
	80	! \| # SISVAT__SnO.OUT \| #aw: OUTPUT/Verification: Albedo Parameteriz. \|
	81	! \| # SISVATe_qSn.OUT \| #em: OUTPUT/Verification: Energy/Water Budget \|
	82	! \| # SISVATu_qSn.OUT \| #su: OUTPUT/Verification: Slush Parameteriz. \|
	83	! \| # SISVATw_qSo.OUT \| #mw: OUTPUT/Verif+Detail: H2O Conservation \|
	84	! \| \| \|
	85	! \| # SISVAT__GSn.OUT \| #VP: OUTPUT/Verification: Snow Properties \|
	86	! \| # SISVAT__wEq.OUT \| #EQ: OUTPUT/Verification: Snow/Ice Water Eqv. \|
	87	! \| \| \|
	88	! \| # \| #VR: VERIFICATION OUTPUT \|
	89	! \| # \| #WR: Additional OUTPUT \|
	90	! \| \|
	91	! +------------------------------------------------------------------------+
	92
	93	USE VAR_SV
	94	USE VARdSV
	95	USE VARxSV
	96	USE VARySV
	97	USE VARtSV
	98
	99	USE VARdCP
	100	USE VARphy, ra_earth=>ra
	101	USE YOMCST
	102
	103	IMPLICIT NONE
	104
	105	! +--INTERFACE Variables
	106	! + ===================
	107
	108	include "dimsoil.h"
	109
	110
	111	! +--Global Variables
	112	! + ================
	113	! Input Variables for SISVAT
	114	INTEGER, INTENT(IN) :: knon
	115	INTEGER, INTENT(IN) :: itime
	116	REAL, INTENT(IN) :: dtime
	117	LOGICAL, INTENT(IN) :: debut ! true if first step
	118	LOGICAL, INTENT(IN) :: lafin ! true if last step
	119
	120	INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i ! Index Decompression
	121	REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat
	122	REAL, DIMENSION(klon), INTENT(IN) :: rmu0 ! cos sol. zenith angle
	123	REAL, DIMENSION(klon), INTENT(IN) :: swdown !
	124	REAL, DIMENSION(klon), INTENT(IN) :: lwdown !
	125	REAL, DIMENSION(klon), INTENT(IN) :: pexner ! Exner potential
	126	REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
	127	REAL, DIMENSION(klon), INTENT(IN) :: precip_snow_adv, snow_adv
	128	!Snow Drift
	129	REAL, DIMENSION(klon), INTENT(IN) :: bl_height, wind_velo
	130	REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum, ps,p1lay
	131	REAL, DIMENSION(klon), INTENT(IN) :: dens_air, tsurf
	132	REAL, DIMENSION(klon), INTENT(IN) :: rugos,snow_cont_air
	133	REAL, DIMENSION(klon), INTENT(IN) :: alb_soil, slope
	134	REAL, DIMENSION(klon), INTENT(IN) :: cloudf
	135	REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ
	136	REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ
	137	REAL, DIMENSION(klon), INTENT(IN) :: cdragh
	138
	139	! Variables exchanged between LMDZ and SISVAT
	140	REAL, DIMENSION(klon,nsoilmx), INTENT(OUT) :: tsoil ! Soil Temperature
	141	REAL, DIMENSION(klon), INTENT(OUT) :: qsol ! Soil Water Content
	142	REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! Tot snow mass [kg/m2]
	143	REAL, DIMENSION(klon), INTENT(IN) :: radsol ! Surface absorbed rad.
	144
	145	! Output Variables from SISVAT
	146	REAL, DIMENSION(klon), INTENT(OUT) :: alb1 ! Albedo SW
	147	REAL, DIMENSION(klon), INTENT(OUT) :: alb2,alb3 ! Albedo LW
	148	REAL, DIMENSION(klon), INTENT(OUT) :: emis_new ! Surface Emissivity
	149	REAL, DIMENSION(klon), INTENT(OUT) :: z0_new ! Momentum Roughn Lgt
	150	REAL, DIMENSION(klon), INTENT(OUT) :: runoff_lic ! Runoff
	151	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s ! d/dT sens. ht flux
	152	REAL, DIMENSION(klon), INTENT(OUT) :: dflux_l ! d/dT latent ht flux
	153	REAL, DIMENSION(klon), INTENT(OUT) :: fluxsens ! Sensible ht flux
	154	REAL, DIMENSION(klon), INTENT(OUT) :: fluxlat ! Latent heat flux
	155	REAL, DIMENSION(klon), INTENT(OUT) :: evap ! Evaporation
	156	REAL, DIMENSION(klon), INTENT(OUT) :: agesno ! Snow age (top layer)
	157	REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! Surface Temperature
	158	REAL, DIMENSION(klon), INTENT(OUT) :: qsurf ! Surface Humidity
	159	REAL, DIMENSION(klon), INTENT(OUT) :: qsnow ! Total H2O snow[kg/m2]
	160	REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt ! Snow height (m)
	161	REAL, DIMENSION(klon), INTENT(OUT) :: to_ice ! Snow passed to ice
	162	REAL, DIMENSION(klon), INTENT(OUT) :: sissnow ! Snow in model (kg/m2)
	163
	164	! +--OUTPUT for NetCDF File
	165	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	166	! REAL, DIMENSION(klonv) :: SOsoKL ! Abs Solar Radiation
	167	! REAL, DIMENSION(klonv) :: IRsoKL ! Abs IR Radiation
	168	! REAL, DIMENSION(klonv) :: HSsoKL ! Abs Sensible Ht Flux
	169	! REAL, DIMENSION(klonv) :: HLsoKL ! Abs Latent Heat Flux
	170	! REAL, DIMENSION(klonv) :: HLs_KL ! Evaporation
	171	! REAL, DIMENSION(klonv) :: HLv_KL ! Transpiration
	172	!
	173	! REAL, DIMENSION(klonv) :: SOsoNC ! Abs Solar Radiation
	174	! REAL, DIMENSION(klonv) :: IRsoNC ! Abs IR Radiation
	175	! REAL, DIMENSION(klonv) :: HSsoNC ! Abs Sensible Ht Flux
	176	! REAL, DIMENSION(klonv) :: HLsoNC ! Abs Latent Heat Flux
	177	! REAL, DIMENSION(klonv) :: HLs_NC ! Evaporation
	178	! REAL, DIMENSION(klonv) :: HLv_NC ! Transpiration
	179	!
	180	! REAL, DIMENSION(klonv,nsoilmx) :: eta_NC ! nsoilmx=nsol+1
	181	! REAL, DIMENSION(klonv,nsoilmx) :: tsolNC !
	182	! REAL, DIMENSION(klonv) :: snowNC !
	183	! INTEGER, DIMENSION(klonv) :: isnoNC !
	184	! INTEGER, DIMENSION(klonv) :: ispiNC !
	185	! INTEGER, DIMENSION(klonv) :: iiceNC !
	186	! REAL, DIMENSION(klonv) :: swaNC !
	187	! REAL, DIMENSION(klonv) :: swsNC !
	188	! INTEGER, DIMENSION(klonv,nsno) :: istoNC !
	189	! REAL, DIMENSION(klonv,nsno) :: dzsnNC !
	190	! REAL, DIMENSION(klonv,nsno) :: rhosnNC !
	191	! REAL, DIMENSION(klonv,nsno) :: etasnNC !
	192	! REAL, DIMENSION(klonv,nsno) :: tsnNC !
	193	! REAL, DIMENSION(klonv,nsno) :: g1snNC !
	194	! REAL, DIMENSION(klonv,nsno) :: g2snNC !
	195	! REAL, DIMENSION(klonv,nsno) :: agsnNC !
	196	! REAL, DIMENSION(klonv) :: meltNC !
	197	! REAL, DIMENSION(klonv) :: refrNC !
	198	! REAL, DIMENSION(klonv) :: alb1NC !
	199	! REAL, DIMENSION(klonv) :: alb2NC !
	200	! REAL, DIMENSION(klonv) :: rnofNC !
	201
	202	! + Optional Variables:
	203	! +--V, dT(a-s) Time Moving Averages
	204	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	205
	206	! #AW real V__mem(klon,ntaver) ! ntaver defined in LMDZ_SL.inc
	207	! #AW real VVmmem(klon) !
	208	! #AW common/SVeSBLmem/V__mem,VVmmem !
	209	! #AH real T__mem(klon,ntaver) !
	210	! #AH real dTmmem(klon) !
	211	! #AH common/STeSBLmem/T__mem,dTmmem !
	212	! +--u, uT, us* Time Moving Averages
	213	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	214	! #AM real u__mem(klon,ntaver) ! ntaver defined in LMDZ_SL.inc
	215	! #AT real uT_mem(klon,ntaver) !
	216	! #AS real us_mem(klon,ntaver) !
	217	! #AM common/S_eSBLmem/u__mem !
	218	! #AT. ,uT_mem !
	219	! #AS. ,us_mem !
	220	!! #AH REAL, DIMENSION(klonv,ntaver) :: T__mem !
	221	!! #AH REAL, DIMENSION(klonv) :: dTmmem !
	222	!! #AW REAL, DIMENSION(klonv,ntaver) :: V__mem !
	223	!! #AM REAL, DIMENSION(klonv) :: u__mem !
	224	!! #AT REAL, DIMENSION(klonv) :: uT_mem !
	225	!! #AS REAL, DIMENSION(klonv) :: us_mem !
	226	!! #AW REAL, DIMENSION(klonv) :: VVmmem !
	227
	228
	229	! +--Internal Variables
	230	! + ===================
	231
	232	CHARACTER(len=20) :: fichnom, fn_outfor ! Name for output file
	233	INTEGER :: i, ig, ikl, isl, isn, nt
	234	INTEGER :: gp_outfor, un_outfor
	235
	236	REAL, PARAMETER :: f1=0.5
	237	REAL, PARAMETER :: sn_upp=5000.,sn_low=500.
	238	REAL, PARAMETER :: sn_add=400.,sn_div=2.
	239	! snow mass upper,lower limit,
	240	! added mass/division lowest layer
	241	REAL, PARAMETER :: c1_zuo=12.960e+4, c2_zuo=2.160e+6
	242	REAL, PARAMETER :: c3_zuo=1.400e+2, czemin=1.e-3
	243	! Parameters for drainage
	244	! c1_zuo/ 2.796e+4/,c2_zuo/2.160e+6/,c3_zuo/1.400e+2/ ! Tuning
	245	! +... Run Off Parameters
	246	! + 864001.5 day ...25 days (Modif. ETH Camp: 86400*0.3day)
	247	! + (Zuo and Oerlemans 1996, J.Glacio. 42, 305--317)
	248
	249	REAL, DIMENSION(klon) :: eps0SL ! surface Emissivity
	250	REAL :: zsigma, Ua_min, Us_min
	251	REAL :: lambda ! Par. soil discret.
	252	REAL, DIMENSION(nsoilmx), SAVE :: dz1,dz2 ! Soil layer thicknesses
	253	!$OMP THREADPRIVATE(dz1)
	254	LOGICAL, SAVE :: firstcall=.TRUE.,SnoMod, ok_outfor=.FALSE.!$OMP THREADPRIVATE(firstcall)
	255
	256
	257
	258	! + Optional:
	259	!c #BW INTEGER :: noUNIT
	260	!c #BW REAL :: BlowST,SnowSB
	261
	262
	263	! +--Internal Variables
	264	! + ==================
	265
	266	INTEGER :: ivg,iso
	267
	268	!========================================================================
	269
	270	SnoMod=.true.
	271	zsigma=1000.
	272	dt__SV=dtime
	273	IF (ok_outfor) THEN
	274	un_outfor=51 ! unit number for point output
	275	gp_outfor=79! 633 !79 ! grid point number for point output
	276	fn_outfor='outfor_SV.dat'
	277	END IF
	278
	279	! write(,)'Start of simulation? ',debut !hj
	280	IF (debut) THEN
	281	firstcall=.TRUE.
	282	INI_SV=.false.
	283	ELSE
	284	firstcall=.false.
	285	INI_SV=.true.
	286	END IF
	287
	288	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	289
	290
	291	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	292	! + INITIALISATION: BEGIN +++
	293	! + -------------------------
	294	! +
	295	! + Compute soil discretization (as for LMDZ)
	296	! + -----------------------------------------
	297	IF (firstcall) THEN
	298
	299	! +--Array size
	300	klonv=klon
	301	knonv=knon
	302	write(,)'klon',klon,'klonv',klonv,'knon',knon
	303
	304
	305	CALL INIT_VARtSV
	306	CALL INIT_VARxSV
	307	CALL INIT_VARySV
	308
	309	eps0SL(:)=0.
	310	! +--Soil layer thickness
	311	! + -----------------------
	312	! write(*,'(/a)') 'Start SISVAT init: soil discretization ', nsoilmx
	313	CALL get_soil_levels(dz1,dz2,lambda)
	314	lambSV=lambda
	315	dz1_SV(1:knon,1:) = 0.
	316	dz2_SV(1:knon,1:) = 0.
	317
	318	DO isl = -nsol,0
	319	dz_dSV(isl) = 0.5e-3*dz2(1-isl) ! Soil layer thickness
	320	DO ikl=1,knon
	321	dz1_SV(ikl,isl) = dz1(1-isl) !1.e-3*
	322	dz2_SV(ikl,isl) = dz2(1-isl) !1.e-3*
	323	END DO
	324	! IF (knon > 0) THEN
	325	! write(,)'level:',dz_dSV(isl),dz1_SV(1,isl),dz2_SV(1,isl)
	326	! END IF
	327	END DO
	328
	329	DO ikl=1,knon
	330	eps0SL(ikl )= 1.
	331	alb0SV(ikl) = alb_soil(ikl) ! Soil Albedo
	332
	333	! + Soil Upward IR Flux, Water Fluxes, roughness length
	334	IRs_SV(ikl) = &
	335	-eps0SL(ikl)* rsigma*temp_air(ikl) & ! Upward IR Flux
	336	temp_air(ikl)temp_air(ikl)*temp_air(ikl)
	337	TvegSV(ikl) = temp_air(ikl)
	338	! + Soil
	339	DO isl = -nsol,0
	340	TsisSV(ikl,isl) = temp_air(ikl) !tsoil(ikl,1-isl) Soil Temperature
	341	eta_SV(ikl,isl) = 0.0001 !etasoil(ikl,1-isl)Soil Water[m3/m3]
	342	ro__SV(ikl,isl) = 50. !rosoil(ikl,1-isl) soil water volumic mass
	343	END DO
	344	END DO
	345
	346	! +--Surface Fall Line Slope
	347	! + -----------------------
	348	IF (SnoMod) THEN
	349	DO ikl=1,knon
	350	slopSV(ikl) = slope(ikl)
	351	SWf_SV(ikl) = & ! Normalized Decay of the
	352	exp(-dt__SV & ! Surficial Water Content
	353	/(c1_zuo & !(Zuo and Oerlemans 1996,
	354	+c2_zuoexp(-c3_zuoabs(slopSV(ikl))))) ! J.Glacio. 42, 305--317)
	355	END DO
	356	END IF
	357
	358	! +--SBL Characteristics history, not stored at restart.
	359	DO ikl=1,knon
	360	DO nt=1,ntaver ! #AW
	361	V__mem(ikl,nt)=wind_velo(ikl) ! #AH
	362	T__mem(ikl,nt)=temp_air(ikl)-tsurf(ikl)
	363	END DO
	364	END DO
	365
	366
	367	! + SISVAT_ini (as for use with MAR, but not computing soil layers)
	368	! + -------------------------------------------------------------
	369	! write(*,'(/a)') 'Start SISVAT initialization: SISVAT_ini'
	370	CALL SISVAT_ini(knon)
	371
	372	! open output file
	373	IF (ok_outfor) THEN
	374	open(unit=un_outfor,status='new',file=fn_outfor)
	375	rewind un_outfor
	376	ikl=gp_outfor ! index sur la grille land ice
	377	ig=611 ! index sur la grille globale
	378	write(un_outfor,501) fn_outfor, ikl, rlon(ig),rlat(ig)
	379	501 format(/,a18,/,'Grid point ',i4,' Long',f9.4,' Lat ',f9.4 &
	380	& ,/,'++++++++++++++++++++++++++++++++++++++++++++++', &
	381	& '++++++++++++++++++++++++++++++++++++++++++++++', &
	382	& '++++++++++++++++++++++++++++++++++++++++++++++', &
	383	& /,' SWdown + IRdown + Wind + Temp. + Humid. ', &
	384	& '+ Press +Precip_l+Precip_s+ Tsrf + Clouds +' &
	385	& '+ Zenith + BLhgt + Densair+ Exner +' &
	386	& ,/,' sol_SV + IRd_SV + VV__SV + TaT_SV + QaT_SV ', &
	387	& '+ ps__SV + drr_SV + dsn_SV + Tsf_SV + cld_SV +' &
	388	& '+ coszSV + za__SV + rhT_SV + ExnrSV+' &
	389	& ,/,' W/m2 + W/m2 + m/s + K + kg/kg ', &
	390	& '+ Pa + kg/m2/s+ kg/m2/s+ K + /1 +' &
	391	& '+ - + m + kg/m3 + +' &
	392	& ,/,'++++++++++++++++++++++++++++++++++++++++++++++', &
	393	& '++++++++++++++++++++++++++++++++++++++++++++++', &
	394	& '++++++++++++++++++++++++++++++++++++++++++++++')
	395	END IF
	396
	397	! +--Read restart file
	398	! + =================================================
	399
	400	! Martin
	401	PRINT*, 'On debranche sisvatetat0'
	402	! Martin
	403
	404	! CALL sisvatetat0("startsis.nc",ikl2i)
	405
	406
	407	END IF ! firstcall
	408	! +
	409	! + +++ INITIALISATION: END +++
	410	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	411
	412
	413
	414	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	415	! + READ FORCINGS
	416	! + ------------------------
	417	! + Update Forcings for SISVAT given by the LMDZ model.
	418	! +
	419	DO ikl=1,knon
	420
	421	! +--Atmospheric Forcing (INPUT)
	422	! + ^^^^^^^^^^^^^^^^^^^ ^^^^^
	423	zSBLSV = 1000. ! [m]
	424	za__SV(ikl) = bl_height(ikl) ! Height boundary layer [m]
	425	Ua_min = epsi !
	426	! #VM Ua_min = 0.2 * sqrt(za__SV(ikl) ) !
	427	Ua_min = 0.2 * sqrt(za__SV(ikl) ) !
	428	VV__SV(ikl) = max(Ua_min, wind_velo(ikl)) ! Wind velocity [m/s]
	429	Us_min = 0.01
	430	us__SV(ikl) = max(Us_min, us__SV(ikl) )
	431	TaT_SV(ikl) = temp_air(ikl) ! BL top Temperature [K]
	432	ExnrSV(ikl) = pexner(ikl) ! Exner potential
	433	rhT_SV(ikl) = dens_air(ikl) ! Air density
	434	QaT_SV(ikl) = spechum(ikl) ! Specific humidity
	435	! #VX dQa_SV(ikl) = 0. !hj dtDiff/zsigma(mz) ! Water Vapor Flux
	436	ps__SV(ikl) = ps(ikl) ! surface pressure [Pa]
	437	p1l_SV(ikl) = p1lay(ikl) ! lowest atm. layer press[Pa]
	438
	439	! +--Energy Fluxes (INPUT)
	440	! + ^^^^^^^^^^^^^ ^^^^^
	441	coszSV(ikl) = max(czemin,rmu0(ikl)) ! cos(zenith.Dist.)
	442	sol_SV(ikl) = swdown(ikl) ! downward Solar
	443	IRd_SV(ikl) = lwdown(ikl) ! downward IR
	444	rsolSV(ikl) = radsol(ikl) ! surface absorbed rad.
	445	!hj 110511
	446	!hj 121011 IRs_SV(ikl) = min(IRs_SV(ikl),-epsi) ! check upward IR
	447	!hj
	448	! +--Water Fluxes (INPUT)
	449	! + ^^^^^^^^^^^^^ ^^^^^
	450	drr_SV(ikl) = precip_rain(ikl) ! Rain fall rate [kg/m2/s]
	451	dsn_SV(ikl) = precip_snow(ikl) ! Snow fall rate [kg/m2/s]
	452	!c #BS dbsnow = -SLussl(i,j,n) ! Erosion
	453	!c #BS. dtPhys rhT_SV(ikl) /ro_Wat
	454	!c #BS dsnbSV(ikl) = snow_adv(ikl) ! min(max(zero,dbsnow)
	455	!c #BS. / max(epsi,d_snow),unun)
	456	!c #BS dbs_SV(ikl) = snow_cont_air(ikl)
	457	!c #BS blowSN(i,j,n) ! [kg/m2]
	458
	459	! +--Soil/Canopy (INPUT)
	460	! + ^^^^^^^^^^^ ^^^^^
	461	alb0SV(ikl) = alb_soil(ikl) ! Soil background Albedo
	462	AcoHSV(ikl) = AcoefH(ikl)
	463	BcoHSV(ikl) = BcoefH(ikl)
	464	AcoQSV(ikl) = AcoefQ(ikl)
	465	BcoQSV(ikl) = BcoefQ(ikl)
	466	cdH_SV(ikl) = cdragh(ikl)
	467	Tsf_SV(ikl) = tsurf(ikl) !hj 12 03 2010
	468
	469	! +--Energy Fluxes (INPUT/OUTPUT)
	470	! + ^^^^^^^^^^^^^ ^^^^^^^^^^^^
	471	IF (.not.firstcall) THEN
	472	! active hj 110411
	473	cld_SV(ikl) = cloudf(ikl) ! Cloudiness
	474	END IF
	475
	476	! +--Time Averages of wind and surface-atmosphere temperature difference
	477	! +--for turbulence calculations ! #AA
	478	! Update stored arrays
	479	DO nt=1,ntaver-1 ! #AW
	480	V__mem(ikl,nt ) = V__mem(ikl,nt+1) ! #AH
	481	T__mem(ikl,nt ) = T__mem(ikl,nt+1) ! #AA
	482	ENDDO ! #AW
	483	V__mem(ikl,ntaver)=wind_velo(ikl) ! #AH
	484	T__mem(ikl,ntaver)=temp_air(ikl)-tsurf(ikl) ! #AW
	485	! Calculate averages
	486	VVmmem(ikl) = 0.0 ! #AH
	487	dTmmem(ikl) = 0.0 ! #AA
	488	DO nt=1,ntaver ! #AW
	489	VVmmem(ikl) = VVmmem(ikl)+V__mem(ikl,nt) ! #AH
	490	dTmmem(ikl) = dTmmem(ikl)+T__mem(ikl,nt) ! #AA
	491	ENDDO ! #AW
	492	VVmmem(ikl) = VVmmem(ikl)/ntaver ! #AH
	493	dTmmem(ikl) = dTmmem(ikl)/ntaver
	494	END DO
	495
	496	!
	497	! + +++ READ FORCINGS: END +++
	498	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	499
	500	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	501	! + OUTPUT FORCING
	502	!
	503	IF (ok_outfor) THEN
	504	ikl=gp_outfor
	505	write(un_outfor,5000) sol_SV(ikl), IRd_SV(ikl), VV__SV(ikl), &
	506	& TaT_SV(ikl), QaT_SV(ikl), ps__SV(ikl), &
	507	& drr_SV(ikl), dsn_SV(ikl), Tsf_SV(ikl), &
	508	& cld_SV(ikl), coszSV(ikl), za__SV(ikl), &
	509	& rhT_SV(ikl), ExnrSV(ikl)
	510
	511	5000 format(f8.3,' ',f8.3,' ',f8.4,' ',f8.4,' ',f10.8,' ',f8.1,' ', &
	512	& f8.6,' ',f8.6,' ',f8.4,' ',f8.6,' ',f8.6,' ',f8.3,' ', &
	513	& f8.4,' ',f8.6)
	514
	515	ENDIF
	516	!
	517	! + OUTPUT FORCINGS: END +++
	518	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	519
	520
	521
	522	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	523	! +--SISVAT EXECUTION
	524	! + ----------------
	525	! +
	526	! write(,) ' Start SISVAT execution!'
	527
	528	call SISVAT(SnoMod,.false.,1)
	529	! BloMod,jjtime)
	530
	531	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	532	! + RETURN RESULTS
	533	! + --------------
	534	! + Return (compressed) SISVAT variables to LMDZ
	535	! +
	536	DO ikl=1,knon ! use only 1:knon (actual ice sheet..)
	537	runoff_lic(ikl) = RnofSV(ikl)dtime ! RunOFF: intensity time step
	538	dflux_s(ikl) = dSdTSV(ikl) ! Sens.H.Flux T-Der.
	539	dflux_l(ikl) = dLdTSV(ikl) ! Latn.H.Flux T-Der.
	540	fluxsens(ikl) = HSs_sv(ikl) ! HS
	541	fluxlat(ikl) = HLs_sv(ikl) ! HL
	542	evap(ikl) = HLs_sv(ikl)/RLVTT ! Evaporation RLVTT=LhvH2O
	543	z0_new(ikl) = Z0h_SV(ikl) ! Moment.Roughn.L.
	544
	545
	546	snow(ikl) = 0.
	547	snowhgt(ikl) = 0.
	548	qsnow(ikl) = 0.
	549	qsol(ikl) = 0.
	550	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	551	! + Check snow thickness, add if too thin, substract if too thick
	552
	553	sissnow(ikl) = 0. !()
	554	DO isn = 1,isnoSV(ikl)
	555	sissnow(ikl) = sissnow(ikl)+dzsnSV(ikl,isn)* ro__SV(ikl,isn)
	556	END DO
	557
	558	IF (sissnow(ikl) .LE. sn_low) THEN !add snow
	559	IF (isnoSV(ikl).GE.1) THEN
	560	dzsnSV(ikl,1) = dzsnSV(ikl,1) + sn_add/max(ro__SV(ikl,1),epsi)
	561	toicSV(ikl) = toicSV(ikl) - sn_add
	562	ELSE
	563	write(,) 'Attention, bare ice... point ',ikl
	564	isnoSV(ikl) = 1
	565	istoSV(ikl,1) = 100
	566	!ym istoSV(ikl) = 100
	567	ro__SV(ikl,1) = 400.
	568	dzsnSV(ikl,1) = sn_add/max(ro__SV(ikl,1),epsi) ! 1.
	569	eta_SV(ikl,1) = epsi
	570	TsisSV(ikl,1) = min(TsisSV(ikl,0),TfSnow-0.2)
	571	G1snSV(ikl,1) = 99.
	572	G2snSV(ikl,1) = 0.3
	573	agsnSV(ikl,1) = 10.
	574	toicSV(ikl) = toicSV(ikl) - sn_add
	575	END IF
	576	END IF
	577
	578	IF (sissnow(ikl) .ge. sn_upp) THEN !thinnen snow layer below
	579	dzsnSV(ikl,1) = dzsnSV(ikl,1)/sn_div
	580	toicSV(ikl) = toicSV(ikl)+dzsnSV(ikl,1)*ro__SV(ikl,1)/sn_div
	581	END IF
	582
	583	sissnow(ikl) = 0. !()
	584	DO isn = 1,isnoSV(ikl)
	585	sissnow(ikl) = sissnow(ikl)+dzsnSV(ikl,isn)* ro__SV(ikl,isn)
	586	! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	587	snowhgt(ikl) = snowhgt(ikl)+dzsnSV(ikl,isn)
	588	qsnow(ikl) = qsnow(ikl)+1e03eta_SV(ikl,isn)dzsnSV(ikl,isn)
	589	END DO
	590	snow(ikl) = sissnow(ikl)+toicSV(ikl)
	591	to_ice(ikl) = toicSV(ikl)
	592	! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	593
	594
	595	DO isl = -nsol,0
	596	tsoil(ikl,1-isl) = TsisSV(ikl,isl) ! Soil Temperature
	597	qsol(ikl) = qsol(ikl) &
	598	+eta_SV(ikl,isl) * dz_dSV(isl)
	599	END DO
	600	agesno(ikl) = agsnSV(ikl,isnoSV(ikl)) ! [day]
	601
	602	alb1(ikl) = alb1sv(ikl) ! Albedo VIS
	603	alb2(ikl) = ((So1dSV-f1)*alb1sv(ikl) &
	604	& +So2dSValb2sv(ikl)+So3dSValb3sv(ikl))/f1
	605	! Albedo NIR
	606	alb3(ikl) = alb3sv(ikl) ! Albedo FIR
	607	tsurf_new(ikl) = TsfnSV(ikl)
	608	!hj220711 Tsrfsv(ikl) ! Surf.Temperature
	609	! TsisSV(ikl,0) *(1-min(1,isnoSV(ikl))) &
	610	! +TsisSV(ikl,max(1,isnoSV(ikl))) * min(1,isnoSV(ikl))
	611	! tsurf_new(ikl) = max(Ts_Min,tsurf_new(ikl) )
	612	! tsurf_new(ikl) = min(Ts_Max,tsurf_new(ikl) )
	613	qsurf(ikl) = QaT_SV(ikl)
	614	emis_new(ikl) = eps0SL(ikl)
	615
	616	!!!hjp 230611 sorties
	617	!! qsnow(ikl)=TsisSV(ikl,isnoSV(ikl))
	618	!! sissnow(ikl)=TsisSV(ikl,isnoSV(ikl)-1)
	619	!! snowhgt(ikl)=TsisSV(ikl,isnoSV(ikl)-2)
	620	!! qsol(ikl)=dzsnSV(ikl,isnoSV(ikl)-1)
	621	!! agesno(ikl)=dzsnSV(ikl,isnoSV(ikl))
	622
	623
	624	END DO
	625
	626
	627	! + -----------------------------
	628	! + END --- RETURN RESULTS
	629	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	630
	631
	632	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	633	! + Uncompress SISVAT output variables & store
	634	! + -----------------------------
	635	! +
	636	! DO ikl = 1,knon
	637	! i = ikl2i(ikl) ! Compression index
	638	! IF (i.GE.0) THEN
	639	! SOsoNC(i) = SOsoKL(ikl) ! Absorb.Sol.Rad.
	640	! IRsoNC(i) = IRsoKL(ikl) ! Absorb.IR Rad.
	641	! HSsoNC(i) = HSsoKL(ikl) ! HS
	642	! HLsoNC(i) = HLsoKL(ikl) ! HL
	643	! HLs_NC(i) = HLs_KL(ikl) ! Evaporation
	644	! HLv_NC(i) = HLv_KL(ikl) ! Transpiration
	645	!
	646	! DO isl = -nsol,0
	647	! eta_NC(i,1-isl) = eta_SV(ikl,isl) ! Soil Humidity
	648	! tsolNC(i,1-isl) = TsisSV(ikl,isl) ! Soil Temperature
	649	! END DO
	650	! snowNC(i) = snow(ikl) ! Snow mass
	651	! isnoNC(i) = isnoSV(ikl) ! Nb Snow/Ice Lay.
	652	! ispiNC(i) = ispiSV(ikl) ! Nb Supr.Ice Lay.
	653	! iiceNC(i) = iiceSV(ikl) ! Nb Ice Lay.
	654	! swaNC(i) = rusnSV(ikl) ! Surficial Water
	655	! swsNC(i) = SWS_SV(ikl) ! Surficial Wat.St.
	656	! DO isn = 1,nsno
	657	! istoNC(i,isn) = istoSV(ikl,isn) ! [-]
	658	! dzsnNC(i,isn) = dzsnSV(ikl,isn) ! [m]
	659	! rhosnNC(i,isn) = ro__SV(ikl,isn) ! [kg/m3]
	660	! etasnNC(i,isn) = eta_SV(ikl,isn) ! [m3/m3]
	661	! tsnNC(i,isn) = TsisSV(ikl,isn) ! [K]
	662	! g1snNC(i,isn) = G1snSV(ikl,isn) ! [-] [-]
	663	! g2snNC(i,isn) = G2snSV(ikl,isn) ! [-] [0.0001 m]
	664	! agsnNC(i,isn) = agsnSV(ikl,isn) ! [day]
	665	! END DO
	666	!!?c #IB depsubNC(i) = wes_SV(ikl) ! Depo. / Subli.
	667	! meltNC(i) = wem_SV(ikl) ! Melting
	668	! refrNC(i) = wer_SV(ikl) ! Refreezing
	669	! alb1NC(i) = alb1sv(ikl) ! Albedo SW
	670	! alb2NC(i) = (alb1sv(ikl)+3*alb2sv(ikl)+alb3sv(ikl))/5
	671	! Albedo LW
	672	! rnofNC(i) = RnofSV(ikl) ! Run OFF Intensity
	673	!
	674	!!? SL_z0(i) = Z0m_SV(ikl) ! Moment.Roughn.L.
	675	!!? SL_r0(i,j,n) = Z0h_SV(ikl) ! Heat Roughn.L.
	676	!!! zWE_NC(i) = zWE_SV(ikl) ! Current *Thick.
	677	!!! zWEcNC(i) = zWEcSV(ikl) ! Non-Erod.*Thick.
	678	!!! hSalNC(i) = hSalSV(ikl) ! Salt.Layer Height
	679	!!! hsenSL(i) = -SLuts(i,j) * rhAir *cp ! Sensible Heat Flux
	680	!!! hlatSL(i) = -SLuqs(i,j) * rhAir *Lv_H2O ! Latent Heat Flux
	681	! END IF
	682	! END DO
	683
	684	!hj + Put storage to Output file here!
	685	IF (lafin) THEN
	686
	687	fichnom = "restartsis.nc"
	688	CALL sisvatredem("restartsis.nc",ikl2i,rlon,rlat)
	689	IF (ok_outfor) THEN
	690	close(unit=un_outfor)
	691	END IF
	692	END IF
	693	! + -----------------------------
	694	! + END --- RETURN RESULTS
	695	! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	696
	697	END SUBROUTINE surf_sisvat
	698
	699
	700
	701
	702
	703
	704
	705	SUBROUTINE get_soil_levels(dz1, dz2, lambda)
	706	! ======================================================================
	707	! Routine to compute the vertical discretization of the soil in analogy
	708	! to LMDZ. In LMDZ it is done in soil.F, which is not used in the case
	709	! of SISVAT, therefore it's needed here.
	710	!
	711	USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root
	712	USE mod_phys_lmdz_para
	713
	714	INCLUDE "dimsoil.h"
	715
	716	REAL, DIMENSION(nsoilmx), INTENT(OUT) :: dz2, dz1
	717	REAL, INTENT(OUT) :: lambda
	718
	719
	720	!-----------------------------------------------------------------------
	721	! Depthts:
	722	! --------
	723	REAL fz,rk,fz1,rk1,rk2
	724	REAL min_period, dalph_soil
	725	INTEGER ierr,jk
	726
	727	fz(rk)=fz1(dalph_soil*rk-1.)/(dalph_soil-1.)
	728
	729	! write(,)'Start soil level computation'
	730	!-----------------------------------------------------------------------
	731	! Calculation of some constants
	732	! NB! These constants do not depend on the sub-surfaces
	733	!-----------------------------------------------------------------------
	734	!-----------------------------------------------------------------------
	735	! ground levels
	736	! grnd=z/l where l is the skin depth of the diurnal cycle:
	737	!-----------------------------------------------------------------------
	738
	739	min_period=1800. ! en secondes
	740	dalph_soil=2. ! rapport entre les epaisseurs de 2 couches succ.
	741	! !$OMP MASTER
	742	! IF (is_mpi_root) THEN
	743	! OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
	744	! IF (ierr == 0) THEN ! Read file only if it exists
	745	! READ(99,*) min_period
	746	! READ(99,*) dalph_soil
	747	! PRINT*,'Discretization for the soil model'
	748	! PRINT*,'First level e-folding depth',min_period, &
	749	! ' dalph',dalph_soil
	750	! CLOSE(99)
	751	! END IF
	752	! ENDIF
	753	! !$OMP END MASTER
	754	! CALL bcast(min_period)
	755	! CALL bcast(dalph_soil)
	756
	757	! la premiere couche represente un dixieme de cycle diurne
	758	fz1=SQRT(min_period/3.14)
	759
	760	DO jk=1,nsoilmx
	761	rk1=jk
	762	rk2=jk-1
	763	dz2(jk)=fz(rk1)-fz(rk2)
	764	ENDDO
	765	DO jk=1,nsoilmx-1
	766	rk1=jk+.5
	767	rk2=jk-.5
	768	dz1(jk)=1./(fz(rk1)-fz(rk2))
	769	ENDDO
	770	lambda=fz(.5)*dz1(1)
	771	PRINT*,'full layers, intermediate layers (seconds)'
	772	DO jk=1,nsoilmx
	773	rk=jk
	774	rk1=jk+.5
	775	rk2=jk-.5
	776	PRINT *,'fz=', &
	777	fz(rk1)fz(rk2)3.14,fz(rk)fz(rk)3.14
	778	ENDDO
	779
	780	END SUBROUTINE get_soil_levels
	781
	782	SUBROUTINE SISVAT_ini(knon)
	783
	784	!C +------------------------------------------------------------------------+
	785	!C \| MAR SISVAT_ini Jd 11-10-2007 MAR \|
	786	!C \| SubRoutine SISVAT_ini generates non time dependant SISVAT parameters \|
	787	!C +------------------------------------------------------------------------+
	788	!C \| PARAMETERS: klonv: Total Number of columns = \|
	789	!C \| ^^^^^^^^^^ = Total Number of continental grid boxes \|
	790	!C \| X Number of Mosaic Cell per grid box \|
	791	!C \| \|
	792	!C \| INPUT: dt__SV : Time Step [s] \|
	793	!C \| ^^^^^ dz_dSV : Layer Thickness [m] \|
	794	!C \| \|
	795	!C \| OUTPUT: RF__SV : Root Fraction in Layer isl [-] \|
	796	!C \| ^^^^^^ rocsSV : Soil Contrib. to (ro c)_s exclud.Water [J/kg/K] \|
	797	!C \| etamSV : Soil Minimum Humidity [m3/m3] \|
	798	!C \| (based on a prescribed Soil Relative Humidity) \|
	799	!C \| s1__SV : Factor of eta**( b+2) in Hydraul.Diffusiv. \|
	800	!C \| s2__SV : Factor of eta**( b+2) in Hydraul.Conduct. \|
	801	!C \| aKdtSV : KHyd: Piecewise Linear Profile: a * dt [m] \|
	802	!C \| bKdtSV : KHyd: Piecewise Linear Profile: b * dt [m/s] \|
	803	!C \| dzsnSV(0): Soil first Layer Thickness [m] \|
	804	!C \| dzmiSV : Distance between two contiguous levels [m] \|
	805	!C \| dz78SV : 7/8 (Layer Thickness) [m] \|
	806	!C \| dz34SV : 3/4 (Layer Thickness) [m] \|
	807	!C \| dz_8SV : 1/8 (Layer Thickness) [m] \|
	808	!C \| dzAvSV : 1/8 dz_(i-1) + 3/4 dz_(i) + 1/8 dz_(i+1) [m] \|
	809	!C \| dtz_SV : dt/dz [s/m] \|
	810	!C \| OcndSV : Swab Ocean / Soil Ratio [-] \|
	811	!C \| Implic : Implicit Parameter (0.5: Crank-Nicholson) \|
	812	!C \| Explic : Explicit Parameter = 1.0 - Implic \|
	813	!C \| \|
	814	!C \| # OPTIONS: #ER: Richards Equation is not smoothed \|
	815	!C \| # ^^^^^^^ #kd: De Ridder Discretization \|
	816	!C \| # #SH: Hapex-Sahel Values !
	817	!C \| \|
	818	!C +------------------------------------------------------------------------+
	819	!
	820	!
	821
	822	!C +--Global Variables
	823	!C + ================
	824
	825	USE VARphy, ra_earth=>ra
	826	USE VAR_SV
	827	USE VARdSV
	828	USE VAR0SV
	829	USE VARxSV
	830	USE VARtSV
	831	USE VARxSV
	832	USE VARySV
	833	IMPLICIT NONE
	834
	835
	836
	837	!C +--Arguments
	838	!C + ==================
	839	INTEGER,INTENT(IN) :: knon
	840
	841	!C +--Internal Variables
	842	!C + ==================
	843
	844	INTEGER :: ivt ,ist ,ivg ,ikl ,isl ,isn ,ikh
	845	INTEGER :: misl_2,nisl_2
	846	REAL :: zDepth
	847	REAL :: d__eta,eta__1,eta__2,Khyd_1,Khyd_2
	848	REAL,PARAMETER :: RHsMin= 0.001 ! Min.Soil Relative Humidity
	849	REAL :: PsiMax ! Max.Soil Water Potential
	850	REAL :: a_Khyd,b_Khyd ! Piecewis.Water Conductivity
	851
	852
	853	!c #WR REAL :: Khyd_x,Khyd_y
	854
	855
	856
	857	!C +--Non Time Dependant SISVAT parameters
	858	!C + ====================================
	859
	860	!C +--Soil Discretization
	861	!C + -------------------
	862
	863	!C +--Numerical Scheme Parameters
	864	!C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^
	865	Implic = 0.75 ! 0.5 <==> Crank-Nicholson
	866	Explic = 1.00 - Implic !
	867
	868	!C +--Soil/Snow Layers Indices
	869	!C + ^^^^^^^^^^^^^^^^^^^^^^^^
	870	DO isl=-nsol,0
	871	islpSV(isl) = isl+1
	872	islpSV(isl) = min( islpSV(isl),0)
	873	islmSV(isl) = isl-1
	874	islmSV(isl) = max(-nsol,islmSV(isl))
	875	END DO
	876
	877	DO isn=1,nsno
	878	isnpSV(isn) = isn+1
	879	isnpSV(isn) = min( isnpSV(isn),nsno)
	880	END DO
	881
	882	!C +--Soil Layers Thicknesses
	883	!C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	884	! Not used here as LMDZ method is applied, see SUBROUTINE get_soil_levels!
	885	!c #kd IF (nsol.gt.4) THEN
	886	!c #kd DO isl=-5,-nsol,-1
	887	!c #kd dz_dSV(isl)= 1.
	888	!c #kd END DO
	889	!c #kd END IF
	890	!
	891	! IF (nsol.ne.4) THEN
	892	! DO isl= 0,-nsol,-1
	893	! misl_2 = -mod(isl,2)
	894	! nisl_2 = -isl/2
	895	! dz_dSV(isl)=(((1-misl_2) * 0.001
	896	! . + misl_2 * 0.003) * 10*(nisl_2)) 4.
	897	!C +... dz_dSV(0) = Hapex-Sahel Calibration: 4 mm
	898	!
	899	!c +SH dz_dSV(isl)=(((1-misl_2) * 0.001
	900	!c +SH. + misl_2 * 0.003) * 10*(nisl_2)) 1.
	901	!
	902	!c #05 dz_dSV(isl)=(((1-misl_2) * 0.001
	903	!c #05. + misl_2 * 0.008) * 10*(nisl_2)) 0.5
	904	! END DO
	905	! dz_dSV(0) = 0.001
	906	! dz_dSV(-1) = dz_dSV(-1) - dz_dSV(0) + 0.004
	907	! END IF
	908
	909	zz_dSV = 0.
	910	DO isl=-nsol,0
	911	dzmiSV(isl) = 0.500*(dz_dSV(isl) +dz_dSV(islmSV(isl)))
	912	dziiSV(isl) = 0.500* dz_dSV(isl) /dzmiSV(isl)
	913	dzi_SV(isl) = 0.500* dz_dSV(islmSV(isl))/dzmiSV(isl)
	914	dtz_SV(isl) = dt__SV /dz_dSV(isl)
	915	dz78SV(isl) = 0.875* dz_dSV(isl)
	916	dz34SV(isl) = 0.750* dz_dSV(isl)
	917	dz_8SV(isl) = 0.125* dz_dSV(isl)
	918	dzAvSV(isl) = 0.125* dz_dSV(islmSV(isl)) &
	919	& + 0.750* dz_dSV(isl) &
	920	& + 0.125* dz_dSV(islpSV(isl))
	921	!c #ER dz78SV(isl) = dz_dSV(isl)
	922	!c #ER dz34SV(isl) = dz_dSV(isl)
	923	!c #ER dz_8SV(isl) = 0.
	924	!c #ER dzAvSV(isl) = dz_dSV(isl)
	925	zz_dSV = zz_dSV+dz_dSV(isl)
	926	END DO
	927	DO ikl=1,knon !v
	928	dzsnSV(ikl,0) = dz_dSV(0)
	929	END DO
	930
	931	!C +--Conversion to a 50 m Swab Ocean Discretization
	932	!C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	933	OcndSV = 0.
	934	DO isl=-nsol,0
	935	OcndSV = OcndSV +dz_dSV(isl)
	936	END DO
	937	OcndSV = 50. /OcndSV
	938
	939
	940	!C +--Secondary Vegetation Parameters
	941	!C + -------------------------------
	942	!
	943	!C +--Minimum Stomatal Resistance (Hapex Sahel Data)
	944	!C + (Taylor et al. 1997, J.Hydrol 188-189, p.1047)
	945	!C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	946	!c #SH DO ivg=1,3 !
	947	!c #SH StodSV(ivg) = 210. ! Millet
	948	!c #SH END DO !
	949	!c #SH StodSV( 4) = 120. ! Sparse Tiger Bush
	950	!c #SH DO ivg=5,6 !
	951	!c #SH StodSV(ivg) = 80. ! Dense Tiger Bush
	952	!c #SH END DO !
	953	!c #SH StodSV( 7) = 80. ! Low Trees (Fallow)
	954	!c #SH StodSV( 10) = 80. !
	955	!
	956	!C +--Minimum Stomatal Resistance (Tropical Forest)
	957	!C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	958	!c #SH StodSV( 8) = 60. ! Medium Trees
	959	!c #SH StodSV( 11) = 60. !
	960	!c #SH StodSV( 9) = 40. ! High Trees
	961	!c #SH StodSV( 12) = 40. !
	962
	963	!C +--Root Fraction
	964	!C + ^^^^^^^^^^^^^
	965	!C + * GENERAL REFERENCE
	966	!C + Jackson et al., 1996: A global analysis of root distributions for
	967	!C + terrestrial biomes. In Oecologia, 108, 389-411.
	968
	969	!C + * ROOT PROFILE
	970	!C + The cumulative root fraction Y is given by
	971	!C + Y = 1 - beta**d with d the depth (in cm),
	972	!C + beta a coefficient (vegetation dependent).
	973
	974	!C + * BETA VALUES (for 11 world biomes)
	975	!C + 1 boreal forest 0.943
	976	!C + 2 crops 0.961
	977	!C + 3 desert 0.975
	978	!C + 4 sclerophyllous shrubs 0.964
	979	!C + 5 temperate coniferous forest 0.976
	980	!C + 6 temperate deciduous forest 0.966
	981	!C + 7 temperate grassland 0.943
	982	!C + 8 tropical deciduous forest 0.961
	983	!C + 9 tropical evergreen forest 0.962
	984	!C + 10 tropical grassland savanna 0.972
	985	!C + 11 tundra 0.914
	986	!
	987	!C + * ADVISED BETA VALUES FOR MAR
	988	!C + (see 'block data SISVAT_dat', variable rbtdSV)
	989	!C +
	990	!C + SVAT veg. type default West Africa
	991	!C + 0 barren soil 0.000 0.000
	992	!C + 1 crops low 0.961 (2) 0.961 (2)
	993	!C + 2 crops medium 0.961 (2) 0.961 (2)
	994	!C + 3 crops high 0.961 (2) 0.961 (2)
	995	!C + 4 grass low 0.943 (7) 0.943 (7)
	996	!C + 5 grass medium 0.943 (7) 0.964 (4)
	997	!C + 6 grass high 0.943 (7) 0.972 (10)
	998	!C + 7 broadleaf low 0.966 (6) 0.968 (4,10)
	999	!C + 8 broadleaf medium 0.966 (6) 0.962 (8,9)
	1000	!C + 9 broadleaf high 0.966 (6) 0.962 (8,9)
	1001	!C + 10 needleleaf low 0.976 (5) 0.971 (5,6)
	1002	!C + 11 needleleaf medium 0.976 (5) 0.976 (5)
	1003	!C + 12 needleleaf high 0.976 (5) 0.976 (5)
	1004
	1005	!C + Numbers between brackets refer to Jackson's biomes. For more details
	1006	!C + about some choices, see the correspondance between the IGBP and SVAT
	1007	!C + vegetation classes (i.e. in NESTOR).
	1008
	1009	!C + * WARNING
	1010	!C + Most of the roots are located in the first 2 m of soil. The root
	1011	!C + fraction per layer depends on the definition of the soil layer
	1012	!C + thickness. It will get wrong if a thick layer is defined around 2 m
	1013	!C + deep.
	1014
	1015	write(*,'(/a)') 'ROOT PROFILES (Jackson, 1996) :'
	1016
	1017	DO ivt = 0, nvgt
	1018	zDepth = 0.
	1019	DO isl = 0, -nsol, -1
	1020	IF (ivt .ne. 0) THEN
	1021	RF__SV(ivt,isl) = rbtdSV(ivt)*zDepth &
	1022	& (1. - rbtdSV(ivt)*(dz_dSV(isl)100) )
	1023	zDepth = zDepth + dz_dSV(isl)*100 !in cm
	1024	ELSE
	1025	RF__SV(ivt,isl) = 0.
	1026	END IF
	1027	END DO
	1028	write(*,'(a,i2,a,i3,a,99f10.5:)') &
	1029	& ' RF__SV(', ivt, ',', -nsol, ':0) =', RF__SV(ivt,:)
	1030	END DO
	1031	! write(6, format( &
	1032	! & ' NOTE: If root fraction is not close to 0 around 2 m deep,', &
	1033	! &/,' Then you should redefine the soil layer thicknesses.', &
	1034	! &/,' See the code for more details.'))
	1035
	1036
	1037	!C +--Secondary Soil Parameters
	1038	!C + -------------------------------
	1039
	1040	DO ist=0,nsot
	1041	rocsSV(ist)=(1.0-etadSV(ist))*1.2E+6 ! Soil Contrib. to (ro c)_s
	1042	s1__SV(ist)= bCHdSV(ist) & ! Factor of (eta)**(b+2)
	1043	& psidSV(ist) Ks_dSV(ist) & ! in DR97, Eqn.(3.36)
	1044	& /(etadSV(ist)**( bCHdSV(ist)+3.)) !
	1045	s2__SV(ist)= Ks_dSV(ist) & ! Factor of (eta)**(2b+3)
	1046	& /(etadSV(ist)*(2.bCHdSV(ist)+3.)) ! in DR97, Eqn.(3.35)
	1047
	1048	!C +--Soil Minimum Humidity (from a prescribed minimum relative Humidity)
	1049	!C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	1050	Psimax = -(log(RHsMin))/7.2E-5 ! DR97, Eqn 3.15 Inversion
	1051	etamSV(ist) = etadSV(ist) &
	1052	& (PsiMax/psidSV(ist))*(-min(10.,1./bCHdSV(ist)))
	1053	END DO
	1054	etamSV(12) = 0.
	1055
	1056	!C +--Piecewise Hydraulic Conductivity Profiles
	1057	!C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	1058	DO ist=0,nsot
	1059	!c #WR write(6, format(' Type \| etaSat \| No \| eta__1 \| eta__2 \|', &
	1060	!c #WR& ' Khyd_1 \| Khyd_x \| Khyd_2 \| Khyd_y \|' &
	1061	!c #WR& /,' -----+-----------+----+-----------+-----------+', &
	1062	!c #WR& '-----------+-----------+-----------+-----------+'))
	1063
	1064	d__eta = etadSV(ist)/nkhy
	1065	eta__1 = 0.
	1066	eta__2 = d__eta
	1067	DO ikh=0,nkhy
	1068	Khyd_1 = s2__SV(ist) & ! DR97, Eqn.(3.35)
	1069	& (eta__1 (2. bCHdSV(ist)+3.)) !
	1070	Khyd_2 = s2__SV(ist) &!
	1071	& (eta__2 (2. bCHdSV(ist)+3.)) !
	1072
	1073	a_Khyd = (Khyd_2-Khyd_1)/d__eta !
	1074	b_Khyd = Khyd_1-a_Khyd *eta__1 !
	1075	!c #WR Khyd_x = a_Khyd*eta__1 +b_Khyd !
	1076	!c #WR Khyd_y = a_Khyd*eta__2 +b_Khyd !
	1077	aKdtSV(ist,ikh) = a_Khyd * dt__SV !
	1078	bKdtSV(ist,ikh) = b_Khyd * dt__SV !
	1079	!c #WR write(6,format(i5,' \|',e10.2,' \|',i3,' \|', 6(e10.2,' \|')) &
	1080	!c #WR& ist,etadSV(ist),ikh,eta__1, &
	1081	!c #WR& eta__2,Khyd_1,Khyd_x,Khyd_2,Khyd_y
	1082
	1083	eta__1 = eta__1 + d__eta
	1084	eta__2 = eta__2 + d__eta
	1085	END DO
	1086	END DO
	1087
	1088
	1089	!c
	1090	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	1091	! + INITIALISATION: ARRAYS on grid +++
	1092	! + -----------------------------------
	1093	! +
	1094	DO ikl=1,knon
	1095
	1096	! + Water Fluxes, roughness length
	1097	us__SV(ikl) = 0.25 ! Frict. Velocity
	1098	! (approx ETH-Camp Lefebre etal CD 2005
	1099	uts_SV(ikl) = 1. ! uT arbitrary
	1100	uqs_SV(ikl) = 0.25 ! turb_vel(ikl,3) ! uq "
	1101	uss_SV(ikl) = 0.25 ! turb_vel(ikl,4) ! us "
	1102	!!c #AE usthSV(ikl) = 0.
	1103	Z0h_SV(ikl) = 0. ! rlength(ikl,9)
	1104	Z0m_SV(ikl) = 0. ! Moment.Roughn.L.
	1105	Z0h_SV(ikl) = 0. ! Heat Roughn.L.
	1106	!!c #OR Z0roSV(ikl) = 0. ! Orogr. Roughn.L.
	1107	LMO_SV(ikl) = 0.
	1108	VVs_SV(ikl) = 0.
	1109	RRs_SV(ikl) = 0.
	1110	DDs_SV(ikl) = 0.
	1111
	1112	! + Snow
	1113	isnoSV(ikl) = 0 ! # snow layers
	1114	BufsSV(ikl) = 0.
	1115	zzsnsv(ikl,:) = 0. ! Snow pack thickness
	1116	qsnoSV(ikl) = 0. ! BL snow content
	1117	zWEcSV(ikl) = 0.
	1118	dbs_SV(ikl) = 0.
	1119	dsnbSV(ikl) = 0.
	1120	esnbSV(ikl) = 0.
	1121	BrosSV(ikl) = 0.
	1122	BG1sSV(ikl) = 0.
	1123	BG2sSV(ikl) = 0.
	1124	SWS_SV(ikl) = 0.
	1125	RnofSV(ikl) = 0. ! RunOFF Intensity
	1126	toicSV(ikl) = 0. !hj preli
	1127	! + Clouds !hj1602
	1128	cld_SV(ikl) = 0.
	1129
	1130	! + Vegetation
	1131	LSmask(ikl) = 1 ! Land/Sea Mask
	1132	isotSV(ikl) = 12 ! Soil Type
	1133	iWaFSV(ikl) = 1 ! Soil Drainage
	1134	Rootsv(ikl,:) = 0.
	1135	ivgtSV(ikl) = 0 ! Vegetation Type
	1136	LAI0SV(ikl) = 0 ! LAI
	1137	glf0SV(ikl) = 0 ! Green Leaf Frac.
	1138	! TVegSV(ikl) = 0.
	1139	SnCaSV(ikl) = 0.
	1140	rrCaSV(ikl) = 0.
	1141	pSivSV(ikl) = 0.
	1142
	1143	! + z0(Sastrugi)
	1144	! #SZ Z0SaSL(ikl) = 0.
	1145	! #ZM DO nt=1,ntavSL
	1146	! #ZM SLn_z0(ikl,nt) = 0.5e-6
	1147	! #ZM SLn_b0(ikl,nt) = 0.5e-6
	1148	! #ZM SLn_r0(ikl,nt) = 0.5e-6
	1149	! #ZM END DO
	1150
	1151
	1152	! +--z0(Orography Roughness)
	1153	! #OR SL_z0 (ikl) = min(SL_z0 (ikl),zsigma/3.)
	1154	! #OR SLzoro(ikl) = min(SLzoro(ikl),zsigma/3.)
	1155
	1156	! +--SBL Characteristics
	1157	! #AA DO nt= 1,ntaver
	1158	! #AW V_0aSL(ikl,nt) = ssvSL(ikl)
	1159	! #AH dT0aSL(ikl,nt) = temp_air(ikl)-tsurf(ikl)
	1160	! #AA END DO
	1161
	1162
	1163	END DO
	1164
	1165
	1166
	1167	return
	1168
	1169	END SUBROUTINE SISVAT_ini
	1170
	1171
	1172
	1173
	1174
	1175
	1176
	1177	!***************************************************************************
	1178
	1179	SUBROUTINE sisvatetat0 (fichnom,ikl2i)
	1180
	1181	USE dimphy
	1182	USE mod_grid_phy_lmdz
	1183	USE mod_phys_lmdz_para
	1184
	1185	USE iostart
	1186	USE VAR_SV
	1187	USE VARdSV
	1188	USE VARxSV
	1189	USE VARtSV
	1190	USE control_mod
	1191	USE indice_sol_mod
	1192
	1193	IMPLICIT none
	1194	!======================================================================
	1195	! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
	1196	! Objet: Ecriture de l'etat de redemarrage pour SISVAT
	1197	!======================================================================
	1198	include "netcdf.inc"
	1199	! include "indicesol.h"
	1200
	1201	include "dimsoil.h"
	1202	include "clesphys.h"
	1203	include "temps.h"
	1204	include "thermcell.h"
	1205	include "compbl.h"
	1206
	1207	!======================================================================
	1208	CHARACTER(LEN=*) :: fichnom
	1209
	1210
	1211	INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
	1212	REAL, DIMENSION(klon) :: rlon
	1213	REAL, DIMENSION(klon) :: rlat
	1214
	1215	! les variables globales ecrites dans le fichier restart
	1216	REAL, DIMENSION(klon) :: isno
	1217	REAL, DIMENSION(klon) :: ispi
	1218	REAL, DIMENSION(klon) :: iice
	1219	REAL, DIMENSION(klon) :: rusn
	1220	REAL, DIMENSION(klon, nsno) :: isto
	1221
	1222	REAL, DIMENSION(klon, nsismx) :: Tsis
	1223	REAL, DIMENSION(klon, nsismx) :: eta
	1224	REAL, DIMENSION(klon, nsismx) :: ro
	1225
	1226	REAL, DIMENSION(klon, nsno) :: dzsn
	1227	REAL, DIMENSION(klon, nsno) :: G1sn
	1228	REAL, DIMENSION(klon, nsno) :: G2sn
	1229	REAL, DIMENSION(klon, nsno) :: agsn
	1230
	1231	REAL, DIMENSION(klon) :: toic
	1232
	1233	! REAL, DIMENSION(klon) :: IRs
	1234	! REAL, DIMENSION(klon) :: LMO
	1235	! REAL, DIMENSION(klon) :: Bufs
	1236	! REAL, DIMENSION(klon, 9) :: rlength
	1237	! REAL, DIMENSION(klon, 5) :: turb_vel
	1238
	1239	INTEGER :: isl, ikl, i, isn , errT, erreta, errro, errdz, snopts
	1240	CHARACTER (len=2) :: str2
	1241	LOGICAL :: found
	1242
	1243	errT=0
	1244	errro=0
	1245	erreta=0
	1246	errdz=0
	1247	snopts=0
	1248	! Ouvrir le fichier contenant l'etat initial:
	1249
	1250	CALL open_startphy(fichnom)
	1251
	1252	! Lecture des latitudes, longitudes (coordonnees):
	1253
	1254	CALL get_field("latitude",rlat,found)
	1255	CALL get_field("longitude",rlon,found)
	1256
	1257	CALL get_field("n_snows", isno,found)
	1258	IF (.NOT. found) THEN
	1259	PRINT*, 'phyetat0: Le champ <n_snows> est absent'
	1260	PRINT *, 'fichier startsisvat non compatible avec sisvatetat0'
	1261	ENDIF
	1262
	1263	CALL get_field("n_ice_top",ispi,found)
	1264	CALL get_field("n_ice",iice,found)
	1265	CALL get_field("surf_water",rusn,found)
	1266	! IF (.NOT. found) THEN
	1267	! PRINT*, 'phyetat0: Le champ <surf_water> est absent'
	1268	! rusn(:)=0.
	1269	! ENDIF
	1270
	1271
	1272	CALL get_field("to_ice",toic,found)
	1273	IF (.NOT. found) THEN
	1274	PRINT*, 'phyetat0: Le champ <to_ice> est absent'
	1275	toic(:)=0.
	1276	ENDIF
	1277
	1278
	1279	! CALL get_field("IR_soil",IRs,found)
	1280	! CALL get_field("LMO",LMO,found)
	1281	! CALL get_field("snow_buffer",Bufs,found)
	1282	! DO i = 1, 5
	1283	! WRITE(str2,'(i2.2)') i
	1284	! CALL get_field("turb_veloc"//str2, &
	1285	! turb_vel(:,i),found)
	1286	! ENDDO
	1287	! DO i = 1, 9
	1288	! WRITE(str2,'(i2.2)') i
	1289	! CALL get_field("rough_length"//str2, &
	1290	! rlength(:,i),found)
	1291	! ENDDO
	1292	DO isn = 1,nsno
	1293	IF (isn.LE.99) THEN
	1294	WRITE(str2,'(i2.2)') isn
	1295	CALL get_field("AGESNOW"//str2, &
	1296	agsn(:,isn),found)
	1297	ELSE
	1298	PRINT*, "Trop de couches"
	1299	CALL abort
	1300	ENDIF
	1301	ENDDO
	1302	DO isn = 1,nsno
	1303	IF (isn.LE.99) THEN
	1304	WRITE(str2,'(i2.2)') isn
	1305	CALL get_field("DZSNOW"//str2, &
	1306	dzsn(:,isn),found)
	1307	ELSE
	1308	PRINT*, "Trop de couches"
	1309	CALL abort
	1310	ENDIF
	1311	ENDDO
	1312	DO isn = 1,nsno
	1313	IF (isn.LE.99) THEN
	1314	WRITE(str2,'(i2.2)') isn
	1315	CALL get_field("G2SNOW"//str2, &
	1316	G2sn(:,isn),found)
	1317	ELSE
	1318	PRINT*, "Trop de couches"
	1319	CALL abort
	1320	ENDIF
	1321	ENDDO
	1322	DO isn = 1,nsno
	1323	IF (isn.LE.99) THEN
	1324	WRITE(str2,'(i2.2)') isn
	1325	CALL get_field("G1SNOW"//str2, &
	1326	G1sn(:,isn),found)
	1327	ELSE
	1328	PRINT*, "Trop de couches"
	1329	CALL abort
	1330	ENDIF
	1331	ENDDO
	1332	DO isn = 1,nsismx
	1333	IF (isn.LE.99) THEN
	1334	WRITE(str2,'(i2.2)') isn
	1335	CALL get_field("ETA"//str2, &
	1336	eta(:,isn),found)
	1337	ELSE
	1338	PRINT*, "Trop de couches"
	1339	CALL abort
	1340	ENDIF
	1341	ENDDO
	1342	DO isn = 1,nsismx
	1343	IF (isn.LE.99) THEN
	1344	WRITE(str2,'(i2.2)') isn
	1345	CALL get_field("RO"//str2, &
	1346	ro(:,isn),found)
	1347	ELSE
	1348	PRINT*, "Trop de couches"
	1349	CALL abort
	1350	ENDIF
	1351	ENDDO
	1352	DO isn = 1,nsismx
	1353	IF (isn.LE.99) THEN
	1354	WRITE(str2,'(i2.2)') isn
	1355	CALL get_field("TSS"//str2, &
	1356	Tsis(:,isn),found)
	1357	ELSE
	1358	PRINT*, "Trop de couches"
	1359	CALL abort
	1360	ENDIF
	1361	ENDDO
	1362	DO isn = 1,nsno
	1363	IF (isn.LE.99) THEN
	1364	WRITE(str2,'(i2.2)') isn
	1365	CALL get_field("HISTORY"//str2, &
	1366	isto(:,isn),found)
	1367	ELSE
	1368	PRINT*, "Trop de couches"
	1369	CALL abort
	1370	ENDIF
	1371	ENDDO
	1372	write(,)'Read ',fichnom,' finished!!'
	1373
	1374	!*********************************************************************************
	1375	! Compress restart file variables for SISVAT
	1376
	1377
	1378	isnoSV(:) = 0
	1379	ispiSV(:) = 0
	1380	iiceSV(:) = 0
	1381
	1382	eta_SV(:,1:nsno) = 0.
	1383	TsisSV(:,1:nsno) = 0.
	1384	istoSV(:,:) = 0
	1385	ro__SV(:,1:nsno) = 0.
	1386	G1snSV(:,:) = 0.
	1387	G2snSV(:,:) = 0.
	1388	dzsnSV(:,1:nsno) = 0.
	1389	agsnSV(:,:) = 0.
	1390	rusnSV(:) = 0.
	1391	toicSV(:) = 0.
	1392
	1393	DO ikl = 1,klon
	1394	i = ikl2i(ikl)
	1395	IF (i > 0) THEN
	1396	isnoSV(ikl) = INT(isno(i)) ! Nb Snow/Ice Lay.
	1397	ispiSV(ikl) = INT(ispi(i)) ! Nb Supr.Ice Lay.
	1398	iiceSV(ikl) = INT(iice(i)) ! Nb Ice Lay.
	1399
	1400	!! IRs_SV(ikl) = IRs(i)
	1401	! LMO_SV(ikl) = LMO(i) !?
	1402	! us__SV(ikl) = turb_vel(i,1)
	1403	! uts_SV(ikl) = turb_vel(i,2)
	1404	! uqs_SV(ikl) = turb_vel(i,3)
	1405	! uss_SV(ikl) = turb_vel(i,4)
	1406	! usthSV(ikl) = turb_vel(i,5)
	1407	! Z0m_SV(ikl) = rlength(i,1) ! Moment.Roughn.L.
	1408	! Z0mmSV(ikl) = rlength(i,2) ! Moment.Roughn.L.
	1409	! Z0mnSV(ikl) = rlength(i,3) ! Moment.Roughn.L.
	1410	! Z0SaSV(ikl) = rlength(i,4) ! Moment.Roughn.L.
	1411	! Z0e_SV(ikl) = rlength(i,5) ! Moment.Roughn.L.
	1412	! Z0emSV(ikl) = rlength(i,6) ! Moment.Roughn.L.
	1413	! Z0enSV(ikl) = rlength(i,7) ! Moment.Roughn.L.
	1414	! Z0roSV(ikl) = rlength(i,8) ! Moment.Roughn.L.
	1415	!! Z0h_SV(ikl) = rlength(i,9) ! Moment.Roughn.L.
	1416
	1417	DO isl = -nsol,0
	1418	ro__SV(ikl,isl) = ro(i,nsno+1-isl) !
	1419	eta_SV(ikl,isl) = eta(i,nsno+1-isl) ! Soil Humidity
	1420	!hjp 15/10/2010
	1421	IF (eta_SV(ikl,isl) <= 1.e-6) THEN !hj check
	1422	eta_SV(ikl,isl) = 1.e-6
	1423	ENDIF
	1424	TsisSV(ikl,isl) = Tsis(i,nsno+1-isl) ! Soil Temperature
	1425	IF (TsisSV(ikl,isl) <= 1.) THEN !hj check
	1426	! errT=errT+1
	1427	TsisSV(ikl,isl) = 273.15
	1428	ENDIF
	1429
	1430	END DO
	1431	write(,)'Copy histo', ikl
	1432
	1433
	1434	istoSV(ikl,0) = 0 ! Snow History
	1435
	1436	G1snSV(ikl,0) = 0.5 ! [-]
	1437	G2snSV(ikl,0) = 3. ! [-] [0.0001 m]
	1438	dzsnSV(ikl,0) = dz_dSV(0) ! [m]
	1439	agsnSV(ikl,0) = 0. ! [day]
	1440
	1441	DO isn = 1,isnoSV(ikl) !nsno
	1442	snopts=snopts+1
	1443	IF (isto(i,isn) > 10.) THEN !hj check
	1444	write(,)'Irregular isto',ikl,i,isn,isto(i,isn)
	1445	isto(i,isn) = 1.
	1446	ENDIF
	1447
	1448	istoSV(ikl,isn) = INT(isto(i,isn)) ! Snow History
	1449	ro__SV(ikl,isn) = ro(i,isn) ! [kg/m3]
	1450	eta_SV(ikl,isn) = eta(i,isn) ! [m3/m3]
	1451	TsisSV(ikl,isn) = Tsis(i,isn) ! [K]
	1452
	1453	IF (TsisSV(ikl,isn) <= 1.) THEN !hj check
	1454	errT=errT+1
	1455	TsisSV(ikl,isn) = TsisSV(ikl,0)
	1456	ENDIF
	1457	IF (TsisSV(ikl,isn) <= 1.) THEN !hj check
	1458	TsisSV(ikl,isn) = 263.15
	1459	ENDIF
	1460	IF (eta_SV(ikl,isn) < 1.e-9) THEN !hj check
	1461	eta_SV(ikl,isn) = 1.e-6
	1462	erreta=erreta+1
	1463	ENDIF
	1464	IF (ro__SV(ikl,isn) <= 10.) THEN !hj check
	1465	ro__SV(ikl,isn) = 11.
	1466	errro=errro+1
	1467	ENDIF
	1468	write(,)ikl,i,isn,Tsis(i,isn),G1sn(i,isn)
	1469	G1snSV(ikl,isn) = G1sn(i,isn) ! [-] [-]
	1470	G2snSV(ikl,isn) = G2sn(i,isn) ! [-] [0.0001 m]
	1471	dzsnSV(ikl,isn) = dzsn(i,isn) ! [m]
	1472	! IF (dzsnSV(ikl,isn) < 5.e-6) THEN !hj check
	1473	! dzsnSV(ikl,isn) = 0.000005
	1474	! errdz=errdz+1
	1475	! ENDIF
	1476	agsnSV(ikl,isn) = agsn(i,isn) ! [day]
	1477	END DO
	1478	rusnSV(ikl) = rusn(i) ! Surficial Water
	1479	toicSV(ikl) = toic(i) ! bilan snow to ice
	1480	!! BufsSV(ikl) = Bufs(i)
	1481	END IF
	1482	END DO
	1483	! write(,)snopts,'snow pts',errT,' T errors,',errro,'ro,'
	1484	! write(,)' ',erreta,'eta,',errdz,'dz'
	1485	END SUBROUTINE sisvatetat0
	1486
	1487
	1488
	1489
	1490	SUBROUTINE sisvatredem (fichnom,ikl2i,rlon,rlat)
	1491	!======================================================================
	1492	! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
	1493	! Objet: Ecriture de l'etat de redemarrage pour SISVAT
	1494	!======================================================================
	1495	USE mod_grid_phy_lmdz
	1496	USE mod_phys_lmdz_para
	1497	USE iostart
	1498	USE VAR_SV
	1499	USE VARxSV
	1500	USE VARySV !hj tmp 12 03 2010
	1501	USE VARtSV
	1502	USE control_mod
	1503	USE indice_sol_mod
	1504
	1505	IMPLICIT none
	1506
	1507	include "netcdf.inc"
	1508	! include "indicesol.h"
	1509	include "dimsoil.h"
	1510	include "clesphys.h"
	1511	include "temps.h"
	1512	include "thermcell.h"
	1513	include "compbl.h"
	1514
	1515	!======================================================================
	1516
	1517	CHARACTER(LEN=*) :: fichnom
	1518	INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
	1519	REAL, DIMENSION(klon), INTENT(IN) :: rlon
	1520	REAL, DIMENSION(klon), INTENT(IN) :: rlat
	1521
	1522	! les variables globales ecrites dans le fichier restart
	1523	REAL, DIMENSION(klon) :: isno
	1524	REAL, DIMENSION(klon) :: ispi
	1525	REAL, DIMENSION(klon) :: iice
	1526	REAL, DIMENSION(klon, nsnowmx) :: isto
	1527
	1528	REAL, DIMENSION(klon, nsismx) :: Tsis
	1529	REAL, DIMENSION(klon, nsismx) :: eta
	1530	REAL, DIMENSION(klon, nsnowmx) :: dzsn
	1531	REAL, DIMENSION(klon, nsismx) :: ro
	1532	REAL, DIMENSION(klon, nsnowmx) :: G1sn
	1533	REAL, DIMENSION(klon, nsnowmx) :: G2sn
	1534	REAL, DIMENSION(klon, nsnowmx) :: agsn
	1535	REAL, DIMENSION(klon) :: IRs
	1536	REAL, DIMENSION(klon) :: LMO
	1537	REAL, DIMENSION(klon) :: rusn
	1538	REAL, DIMENSION(klon) :: toic
	1539	REAL, DIMENSION(klon) :: Bufs
	1540	REAL, DIMENSION(klon) :: alb1,alb2,alb3
	1541	REAL, DIMENSION(klon, 9) :: rlength
	1542	REAL, DIMENSION(klon, 5) :: turb_vel
	1543
	1544	INTEGER isl, ikl, i, isn
	1545	CHARACTER (len=2) :: str2
	1546
	1547	isno(:) = 0
	1548	ispi(:) = 0
	1549	iice(:) = 0
	1550	IRs(:) = 0.
	1551	LMO(:) = 0.
	1552	turb_vel(:,:) = 0.
	1553	rlength(:,:) = 0.
	1554	eta(:,:) = 0.
	1555	Tsis(:,:) = 0.
	1556	isto(:,:) = 0
	1557	ro(:,:) = 0.
	1558	G1sn(:,:) = 0.
	1559	G2sn(:,:) = 0.
	1560	dzsn(:,:) = 0.
	1561	agsn(:,:) = 0.
	1562	rusn(:) = 0.
	1563	toic(:) = 0.
	1564	Bufs(:) = 0.
	1565	alb1(:) = 0.
	1566	alb2(:) = 0.
	1567	alb3(:) = 0.
	1568
	1569	!***************************************************************************
	1570	! Uncompress SISVAT output variables for storage
	1571
	1572	DO ikl = 1,klon
	1573	i = ikl2i(ikl)
	1574	IF (i > 0) THEN
	1575	isno(i) = 1.*isnoSV(ikl) ! Nb Snow/Ice Lay.
	1576	ispi(i) = 1.*ispiSV(ikl) ! Nb Supr.Ice Lay.
	1577	iice(i) = 1.*iiceSV(ikl) ! Nb Ice Lay.
	1578
	1579	! IRs(i) = IRs_SV(ikl)
	1580	! LMO(i) = LMO_SV(ikl)
	1581	! turb_vel(i,1) = us__SV(ikl)
	1582	! turb_vel(i,2) = uts_SV(ikl)
	1583	! turb_vel(i,3) = uqs_SV(ikl)
	1584	! turb_vel(i,4) = uss_SV(ikl)
	1585	! turb_vel(i,5) = usthSV(ikl)
	1586	! rlength(i,1) = Z0m_SV(ikl) ! Moment.Roughn.L.
	1587	! rlength(i,2) = Z0mmSV(ikl) !
	1588	! rlength(i,3) = Z0mnSV(ikl) !
	1589	! rlength(i,4) = Z0SaSV(ikl) !
	1590	! rlength(i,5) = Z0e_SV(ikl) !
	1591	! rlength(i,6) = Z0emSV(ikl) !
	1592	! rlength(i,7) = Z0enSV(ikl) !
	1593	! rlength(i,8) = Z0roSV(ikl) !
	1594	! rlength(i,9) = Z0h_SV(ikl) !
	1595
	1596	DO isl = -nsol,0 !
	1597	eta(i,nsno+1-isl) = eta_SV(ikl,isl) ! Soil Humidity
	1598	Tsis(i,nsno+1-isl) = TsisSV(ikl,isl) ! Soil Temperature
	1599	ro(i,nsno+1-isl) = ro__SV(ikl,isl) ! [kg/m3]
	1600	END DO
	1601
	1602
	1603	DO isn = 1,nsno
	1604	isto(i,isn) = 1.*istoSV(ikl,isn) ! Snow History
	1605	ro(i,isn) = ro__SV(ikl,isn) ! [kg/m3]
	1606	eta(i,isn) = eta_SV(ikl,isn) ! [m3/m3]
	1607	Tsis(i,isn) = TsisSV(ikl,isn) ! [K]
	1608	G1sn(i,isn) = G1snSV(ikl,isn) ! [-] [-]
	1609	G2sn(i,isn) = G2snSV(ikl,isn) ! [-] [0.0001 m]
	1610	dzsn(i,isn) = dzsnSV(ikl,isn) ! [m]
	1611	agsn(i,isn) = agsnSV(ikl,isn) ! [day]
	1612	END DO
	1613	rusn(i) = rusnSV(ikl) ! Surficial Water
	1614	toic(i) = toicSV(ikl) ! Surficial Water
	1615	alb1(i) = alb1sv(ikl)
	1616	alb2(i) = alb2sv(ikl)
	1617	alb3(i) = alb3sv(ikl)
	1618	! Bufs(i) = BufsSV(ikl)
	1619	END IF
	1620	END DO
	1621
	1622	CALL open_restartphy(fichnom)
	1623	CALL put_field("longitude", &
	1624	"Longitudes de la grille physique",rlon)
	1625	CALL put_field("latitude","Latitudes de la grille physique",rlat)
	1626
	1627	CALL put_field("n_snows", "number of snow/ice layers",isno)
	1628	CALL put_field("n_ice_top", "number of top ice layers",ispi)
	1629	CALL put_field("n_ice", "number of ice layers",iice)
	1630	CALL put_field("IR_soil", "Soil IR flux",IRs)
	1631	CALL put_field("LMO", "Monin-Obukhov Scale",LMO)
	1632	CALL put_field("surf_water", "Surficial water",rusn)
	1633	CALL put_field("snow_buffer", "Snow buffer layer",Bufs)
	1634	CALL put_field("alb_1", "albedo sw",alb1)
	1635	CALL put_field("alb_2", "albedo nIR",alb2)
	1636	CALL put_field("alb_3", "albedo fIR",alb3)
	1637	CALL put_field("to_ice", "Snow passed to ice",toic)
	1638
	1639	! DO i = 1, 5
	1640	! WRITE(str2,'(i2.2)') i
	1641	! CALL put_field("turb_veloc"//str2, &
	1642	! "various turbulent velocities"//str2, &
	1643	! turb_vel(:,i))
	1644	! ENDDO
	1645	! DO i = 1, 9
	1646	! WRITE(str2,'(i2.2)') i
	1647	! CALL put_field("rough_length"//str2, &
	1648	! "various roughness lengths"//str2, &
	1649	! rlength(:,i))
	1650	! ENDDO
	1651	DO isn = 1,nsno
	1652	IF (isn.LE.99) THEN
	1653	WRITE(str2,'(i2.2)') isn
	1654	CALL put_field("AGESNOW"//str2, &
	1655	"Age de la neige layer No."//str2, &
	1656	agsn(:,isn))
	1657	ELSE
	1658	PRINT*, "Trop de couches"
	1659	CALL abort
	1660	ENDIF
	1661	ENDDO
	1662	DO isn = 1,nsno
	1663	IF (isn.LE.99) THEN
	1664	WRITE(str2,'(i2.2)') isn
	1665	CALL put_field("DZSNOW"//str2, &
	1666	"Snow/ice thickness layer No."//str2, &
	1667	dzsn(:,isn))
	1668	ELSE
	1669	PRINT*, "Trop de couches"
	1670	CALL abort
	1671	ENDIF
	1672	ENDDO
	1673	DO isn = 1,nsno
	1674	IF (isn.LE.99) THEN
	1675	WRITE(str2,'(i2.2)') isn
	1676	CALL put_field("G2SNOW"//str2, &
	1677	"Snow Property 2, layer No."//str2, &
	1678	G2sn(:,isn))
	1679	ELSE
	1680	PRINT*, "Trop de couches"
	1681	CALL abort
	1682	ENDIF
	1683	ENDDO
	1684	DO isn = 1,nsno
	1685	IF (isn.LE.99) THEN
	1686	WRITE(str2,'(i2.2)') isn
	1687	CALL put_field("G1SNOW"//str2, &
	1688	"Snow Property 1, layer No."//str2, &
	1689	G1sn(:,isn))
	1690	ELSE
	1691	PRINT*, "Trop de couches"
	1692	CALL abort
	1693	ENDIF
	1694	ENDDO
	1695	DO isn = 1,nsismx
	1696	IF (isn.LE.99) THEN
	1697	WRITE(str2,'(i2.2)') isn
	1698	CALL put_field("ETA"//str2, &
	1699	"Soil/snow water content layer No."//str2, &
	1700	eta(:,isn))
	1701	ELSE
	1702	PRINT*, "Trop de couches"
	1703	CALL abort
	1704	ENDIF
	1705	ENDDO
	1706	DO isn = 1,nsismx !nsno
	1707	IF (isn.LE.99) THEN
	1708	WRITE(str2,'(i2.2)') isn
	1709	CALL put_field("RO"//str2, &
	1710	"Snow density layer No."//str2, &
	1711	ro(:,isn))
	1712	ELSE
	1713	PRINT*, "Trop de couches"
	1714	CALL abort
	1715	ENDIF
	1716	ENDDO
	1717	DO isn = 1,nsismx
	1718	IF (isn.LE.99) THEN
	1719	WRITE(str2,'(i2.2)') isn
	1720	CALL put_field("TSS"//str2, &
	1721	"Soil/snow temperature layer No."//str2, &
	1722	Tsis(:,isn))
	1723	ELSE
	1724	PRINT*, "Trop de couches"
	1725	CALL abort
	1726	ENDIF
	1727	ENDDO
	1728	DO isn = 1,nsno
	1729	IF (isn.LE.99) THEN
	1730	WRITE(str2,'(i2.2)') isn
	1731	CALL put_field("HISTORY"//str2, &
	1732	"Snow history layer No."//str2, &
	1733	isto(:,isn))
	1734	ELSE
	1735	PRINT*, "Trop de couches"
	1736	CALL abort
	1737	ENDIF
	1738	ENDDO
	1739
	1740	END SUBROUTINE sisvatredem
	1741
	1742	END MODULE surf_sisvat_mod

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

Download in other formats:

Original Format