Context Navigation

explicitsnow.f90 @ 7595

Last change on this file since 7595 was 7595, checked in by josefine.ghattas, 2 years ago
Correction on water balance, see ticket #844 and also ticket #642. Correction given by Fabienne Maignan and Catherine Ottle.
Property svn:keywords set to `Date Revision HeadURL`
File size: 73.1 KB

Line
1	! ================================================================================================================================
2	! MODULE : explicitsnow
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Computes hydrologic snow processes on continental points.
10	!!
11	!!\n DESCRIPTION: This module computes hydrologic snow processes on continental points.
12	!!
13	!! RECENT CHANGES:
14	!!
15	!! REFERENCES : None
16	!!
17	!! SVN :
18	!! $HeadURL$
19	!! $Date$
20	!! $Revision$
21	!! \n
22	!_ ================================================================================================================================
23
24	MODULE explicitsnow
25	USE ioipsl_para
26	USE constantes_soil
27	USE constantes
28	USE time, ONLY : one_day, dt_sechiba
29	USE pft_parameters
30	USE qsat_moisture
31	USE sechiba_io_p
32	USE xios_orchidee
33
34	IMPLICIT NONE
35
36	! Public routines :
37	PRIVATE
38	PUBLIC :: explicitsnow_main, explicitsnow_initialize, explicitsnow_finalize
39
40	CONTAINS
41
42	!================================================================================================================================
43	!! SUBROUTINE : explicitsnow_initialize
44	!!
45	!>\BRIEF Read variables for explictsnow module from restart file
46	!!
47	!! DESCRIPTION : Read variables for explictsnow module from restart file
48	!!
49	!! \n
50	!_
51	!================================================================================================================================
52	SUBROUTINE explicitsnow_initialize( kjit, kjpindex, rest_id, snowrho, &
53	snowtemp, snowdz, snowheat, snowgrain)
54
55	!! 0.1 Input variables
56	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
57	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
58	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
59
60	!! 0.2 Output variables
61	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowrho !! Snow density (Kg/m^3)
62	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowtemp !! Snow temperature
63	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowdz !! Snow layer thickness
64	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowheat !! Snow heat content (J/m^2)
65	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowgrain !! Snow grainsize (m)
66
67
68	!! 1. Read from restart file
69	CALL restget_p (rest_id, 'snowrho', nbp_glo, nsnow, 1, kjit,.TRUE.,snowrho, "gather", nbp_glo, index_g)
70	CALL setvar_p (snowrho, val_exp, 'Snow Density profile', xrhosmin)
71
72	CALL restget_p (rest_id, 'snowtemp', nbp_glo, nsnow, 1, kjit,.TRUE.,snowtemp, "gather", nbp_glo, index_g)
73	CALL setvar_p (snowtemp, val_exp, 'Snow Temperature profile', tp_00)
74
75	CALL restget_p (rest_id, 'snowdz', nbp_glo, nsnow, 1, kjit,.TRUE.,snowdz, "gather", nbp_glo, index_g)
76	CALL setvar_p (snowdz, val_exp, 'Snow depth profile', 0.0)
77
78	CALL restget_p (rest_id, 'snowheat', nbp_glo, nsnow, 1, kjit,.TRUE.,snowheat, "gather", nbp_glo, index_g)
79	CALL setvar_p (snowheat, val_exp, 'Snow Heat profile', 0.0)
80
81	CALL restget_p (rest_id, 'snowgrain', nbp_glo, nsnow, 1, kjit,.TRUE.,snowgrain, "gather", nbp_glo, index_g)
82	CALL setvar_p (snowgrain, val_exp, 'Snow grain profile', 0.0)
83
84	END SUBROUTINE explicitsnow_initialize
85
86
87	!================================================================================================================================
88	!! SUBROUTINE : explicitsnow_main
89	!!
90	!>\BRIEF Main call for snow calculations
91	!!
92	!! DESCRIPTION : Main routine for calculation of the snow processes with explictsnow module.
93	!!
94	!! RECENT CHANGE(S) : None
95	!!
96	!! MAIN OUTPUT VARIABLE(S): None
97	!!
98	!! REFERENCE(S) :
99	!!
100	!! FLOWCHART : None
101	!! \n
102	!_
103	!================================================================================================================================
104	SUBROUTINE explicitsnow_main(kjpindex, precip_rain, precip_snow, temp_air, pb, & ! in
105	u, v, temp_sol_new, soilcap, pgflux, & ! in
106	frac_nobio, totfrac_nobio,gtemp, & ! in
107	lambda_snow, cgrnd_snow, dgrnd_snow, contfrac, & ! in
108	vevapsno, snow_age, snow_nobio_age,snow_nobio, snowrho, & ! inout
109	snowgrain, snowdz, snowtemp, snowheat, snow, & ! inout
110	temp_sol_add, & ! inout
111	snowliq, subsnownobio, grndflux, snowmelt, tot_melt, & ! output
112	subsinksoil ) ! output
113
114
115	!! 0. Variable and parameter declaration
116
117	!! 0.1 Input variables
118	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
119	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall
120	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_snow !! Snowfall
121	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_air !! Air temperature
122	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: pb !! Surface pressure
123	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: u,v !! Horizontal wind speed
124	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_sol_new !! Surface temperature
125	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: soilcap !! Soil capacity
126	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: pgflux !! Net energy into snowpack
127	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
128	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+ ...
129	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: gtemp !! First soil layer temperature
130	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
131	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
132	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
133	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid
134
135	!! 0.2 Output fields
136	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowliq !! Snow liquid content (m)
137	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(out) :: subsnownobio !! Sublimation of snow on other surface types (ice, lakes, ...)
138	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: grndflux !! Net flux into soil [W/m2]
139	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: snowmelt !! Snow melt [mm/dt_sechiba]
140	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_melt !! Total melt from ice and snow
141	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: subsinksoil !! Excess of sublimation as a sink for the soil
142
143	!! 0.3 Modified fields
144	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapsno !! Snow evaporation @tex ($kg m^{-2}$) @endtex
145	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow_age !! Snow age
146	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio !! Ice water balance
147	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio_age !! Snow age on ice, lakes, ...
148	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
149	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowgrain !! Snow grainsize (m)
150	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow layer thickness [m]
151	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature
152	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowheat !! Snow heat content (J/m^2)
153	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow !! Snow mass [Kg/m^2]
154	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
155
156
157	!! 0.4 Local declaration
158	INTEGER(i_std) :: ji, iv, jj,m,jv
159	REAL(r_std),DIMENSION (kjpindex) :: snow_depth_tmp
160	REAL(r_std),DIMENSION (kjpindex) :: snowmelt_from_maxmass
161	REAL(r_std) :: snowdzm1
162	REAL(r_std), DIMENSION (kjpindex) :: thrufal !! Water leaving snow pack [kg/m2/s]
163	REAL(r_std), DIMENSION (kjpindex) :: d_age !! Snow age change
164	REAL(r_std), DIMENSION (kjpindex) :: xx !! Temporary
165	REAL(r_std), DIMENSION (kjpindex) :: snowmelt_tmp,snowmelt_ice,icemelt,temp_sol_new_old
166	REAL(r_std), DIMENSION (kjpindex,nsnow) :: snowdz_old
167	REAL(r_std), DIMENSION (kjpindex) :: ZLIQHEATXS
168	REAL(r_std), DIMENSION (kjpindex) :: ZSNOWEVAPS, ZSNOWDZ,subsnowveg,ZSNOWMELT_XS
169	REAL(r_std) :: maxmass_snowdepth, snow_d1k
170	REAL(r_std), DIMENSION (kjpindex,nsnow) :: WSNOWDZ,SMASS
171	REAL(r_std), DIMENSION (kjpindex) :: SMASSC,snowacc,snow_remove
172	INTEGER(i_std) :: locjj
173	REAL(r_std) :: grndflux_tmp
174	REAL(r_std), DIMENSION (nsnow) :: snowtemp_tmp
175	REAL(r_std) :: s2flux_tmp,fromsoilflux
176	REAL(r_std), DIMENSION (kjpindex,nsnow) :: pcapa_snow
177	REAL(r_std), DIMENSION (kjpindex) :: psnowhmass
178	REAL(r_std), PARAMETER :: XP00 = 1.E5
179	REAL(r_std), DIMENSION (kjpindex) :: snowheattot_begin,snowheattot_end,del_snowheattot !![J/m2]
180	REAL(r_std), DIMENSION (kjpindex,nsnow) :: snowliq_diag ! Snow liquid content per snow layer [kg/m2]
181	REAL(r_std), DIMENSION (kjpindex) :: snowliqtot_diag ! Snow liquid content integrated over snow depth [kg/m2]
182	REAL(r_std), DIMENSION (kjpindex,nsnow) :: snowrho_diag
183	REAL(r_std), DIMENSION (kjpindex,nsnow) :: snowheat_diag
184	REAL(r_std), DIMENSION (kjpindex,nsnow) :: snowgrain_diag
185
186	!! 1. Initialization
187
188	temp_sol_new_old = temp_sol_new
189	snowmelt_ice(:) = zero
190	icemelt(:) = zero
191	tot_melt(:) = zero
192	snowmelt(:) = zero
193	snowmelt_from_maxmass(:) = zero
194
195
196	!! 2. on Vegetation
197	! 2.1 Snow fall
198	snowheattot_begin(:)=SUM(snowheat(:,:),2)
199	CALL explicitsnow_fall(kjpindex,precip_snow,temp_air,u,v,totfrac_nobio,snowrho,snowdz,&
200	& snowheat,snowgrain,snowtemp,psnowhmass)
201
202	! 2.2 calculate the new snow discretization
203	snow_depth_tmp(:) = SUM(snowdz(:,:),2)
204
205	snowdz_old = snowdz
206
207	CALL explicitsnow_levels(kjpindex,snow_depth_tmp, snowdz)
208
209	! 2.3 Snow heat redistribution
210	CALL explicitsnow_transf(kjpindex,snowdz_old,snowdz,snowrho,snowheat,snowgrain)
211
212	! 2.4 Diagonize water portion of the snow from snow heat content:
213	DO ji=1, kjpindex
214	IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
215	snowtemp(ji,:) = snow3ltemp_1d(snowheat(ji,:),snowrho(ji,:),snowdz(ji,:))
216	snowliq(ji,:) = snow3lliq_1d(snowheat(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:))
217	ELSE
218	snowliq(ji,:) = zero
219	snowtemp(ji,:) = tp_00
220	ENDIF
221	END DO
222
223	! 2.5 snow compaction
224	CALL explicitsnow_compactn(kjpindex,snowtemp,snowrho,snowdz)
225	! Update snow heat
226	DO ji = 1, kjpindex
227	snowheat(ji,:) = snow3lheat_1d(snowliq(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:))
228	ENDDO
229
230	! 2.6 Calculate the snow temperature profile based on heat diffusion
231	CALL explicitsnow_profile (kjpindex,cgrnd_snow,dgrnd_snow,lambda_snow,temp_sol_new, snowtemp,snowdz,temp_sol_add)
232
233	! 2.7 Test whether snow is existed on the ground or not
234	grndflux(:)=0.0
235	CALL explicitsnow_gone(kjpindex,pgflux,&
236	& snowheat,snowtemp,snowdz,snowrho,snowliq,grndflux,snowmelt)
237
238	! 2.8 Calculate snow melt/refreezing processes
239	CALL explicitsnow_melt_refrz(kjpindex,precip_rain,pgflux,soilcap,&
240	& snowtemp,snowdz,snowrho,snowliq,snowmelt,grndflux,temp_air)
241
242
243	! 2.9 Snow sublimation changing snow thickness
244	snow(:) = 0.0
245	DO ji=1,kjpindex !domain size
246	snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
247	ENDDO
248
249	subsnownobio(:,:) = zero
250	subsinksoil(:) = zero
251
252	DO ji=1, kjpindex ! domain size
253	IF ( snow(ji) > snowcri ) THEN
254	! There is enough snow to split the sublimation demand between the nobio and vegetation parts.
255	IF (snow_nobio(ji,iice) .GE. vevapsno(ji)) THEN
256	! There is enough snow on the ice fraction to sustain the sublimation demand.
257	subsnownobio(ji,iice) = frac_nobio(ji,iice)*vevapsno(ji)
258	ELSE
259	! There is not enough snow on the ice fraction to sustain the full sublimation demand.
260	! We take all the nobio snow and set the remaining sublimation demand on the vegetation fraction.
261	! We do this because the nobio snow cannot deal with a too high sublimation demand, whereas the
262	! vegetation fraction can, through the subsinksoil variable (see below 2.8.1), which is taken into
263	! account in the water budget (see TWBR in hydrol.f90).
264	subsnownobio(ji,iice) = snow_nobio(ji,iice)
265	ENDIF
266	subsnowveg(ji) = vevapsno(ji) - subsnownobio(ji,iice)
267	ELSE
268	! There is a small amount of snow.
269	IF ( frac_nobio(ji,iice) .GT. min_sechiba) THEN
270	! The ice fraction is not too small. In this case, we make the hypothesis that the snow is mainly on the ice fraction.
271	IF (snow_nobio(ji,iice) .GE. vevapsno(ji)) THEN
272	! There is enough snow on the ice fraction to sustain the sublimation demand.
273	subsnownobio(ji,iice) = vevapsno(ji)
274	subsnowveg(ji) = zero
275	ELSE
276	! There is not enough snow on the ice fraction to sustain the sublimation demand.
277	! We take all the nobio snow and set the remaining sublimation demand on the vegetation fraction.
278	! We do this beacuse the nobio snow cannot deal with a too high sublimation demand, whereas the
279	! vegetation fraction can, through the subsinksoil variable (see below 2.8.1), which is taken into
280	! account in the water budget (see TWBR in hydrol.f90).
281	subsnownobio(ji,iice) = snow_nobio(ji,iice)
282	subsnowveg(ji) = vevapsno(ji) - subsnownobio(ji,iice)
283	ENDIF
284	ELSE
285	! The nobio snow fraction is too small, the vegetation fraction deals with the whole sublimation demand.
286	subsnownobio(ji,iice) = zero
287	subsnowveg(ji) = vevapsno(ji)
288	ENDIF
289	ENDIF
290
291	!! 2.8.1 Check that sublimation on the vegetated fraction is possible.
292	IF (subsnowveg(ji) .GT. snow(ji)) THEN
293	IF( (un - totfrac_nobio(ji)).GT.min_sechiba) THEN
294	subsinksoil (ji) = (subsnowveg(ji) - snow(ji))/ (un - totfrac_nobio(ji))
295	END IF
296	! Sublimation is thus limited to what is available
297	subsnowveg(ji) = snow(ji)
298	snow(ji) = zero
299	snowdz(ji,:) = 0
300	snowliq(ji,:) = 0
301	snowtemp(ji,:) = tp_00
302	vevapsno(ji) = subsnowveg(ji) + subsnownobio(ji,iice)
303	ELSE
304	! Calculating the snow accumulation
305	WSNOWDZ(ji,:)= snowdz(ji,:)*snowrho(ji,:)
306	SMASSC (ji)= 0.0
307	DO jj=1,nsnow
308	SMASS(ji,jj) = SMASSC(ji) + WSNOWDZ(ji,jj)
309	SMASSC(ji) = SMASSC(ji) + WSNOWDZ(ji,jj)
310	ENDDO
311	! Finding the layer
312	locjj=1
313	DO jj=1,nsnow-1
314	IF ((SMASS(ji,jj) .LE. subsnowveg(ji)) .AND. (SMASS(ji,jj+1) .GE. subsnowveg(ji)) ) THEN
315	locjj=jj+1
316	ENDIF
317	ENDDO
318
319	! Calculating the removal of snow depth
320	IF (locjj .EQ. 1) THEN
321	ZSNOWEVAPS(ji) = subsnowveg(ji)/snowrho(ji,1)
322	ZSNOWDZ(ji) = snowdz(ji,1) - ZSNOWEVAPS(ji)
323	snowdz(ji,1) = MAX(0.0, ZSNOWDZ(ji))
324	ELSE IF (locjj .GT. 1) THEN
325	snowacc(ji)=0
326	DO jj=1,locjj-1
327	snowacc(ji)=snowacc(ji)+snowdz(ji,jj)*snowrho(ji,jj)
328	snowdz(ji,jj)=0
329	ENDDO
330	ZSNOWEVAPS(ji) = (subsnowveg(ji)-snowacc(ji))/snowrho(ji,locjj)
331	ZSNOWDZ(ji) = snowdz(ji,locjj) - ZSNOWEVAPS(ji)
332	snowdz(ji,locjj) = MAX(0.0, ZSNOWDZ(ji))
333	! ELSE
334	! ZSNOWEVAPS(ji) = subsnowveg(ji)/snowrho(ji,1)
335	! ZSNOWDZ(ji) = snowdz(ji,1) - ZSNOWEVAPS(ji)
336	! snowdz(ji,1) = MAX(0.0, ZSNOWDZ(ji))
337	ENDIF
338
339	ENDIF
340	ENDDO
341
342	! 2.9.1 Snowmelt_from_maxmass
343
344	DO ji=1,kjpindex !domain size
345	snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
346	IF(snow(ji).GT.maxmass_snow)THEN
347	IF (snow(ji).GT.(1.2*maxmass_snow)) THEN
348	! melt faster for points where accumulation is too high
349	snow_d1k = trois * soilcap(ji) / chalfu0
350	ELSE
351	! standard melting
352	snow_d1k = un * soilcap(ji) / chalfu0
353	ENDIF
354	snowmelt_from_maxmass(ji) = MIN((snow(ji) - maxmass_snow),snow_d1k)
355	! Calculating the snow accumulation
356	WSNOWDZ(ji,:)= snowdz(ji,:)*snowrho(ji,:) ! WSNOWDZ en kg/m2
357	SMASSC (ji)= 0.0
358	DO jj=nsnow,1,-1
359	SMASS(ji,jj) = SMASSC(ji) + WSNOWDZ(ji,jj)
360	SMASSC(ji) = SMASSC(ji) + WSNOWDZ(ji,jj)
361	ENDDO
362
363	! Finding the layer
364	locjj = nsnow
365	DO jj=nsnow,2,-1
366	IF ((SMASS(ji,jj) .LE. snowmelt_from_maxmass(ji)) .AND. (SMASS(ji,jj-1) .GE. snowmelt_from_maxmass(ji)) ) THEN
367	locjj=jj-1
368	ENDIF
369	ENDDO
370
371	! Calculating the removal of snow depth
372	IF (locjj .EQ. 3) THEN
373	! ZSNOWMELT_XS : Epaisseur de la couche qui a disparu partiellement
374	ZSNOWMELT_XS(ji) = snowmelt_from_maxmass(ji)/snowrho(ji,3)
375	ZSNOWDZ(ji) = snowdz(ji,3) - ZSNOWMELT_XS(ji)
376	snowdz(ji,3) = MAX(0.0, ZSNOWDZ(ji))
377	ELSE IF (locjj .LT. 3) THEN
378	snow_remove(ji)=0 ! masse de neige des couches qui disparaissent totalement
379	DO jj=3,locjj+1,-1
380	snow_remove(ji)=snow_remove(ji)+snowdz(ji,jj)*snowrho(ji,jj)
381	snowdz(ji,jj)=0
382	ENDDO
383	ZSNOWMELT_XS(ji) = (snowmelt_from_maxmass(ji) - snow_remove(ji))/snowrho(ji,locjj)
384	ZSNOWDZ(ji) = snowdz(ji,locjj) - ZSNOWMELT_XS(ji)
385	snowdz(ji,locjj) = MAX(0.0, ZSNOWDZ(ji))
386	ENDIF
387	ENDIF
388
389	ENDDO
390
391	!2.10 Calculate snow grain size using the updated thermal gradient
392
393	CALL explicitsnow_grain(kjpindex,snowliq,snowdz,gtemp,snowtemp,pb,snowgrain)
394
395	!2.11 Update snow heat
396	! Update the heat content (variable stored each time step)
397	! using current snow temperature and liquid water content:
398	!
399	! First, make check to make sure heat content not too large
400	! (this can result due to signifigant heating of thin snowpacks):
401	! add any excess heat to ground flux:
402	!
403	DO ji=1,kjpindex
404	DO jj=1,nsnow
405	ZLIQHEATXS(ji) = MAX(0.0, snowliq(ji,jj)ph2o - 0.10snowdz(ji,jj)snowrho(ji,jj))chalfu0/dt_sechiba
406	snowliq(ji,jj) = snowliq(ji,jj) - ZLIQHEATXS(ji)dt_sechiba/(ph2ochalfu0)
407	snowliq(ji,jj) = MAX(0.0, snowliq(ji,jj))
408	grndflux(ji) = grndflux(ji) + ZLIQHEATXS(ji)
409	ENDDO
410	ENDDO
411
412	snow(:) = 0.0
413	DO ji=1,kjpindex !domain size
414	snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
415	ENDDO
416
417	DO ji = 1, kjpindex
418	snowheat(ji,:) = snow3lheat_1d(snowliq(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:))
419	ENDDO
420	snowheattot_end(:)=SUM(snowheat(:,:),2)
421	del_snowheattot(:)=snowheattot_end(:)-snowheattot_begin(:)
422
423	!3. on land ice (using the default ORCHIDEE snow scheme)
424
425	DO ji = 1,kjpindex
426
427	!! 3.1. It is snowing
428
429	snow_nobio(ji,iice) = snow_nobio(ji,iice) + frac_nobio(ji,iice)*precip_snow(ji) + &
430	& frac_nobio(ji,iice)*precip_rain(ji)
431
432	!! 3.2. Sublimation - was calculated before, it cannot give us negative snow_nobio (see 2.9).
433
434	snow_nobio(ji,iice) = snow_nobio(ji,iice) - subsnownobio(ji,iice)
435
436	!! 3.3. ice melt only for continental ice fraction
437
438	snowmelt_tmp(ji) = zero
439	IF (temp_sol_new_old(ji) .GT. tp_00) THEN
440
441	!! 3.3.1 If there is snow on the ice-fraction it can melt
442
443	snowmelt_tmp(ji) = frac_nobio(ji,iice)(temp_sol_new_old(ji) - tp_00) soilcap(ji) / chalfu0
444
445	IF ( snowmelt_tmp(ji) .GT. snow_nobio(ji,iice) ) THEN
446	snowmelt_tmp(ji) = MAX( 0., snow_nobio(ji,iice))
447	ENDIF
448	snowmelt_ice(ji) = snowmelt_ice(ji) + snowmelt_tmp(ji)
449	snow_nobio(ji,iice) = snow_nobio(ji,iice) - snowmelt_tmp(ji)
450
451	ENDIF
452
453	!! Ice melt only if there is more than a given mass : maxmass_snow, i.e. only weight melts glaciers !
454
455	IF ( snow_nobio(ji,iice) .GE. maxmass_snow ) THEN
456	icemelt(ji) = snow_nobio(ji,iice) - maxmass_snow
457	snow_nobio(ji,iice) = maxmass_snow
458	ENDIF
459
460	END DO
461
462
463	!! 4. On other surface types - not done yet
464
465	IF ( nnobio .GT. 1 ) THEN
466	WRITE(numout,*) 'WE HAVE',nnobio-1,' SURFACE TYPES I DO NOT KNOW'
467	WRITE(numout,*) 'CANNOT TREAT SNOW ON THESE SURFACE TYPES'
468	CALL ipslerr_p(3,'explicitsnow_main','Cannot treat snow on these surface types.','','')
469	ENDIF
470
471	!! 5. Computes snow age on land and land ice (for albedo)
472
473	DO ji = 1, kjpindex
474
475	!! 5.1. Snow age on land
476
477	IF (snow(ji) .LE. zero) THEN
478	snow_age(ji) = zero
479	ELSE
480	snow_age(ji) =(snow_age(ji) + (un - snow_age(ji)/max_snow_age) * dt_sechiba/one_day) &
481	& * EXP(-precip_snow(ji) / snow_trans)
482	ENDIF
483
484	!! 5.2. Snow age on land ice
485
486	!! Age of snow on ice: a little bit different because in cold regions, we really
487	!! cannot negect the effect of cold temperatures on snow metamorphism any more.
488
489	IF (snow_nobio(ji,iice) .LE. zero) THEN
490	snow_nobio_age(ji,iice) = zero
491	ELSE
492
493	d_age(ji) = ( snow_nobio_age(ji,iice) + &
494	& (un - snow_nobio_age(ji,iice)/max_snow_age) * dt_sechiba/one_day ) * &
495	& EXP(-precip_snow(ji) / snow_trans) - snow_nobio_age(ji,iice)
496	IF (d_age(ji) .GT. 0. ) THEN
497	xx(ji) = MAX( tp_00 - temp_sol_new(ji), zero )
498	xx(ji) = ( xx(ji) / 7._r_std ) ** 4._r_std
499	d_age(ji) = d_age(ji) / (un+xx(ji))
500	ENDIF
501	snow_nobio_age(ji,iice) = MAX( snow_nobio_age(ji,iice) + d_age(ji), zero )
502
503	ENDIF
504
505	ENDDO
506
507
508	!! 6. Check the snow on land
509	DO ji=1,kjpindex
510	IF (snow(ji) .EQ. 0) THEN
511	snowrho(ji,:)=50.0
512	snowgrain(ji,:)=0.0
513	snowdz(ji,:)=0.0
514	snowliq(ji,:)=0.0
515	ENDIF
516	ENDDO
517
518
519	! Snow melt only if there is more than a given mass : maxmass_snow
520	! Here I suggest to remove the snow based on a certain threshold of snow
521	! depth instead of snow mass because it is quite difficult for
522	! explictsnow module to calculate other snow properties following the
523	! removal of snow mass
524	! to define the threshold of snow depth based on old snow density (330
525	! kg/m3)
526	! maxmass_snowdepth=maxmass_snow/sn_dens
527	! snowmelt_from_maxmass(:) = 0.0
528
529	!! 7. compute total melt
530
531	DO ji=1,kjpindex
532	tot_melt(ji) = icemelt(ji) + snowmelt(ji) + snowmelt_ice(ji) + snowmelt_from_maxmass(ji)
533	ENDDO
534	IF (printlev>=3) WRITE(numout,*) 'explicitsnow_main done'
535
536
537	!! 8. Recalculate the new snow discretization
538	! This is done to make sure that all levels have snow
539	snow_depth_tmp(:) = SUM(snowdz(:,:),2)
540	snowdz_old = snowdz
541	CALL explicitsnow_levels(kjpindex,snow_depth_tmp, snowdz)
542	CALL explicitsnow_transf(kjpindex,snowdz_old,snowdz,snowrho,snowheat,snowgrain)
543
544	! Write output variables
545
546	! Add XIOS default value where no snow
547	DO ji=1,kjpindex
548	IF (snow(ji) .GT. zero) THEN
549	! Output for snowliq and snowliqtot change unit from m to kg/m2 and are multiplied by contfrac
550	! to take into acount the portion form land fraction only
551	snowliq_diag(ji,:) = snowliq(ji,:) * 1000 * contfrac(ji)
552	snowliqtot_diag(ji) = SUM(snowliq_diag(ji,:)) * 1000 * contfrac(ji)
553	snowrho_diag(ji,:) = snowrho(ji,:)
554	snowheat_diag(ji,:) = snowheat(ji,:)
555	snowgrain_diag(ji,:) = snowgrain(ji,:)
556	ELSE
557	snowliq_diag(ji,:) = xios_default_val
558	snowliqtot_diag(ji) = xios_default_val
559	snowrho_diag(ji,:) = xios_default_val
560	snowheat_diag(ji,:) = xios_default_val
561	snowgrain_diag(ji,:) = xios_default_val
562	END IF
563	END DO
564
565	CALL xios_orchidee_send_field("snowliq",snowliq_diag)
566	CALL xios_orchidee_send_field("snowliqtot", snowliqtot_diag)
567	CALL xios_orchidee_send_field("snowrho",snowrho_diag)
568	CALL xios_orchidee_send_field("snowheat",snowheat_diag)
569	CALL xios_orchidee_send_field("snowgrain",snowgrain_diag)
570	CALL xios_orchidee_send_field("snowmelt_from_maxmass",snowmelt_from_maxmass)
571	CALL xios_orchidee_send_field("soilcap",soilcap)
572	CALL xios_orchidee_send_field("del_snowheattot",del_snowheattot)
573
574	END SUBROUTINE explicitsnow_main
575
576
577	!================================================================================================================================
578	!! SUBROUTINE : explicitsnow_finalize
579	!!
580	!>\BRIEF Write variables for explictsnow module to restart file
581	!!
582	!! DESCRIPTION : Write variables for explictsnow module to restart file
583	!!
584	!! \n
585	!_
586	!================================================================================================================================
587	SUBROUTINE explicitsnow_finalize ( kjit, kjpindex, rest_id, snowrho, &
588	snowtemp, snowdz, snowheat, snowgrain)
589
590	!! 0.1 Input variables
591	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
592	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
593	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
594	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density (Kg/m^3)
595	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature
596	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow layer thickness
597	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowheat !! Snow heat content (J/m^2)
598	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowgrain !! Snow grainsize (m)
599
600
601	!! 1. Write to restart file
602	CALL restput_p(rest_id, 'snowrho', nbp_glo, nsnow, 1, kjit, snowrho, 'scatter', nbp_glo, index_g)
603	CALL restput_p(rest_id, 'snowtemp', nbp_glo, nsnow, 1, kjit, snowtemp, 'scatter', nbp_glo, index_g)
604	CALL restput_p(rest_id, 'snowdz', nbp_glo, nsnow, 1, kjit, snowdz, 'scatter', nbp_glo, index_g)
605	CALL restput_p(rest_id, 'snowheat', nbp_glo, nsnow, 1, kjit, snowheat, 'scatter', nbp_glo, index_g)
606	CALL restput_p(rest_id, 'snowgrain', nbp_glo, nsnow, 1, kjit, snowgrain, 'scatter', nbp_glo, index_g)
607
608	END SUBROUTINE explicitsnow_finalize
609
610
611	!================================================================================================================================
612	!! SUBROUTINE : explicitsnow_grain
613	!!
614	!>\BRIEF Compute evolution of snow grain size
615	!!
616	!! DESCRIPTION :
617	!!
618	!! RECENT CHANGE(S) : None
619	!!
620	!! MAIN OUTPUT VARIABLE(S): None
621	!!
622	!! REFERENCE(S) : R. Jordan (1991)
623	!!
624	!! FLOWCHART : None
625	!! \n
626	!_
627	!================================================================================================================================
628
629
630	SUBROUTINE explicitsnow_grain(kjpindex,snowliq,snowdz,gtemp,snowtemp,pb,snowgrain)
631
632	!! 0.1 Input variables
633	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
634	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowliq !! Liquid water content
635	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth (m)
636	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: gtemp !! First soil layer temperature
637	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature
638	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: pb !! Surface pressure (hpa)
639
640	!! 0.2 Output variables
641
642	!! 0.3 Modified variables
643
644	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m)
645
646	!! 0.4 Local variables
647	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowdz,zdz,ztheta
648	REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: ztemp,zdiff,ztgrad,zthetaa,zfrac,&
649	zexpo,zckt_liq,zckt_ice,zckt
650	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zrhomin,zgrainmin
651	INTEGER(i_std) :: ji,jj
652
653	!! 0.5 Local parameters
654	REAL(r_std), PARAMETER :: ztheta_crit = 0.02 !! m3 m-3
655	REAL(r_std), PARAMETER :: zc1_ice = 8.047E+9 !! kg m-3 K
656	REAL(r_std), PARAMETER :: zc1_liq = 5.726E+8 !! kg m-3 K
657	REAL(r_std), PARAMETER :: zdeos = 0.92E-4 !! effective diffusion
658	!! coef for water vapor in snow
659	!! at 0C and 1000 mb (m2 s-1)
660	REAL(r_std), PARAMETER :: zg1 = 5.0E-7 !! m4 kg-1
661	REAL(r_std), PARAMETER :: zg2 = 4.0E-12 !! m2 s-1
662	REAL(r_std), PARAMETER :: ztheta_w = 0.05 !! m3 m-3
663	REAL(r_std), PARAMETER :: ztheta_crit_w = 0.14 !! m3 m-3
664	REAL(r_std), PARAMETER :: zdzmin = 0.01 !! m : minimum thickness
665	!! for thermal gradient evaluation:
666	!! to prevent excessive gradients
667	!! for vanishingly thin snowpacks.
668	REAL(r_std), PARAMETER :: xp00=1.E5
669	!! 1. initialize
670
671	DO ji=1,kjpindex
672
673
674	zsnowdz(ji,:) = MAX(xsnowdmin/nsnow, snowdz(ji,:))
675
676	DO jj=1,nsnow-1
677	zdz(ji,jj) = zsnowdz(ji,jj) + zsnowdz(ji,jj+1)
678	ENDDO
679	zdz(ji,nsnow) = zsnowdz(ji,nsnow)
680
681	! compute interface average volumetric water content (m3 m-3):
682	! first, layer avg VWC:
683	!
684	ztheta(ji,:) = snowliq(ji,:)/MAX(xsnowdmin, zsnowdz(ji,:))
685
686	! at interfaces:
687	zthetaa(ji,0) = ztheta(ji,1)
688	DO jj=1,nsnow-1
689	zthetaa(ji,jj) = (zsnowdz(ji,jj) *ztheta(ji,jj) + &
690	zsnowdz(ji,jj+1)*ztheta(ji,jj+1))/zdz(ji,jj)
691	ENDDO
692	zthetaa(ji,nsnow) = ztheta(ji,nsnow)
693	! compute interface average temperatures (K):
694	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
695	!
696	ztemp(ji,0) = snowtemp(ji,1)
697	DO jj=1,nsnow-1
698	ztemp(ji,jj) = (zsnowdz(ji,jj) *snowtemp(ji,jj) + &
699	zsnowdz(ji,jj+1)*snowtemp(ji,jj+1))/zdz(ji,jj)
700	ENDDO
701	ztemp(ji,nsnow) = snowtemp(ji,nsnow)
702	!
703	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
704	! compute variation of saturation vapor pressure with temperature
705	! for solid and liquid phases:
706	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
707	zexpo(ji,:) = chalsu0/(xrv*ztemp(ji,:))
708	zckt_ice(ji,:) = (zc1_ice/ztemp(ji,:)*2)(zexpo(ji,:) - 1.0)*EXP(-zexpo(ji,:))
709	!
710	zexpo(ji,:) = chalev0/(xrv*ztemp(ji,:))
711	zckt_liq(ji,:) = (zc1_liq/ztemp(ji,:)*2)(zexpo(ji,:) - 1.0)*EXP(-zexpo(ji,:))
712	!
713	! compute the weighted ice/liquid total variation (N m-2 K):
714	!
715	zfrac(ji,:) = MIN(1.0, zthetaa(ji,:)/ztheta_crit)
716	zckt(ji,:) = zfrac(ji,:)zckt_liq(ji,:) + (1.0 - zfrac(ji,:))zckt_ice(ji,:)
717
718	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
719	! Compute effective diffusion coefficient (m2 s-1):
720	! -diffudivity relative to a reference diffusion at 1000 mb and freezing point
721	! multiplied by phase energy coefficient
722	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
723	!
724	DO jj=0,nsnow
725	zdiff(ji,jj) = zdeos(xp00/(pb(ji)100.))((ztemp(ji,jj)/tp_00)6)zckt(ji,jj)
726	ENDDO
727
728	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
729	! Temperature gradient (K m-1):
730
731	ztgrad(ji,0) = 0.0 ! uppermost layer-mean and surface T's are assumed to be equal
732	DO jj=1,nsnow-1
733	ztgrad(ji,jj) = 2*(snowtemp(ji,jj)-snowtemp(ji,jj+1))/MAX(zdzmin, zdz(ji,jj))
734	ENDDO
735	!
736	! assume at base of snow, temperature is in equilibrium with soil
737	! (but obviously must be at or below freezing point!)
738	!
739	ztgrad(ji,nsnow) = 2*(snowtemp(ji,nsnow) - MIN(tp_00, gtemp(ji)))/MAX(zdzmin, zdz(ji,nsnow))
740	! prognostic grain size (m) equation:
741	!-------------------------------------------------------------------
742	! first compute the minimum grain size (m):
743	!
744	zrhomin(ji,:) = xrhosmin
745	zgrainmin(ji,:) = snow3lgrain_1d(zrhomin(ji,:))
746
747	! dry snow:
748	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
749	!
750	DO jj=1,nsnow
751
752	IF(ztheta(ji,jj) == 0.0) THEN
753
754	! only allow growth due to vapor flux INTO layer:
755	! aab add sublimation(only condensation) as upper BC...?
756
757	snowgrain(ji,jj) = snowgrain(ji,jj) + &
758	(dt_sechibazg1/MAX(zgrainmin(ji,jj),snowgrain(ji,jj))) &
759	( zdiff(ji,jj-1)*MAX(0.0,ztgrad(ji,jj-1)) - &
760	zdiff(ji,jj) *MIN(0.0,ztgrad(ji,jj)) )
761	ELSE
762
763	! wet snow
764	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
765	!
766	snowgrain(ji,jj) = snowgrain(ji,jj) + &
767	(dt_sechibazg2/MAX(zgrainmin(ji,jj),snowgrain(ji,jj))) &
768	MIN(ztheta_crit_w, ztheta(ji,jj) + ztheta_w)
769	END IF
770
771	ENDDO
772
773
774	ENDDO
775
776
777	END SUBROUTINE explicitsnow_grain
778
779	!!
780	!================================================================================================================================
781	!! SUBROUTINE : explicitsnow_compactn
782	!!
783	!>\BRIEF Compute Compaction/Settling
784	!!
785	!! DESCRIPTION :
786	!! Snow compaction due to overburden and settling.
787	!! Mass is unchanged: layer thickness is reduced
788	!! in proportion to density increases. Method
789	!! of Anderson (1976): see Loth and Graf, 1993,
790	!! J. of Geophys. Res., 98, 10,451-10,464.
791	!!
792	!! RECENT CHANGE(S) : None
793	!!
794	!! MAIN OUTPUT VARIABLE(S): None
795	!!
796	!! REFERENCE(S) : Loth and Graf (1993), Mellor (1964) and Anderson (1976)
797	!!
798	!! FLOWCHART : None
799	!! \n
800	!_
801	!================================================================================================================================
802
803
804	SUBROUTINE explicitsnow_compactn(kjpindex,snowtemp,snowrho,snowdz)
805
806	!! 0.1 Input variables
807
808	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
809	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature
810
811	!! 0.2 Output variables
812
813	!! 0.3 Modified variables
814
815	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
816	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth
817
818	!! 0.4 Local variables
819
820	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zwsnowdz,zsmass,zsnowrho2,zviscocity,zsettle
821	REAL(r_std),DIMENSION(kjpindex) :: zsmassc !! cummulative snow mass (kg/m2)
822	REAL(r_std),DIMENSION(kjpindex) :: snowdepth_crit
823	INTEGER(i_std) :: ji,jj
824
825	!! 1. initialize
826
827	zsnowrho2 = snowrho
828	zsettle(:,:) = ZSNOWCMPCT_ACM
829	zviscocity(:,:) = ZSNOWCMPCT_V0
830
831	!! 2. Calculating Cumulative snow mass (kg/m2):
832
833	DO ji=1, kjpindex
834
835
836	IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
837
838	zwsnowdz(ji,:)= snowdz(ji,:)*snowrho(ji,:)
839
840	zsmassc (ji)= 0.0
841
842	DO jj=1,nsnow
843	zsmass(ji,jj) = zsmassc(ji) + zwsnowdz(ji,jj)
844	zsmassc(ji) = zsmassc(ji) + zwsnowdz(ji,jj)
845	ENDDO
846
847
848	!! 3. Computing compaction/Settling
849	! ----------------------
850	! Compaction/settling if density below upper limit
851	! (compaction is generally quite small above ~ 500 kg m-3):
852	!
853	DO jj=1,nsnow
854	IF (snowrho(ji,jj) .LT. xrhosmax) THEN
855
856	!
857	! First calculate settling due to freshly fallen snow: (NOTE:bug here for the snow temperature profile)
858	!
859
860	zsettle(ji,jj) = ZSNOWCMPCT_ACM*EXP( &
861	-ZSNOWCMPCT_BCM*(tp_00-MIN(tp_00,snowtemp(ji,jj))) &
862	-ZSNOWCMPCT_CCM*MAX(0.0, &
863	snowrho(ji,jj)-ZSNOWCMPCT_RHOD))
864	!
865	! Snow viscocity:
866	!
867	zviscocity(ji,jj) = ZSNOWCMPCT_V0EXP( ZSNOWCMPCT_VT(tp_00-MIN(tp_00,snowtemp(ji,jj))) + &
868	ZSNOWCMPCT_VR*snowrho(ji,jj) )
869
870	! Calculate snow density: compaction from weight/over-burden
871	! Anderson 1976 method:
872	zsnowrho2(ji,jj) = snowrho(ji,jj) + snowrho(ji,jj)dt_sechiba( &
873	(cte_grav*zsmass(ji,jj)/zviscocity(ji,jj)) &
874	+ zsettle(ji,jj) )
875	! Conserve mass by decreasing grid thicknesses in response
876	! to density increases
877	!
878	snowdz(ji,jj) = snowdz(ji,jj)*(snowrho(ji,jj)/zsnowrho2(ji,jj))
879
880	ENDIF
881	ENDDO
882
883	! Update density (kg m-3):
884	snowrho(ji,:) = zsnowrho2(ji,:)
885
886	ENDIF
887
888	ENDDO
889
890	END SUBROUTINE explicitsnow_compactn
891
892	!!
893	!================================================================================================================================
894	!! SUBROUTINE : explicitsnow_transf
895	!!
896	!>\BRIEF Computing snow mass and heat redistribution due to grid thickness configuration resetting
897	!!
898	!! DESCRIPTION : Snow mass and heat redistibution due to grid thickness
899	!! configuration resetting. Total mass and heat content
900	!! of the overall snowpack unchanged/conserved within this routine.
901	!! RECENT CHANGE(S) : None
902	!!
903	!! MAIN OUTPUT VARIABLE(S): None
904	!!
905	!! REFERENCE(S) :
906	!!
907	!! FLOWCHART : None
908	!! \n
909	!_
910	!================================================================================================================================
911
912	SUBROUTINE explicitsnow_transf(kjpindex,snowdz_old,snowdz,snowrho,snowheat,snowgrain)
913
914	!! 0.1 Input variables
915
916	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
917	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowdz_old !! Snow depth at the previous time step
918
919	!! 0.2 Output variables
920
921	!! 0.3 Modified variables
922
923	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
924	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m)
925	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth (m)
926	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowheat !! Snow heat content/enthalpy (J/m2)
927
928	!! 0.4 Local varibles
929
930	REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: zsnowzo
931	REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: zsnowzn
932	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowddz
933	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zdelta
934	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowrhon,zsnowheatn,zsnowgrainn
935	REAL(r_std),DIMENSION(kjpindex) :: zsnowmix_delta
936	REAL(r_std),DIMENSION(kjpindex) :: zsumheat, zsumswe, zsumgrain
937	INTEGER(i_std),DIMENSION(nsnow,2) :: locflag
938	REAL(r_std) :: psnow
939	INTEGER(i_std) :: ji,jj,jjj
940
941	! Initialization
942	zsumheat(:) = 0.0
943	zsumswe(:) = 0.0
944	zsumgrain(:) = 0.0
945	zsnowmix_delta(:) = 0.0
946	locflag(:,:) = 0
947
948	DO ji=1, kjpindex
949
950
951	psnow = SUM(snowdz(ji,:))
952
953	IF (psnow .GE. xsnowcritd .AND. snowdz_old(ji,1) .NE. 0 .AND. snowdz_old(ji,2) .NE. 0 .AND. snowdz_old(ji,3) .NE. 0) THEN
954	!
955	zsnowzo(ji,0) = 0.
956	zsnowzn(ji,0) = 0.
957	zsnowzo(ji,1) = snowdz_old(ji,1)
958	zsnowzn(ji,1) = snowdz(ji,1)
959	!
960	DO jj=2,nsnow
961	zsnowzo(ji,jj) = zsnowzo(ji,jj-1) + snowdz_old(ji,jj)
962	zsnowzn(ji,jj) = zsnowzn(ji,jj-1) + snowdz(ji,jj)
963	ENDDO
964
965	DO jj=1,nsnow
966
967	IF (jj .EQ. 1) THEN
968	locflag(jj,1)=1
969	ELSE
970	DO jjj=nsnow,1,-1
971
972	!upper bound of the snow layer
973	IF (zsnowzn(ji,jj-1) .LE. zsnowzo(ji,jjj)) THEN
974	locflag(jj,1) = jjj
975	ENDIF
976	ENDDO
977	ENDIF
978
979	IF (jj .EQ. nsnow) THEN
980	locflag(jj,2)=nsnow
981	ELSE
982	DO jjj=nsnow,1,-1
983
984	!lower bound of the snow layer
985	IF (zsnowzn(ji,jj) .LE. zsnowzo(ji,jjj)) THEN
986	locflag(jj,2) = jjj
987	ENDIF
988	ENDDO
989	ENDIF
990
991	!to interpolate
992	! when heavy snow occurred
993	IF (locflag(jj,1) .EQ. locflag(jj,2)) THEN
994
995	zsnowrhon(ji,jj) = snowrho(ji,locflag(jj,1))
996
997	zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,1))*snowdz(ji,jj)/snowdz_old(ji,locflag(jj,1))
998
999	zsnowgrainn(ji,jj) = snowgrain(ji,locflag(jj,1))
1000
1001	ELSE
1002
1003	!snow density
1004	zsnowrhon(ji,jj) = snowrho(ji,locflag(jj,1)) * &
1005	(zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1)) + &
1006	snowrho(ji,locflag(jj,2)) * &
1007	(zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))
1008
1009	DO jjj=locflag(jj,1),locflag(jj,2)-1
1010	zsnowrhon(ji,jj)=zsnowrhon(ji,jj) + &
1011	(jjj-locflag(jj,1))snowrho(ji,jjj)snowdz_old(ji,jjj)
1012	ENDDO
1013
1014	zsnowrhon(ji,jj) = zsnowrhon(ji,jj) / snowdz(ji,jj)
1015
1016
1017	!snow heat
1018	IF (snowdz_old(ji,locflag(jj,1)) .GT. 0.0) THEN
1019
1020	zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,1)) * &
1021	(zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1))/snowdz_old(ji,locflag(jj,1)) + &
1022	snowheat(ji,locflag(jj,2)) * &
1023	(zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))/snowdz_old(ji,locflag(jj,2))
1024	ELSE
1025	zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,2))/snowdz_old(ji,locflag(jj,2)) * &
1026	(zsnowzn(ji,locflag(jj,1))-zsnowzo(ji,locflag(jj,1)))
1027	ENDIF
1028
1029	DO jjj=locflag(jj,1),locflag(jj,2)-1
1030	zsnowheatn(ji,jj)=zsnowheatn(ji,jj) + &
1031	(jjj-locflag(jj,1))*snowheat(ji,jjj)
1032	ENDDO
1033
1034
1035	!snow grain
1036	zsnowgrainn(ji,jj) = snowgrain(ji,locflag(jj,1)) * &
1037	(zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1)) + &
1038	snowgrain(ji,locflag(jj,2)) * &
1039	(zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))
1040
1041	DO jjj=locflag(jj,1),locflag(jj,2)-1
1042	zsnowgrainn(ji,jj)=zsnowgrainn(ji,jj) + &
1043	(jjj-locflag(jj,1))snowgrain(ji,jjj)snowdz_old(ji,jjj)
1044	ENDDO
1045
1046	zsnowgrainn(ji,jj) = zsnowgrainn(ji,jj) / snowdz(ji,jj)
1047
1048	ENDIF
1049
1050	ENDDO
1051
1052	snowrho(ji,:) = zsnowrhon(ji,:)
1053	snowheat(ji,:) = zsnowheatn(ji,:)
1054	snowgrain(ji,:) = zsnowgrainn(ji,:)
1055
1056	ENDIF
1057
1058	! Vanishing or very thin snowpack check:
1059	! -----------------------------------------
1060	!
1061	! NOTE: ONLY for very shallow snowpacks, mix properties (homogeneous):
1062	! this avoids problems related to heat and mass exchange for
1063	! thin layers during heavy snowfall or signifigant melt: one
1064	! new/old layer can exceed the thickness of several old/new layers.
1065	! Therefore, mix (conservative):
1066	!
1067	! modified by Tao Wang
1068	IF (psnow > 0 .AND. (psnow < xsnowcritd .OR. snowdz_old(ji,1) &
1069	& .eq. 0 .OR. snowdz_old(ji,2) .eq. 0 .OR. snowdz_old(ji,3) .eq. 0)) THEN
1070	zsumheat(ji) = SUM(snowheat(ji,:))
1071	zsumswe(ji) = SUM(snowrho(ji,:)*snowdz_old(ji,:))
1072	zsumgrain(ji)= SUM(snowgrain(ji,:)*snowdz_old(ji,:))
1073	zsnowmix_delta(ji) = 1.0
1074	DO jj=1,nsnow
1075	zsnowheatn(ji,jj) = zsnowmix_delta(ji)*(zsumheat(ji)/nsnow)
1076	snowdz(ji,jj) = zsnowmix_delta(ji)*(psnow/nsnow)
1077	zsnowrhon(ji,jj) = zsnowmix_delta(ji)*(zsumswe(ji)/psnow)
1078	zsnowgrainn(ji,jj) = zsnowmix_delta(ji)*(zsumgrain(ji)/psnow)
1079	ENDDO
1080	! Update mass (density and thickness), heat and grain size:
1081	! ------------------------------------------------------------
1082	!
1083	snowrho(ji,:) = zsnowrhon(ji,:)
1084	snowheat(ji,:) = zsnowheatn(ji,:)
1085	snowgrain(ji,:) = zsnowgrainn(ji,:)
1086
1087	ENDIF
1088
1089	ENDDO
1090
1091
1092	END SUBROUTINE explicitsnow_transf
1093
1094
1095	!!
1096	!================================================================================================================================
1097	!! SUBROUTINE : explicitsnow_fall
1098	!!
1099	!>\BRIEF Computes snowfall
1100	!!
1101	!! DESCRIPTION : Computes snowfall
1102	!routine.
1103	!! RECENT CHANGE(S) : None
1104	!!
1105	!! MAIN OUTPUT VARIABLE(S): None
1106	!!
1107	!! REFERENCE(S) :
1108	!!
1109	!! FLOWCHART : None
1110	!! \n
1111	!_
1112	!================================================================================================================================
1113
1114	SUBROUTINE explicitsnow_fall(kjpindex,precip_snow,temp_air,u,v,totfrac_nobio,snowrho,snowdz,&
1115	& snowheat,snowgrain,snowtemp,psnowhmass)
1116
1117	!! 0.1 Input variables
1118
1119	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
1120	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: precip_snow !! Snow rate (SWE) (kg/m2 per dt_sechiba)
1121	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_air !! Air temperature
1122	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u,v !! Horizontal wind speed
1123	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature
1124	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: totfrac_nobio
1125	!! 0.2 Output variables
1126
1127	REAL(r_std), DIMENSION(kjpindex),INTENT(out) :: psnowhmass !! Heat content of snowfall (J/m2)
1128
1129	!! 0.3 Modified variables
1130
1131	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density profile (kg/m3)
1132	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow layer thickness profile (m)
1133	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowheat !! Snow heat content/enthalpy (J/m2)
1134	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m)
1135
1136	!! 0.4 Local variables
1137
1138	REAL(r_std), DIMENSION(kjpindex) :: rhosnew !! Snowfall density
1139	REAL(r_std), DIMENSION(kjpindex) :: dsnowfall !! Snowfall thickness (m)
1140	REAL(r_std), DIMENSION(kjpindex,nsnow) :: snowdz_old
1141	REAL(r_std), DIMENSION(kjpindex) :: snow_depth,snow_depth_old,newgrain
1142	REAL(r_std) :: snowfall_delta,speed
1143	INTEGER(i_std) :: ji,jj
1144
1145	!! 1. initialize the variables
1146
1147	snowdz_old = snowdz
1148	DO ji=1,kjpindex
1149	snow_depth(ji) = SUM(snowdz(ji,:))
1150	ENDDO
1151
1152	snow_depth_old = snow_depth
1153
1154	snowfall_delta = 0.0
1155
1156	!! 2. incorporate snowfall into snowpack
1157	DO ji = 1, kjpindex
1158
1159	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
1160
1161	! new snow fall on snowpack
1162	! NOTE: when the surface temperature is zero, it means that snowfall has no
1163	! heat content but it can be used for increasing the thickness and changing the density (maybe it is a bug)
1164	psnowhmass(ji) = 0.0
1165	IF ( (precip_snow(ji) .GT. 0.0) ) THEN
1166
1167	!calculate
1168
1169	psnowhmass(ji) = precip_snow(ji)(un-totfrac_nobio(ji)) &
1170	(xci*(snowtemp(ji,1)-tp_00)-chalfu0)
1171
1172	! Snowfall density: Following CROCUS (Pahaut 1976)
1173	!
1174	rhosnew(ji) = MAX(xrhosmin, snowfall_a_sn + snowfall_b_sn*(temp_air(ji)-tp_00)+ &
1175	snowfall_c_sn*SQRT(speed))
1176
1177	! Augment total pack depth:
1178	! Note that dsnowfall(ji) can be equal to zero if totfrac_nobio(ji)=1
1179	dsnowfall(ji) = (precip_snow(ji)*(un-totfrac_nobio(ji)))/rhosnew(ji) !snowfall thickness (m)
1180
1181	snow_depth(ji) = snow_depth(ji) + dsnowfall(ji)
1182
1183
1184	! Fresh snowfall changes the snowpack density and liquid content in uppermost layer
1185
1186	IF (dsnowfall(ji) .GT. zero) THEN
1187	snowrho(ji,1) = (snowdz(ji,1)snowrho(ji,1) + dsnowfall(ji)rhosnew(ji))/ &
1188	(snowdz(ji,1)+dsnowfall(ji))
1189
1190	snowdz(ji,1) = snowdz(ji,1) + dsnowfall(ji)
1191
1192	! Add energy of snowfall to snowpack:
1193	! Update heat content (J/m2) (therefore the snow temperature
1194	! and liquid content):
1195	!
1196	snowheat(ji,1) = snowheat(ji,1) + psnowhmass(ji)
1197	!
1198	! Incorporate snowfall grain size:
1199	!
1200	newgrain(ji) = MIN(dgrain_new_max, snow3lgrain_0d(rhosnew(ji)))
1201
1202	snowgrain(ji,1) = (snowdz_old(ji,1)snowgrain(ji,1) + dsnowfall(ji)newgrain(ji))/ &
1203	snowdz(ji,1)
1204	ENDIF
1205	ELSE
1206	dsnowfall(ji) = 0.
1207	ENDIF
1208
1209	! new snow fall on snow free surface.
1210	! we use the linearization for the new snow fall on snow-free ground
1211
1212	IF ( (dsnowfall(ji) .GT. zero) .AND. (snow_depth_old(ji) .EQ. zero) ) THEN
1213
1214	snowfall_delta = 1.0
1215
1216	DO jj=1,nsnow
1217
1218	snowdz(ji,jj) = snowfall_delta*(dsnowfall(ji)/nsnow) + &
1219	(1.0-snowfall_delta)*snowdz(ji,jj)
1220
1221	snowheat(ji,jj) = snowfall_delta*(psnowhmass(ji)/nsnow) + &
1222	(1.0-snowfall_delta)*snowheat(ji,jj)
1223
1224	snowrho(ji,jj) = snowfall_delta*rhosnew(ji) + &
1225	(1.0-snowfall_delta)*snowrho(ji,jj)
1226
1227	snowgrain(ji,jj)= snowfall_delta*newgrain(ji) + &
1228	(1.0-snowfall_delta)*snowgrain(ji,jj)
1229
1230	ENDDO
1231
1232
1233	ENDIF
1234
1235
1236	ENDDO
1237
1238	END SUBROUTINE explicitsnow_fall
1239
1240	!!
1241	!================================================================================================================================
1242	!! SUBROUTINE : explicitsnow_gone
1243	!!
1244	!>\BRIEF Check whether snow is gone
1245	!!
1246	!! DESCRIPTION : If so, set thickness (and therefore mass and heat) and liquid
1247	!! content to zero, and adjust fluxes of water, evaporation and
1248	!! heat into underlying surface.
1249	!! RECENT CHANGE(S) : None
1250	!!
1251	!! MAIN OUTPUT VARIABLE(S): None
1252	!!
1253	!! REFERENCE(S) :
1254	!!
1255	!! FLOWCHART : None
1256	!! \n
1257	!_
1258	!================================================================================================================================
1259
1260	SUBROUTINE explicitsnow_gone(kjpindex,pgflux,&
1261	snowheat,snowtemp,snowdz,snowrho,snowliq,grndflux,snowmelt)
1262
1263	!! 0.1 Input variables
1264
1265	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1266	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pgflux !! Net energy into snow pack(w/m2)
1267	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowheat !! Snow heat content (J/m^2)
1268
1269	!! 0.2 Output variables
1270
1271	!! 0.3 Modified variables
1272
1273	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature
1274	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth
1275	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
1276	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowliq !! Liquid water content
1277	REAL(r_std),DIMENSION(kjpindex), INTENT(inout) :: grndflux !! Soil/snow interface heat flux (W/m2)
1278	REAL(r_std),DIMENSION(kjpindex),INTENT(inout) :: snowmelt !! Snow melt
1279	REAL(r_std),DIMENSION(kjpindex) :: thrufal !! Water leaving snowpack(kg/m2/s)
1280
1281	!! 0.4 Local variables
1282
1283	INTEGER(i_std) :: ji,jj
1284	REAL(r_std),DIMENSION(kjpindex) :: snowgone_delta
1285	REAL(r_std),DIMENSION (kjpindex) :: totsnowheat !!snow heat content at each layer
1286	REAL(r_std),DIMENSION(kjpindex) :: snowdepth_crit
1287
1288	! first caculate total snowpack snow heat content
1289	snowgone_delta(:) = un
1290	thrufal(:)=0.0
1291	snowmelt(:)=0
1292	totsnowheat(:) = SUM(snowheat(:,:),2)
1293
1294	DO ji = 1, kjpindex
1295
1296	IF ( (SUM(snowdz(ji,:)) .GT. min_sechiba)) THEN
1297
1298	IF( pgflux(ji) >= (-totsnowheat(ji)/dt_sechiba) ) THEN
1299
1300	grndflux(ji) = pgflux(ji) + (totsnowheat(ji)/dt_sechiba)
1301
1302	thrufal(ji)=SUM(snowrho(ji,:)*snowdz(ji,:))
1303
1304	snowgone_delta(ji) = 0.0
1305
1306	snowmelt(ji) = snowmelt(ji)+thrufal(ji)
1307
1308	ENDIF
1309
1310	! update of snow state (either still present or not)
1311
1312	DO jj=1,nsnow
1313	snowdz(ji,jj) = snowdz(ji,jj) *snowgone_delta(ji)
1314	snowliq(ji,jj) = snowliq(ji,jj) *snowgone_delta(ji)
1315	snowtemp(ji,jj) = (1.0-snowgone_delta(ji))tp_00 + snowtemp(ji,jj)snowgone_delta(ji)
1316	ENDDO
1317
1318	ELSE
1319	snowdz(ji,:)=0.
1320	snowliq(ji,:)=0.
1321	snowmelt(ji)=snowmelt(ji)+SUM(snowrho(ji,:)*snowdz(ji,:))
1322	ENDIF
1323	ENDDO
1324
1325	END SUBROUTINE explicitsnow_gone
1326
1327	!================================================================================================================================
1328	!! SUBROUTINE : explicitsnow_melt_refrz
1329	!!
1330	!>\BRIEF Computes snow melt and refreezing processes within snowpack
1331	!!
1332	!! DESCRIPTION :
1333	!! RECENT CHANGE(S) : None
1334	!!
1335	!! MAIN OUTPUT VARIABLE(S): None
1336	!!
1337	!! REFERENCE(S) :
1338	!!
1339	!! FLOWCHART : None
1340	!! \n
1341	!_
1342	!================================================================================================================================
1343
1344	SUBROUTINE explicitsnow_melt_refrz(kjpindex,precip_rain,pgflux,soilcap,&
1345	snowtemp,snowdz,snowrho,snowliq,snowmelt,grndflux,temp_air)
1346
1347	!! 0.1 Input variables
1348
1349	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
1350	REAL(r_std),DIMENSION (kjpindex,nsnow) :: pcapa_snow !! Heat capacity for snow
1351	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall
1352	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: temp_air !! Air temperature
1353	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: pgflux !! Net energy into snowpack(w/m2)
1354	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilcap !! Soil heat capacity
1355
1356	!! 0.2 Output variables
1357
1358	!! 0.3 Modified variables
1359
1360	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature
1361	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth
1362	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow layer density (Kg/m^3)
1363	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowliq !! Liquid water content
1364	REAL(r_std),DIMENSION(kjpindex),INTENT(inout) :: grndflux !! Net energy input to soil
1365	REAL(r_std),DIMENSION (kjpindex), INTENT(inout) :: snowmelt !! Snowmelt
1366
1367	!! 0.4 Local variables
1368
1369	REAL(r_std),DIMENSION (kjpindex) :: meltxs !! Residual snowmelt energy applied to underlying soil
1370	REAL(r_std) :: enerin,melttot,hrain
1371	REAL(r_std),DIMENSION (nsnow) :: zsnowlwe
1372	REAL(r_std),DIMENSION (nsnow) :: flowliq
1373	REAL(r_std),DIMENSION (kjpindex) :: snowmass
1374	REAL(r_std),DIMENSION (nsnow) :: zphase !! Phase change (from ice to water) (J/m2)
1375	REAL(r_std),DIMENSION (nsnow) :: zphase2 !! Phase change (from water to ice)
1376	REAL(r_std),DIMENSION (nsnow) :: zphase3 !! Phase change related with net energy input to snowpack
1377	REAL(r_std),DIMENSION (nsnow) :: zsnowdz !! Snow layer depth
1378	REAL(r_std),DIMENSION (nsnow) :: zsnowmelt !! Snow melt (liquid water) (m)
1379	REAL(r_std),DIMENSION (nsnow) :: zsnowtemp
1380	REAL(r_std),DIMENSION (nsnow) :: zmeltxs !! Excess melt
1381	REAL(r_std),DIMENSION (nsnow) :: zwholdmax !! Maximum liquid water holding (m)
1382	REAL(r_std),DIMENSION (nsnow) :: zcmprsfact !! Compression factor due to densification from melting
1383	REAL(r_std),DIMENSION (nsnow) :: zscap !! Snow heat capacity (J/m3 K)
1384	REAL(r_std),DIMENSION (nsnow) :: zsnowliq !! (m)
1385	REAL(r_std),DIMENSION (nsnow) :: snowtemp_old
1386	REAL(r_std),DIMENSION (0:nsnow) :: zflowliqt !!(m)
1387	REAL(r_std) :: zrainfall,zpcpxs
1388	REAL(r_std) :: ztotwcap
1389	REAL(r_std),DIMENSION(kjpindex,nsnow) :: snowdz_old,snowliq_old
1390	INTEGER(i_std) :: jj,ji, iv
1391	REAL(r_std),DIMENSION(nsnow) :: snowdz_old2
1392	REAL(r_std),DIMENSION(nsnow) :: zsnowrho
1393
1394	!initialize
1395	snowdz_old = snowdz
1396	snowliq_old = snowliq
1397
1398	DO ji = 1, kjpindex
1399
1400
1401	snowmass(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
1402	IF ((snowmass(ji) .GT. min_sechiba)) THEN
1403
1404	!! 1 snow melting due to positive snowpack snow temperature
1405
1406	!! 1.0 total liquid equivalent water content of each snow layer
1407
1408	zsnowlwe(:) = snowrho(ji,:) * snowdz(ji,:)/ ph2o
1409
1410	!! 1.1 phase change (J/m2)
1411
1412	pcapa_snow(ji,:) = snowrho(ji,:)*xci
1413
1414	zphase(:) = MIN(pcapa_snow(ji,:)MAX(0.0, snowtemp(ji,:)-tp_00) &
1415	snowdz(ji,:), &
1416	MAX(0.0,zsnowlwe(:)-snowliq(ji,:))chalfu0ph2o)
1417
1418	!! 1.2 update snow liq water and temperature if melting
1419
1420	zsnowmelt(:) = zphase(:)/(chalfu0*ph2o)
1421
1422	!! 1.3 cool off the snow layer temperature due to melt
1423
1424	zsnowtemp(:) = snowtemp(ji,:) - zphase(:)/(pcapa_snow(ji,:)* snowdz(ji,:))
1425
1426	snowtemp(ji,:) = MIN(tp_00, zsnowtemp(:))
1427
1428	zmeltxs(:) = (zsnowtemp(:)-snowtemp(ji,:))pcapa_snow(ji,:)snowdz(ji,:)
1429
1430	!! 1.4 loss of snowpack depth and liquid equivalent water
1431
1432	zwholdmax(:) = snow3lhold_1d(snowrho(ji,:),snowdz(ji,:)) ! 1 dimension
1433
1434	zcmprsfact(:) = (zsnowlwe(:)-MIN(snowliq(ji,:)+zsnowmelt(:),zwholdmax(:)))/ &
1435	(zsnowlwe(:)-MIN(snowliq(ji,:),zwholdmax(:)))
1436
1437	snowdz(ji,:) = snowdz(ji,:)*zcmprsfact(:)
1438
1439	snowrho(ji,:) = zsnowlwe(:)*ph2o/snowdz(ji,:)
1440
1441	snowliq(ji,:) = snowliq(ji,:) + zsnowmelt(:)
1442
1443
1444	!! 2 snow refreezing process
1445	!! 2.1 freeze liquid water in any layer
1446	zscap(:) = snowrho(ji,:)*xci !J K-1 m-3
1447	zphase2(:) = MIN(zscap(:)* &
1448	MAX(0.0, tp_00 - snowtemp(ji,:))*snowdz(ji,:), &
1449	snowliq(ji,:)chalfu0ph2o)
1450
1451	! warm layer and reduce liquid if freezing occurs
1452	zsnowdz(:) = MAX(xsnowdmin/nsnow, snowdz(ji,:))
1453	snowtemp_old(:) = snowtemp(ji,:)
1454	snowtemp(ji,:) = snowtemp(ji,:) + zphase2(:)/(zscap(:)*zsnowdz(:))
1455
1456	! Reduce liquid portion if freezing occurs:
1457	snowliq(ji,:) = snowliq(ji,:) - ( (snowtemp(ji,:)-snowtemp_old(:))* &
1458	zscap(:)zsnowdz(:)/(chalfu0ph2o) )
1459	snowliq(ji,:) = MAX(snowliq(ji,:), 0.0)
1460	!! 3. thickness change due to snowmelt in excess of holding capacity
1461	zwholdmax(:) = snow3lhold_1d(snowrho(ji,:),snowdz(ji,:)) ! 1 dimension
1462	flowliq(:) = MAX(0.,(snowliq(ji,:)-zwholdmax(:)))
1463	snowliq(ji,:) = snowliq(ji,:) - flowliq(:)
1464	snowdz(ji,:) = snowdz(ji,:) - flowliq(:)*ph2o/snowrho(ji,:)
1465	! to prevent possible very small negative values (machine
1466	! prescision as snow vanishes
1467	snowdz(ji,:) = MAX(0.0, snowdz(ji,:))
1468
1469	!! 4. liquid water flow
1470	ztotwcap = SUM(zwholdmax(:))
1471	! Rain entering snow (m):
1472	!ORCHIDEE has assumed that all precipitation entering snow has
1473	!left the snowpack and finally become runoff in hydrol_soil. Here we just put zrainfall as 0.0
1474	!!zrainfall = precip_rain(ji)/ph2o ! rainfall (m)
1475	zrainfall = 0.0
1476
1477	zflowliqt(0) = MIN(zrainfall,ztotwcap)
1478	! Rain assumed to directly pass through the pack to runoff (m):
1479	zpcpxs = zrainfall - zflowliqt(0)
1480	!
1481	DO jj=1,nsnow
1482	zflowliqt(jj) = flowliq(jj)
1483	ENDDO
1484
1485	! translated into a density increase:
1486	flowliq(:) = 0.0 ! clear this array for work
1487	zsnowliq(:) = snowliq(ji,:) ! reset liquid water content
1488	!
1489
1490	DO jj=1,nsnow
1491	snowliq(ji,jj) = snowliq(ji,jj) + zflowliqt(jj-1)
1492	flowliq(jj) = MAX(0.0, (snowliq(ji,jj)-zwholdmax(jj)))
1493	snowliq(ji,jj) = snowliq(ji,jj) - flowliq(jj)
1494
1495	!Modified by Tao Wang based on previous ISBA-ES scheme
1496	snowrho(ji,jj) = snowrho(ji,jj) + (snowliq(ji,jj) - zsnowliq(jj))* &
1497	& ph2o/MAX(xsnowdmin/nsnow,snowdz(ji,jj))
1498
1499	zflowliqt(jj) = zflowliqt(jj) + flowliq(jj)
1500	ENDDO
1501
1502	snowmelt(ji) = snowmelt(ji) + (zflowliqt(nsnow) + zpcpxs) * ph2o
1503
1504	! excess heat from melting, using it to warm underlying ground to conserve energy
1505	meltxs(ji) = SUM(zmeltxs(:))/dt_sechiba ! (W/m2)
1506
1507	! energy flux into the soil
1508	grndflux(ji) = grndflux(ji) + meltxs(ji)
1509
1510	ELSE
1511	snowdz(ji,:)=0.
1512	snowliq(ji,:)=0.
1513	snowmelt(ji)=snowmelt(ji)+SUM(snowrho(ji,:)*snowdz(ji,:))
1514	ENDIF
1515
1516	ENDDO
1517
1518	END SUBROUTINE explicitsnow_melt_refrz
1519
1520	!================================================================================================================================
1521	!! SUBROUTINE : explicitsnow_levels
1522	!!
1523	!>\BRIEF Computes snow discretization based on given total snow depth
1524	!!
1525	!! DESCRIPTION :
1526	!! RECENT CHANGE(S) : None
1527	!!
1528	!! MAIN OUTPUT VARIABLE(S): None
1529	!!
1530	!! REFERENCE(S) :
1531	!!
1532	!! FLOWCHART : None
1533	!! \n
1534	!_
1535	!================================================================================================================================
1536	SUBROUTINE explicitsnow_levels( kjpindex,snow_thick, snowdz)
1537
1538	!! 0.1 Input variables
1539	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1540	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow_thick !! Total snow depth
1541
1542	!! 0.2 Output variables
1543
1544	!! 0.3 Modified variables
1545
1546	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth
1547
1548	!! 0.4 Local variables
1549
1550	INTEGER(i_std) :: il,it,ji
1551	INTEGER(i_std) :: i,j, ix
1552	REAL(r_std), DIMENSION(kjpindex) :: xi, xf
1553	REAL(r_std), PARAMETER, DIMENSION(3) :: ZSGCOEF1 = (/0.25, 0.50, 0.25/) !! Snow grid parameters
1554	REAL(r_std), PARAMETER, DIMENSION(2) :: ZSGCOEF2 = (/0.05, 0.34/) !! Snow grid parameters
1555	REAL(r_std), PARAMETER :: ZSNOWTRANS = 0.20 !! Minimum total snow depth at which surface layer thickness is constant (m)
1556	REAL(r_std), PARAMETER :: XSNOWCRITD = 0.03 !! (m)
1557
1558	DO ji=1,kjpindex
1559	IF ( snow_thick(ji) .LE. (XSNOWCRITD+0.01)) THEN
1560
1561	snowdz(ji,1) = MIN(0.01, snow_thick(ji)/nsnow)
1562	snowdz(ji,3) = MIN(0.01, snow_thick(ji)/nsnow)
1563	snowdz(ji,2) = snow_thick(ji) - snowdz(ji,1) - snowdz(ji,3)
1564
1565	ENDIF
1566	ENDDO
1567
1568	WHERE ( snow_thick(:) .LE. ZSNOWTRANS .AND. &
1569	snow_thick(:) .GT. (XSNOWCRITD+0.01) )
1570	!
1571	snowdz(:,1) = snow_thick(:)*ZSGCOEF1(1)
1572	snowdz(:,2) = snow_thick(:)*ZSGCOEF1(2)
1573	snowdz(:,3) = snow_thick(:)*ZSGCOEF1(3)
1574	!
1575	END WHERE
1576
1577	DO ji = 1,kjpindex
1578	IF (snow_thick(ji) .GT. ZSNOWTRANS) THEN
1579	snowdz(ji,1) = ZSGCOEF2(1)
1580	snowdz(ji,2) = (snow_thick(ji)-ZSGCOEF2(1))*ZSGCOEF2(2) + ZSGCOEF2(1)
1581	!When using simple finite differences, limit the thickness
1582	!factor between the top and 2nd layers to at most 10
1583	snowdz(ji,2) = MIN(10*ZSGCOEF2(1), snowdz(ji,2) )
1584	snowdz(ji,3) = snow_thick(ji) - snowdz(ji,2) - snowdz(ji,1)
1585
1586	ENDIF
1587	ENDDO
1588
1589
1590	END SUBROUTINE explicitsnow_levels
1591
1592	!!
1593	!================================================================================================================================
1594	!! SUBROUTINE : explicitsnow_profile
1595	!!
1596	!>\BRIEF
1597	!!
1598	!! DESCRIPTION : In this routine solves the numerical soil thermal scheme, ie calculates the new soil temperature profile.
1599	!!
1600	!! RECENT CHANGE(S) : None
1601	!!
1602	!! MAIN OUTPUT VARIABLE(S):
1603	!!
1604	!! REFERENCE(S) :
1605	!!
1606	!! FLOWCHART : None
1607	!! \n
1608	!_
1609	!================================================================================================================================
1610	SUBROUTINE explicitsnow_profile (kjpindex, cgrnd_snow,dgrnd_snow,lambda_snow,temp_sol_new, snowtemp,snowdz,temp_sol_add)
1611
1612	!! 0. Variables and parameter declaration
1613
1614	!! 0.1 Input variables
1615
1616	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1617	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! skin temperature
1618	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
1619	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
1620	REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
1621	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow layer thickness
1622	REAL(r_std), DIMENSION (kjpindex),INTENT(inout) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
1623
1624	!! 0.3 Modified variable
1625
1626	!! 0.2 Output variables
1627
1628	!! 0.3 Modified variables
1629	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout) :: snowtemp
1630
1631	!! 0.4 Local variables
1632
1633	INTEGER(i_std) :: ji, jg
1634	!_
1635	!================================================================================================================================
1636	!! 1. Computes the snow temperatures ptn.
1637	DO ji = 1,kjpindex
1638	IF (SUM(snowdz(ji,:)) .GT. 0) THEN
1639	snowtemp(ji,1) = (lambda_snow(ji) * cgrnd_snow(ji,1) + (temp_sol_new(ji)+temp_sol_add(ji))) / &
1640	(lambda_snow(ji) * (un - dgrnd_snow(ji,1)) + un)
1641	temp_sol_add(ji) = zero
1642	DO jg = 1,nsnow-1
1643	snowtemp(ji,jg+1) = cgrnd_snow(ji,jg) + dgrnd_snow(ji,jg) * snowtemp(ji,jg)
1644	ENDDO
1645	ENDIF
1646	ENDDO
1647	IF (printlev>=3) WRITE (numout,*) ' explicitsnow_profile done '
1648
1649	END SUBROUTINE explicitsnow_profile
1650
1651	END MODULE explicitsnow

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

Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/explicitsnow.f90 @ 7595

Download in other formats: