Context Navigation

surf_sisvat_mod.F90 @ 5442

Last change on this file since 5442 was 1872, checked in by Laurent Fairhead, 11 years ago

Correction de bugs suite à l'intégration de SISVAT

Bug corrections following integration of SISVAT

File size: 84.4 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
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