Context Navigation

surf_sisvat_mod.F90 @ 4799

Last change on this file since 4799 was 3518, checked in by Laurent Fairhead, 5 years ago
Modifications needed as put_field call has changed
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.7 KB

Line
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_SISVAT
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 indice_sol_mod
1191
1192	IMPLICIT none
1193	!======================================================================
1194	! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
1195	! Objet: Ecriture de l'etat de redemarrage pour SISVAT
1196	!======================================================================
1197	include "netcdf.inc"
1198	! include "indicesol.h"
1199
1200	include "dimsoil.h"
1201	include "clesphys.h"
1202	include "thermcell.h"
1203	include "compbl.h"
1204
1205	!======================================================================
1206	CHARACTER(LEN=*) :: fichnom
1207
1208
1209	INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
1210	REAL, DIMENSION(klon) :: rlon
1211	REAL, DIMENSION(klon) :: rlat
1212
1213	! les variables globales ecrites dans le fichier restart
1214	REAL, DIMENSION(klon) :: isno
1215	REAL, DIMENSION(klon) :: ispi
1216	REAL, DIMENSION(klon) :: iice
1217	REAL, DIMENSION(klon) :: rusn
1218	REAL, DIMENSION(klon, nsno) :: isto
1219
1220	REAL, DIMENSION(klon, nsismx) :: Tsis
1221	REAL, DIMENSION(klon, nsismx) :: eta
1222	REAL, DIMENSION(klon, nsismx) :: ro
1223
1224	REAL, DIMENSION(klon, nsno) :: dzsn
1225	REAL, DIMENSION(klon, nsno) :: G1sn
1226	REAL, DIMENSION(klon, nsno) :: G2sn
1227	REAL, DIMENSION(klon, nsno) :: agsn
1228
1229	REAL, DIMENSION(klon) :: toic
1230
1231	! REAL, DIMENSION(klon) :: IRs
1232	! REAL, DIMENSION(klon) :: LMO
1233	! REAL, DIMENSION(klon) :: Bufs
1234	! REAL, DIMENSION(klon, 9) :: rlength
1235	! REAL, DIMENSION(klon, 5) :: turb_vel
1236
1237	INTEGER :: isl, ikl, i, isn , errT, erreta, errro, errdz, snopts
1238	CHARACTER (len=2) :: str2
1239	LOGICAL :: found
1240
1241	errT=0
1242	errro=0
1243	erreta=0
1244	errdz=0
1245	snopts=0
1246	! Ouvrir le fichier contenant l'etat initial:
1247
1248	CALL open_startphy(fichnom)
1249
1250	! Lecture des latitudes, longitudes (coordonnees):
1251
1252	CALL get_field("latitude",rlat,found)
1253	CALL get_field("longitude",rlon,found)
1254
1255	CALL get_field("n_snows", isno,found)
1256	IF (.NOT. found) THEN
1257	PRINT*, 'phyetat0: Le champ <n_snows> est absent'
1258	PRINT *, 'fichier startsisvat non compatible avec sisvatetat0'
1259	ENDIF
1260
1261	CALL get_field("n_ice_top",ispi,found)
1262	CALL get_field("n_ice",iice,found)
1263	CALL get_field("surf_water",rusn,found)
1264	! IF (.NOT. found) THEN
1265	! PRINT*, 'phyetat0: Le champ <surf_water> est absent'
1266	! rusn(:)=0.
1267	! ENDIF
1268
1269
1270	CALL get_field("to_ice",toic,found)
1271	IF (.NOT. found) THEN
1272	PRINT*, 'phyetat0: Le champ <to_ice> est absent'
1273	toic(:)=0.
1274	ENDIF
1275
1276
1277	! CALL get_field("IR_soil",IRs,found)
1278	! CALL get_field("LMO",LMO,found)
1279	! CALL get_field("snow_buffer",Bufs,found)
1280	! DO i = 1, 5
1281	! WRITE(str2,'(i2.2)') i
1282	! CALL get_field("turb_veloc"//str2, &
1283	! turb_vel(:,i),found)
1284	! ENDDO
1285	! DO i = 1, 9
1286	! WRITE(str2,'(i2.2)') i
1287	! CALL get_field("rough_length"//str2, &
1288	! rlength(:,i),found)
1289	! ENDDO
1290	DO isn = 1,nsno
1291	IF (isn.LE.99) THEN
1292	WRITE(str2,'(i2.2)') isn
1293	CALL get_field("AGESNOW"//str2, &
1294	agsn(:,isn),found)
1295	ELSE
1296	PRINT*, "Trop de couches"
1297	CALL abort
1298	ENDIF
1299	ENDDO
1300	DO isn = 1,nsno
1301	IF (isn.LE.99) THEN
1302	WRITE(str2,'(i2.2)') isn
1303	CALL get_field("DZSNOW"//str2, &
1304	dzsn(:,isn),found)
1305	ELSE
1306	PRINT*, "Trop de couches"
1307	CALL abort
1308	ENDIF
1309	ENDDO
1310	DO isn = 1,nsno
1311	IF (isn.LE.99) THEN
1312	WRITE(str2,'(i2.2)') isn
1313	CALL get_field("G2SNOW"//str2, &
1314	G2sn(:,isn),found)
1315	ELSE
1316	PRINT*, "Trop de couches"
1317	CALL abort
1318	ENDIF
1319	ENDDO
1320	DO isn = 1,nsno
1321	IF (isn.LE.99) THEN
1322	WRITE(str2,'(i2.2)') isn
1323	CALL get_field("G1SNOW"//str2, &
1324	G1sn(:,isn),found)
1325	ELSE
1326	PRINT*, "Trop de couches"
1327	CALL abort
1328	ENDIF
1329	ENDDO
1330	DO isn = 1,nsismx
1331	IF (isn.LE.99) THEN
1332	WRITE(str2,'(i2.2)') isn
1333	CALL get_field("ETA"//str2, &
1334	eta(:,isn),found)
1335	ELSE
1336	PRINT*, "Trop de couches"
1337	CALL abort
1338	ENDIF
1339	ENDDO
1340	DO isn = 1,nsismx
1341	IF (isn.LE.99) THEN
1342	WRITE(str2,'(i2.2)') isn
1343	CALL get_field("RO"//str2, &
1344	ro(:,isn),found)
1345	ELSE
1346	PRINT*, "Trop de couches"
1347	CALL abort
1348	ENDIF
1349	ENDDO
1350	DO isn = 1,nsismx
1351	IF (isn.LE.99) THEN
1352	WRITE(str2,'(i2.2)') isn
1353	CALL get_field("TSS"//str2, &
1354	Tsis(:,isn),found)
1355	ELSE
1356	PRINT*, "Trop de couches"
1357	CALL abort
1358	ENDIF
1359	ENDDO
1360	DO isn = 1,nsno
1361	IF (isn.LE.99) THEN
1362	WRITE(str2,'(i2.2)') isn
1363	CALL get_field("HISTORY"//str2, &
1364	isto(:,isn),found)
1365	ELSE
1366	PRINT*, "Trop de couches"
1367	CALL abort
1368	ENDIF
1369	ENDDO
1370	write(,)'Read ',fichnom,' finished!!'
1371
1372	!*********************************************************************************
1373	! Compress restart file variables for SISVAT
1374
1375
1376	isnoSV(:) = 0
1377	ispiSV(:) = 0
1378	iiceSV(:) = 0
1379
1380	eta_SV(:,1:nsno) = 0.
1381	TsisSV(:,1:nsno) = 0.
1382	istoSV(:,:) = 0
1383	ro__SV(:,1:nsno) = 0.
1384	G1snSV(:,:) = 0.
1385	G2snSV(:,:) = 0.
1386	dzsnSV(:,1:nsno) = 0.
1387	agsnSV(:,:) = 0.
1388	rusnSV(:) = 0.
1389	toicSV(:) = 0.
1390
1391	DO ikl = 1,klon
1392	i = ikl2i(ikl)
1393	IF (i > 0) THEN
1394	isnoSV(ikl) = INT(isno(i)) ! Nb Snow/Ice Lay.
1395	ispiSV(ikl) = INT(ispi(i)) ! Nb Supr.Ice Lay.
1396	iiceSV(ikl) = INT(iice(i)) ! Nb Ice Lay.
1397
1398	!! IRs_SV(ikl) = IRs(i)
1399	! LMO_SV(ikl) = LMO(i) !?
1400	! us__SV(ikl) = turb_vel(i,1)
1401	! uts_SV(ikl) = turb_vel(i,2)
1402	! uqs_SV(ikl) = turb_vel(i,3)
1403	! uss_SV(ikl) = turb_vel(i,4)
1404	! usthSV(ikl) = turb_vel(i,5)
1405	! Z0m_SV(ikl) = rlength(i,1) ! Moment.Roughn.L.
1406	! Z0mmSV(ikl) = rlength(i,2) ! Moment.Roughn.L.
1407	! Z0mnSV(ikl) = rlength(i,3) ! Moment.Roughn.L.
1408	! Z0SaSV(ikl) = rlength(i,4) ! Moment.Roughn.L.
1409	! Z0e_SV(ikl) = rlength(i,5) ! Moment.Roughn.L.
1410	! Z0emSV(ikl) = rlength(i,6) ! Moment.Roughn.L.
1411	! Z0enSV(ikl) = rlength(i,7) ! Moment.Roughn.L.
1412	! Z0roSV(ikl) = rlength(i,8) ! Moment.Roughn.L.
1413	!! Z0h_SV(ikl) = rlength(i,9) ! Moment.Roughn.L.
1414
1415	DO isl = -nsol,0
1416	ro__SV(ikl,isl) = ro(i,nsno+1-isl) !
1417	eta_SV(ikl,isl) = eta(i,nsno+1-isl) ! Soil Humidity
1418	!hjp 15/10/2010
1419	IF (eta_SV(ikl,isl) <= 1.e-6) THEN !hj check
1420	eta_SV(ikl,isl) = 1.e-6
1421	ENDIF
1422	TsisSV(ikl,isl) = Tsis(i,nsno+1-isl) ! Soil Temperature
1423	IF (TsisSV(ikl,isl) <= 1.) THEN !hj check
1424	! errT=errT+1
1425	TsisSV(ikl,isl) = 273.15
1426	ENDIF
1427
1428	END DO
1429	write(,)'Copy histo', ikl
1430
1431
1432	istoSV(ikl,0) = 0 ! Snow History
1433
1434	G1snSV(ikl,0) = 0.5 ! [-]
1435	G2snSV(ikl,0) = 3. ! [-] [0.0001 m]
1436	dzsnSV(ikl,0) = dz_dSV(0) ! [m]
1437	agsnSV(ikl,0) = 0. ! [day]
1438
1439	DO isn = 1,isnoSV(ikl) !nsno
1440	snopts=snopts+1
1441	IF (isto(i,isn) > 10.) THEN !hj check
1442	write(,)'Irregular isto',ikl,i,isn,isto(i,isn)
1443	isto(i,isn) = 1.
1444	ENDIF
1445
1446	istoSV(ikl,isn) = INT(isto(i,isn)) ! Snow History
1447	ro__SV(ikl,isn) = ro(i,isn) ! [kg/m3]
1448	eta_SV(ikl,isn) = eta(i,isn) ! [m3/m3]
1449	TsisSV(ikl,isn) = Tsis(i,isn) ! [K]
1450
1451	IF (TsisSV(ikl,isn) <= 1.) THEN !hj check
1452	errT=errT+1
1453	TsisSV(ikl,isn) = TsisSV(ikl,0)
1454	ENDIF
1455	IF (TsisSV(ikl,isn) <= 1.) THEN !hj check
1456	TsisSV(ikl,isn) = 263.15
1457	ENDIF
1458	IF (eta_SV(ikl,isn) < 1.e-9) THEN !hj check
1459	eta_SV(ikl,isn) = 1.e-6
1460	erreta=erreta+1
1461	ENDIF
1462	IF (ro__SV(ikl,isn) <= 10.) THEN !hj check
1463	ro__SV(ikl,isn) = 11.
1464	errro=errro+1
1465	ENDIF
1466	write(,)ikl,i,isn,Tsis(i,isn),G1sn(i,isn)
1467	G1snSV(ikl,isn) = G1sn(i,isn) ! [-] [-]
1468	G2snSV(ikl,isn) = G2sn(i,isn) ! [-] [0.0001 m]
1469	dzsnSV(ikl,isn) = dzsn(i,isn) ! [m]
1470	! IF (dzsnSV(ikl,isn) < 5.e-6) THEN !hj check
1471	! dzsnSV(ikl,isn) = 0.000005
1472	! errdz=errdz+1
1473	! ENDIF
1474	agsnSV(ikl,isn) = agsn(i,isn) ! [day]
1475	END DO
1476	rusnSV(ikl) = rusn(i) ! Surficial Water
1477	toicSV(ikl) = toic(i) ! bilan snow to ice
1478	!! BufsSV(ikl) = Bufs(i)
1479	END IF
1480	END DO
1481	! write(,)snopts,'snow pts',errT,' T errors,',errro,'ro,'
1482	! write(,)' ',erreta,'eta,',errdz,'dz'
1483	END SUBROUTINE sisvatetat0
1484
1485
1486
1487
1488	SUBROUTINE sisvatredem (fichnom,ikl2i,rlon,rlat)
1489	!======================================================================
1490	! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
1491	! Objet: Ecriture de l'etat de redemarrage pour SISVAT
1492	!======================================================================
1493	USE mod_grid_phy_lmdz
1494	USE mod_phys_lmdz_para
1495	USE iostart
1496	USE VAR_SV
1497	USE VARxSV
1498	USE VARySV !hj tmp 12 03 2010
1499	USE VARtSV
1500	USE indice_sol_mod
1501
1502	IMPLICIT none
1503
1504	include "netcdf.inc"
1505	! include "indicesol.h"
1506	include "dimsoil.h"
1507	include "clesphys.h"
1508	include "thermcell.h"
1509	include "compbl.h"
1510
1511	!======================================================================
1512
1513	CHARACTER(LEN=*) :: fichnom
1514	INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
1515	REAL, DIMENSION(klon), INTENT(IN) :: rlon
1516	REAL, DIMENSION(klon), INTENT(IN) :: rlat
1517
1518	! les variables globales ecrites dans le fichier restart
1519	REAL, DIMENSION(klon) :: isno
1520	REAL, DIMENSION(klon) :: ispi
1521	REAL, DIMENSION(klon) :: iice
1522	REAL, DIMENSION(klon, nsnowmx) :: isto
1523
1524	REAL, DIMENSION(klon, nsismx) :: Tsis
1525	REAL, DIMENSION(klon, nsismx) :: eta
1526	REAL, DIMENSION(klon, nsnowmx) :: dzsn
1527	REAL, DIMENSION(klon, nsismx) :: ro
1528	REAL, DIMENSION(klon, nsnowmx) :: G1sn
1529	REAL, DIMENSION(klon, nsnowmx) :: G2sn
1530	REAL, DIMENSION(klon, nsnowmx) :: agsn
1531	REAL, DIMENSION(klon) :: IRs
1532	REAL, DIMENSION(klon) :: LMO
1533	REAL, DIMENSION(klon) :: rusn
1534	REAL, DIMENSION(klon) :: toic
1535	REAL, DIMENSION(klon) :: Bufs
1536	REAL, DIMENSION(klon) :: alb1,alb2,alb3
1537	REAL, DIMENSION(klon, 9) :: rlength
1538	REAL, DIMENSION(klon, 5) :: turb_vel
1539
1540	INTEGER isl, ikl, i, isn
1541	CHARACTER (len=2) :: str2
1542	INTEGER :: pass
1543
1544	isno(:) = 0
1545	ispi(:) = 0
1546	iice(:) = 0
1547	IRs(:) = 0.
1548	LMO(:) = 0.
1549	turb_vel(:,:) = 0.
1550	rlength(:,:) = 0.
1551	eta(:,:) = 0.
1552	Tsis(:,:) = 0.
1553	isto(:,:) = 0
1554	ro(:,:) = 0.
1555	G1sn(:,:) = 0.
1556	G2sn(:,:) = 0.
1557	dzsn(:,:) = 0.
1558	agsn(:,:) = 0.
1559	rusn(:) = 0.
1560	toic(:) = 0.
1561	Bufs(:) = 0.
1562	alb1(:) = 0.
1563	alb2(:) = 0.
1564	alb3(:) = 0.
1565
1566	!***************************************************************************
1567	! Uncompress SISVAT output variables for storage
1568
1569	DO ikl = 1,klon
1570	i = ikl2i(ikl)
1571	IF (i > 0) THEN
1572	isno(i) = 1.*isnoSV(ikl) ! Nb Snow/Ice Lay.
1573	ispi(i) = 1.*ispiSV(ikl) ! Nb Supr.Ice Lay.
1574	iice(i) = 1.*iiceSV(ikl) ! Nb Ice Lay.
1575
1576	! IRs(i) = IRs_SV(ikl)
1577	! LMO(i) = LMO_SV(ikl)
1578	! turb_vel(i,1) = us__SV(ikl)
1579	! turb_vel(i,2) = uts_SV(ikl)
1580	! turb_vel(i,3) = uqs_SV(ikl)
1581	! turb_vel(i,4) = uss_SV(ikl)
1582	! turb_vel(i,5) = usthSV(ikl)
1583	! rlength(i,1) = Z0m_SV(ikl) ! Moment.Roughn.L.
1584	! rlength(i,2) = Z0mmSV(ikl) !
1585	! rlength(i,3) = Z0mnSV(ikl) !
1586	! rlength(i,4) = Z0SaSV(ikl) !
1587	! rlength(i,5) = Z0e_SV(ikl) !
1588	! rlength(i,6) = Z0emSV(ikl) !
1589	! rlength(i,7) = Z0enSV(ikl) !
1590	! rlength(i,8) = Z0roSV(ikl) !
1591	! rlength(i,9) = Z0h_SV(ikl) !
1592
1593	DO isl = -nsol,0 !
1594	eta(i,nsno+1-isl) = eta_SV(ikl,isl) ! Soil Humidity
1595	Tsis(i,nsno+1-isl) = TsisSV(ikl,isl) ! Soil Temperature
1596	ro(i,nsno+1-isl) = ro__SV(ikl,isl) ! [kg/m3]
1597	END DO
1598
1599
1600	DO isn = 1,nsno
1601	isto(i,isn) = 1.*istoSV(ikl,isn) ! Snow History
1602	ro(i,isn) = ro__SV(ikl,isn) ! [kg/m3]
1603	eta(i,isn) = eta_SV(ikl,isn) ! [m3/m3]
1604	Tsis(i,isn) = TsisSV(ikl,isn) ! [K]
1605	G1sn(i,isn) = G1snSV(ikl,isn) ! [-] [-]
1606	G2sn(i,isn) = G2snSV(ikl,isn) ! [-] [0.0001 m]
1607	dzsn(i,isn) = dzsnSV(ikl,isn) ! [m]
1608	agsn(i,isn) = agsnSV(ikl,isn) ! [day]
1609	END DO
1610	rusn(i) = rusnSV(ikl) ! Surficial Water
1611	toic(i) = toicSV(ikl) ! Surficial Water
1612	alb1(i) = alb1sv(ikl)
1613	alb2(i) = alb2sv(ikl)
1614	alb3(i) = alb3sv(ikl)
1615	! Bufs(i) = BufsSV(ikl)
1616	END IF
1617	END DO
1618
1619	CALL open_restartphy(fichnom)
1620	DO pass = 1, 2
1621	CALL put_field(pass,"longitude", &
1622	"Longitudes de la grille physique",rlon)
1623	CALL put_field(pass,"latitude","Latitudes de la grille physique",rlat)
1624
1625	CALL put_field(pass,"n_snows", "number of snow/ice layers",isno)
1626	CALL put_field(pass,"n_ice_top", "number of top ice layers",ispi)
1627	CALL put_field(pass,"n_ice", "number of ice layers",iice)
1628	CALL put_field(pass,"IR_soil", "Soil IR flux",IRs)
1629	CALL put_field(pass,"LMO", "Monin-Obukhov Scale",LMO)
1630	CALL put_field(pass,"surf_water", "Surficial water",rusn)
1631	CALL put_field(pass,"snow_buffer", "Snow buffer layer",Bufs)
1632	CALL put_field(pass,"alb_1", "albedo sw",alb1)
1633	CALL put_field(pass,"alb_2", "albedo nIR",alb2)
1634	CALL put_field(pass,"alb_3", "albedo fIR",alb3)
1635	CALL put_field(pass,"to_ice", "Snow passed to ice",toic)
1636
1637	! DO i = 1, 5
1638	! WRITE(str2,'(i2.2)') i
1639	! CALL put_field(pass,"turb_veloc"//str2, &
1640	! "various turbulent velocities"//str2, &
1641	! turb_vel(:,i))
1642	! ENDDO
1643	! DO i = 1, 9
1644	! WRITE(str2,'(i2.2)') i
1645	! CALL put_field(pass,"rough_length"//str2, &
1646	! "various roughness lengths"//str2, &
1647	! rlength(:,i))
1648	! ENDDO
1649	DO isn = 1,nsno
1650	IF (isn.LE.99) THEN
1651	WRITE(str2,'(i2.2)') isn
1652	CALL put_field(pass,"AGESNOW"//str2, &
1653	"Age de la neige layer No."//str2, &
1654	agsn(:,isn))
1655	ELSE
1656	PRINT*, "Trop de couches"
1657	CALL abort
1658	ENDIF
1659	ENDDO
1660	DO isn = 1,nsno
1661	IF (isn.LE.99) THEN
1662	WRITE(str2,'(i2.2)') isn
1663	CALL put_field(pass,"DZSNOW"//str2, &
1664	"Snow/ice thickness layer No."//str2, &
1665	dzsn(:,isn))
1666	ELSE
1667	PRINT*, "Trop de couches"
1668	CALL abort
1669	ENDIF
1670	ENDDO
1671	DO isn = 1,nsno
1672	IF (isn.LE.99) THEN
1673	WRITE(str2,'(i2.2)') isn
1674	CALL put_field(pass,"G2SNOW"//str2, &
1675	"Snow Property 2, layer No."//str2, &
1676	G2sn(:,isn))
1677	ELSE
1678	PRINT*, "Trop de couches"
1679	CALL abort
1680	ENDIF
1681	ENDDO
1682	DO isn = 1,nsno
1683	IF (isn.LE.99) THEN
1684	WRITE(str2,'(i2.2)') isn
1685	CALL put_field(pass,"G1SNOW"//str2, &
1686	"Snow Property 1, layer No."//str2, &
1687	G1sn(:,isn))
1688	ELSE
1689	PRINT*, "Trop de couches"
1690	CALL abort
1691	ENDIF
1692	ENDDO
1693	DO isn = 1,nsismx
1694	IF (isn.LE.99) THEN
1695	WRITE(str2,'(i2.2)') isn
1696	CALL put_field(pass,"ETA"//str2, &
1697	"Soil/snow water content layer No."//str2, &
1698	eta(:,isn))
1699	ELSE
1700	PRINT*, "Trop de couches"
1701	CALL abort
1702	ENDIF
1703	ENDDO
1704	DO isn = 1,nsismx !nsno
1705	IF (isn.LE.99) THEN
1706	WRITE(str2,'(i2.2)') isn
1707	CALL put_field(pass,"RO"//str2, &
1708	"Snow density layer No."//str2, &
1709	ro(:,isn))
1710	ELSE
1711	PRINT*, "Trop de couches"
1712	CALL abort
1713	ENDIF
1714	ENDDO
1715	DO isn = 1,nsismx
1716	IF (isn.LE.99) THEN
1717	WRITE(str2,'(i2.2)') isn
1718	CALL put_field(pass,"TSS"//str2, &
1719	"Soil/snow temperature layer No."//str2, &
1720	Tsis(:,isn))
1721	ELSE
1722	PRINT*, "Trop de couches"
1723	CALL abort
1724	ENDIF
1725	ENDDO
1726	DO isn = 1,nsno
1727	IF (isn.LE.99) THEN
1728	WRITE(str2,'(i2.2)') isn
1729	CALL put_field(pass,"HISTORY"//str2, &
1730	"Snow history layer No."//str2, &
1731	isto(:,isn))
1732	ELSE
1733	PRINT*, "Trop de couches"
1734	CALL abort
1735	ENDIF
1736	ENDDO
1737	ENDDO
1738
1739	END SUBROUTINE sisvatredem
1740
1741	END MODULE surf_sisvat_mod

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

Download in other formats:

Original Format