Context Navigation

surf_inlandsis_mod.F90 @ 3871

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

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

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

Download in other formats:

Original Format