Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

explicitsnow.f90 @ 7693

Last change on this file since 7693 was 6830, checked in by nicolas.vuichard, 4 years ago
update growth respiration and bud fix ticket 642 on water balance
Property svn:keywords set to `Date Revision HeadURL`
File size: 73.4 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 [m]
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	! CALL explicitsnow_levels(kjpindex,snow_depth_tmp, snowdz)
541	DO ji=1,kjpindex
542	IF (snowdz(ji,1) .NE. 0 .AND. snowdz(ji,2) .EQ. 0) THEN
543	!Condensation on first layer, see Ticket #642.
544	snow_depth_tmp(ji) = SUM(snowdz(ji,:))
545	snowdz(ji,:)=snow_depth_tmp(ji)/3
546	snowrho(ji,2:3)=snowrho(ji,1)
547	ENDIF
548	ENDDO
549
550	! Write output variables
551
552	! Add XIOS default value where no snow
553	DO ji=1,kjpindex
554	IF (snow(ji) .GT. zero) THEN
555	! Output for snowliq and snowliqtot change unit from m to kg/m2 and are multiplied by contfrac
556	! to take into acount the portion form land fraction only
557	snowliq_diag(ji,:) = snowliq(ji,:) * 1000 * contfrac(ji)
558	snowliqtot_diag(ji) = SUM(snowliq_diag(ji,:)) * 1000 * contfrac(ji)
559	snowrho_diag(ji,:) = snowrho(ji,:)
560	snowheat_diag(ji,:) = snowheat(ji,:)
561	snowgrain_diag(ji,:) = snowgrain(ji,:)
562	ELSE
563	snowliq_diag(ji,:) = xios_default_val
564	snowliqtot_diag(ji) = xios_default_val
565	snowrho_diag(ji,:) = xios_default_val
566	snowheat_diag(ji,:) = xios_default_val
567	snowgrain_diag(ji,:) = xios_default_val
568	END IF
569	END DO
570
571	CALL xios_orchidee_send_field("snowliq",snowliq_diag)
572	CALL xios_orchidee_send_field("snowliqtot", snowliqtot_diag)
573	CALL xios_orchidee_send_field("snowrho",snowrho_diag)
574	CALL xios_orchidee_send_field("snowheat",snowheat_diag)
575	CALL xios_orchidee_send_field("snowgrain",snowgrain_diag)
576	CALL xios_orchidee_send_field("snowmelt_from_maxmass",snowmelt_from_maxmass)
577	CALL xios_orchidee_send_field("soilcap",soilcap)
578	CALL xios_orchidee_send_field("del_snowheattot",del_snowheattot)
579
580	END SUBROUTINE explicitsnow_main
581
582
583	!================================================================================================================================
584	!! SUBROUTINE : explicitsnow_finalize
585	!!
586	!>\BRIEF Write variables for explictsnow module to restart file
587	!!
588	!! DESCRIPTION : Write variables for explictsnow module to restart file
589	!!
590	!! \n
591	!_
592	!================================================================================================================================
593	SUBROUTINE explicitsnow_finalize ( kjit, kjpindex, rest_id, snowrho, &
594	snowtemp, snowdz, snowheat, snowgrain)
595
596	!! 0.1 Input variables
597	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
598	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
599	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
600	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density (Kg/m^3)
601	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K)
602	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow layer thickness [m]
603	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowheat !! Snow heat content (J/m^2)
604	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowgrain !! Snow grainsize (m)
605
606
607	!! 1. Write to restart file
608	CALL restput_p(rest_id, 'snowrho', nbp_glo, nsnow, 1, kjit, snowrho, 'scatter', nbp_glo, index_g)
609	CALL restput_p(rest_id, 'snowtemp', nbp_glo, nsnow, 1, kjit, snowtemp, 'scatter', nbp_glo, index_g)
610	CALL restput_p(rest_id, 'snowdz', nbp_glo, nsnow, 1, kjit, snowdz, 'scatter', nbp_glo, index_g)
611	CALL restput_p(rest_id, 'snowheat', nbp_glo, nsnow, 1, kjit, snowheat, 'scatter', nbp_glo, index_g)
612	CALL restput_p(rest_id, 'snowgrain', nbp_glo, nsnow, 1, kjit, snowgrain, 'scatter', nbp_glo, index_g)
613
614	END SUBROUTINE explicitsnow_finalize
615
616
617	!================================================================================================================================
618	!! SUBROUTINE : explicitsnow_grain
619	!!
620	!>\BRIEF Compute evolution of snow grain size
621	!!
622	!! DESCRIPTION :
623	!!
624	!! RECENT CHANGE(S) : None
625	!!
626	!! MAIN OUTPUT VARIABLE(S): None
627	!!
628	!! REFERENCE(S) : R. Jordan (1991)
629	!!
630	!! FLOWCHART : None
631	!! \n
632	!_
633	!================================================================================================================================
634
635
636	SUBROUTINE explicitsnow_grain(kjpindex,snowliq,snowdz,gtemp,snowtemp,pb,snowgrain)
637
638	!! 0.1 Input variables
639	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
640	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowliq !! Liquid water content
641	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth (m)
642	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: gtemp !! First soil layer temperature
643	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature (K)
644	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: pb !! Surface pressure (hpa)
645
646	!! 0.2 Output variables
647
648	!! 0.3 Modified variables
649
650	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m)
651
652	!! 0.4 Local variables
653	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowdz,zdz,ztheta
654	REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: ztemp,zdiff,ztgrad,zthetaa,zfrac,&
655	zexpo,zckt_liq,zckt_ice,zckt
656	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zrhomin,zgrainmin
657	INTEGER(i_std) :: ji,jj
658
659	!! 0.5 Local parameters
660	REAL(r_std), PARAMETER :: ztheta_crit = 0.02 !! m3 m-3
661	REAL(r_std), PARAMETER :: zc1_ice = 8.047E+9 !! kg m-3 K
662	REAL(r_std), PARAMETER :: zc1_liq = 5.726E+8 !! kg m-3 K
663	REAL(r_std), PARAMETER :: zdeos = 0.92E-4 !! effective diffusion
664	!! coef for water vapor in snow
665	!! at 0C and 1000 mb (m2 s-1)
666	REAL(r_std), PARAMETER :: zg1 = 5.0E-7 !! m4 kg-1
667	REAL(r_std), PARAMETER :: zg2 = 4.0E-12 !! m2 s-1
668	REAL(r_std), PARAMETER :: ztheta_w = 0.05 !! m3 m-3
669	REAL(r_std), PARAMETER :: ztheta_crit_w = 0.14 !! m3 m-3
670	REAL(r_std), PARAMETER :: zdzmin = 0.01 !! m : minimum thickness
671	!! for thermal gradient evaluation:
672	!! to prevent excessive gradients
673	!! for vanishingly thin snowpacks.
674	REAL(r_std), PARAMETER :: xp00=1.E5
675	!! 1. initialize
676
677	DO ji=1,kjpindex
678
679
680	zsnowdz(ji,:) = MAX(xsnowdmin/nsnow, snowdz(ji,:))
681
682	DO jj=1,nsnow-1
683	zdz(ji,jj) = zsnowdz(ji,jj) + zsnowdz(ji,jj+1)
684	ENDDO
685	zdz(ji,nsnow) = zsnowdz(ji,nsnow)
686
687	! compute interface average volumetric water content (m3 m-3):
688	! first, layer avg VWC:
689	!
690	ztheta(ji,:) = snowliq(ji,:)/MAX(xsnowdmin, zsnowdz(ji,:))
691
692	! at interfaces:
693	zthetaa(ji,0) = ztheta(ji,1)
694	DO jj=1,nsnow-1
695	zthetaa(ji,jj) = (zsnowdz(ji,jj) *ztheta(ji,jj) + &
696	zsnowdz(ji,jj+1)*ztheta(ji,jj+1))/zdz(ji,jj)
697	ENDDO
698	zthetaa(ji,nsnow) = ztheta(ji,nsnow)
699	! compute interface average temperatures (K):
700	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
701	!
702	ztemp(ji,0) = snowtemp(ji,1)
703	DO jj=1,nsnow-1
704	ztemp(ji,jj) = (zsnowdz(ji,jj) *snowtemp(ji,jj) + &
705	zsnowdz(ji,jj+1)*snowtemp(ji,jj+1))/zdz(ji,jj)
706	ENDDO
707	ztemp(ji,nsnow) = snowtemp(ji,nsnow)
708	!
709	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
710	! compute variation of saturation vapor pressure with temperature
711	! for solid and liquid phases:
712	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
713	zexpo(ji,:) = chalsu0/(xrv*ztemp(ji,:))
714	zckt_ice(ji,:) = (zc1_ice/ztemp(ji,:)*2)(zexpo(ji,:) - 1.0)*EXP(-zexpo(ji,:))
715	!
716	zexpo(ji,:) = chalev0/(xrv*ztemp(ji,:))
717	zckt_liq(ji,:) = (zc1_liq/ztemp(ji,:)*2)(zexpo(ji,:) - 1.0)*EXP(-zexpo(ji,:))
718	!
719	! compute the weighted ice/liquid total variation (N m-2 K):
720	!
721	zfrac(ji,:) = MIN(1.0, zthetaa(ji,:)/ztheta_crit)
722	zckt(ji,:) = zfrac(ji,:)zckt_liq(ji,:) + (1.0 - zfrac(ji,:))zckt_ice(ji,:)
723
724	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
725	! Compute effective diffusion coefficient (m2 s-1):
726	! -diffudivity relative to a reference diffusion at 1000 mb and freezing point
727	! multiplied by phase energy coefficient
728	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
729	!
730	DO jj=0,nsnow
731	zdiff(ji,jj) = zdeos(xp00/(pb(ji)100.))((ztemp(ji,jj)/tp_00)6)zckt(ji,jj)
732	ENDDO
733
734	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
735	! Temperature gradient (K m-1):
736
737	ztgrad(ji,0) = 0.0 ! uppermost layer-mean and surface T's are assumed to be equal
738	DO jj=1,nsnow-1
739	ztgrad(ji,jj) = 2*(snowtemp(ji,jj)-snowtemp(ji,jj+1))/MAX(zdzmin, zdz(ji,jj))
740	ENDDO
741	!
742	! assume at base of snow, temperature is in equilibrium with soil
743	! (but obviously must be at or below freezing point!)
744	!
745	ztgrad(ji,nsnow) = 2*(snowtemp(ji,nsnow) - MIN(tp_00, gtemp(ji)))/MAX(zdzmin, zdz(ji,nsnow))
746	! prognostic grain size (m) equation:
747	!-------------------------------------------------------------------
748	! first compute the minimum grain size (m):
749	!
750	zrhomin(ji,:) = xrhosmin
751	zgrainmin(ji,:) = snow3lgrain_1d(zrhomin(ji,:))
752
753	! dry snow:
754	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
755	!
756	DO jj=1,nsnow
757
758	IF(ztheta(ji,jj) == 0.0) THEN
759
760	! only allow growth due to vapor flux INTO layer:
761	! aab add sublimation(only condensation) as upper BC...?
762
763	snowgrain(ji,jj) = snowgrain(ji,jj) + &
764	(dt_sechibazg1/MAX(zgrainmin(ji,jj),snowgrain(ji,jj))) &
765	( zdiff(ji,jj-1)*MAX(0.0,ztgrad(ji,jj-1)) - &
766	zdiff(ji,jj) *MIN(0.0,ztgrad(ji,jj)) )
767	ELSE
768
769	! wet snow
770	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
771	!
772	snowgrain(ji,jj) = snowgrain(ji,jj) + &
773	(dt_sechibazg2/MAX(zgrainmin(ji,jj),snowgrain(ji,jj))) &
774	MIN(ztheta_crit_w, ztheta(ji,jj) + ztheta_w)
775	END IF
776
777	ENDDO
778
779
780	ENDDO
781
782
783	END SUBROUTINE explicitsnow_grain
784
785	!!
786	!================================================================================================================================
787	!! SUBROUTINE : explicitsnow_compactn
788	!!
789	!>\BRIEF Compute Compaction/Settling
790	!!
791	!! DESCRIPTION :
792	!! Snow compaction due to overburden and settling.
793	!! Mass is unchanged: layer thickness is reduced
794	!! in proportion to density increases. Method
795	!! of Anderson (1976): see Loth and Graf, 1993,
796	!! J. of Geophys. Res., 98, 10,451-10,464.
797	!!
798	!! RECENT CHANGE(S) : None
799	!!
800	!! MAIN OUTPUT VARIABLE(S): None
801	!!
802	!! REFERENCE(S) : Loth and Graf (1993), Mellor (1964) and Anderson (1976)
803	!!
804	!! FLOWCHART : None
805	!! \n
806	!_
807	!================================================================================================================================
808
809
810	SUBROUTINE explicitsnow_compactn(kjpindex,snowtemp,snowrho,snowdz)
811
812	!! 0.1 Input variables
813
814	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
815	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature
816
817	!! 0.2 Output variables
818
819	!! 0.3 Modified variables
820
821	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
822	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth [m]
823
824	!! 0.4 Local variables
825
826	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zwsnowdz,zsmass,zsnowrho2,zviscocity,zsettle
827	REAL(r_std),DIMENSION(kjpindex) :: zsmassc !! cummulative snow mass (kg/m2)
828	REAL(r_std),DIMENSION(kjpindex) :: snowdepth_crit
829	INTEGER(i_std) :: ji,jj
830
831	!! 1. initialize
832
833	zsnowrho2 = snowrho
834	zsettle(:,:) = ZSNOWCMPCT_ACM
835	zviscocity(:,:) = ZSNOWCMPCT_V0
836
837	!! 2. Calculating Cumulative snow mass (kg/m2):
838
839	DO ji=1, kjpindex
840
841
842	IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
843
844	zwsnowdz(ji,:)= snowdz(ji,:)*snowrho(ji,:)
845
846	zsmassc (ji)= 0.0
847
848	DO jj=1,nsnow
849	zsmass(ji,jj) = zsmassc(ji) + zwsnowdz(ji,jj)
850	zsmassc(ji) = zsmassc(ji) + zwsnowdz(ji,jj)
851	ENDDO
852
853
854	!! 3. Computing compaction/Settling
855	! ----------------------
856	! Compaction/settling if density below upper limit
857	! (compaction is generally quite small above ~ 500 kg m-3):
858	!
859	DO jj=1,nsnow
860	IF (snowrho(ji,jj) .LT. xrhosmax) THEN
861
862	!
863	! First calculate settling due to freshly fallen snow: (NOTE:bug here for the snow temperature profile)
864	!
865
866	zsettle(ji,jj) = ZSNOWCMPCT_ACM*EXP( &
867	-ZSNOWCMPCT_BCM*(tp_00-MIN(tp_00,snowtemp(ji,jj))) &
868	-ZSNOWCMPCT_CCM*MAX(0.0, &
869	snowrho(ji,jj)-ZSNOWCMPCT_RHOD))
870	!
871	! Snow viscocity:
872	!
873	zviscocity(ji,jj) = ZSNOWCMPCT_V0EXP( ZSNOWCMPCT_VT(tp_00-MIN(tp_00,snowtemp(ji,jj))) + &
874	ZSNOWCMPCT_VR*snowrho(ji,jj) )
875
876	! Calculate snow density: compaction from weight/over-burden
877	! Anderson 1976 method:
878	zsnowrho2(ji,jj) = snowrho(ji,jj) + snowrho(ji,jj)dt_sechiba( &
879	(cte_grav*zsmass(ji,jj)/zviscocity(ji,jj)) &
880	+ zsettle(ji,jj) )
881	! Conserve mass by decreasing grid thicknesses in response
882	! to density increases
883	!
884	snowdz(ji,jj) = snowdz(ji,jj)*(snowrho(ji,jj)/zsnowrho2(ji,jj))
885
886	ENDIF
887	ENDDO
888
889	! Update density (kg m-3):
890	snowrho(ji,:) = zsnowrho2(ji,:)
891
892	ENDIF
893
894	ENDDO
895
896	END SUBROUTINE explicitsnow_compactn
897
898	!!
899	!================================================================================================================================
900	!! SUBROUTINE : explicitsnow_transf
901	!!
902	!>\BRIEF Computing snow mass and heat redistribution due to grid thickness configuration resetting
903	!!
904	!! DESCRIPTION : Snow mass and heat redistibution due to grid thickness
905	!! configuration resetting. Total mass and heat content
906	!! of the overall snowpack unchanged/conserved within this routine.
907	!! RECENT CHANGE(S) : None
908	!!
909	!! MAIN OUTPUT VARIABLE(S): None
910	!!
911	!! REFERENCE(S) :
912	!!
913	!! FLOWCHART : None
914	!! \n
915	!_
916	!================================================================================================================================
917
918	SUBROUTINE explicitsnow_transf(kjpindex,snowdz_old,snowdz,snowrho,snowheat,snowgrain)
919
920	!! 0.1 Input variables
921
922	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
923	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowdz_old !! Snow depth at the previous time step [m]
924
925	!! 0.2 Output variables
926
927	!! 0.3 Modified variables
928
929	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
930	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m)
931	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth (m)
932	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowheat !! Snow heat content/enthalpy (J/m2)
933
934	!! 0.4 Local varibles
935
936	REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: zsnowzo
937	REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: zsnowzn
938	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowddz
939	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zdelta
940	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowrhon,zsnowheatn,zsnowgrainn
941	REAL(r_std),DIMENSION(kjpindex) :: zsnowmix_delta
942	REAL(r_std),DIMENSION(kjpindex) :: zsumheat, zsumswe, zsumgrain
943	INTEGER(i_std),DIMENSION(nsnow,2) :: locflag
944	REAL(r_std) :: psnow
945	INTEGER(i_std) :: ji,jj,jjj
946
947	! Initialization
948	zsumheat(:) = 0.0
949	zsumswe(:) = 0.0
950	zsumgrain(:) = 0.0
951	zsnowmix_delta(:) = 0.0
952	locflag(:,:) = 0
953
954	DO ji=1, kjpindex
955
956
957	psnow = SUM(snowdz(ji,:))
958
959	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
960	!
961	zsnowzo(ji,0) = 0.
962	zsnowzn(ji,0) = 0.
963	zsnowzo(ji,1) = snowdz_old(ji,1)
964	zsnowzn(ji,1) = snowdz(ji,1)
965	!
966	DO jj=2,nsnow
967	zsnowzo(ji,jj) = zsnowzo(ji,jj-1) + snowdz_old(ji,jj)
968	zsnowzn(ji,jj) = zsnowzn(ji,jj-1) + snowdz(ji,jj)
969	ENDDO
970
971	DO jj=1,nsnow
972
973	IF (jj .EQ. 1) THEN
974	locflag(jj,1)=1
975	ELSE
976	DO jjj=nsnow,1,-1
977
978	!upper bound of the snow layer
979	IF (zsnowzn(ji,jj-1) .LE. zsnowzo(ji,jjj)) THEN
980	locflag(jj,1) = jjj
981	ENDIF
982	ENDDO
983	ENDIF
984
985	IF (jj .EQ. nsnow) THEN
986	locflag(jj,2)=nsnow
987	ELSE
988	DO jjj=nsnow,1,-1
989
990	!lower bound of the snow layer
991	IF (zsnowzn(ji,jj) .LE. zsnowzo(ji,jjj)) THEN
992	locflag(jj,2) = jjj
993	ENDIF
994	ENDDO
995	ENDIF
996
997	!to interpolate
998	! when heavy snow occurred
999	IF (locflag(jj,1) .EQ. locflag(jj,2)) THEN
1000
1001	zsnowrhon(ji,jj) = snowrho(ji,locflag(jj,1))
1002
1003	zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,1))*snowdz(ji,jj)/snowdz_old(ji,locflag(jj,1))
1004
1005	zsnowgrainn(ji,jj) = snowgrain(ji,locflag(jj,1))
1006
1007	ELSE
1008
1009	!snow density
1010	zsnowrhon(ji,jj) = snowrho(ji,locflag(jj,1)) * &
1011	(zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1)) + &
1012	snowrho(ji,locflag(jj,2)) * &
1013	(zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))
1014
1015	DO jjj=locflag(jj,1),locflag(jj,2)-1
1016	zsnowrhon(ji,jj)=zsnowrhon(ji,jj) + &
1017	(jjj-locflag(jj,1))snowrho(ji,jjj)snowdz_old(ji,jjj)
1018	ENDDO
1019
1020	zsnowrhon(ji,jj) = zsnowrhon(ji,jj) / snowdz(ji,jj)
1021
1022
1023	!snow heat
1024	IF (snowdz_old(ji,locflag(jj,1)) .GT. 0.0) THEN
1025
1026	zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,1)) * &
1027	(zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1))/snowdz_old(ji,locflag(jj,1)) + &
1028	snowheat(ji,locflag(jj,2)) * &
1029	(zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))/snowdz_old(ji,locflag(jj,2))
1030	ELSE
1031	zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,2))/snowdz_old(ji,locflag(jj,2)) * &
1032	(zsnowzn(ji,locflag(jj,1))-zsnowzo(ji,locflag(jj,1)))
1033	ENDIF
1034
1035	DO jjj=locflag(jj,1),locflag(jj,2)-1
1036	zsnowheatn(ji,jj)=zsnowheatn(ji,jj) + &
1037	(jjj-locflag(jj,1))*snowheat(ji,jjj)
1038	ENDDO
1039
1040
1041	!snow grain
1042	zsnowgrainn(ji,jj) = snowgrain(ji,locflag(jj,1)) * &
1043	(zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1)) + &
1044	snowgrain(ji,locflag(jj,2)) * &
1045	(zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))
1046
1047	DO jjj=locflag(jj,1),locflag(jj,2)-1
1048	zsnowgrainn(ji,jj)=zsnowgrainn(ji,jj) + &
1049	(jjj-locflag(jj,1))snowgrain(ji,jjj)snowdz_old(ji,jjj)
1050	ENDDO
1051
1052	zsnowgrainn(ji,jj) = zsnowgrainn(ji,jj) / snowdz(ji,jj)
1053
1054	ENDIF
1055
1056	ENDDO
1057
1058	snowrho(ji,:) = zsnowrhon(ji,:)
1059	snowheat(ji,:) = zsnowheatn(ji,:)
1060	snowgrain(ji,:) = zsnowgrainn(ji,:)
1061
1062	ENDIF
1063
1064	! Vanishing or very thin snowpack check:
1065	! -----------------------------------------
1066	!
1067	! NOTE: ONLY for very shallow snowpacks, mix properties (homogeneous):
1068	! this avoids problems related to heat and mass exchange for
1069	! thin layers during heavy snowfall or signifigant melt: one
1070	! new/old layer can exceed the thickness of several old/new layers.
1071	! Therefore, mix (conservative):
1072	!
1073	! modified by Tao Wang
1074	IF (psnow > 0 .AND. (psnow < xsnowcritd .OR. snowdz_old(ji,1) &
1075	& .eq. 0 .OR. snowdz_old(ji,2) .eq. 0 .OR. snowdz_old(ji,3) .eq. 0)) THEN
1076	zsumheat(ji) = SUM(snowheat(ji,:))
1077	zsumswe(ji) = SUM(snowrho(ji,:)*snowdz_old(ji,:))
1078	zsumgrain(ji)= SUM(snowgrain(ji,:)*snowdz_old(ji,:))
1079	zsnowmix_delta(ji) = 1.0
1080	DO jj=1,nsnow
1081	zsnowheatn(ji,jj) = zsnowmix_delta(ji)*(zsumheat(ji)/nsnow)
1082	snowdz(ji,jj) = zsnowmix_delta(ji)*(psnow/nsnow)
1083	zsnowrhon(ji,jj) = zsnowmix_delta(ji)*(zsumswe(ji)/psnow)
1084	zsnowgrainn(ji,jj) = zsnowmix_delta(ji)*(zsumgrain(ji)/psnow)
1085	ENDDO
1086	! Update mass (density and thickness), heat and grain size:
1087	! ------------------------------------------------------------
1088	!
1089	snowrho(ji,:) = zsnowrhon(ji,:)
1090	snowheat(ji,:) = zsnowheatn(ji,:)
1091	snowgrain(ji,:) = zsnowgrainn(ji,:)
1092
1093	ENDIF
1094
1095	ENDDO
1096
1097
1098	END SUBROUTINE explicitsnow_transf
1099
1100
1101	!!
1102	!================================================================================================================================
1103	!! SUBROUTINE : explicitsnow_fall
1104	!!
1105	!>\BRIEF Computes snowfall
1106	!!
1107	!! DESCRIPTION : Computes snowfall
1108	!routine.
1109	!! RECENT CHANGE(S) : None
1110	!!
1111	!! MAIN OUTPUT VARIABLE(S): None
1112	!!
1113	!! REFERENCE(S) :
1114	!!
1115	!! FLOWCHART : None
1116	!! \n
1117	!_
1118	!================================================================================================================================
1119
1120	SUBROUTINE explicitsnow_fall(kjpindex,precip_snow,temp_air,u,v,totfrac_nobio,snowrho,snowdz,&
1121	& snowheat,snowgrain,snowtemp,psnowhmass)
1122
1123	!! 0.1 Input variables
1124
1125	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
1126	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: precip_snow !! Snow rate (SWE) (kg/m2 per dt_sechiba)
1127	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_air !! Air temperature
1128	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u,v !! Horizontal wind speed
1129	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature
1130	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: totfrac_nobio
1131	!! 0.2 Output variables
1132
1133	REAL(r_std), DIMENSION(kjpindex),INTENT(out) :: psnowhmass !! Heat content of snowfall (J/m2)
1134
1135	!! 0.3 Modified variables
1136
1137	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density profile (kg/m3)
1138	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow layer thickness profile (m)
1139	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowheat !! Snow heat content/enthalpy (J/m2)
1140	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m)
1141
1142	!! 0.4 Local variables
1143
1144	REAL(r_std), DIMENSION(kjpindex) :: rhosnew !! Snowfall density
1145	REAL(r_std), DIMENSION(kjpindex) :: dsnowfall !! Snowfall thickness (m)
1146	REAL(r_std), DIMENSION(kjpindex,nsnow) :: snowdz_old
1147	REAL(r_std), DIMENSION(kjpindex) :: snow_depth,snow_depth_old,newgrain
1148	REAL(r_std) :: snowfall_delta,speed
1149	INTEGER(i_std) :: ji,jj
1150
1151	!! 1. initialize the variables
1152
1153	snowdz_old = snowdz
1154	DO ji=1,kjpindex
1155	snow_depth(ji) = SUM(snowdz(ji,:))
1156	ENDDO
1157
1158	snow_depth_old = snow_depth
1159
1160	snowfall_delta = 0.0
1161
1162	!! 2. incorporate snowfall into snowpack
1163	DO ji = 1, kjpindex
1164
1165	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
1166
1167	! new snow fall on snowpack
1168	! NOTE: when the surface temperature is zero, it means that snowfall has no
1169	! heat content but it can be used for increasing the thickness and changing the density (maybe it is a bug)
1170	psnowhmass(ji) = 0.0
1171	IF ( (precip_snow(ji) .GT. 0.0) ) THEN
1172
1173	!calculate
1174
1175	psnowhmass(ji) = precip_snow(ji)(un-totfrac_nobio(ji)) &
1176	(xci*(snowtemp(ji,1)-tp_00)-chalfu0)
1177
1178	! Snowfall density: Following CROCUS (Pahaut 1976)
1179	!
1180	rhosnew(ji) = MAX(xrhosmin, snowfall_a_sn + snowfall_b_sn*(temp_air(ji)-tp_00)+ &
1181	snowfall_c_sn*SQRT(speed))
1182
1183	! Augment total pack depth:
1184	! Note that dsnowfall(ji) can be equal to zero if totfrac_nobio(ji)=1
1185	dsnowfall(ji) = (precip_snow(ji)*(un-totfrac_nobio(ji)))/rhosnew(ji) !snowfall thickness (m)
1186
1187	snow_depth(ji) = snow_depth(ji) + dsnowfall(ji)
1188
1189
1190	! Fresh snowfall changes the snowpack density and liquid content in uppermost layer
1191
1192	IF (dsnowfall(ji) .GT. zero) THEN
1193	snowrho(ji,1) = (snowdz(ji,1)snowrho(ji,1) + dsnowfall(ji)rhosnew(ji))/ &
1194	(snowdz(ji,1)+dsnowfall(ji))
1195
1196	snowdz(ji,1) = snowdz(ji,1) + dsnowfall(ji)
1197
1198	! Add energy of snowfall to snowpack:
1199	! Update heat content (J/m2) (therefore the snow temperature
1200	! and liquid content):
1201	!
1202	snowheat(ji,1) = snowheat(ji,1) + psnowhmass(ji)
1203	!
1204	! Incorporate snowfall grain size:
1205	!
1206	newgrain(ji) = MIN(dgrain_new_max, snow3lgrain_0d(rhosnew(ji)))
1207
1208	snowgrain(ji,1) = (snowdz_old(ji,1)snowgrain(ji,1) + dsnowfall(ji)newgrain(ji))/ &
1209	snowdz(ji,1)
1210	ENDIF
1211	ELSE
1212	dsnowfall(ji) = 0.
1213	ENDIF
1214
1215	! new snow fall on snow free surface.
1216	! we use the linearization for the new snow fall on snow-free ground
1217
1218	IF ( (dsnowfall(ji) .GT. zero) .AND. (snow_depth_old(ji) .EQ. zero) ) THEN
1219
1220	snowfall_delta = 1.0
1221
1222	DO jj=1,nsnow
1223
1224	snowdz(ji,jj) = snowfall_delta*(dsnowfall(ji)/nsnow) + &
1225	(1.0-snowfall_delta)*snowdz(ji,jj)
1226
1227	snowheat(ji,jj) = snowfall_delta*(psnowhmass(ji)/nsnow) + &
1228	(1.0-snowfall_delta)*snowheat(ji,jj)
1229
1230	snowrho(ji,jj) = snowfall_delta*rhosnew(ji) + &
1231	(1.0-snowfall_delta)*snowrho(ji,jj)
1232
1233	snowgrain(ji,jj)= snowfall_delta*newgrain(ji) + &
1234	(1.0-snowfall_delta)*snowgrain(ji,jj)
1235
1236	ENDDO
1237
1238
1239	ENDIF
1240
1241
1242	ENDDO
1243
1244	END SUBROUTINE explicitsnow_fall
1245
1246	!!
1247	!================================================================================================================================
1248	!! SUBROUTINE : explicitsnow_gone
1249	!!
1250	!>\BRIEF Check whether snow is gone
1251	!!
1252	!! DESCRIPTION : If so, set thickness (and therefore mass and heat) and liquid
1253	!! content to zero, and adjust fluxes of water, evaporation and
1254	!! heat into underlying surface.
1255	!! RECENT CHANGE(S) : None
1256	!!
1257	!! MAIN OUTPUT VARIABLE(S): None
1258	!!
1259	!! REFERENCE(S) :
1260	!!
1261	!! FLOWCHART : None
1262	!! \n
1263	!_
1264	!================================================================================================================================
1265
1266	SUBROUTINE explicitsnow_gone(kjpindex,pgflux,&
1267	snowheat,snowtemp,snowdz,snowrho,snowliq,grndflux,snowmelt)
1268
1269	!! 0.1 Input variables
1270
1271	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1272	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pgflux !! Net energy into snow pack(w/m2)
1273	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowheat !! Snow heat content (J/m^2)
1274
1275	!! 0.2 Output variables
1276
1277	!! 0.3 Modified variables
1278
1279	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature
1280	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth [m]
1281	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
1282	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowliq !! Liquid water content
1283	REAL(r_std),DIMENSION(kjpindex), INTENT(inout) :: grndflux !! Soil/snow interface heat flux (W/m2)
1284	REAL(r_std),DIMENSION(kjpindex),INTENT(inout) :: snowmelt !! Snow melt
1285	REAL(r_std),DIMENSION(kjpindex) :: thrufal !! Water leaving snowpack(kg/m2/s)
1286
1287	!! 0.4 Local variables
1288
1289	INTEGER(i_std) :: ji,jj
1290	REAL(r_std),DIMENSION(kjpindex) :: snowgone_delta
1291	REAL(r_std),DIMENSION (kjpindex) :: totsnowheat !!snow heat content at each layer
1292	REAL(r_std),DIMENSION(kjpindex) :: snowdepth_crit
1293
1294	! first caculate total snowpack snow heat content
1295	snowgone_delta(:) = un
1296	thrufal(:)=0.0
1297	snowmelt(:)=0
1298	totsnowheat(:) = SUM(snowheat(:,:),2)
1299
1300	DO ji = 1, kjpindex
1301
1302	IF ( (SUM(snowdz(ji,:)) .GT. min_sechiba)) THEN
1303
1304	IF( pgflux(ji) >= (-totsnowheat(ji)/dt_sechiba) ) THEN
1305
1306	grndflux(ji) = pgflux(ji) + (totsnowheat(ji)/dt_sechiba)
1307
1308	thrufal(ji)=SUM(snowrho(ji,:)*snowdz(ji,:))
1309
1310	snowgone_delta(ji) = 0.0
1311
1312	snowmelt(ji) = snowmelt(ji)+thrufal(ji)
1313
1314	ENDIF
1315
1316	! update of snow state (either still present or not)
1317
1318	DO jj=1,nsnow
1319	snowdz(ji,jj) = snowdz(ji,jj) *snowgone_delta(ji)
1320	snowliq(ji,jj) = snowliq(ji,jj) *snowgone_delta(ji)
1321	snowtemp(ji,jj) = (1.0-snowgone_delta(ji))tp_00 + snowtemp(ji,jj)snowgone_delta(ji)
1322	ENDDO
1323
1324	ELSE
1325	snowdz(ji,:)=0.
1326	snowliq(ji,:)=0.
1327	snowmelt(ji)=snowmelt(ji)+SUM(snowrho(ji,:)*snowdz(ji,:))
1328	ENDIF
1329	ENDDO
1330
1331	END SUBROUTINE explicitsnow_gone
1332
1333	!================================================================================================================================
1334	!! SUBROUTINE : explicitsnow_melt_refrz
1335	!!
1336	!>\BRIEF Computes snow melt and refreezing processes within snowpack
1337	!!
1338	!! DESCRIPTION :
1339	!! RECENT CHANGE(S) : None
1340	!!
1341	!! MAIN OUTPUT VARIABLE(S): None
1342	!!
1343	!! REFERENCE(S) :
1344	!!
1345	!! FLOWCHART : None
1346	!! \n
1347	!_
1348	!================================================================================================================================
1349
1350	SUBROUTINE explicitsnow_melt_refrz(kjpindex,precip_rain,pgflux,soilcap,&
1351	snowtemp,snowdz,snowrho,snowliq,snowmelt,grndflux,temp_air)
1352
1353	!! 0.1 Input variables
1354
1355	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
1356	REAL(r_std),DIMENSION (kjpindex,nsnow) :: pcapa_snow !! Heat capacity for snow
1357	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall
1358	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: temp_air !! Air temperature
1359	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: pgflux !! Net energy into snowpack(w/m2)
1360	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilcap !! Soil heat capacity
1361
1362	!! 0.2 Output variables
1363
1364	!! 0.3 Modified variables
1365
1366	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature
1367	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth [m]
1368	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow layer density (Kg/m^3)
1369	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowliq !! Liquid water content
1370	REAL(r_std),DIMENSION(kjpindex),INTENT(inout) :: grndflux !! Net energy input to soil
1371	REAL(r_std),DIMENSION (kjpindex), INTENT(inout) :: snowmelt !! Snowmelt
1372
1373	!! 0.4 Local variables
1374
1375	REAL(r_std),DIMENSION (kjpindex) :: meltxs !! Residual snowmelt energy applied to underlying soil
1376	REAL(r_std) :: enerin,melttot,hrain
1377	REAL(r_std),DIMENSION (nsnow) :: zsnowlwe
1378	REAL(r_std),DIMENSION (nsnow) :: flowliq
1379	REAL(r_std),DIMENSION (kjpindex) :: snowmass
1380	REAL(r_std),DIMENSION (nsnow) :: zphase !! Phase change (from ice to water) (J/m2)
1381	REAL(r_std),DIMENSION (nsnow) :: zphase2 !! Phase change (from water to ice)
1382	REAL(r_std),DIMENSION (nsnow) :: zphase3 !! Phase change related with net energy input to snowpack
1383	REAL(r_std),DIMENSION (nsnow) :: zsnowdz !! Snow layer depth [m]
1384	REAL(r_std),DIMENSION (nsnow) :: zsnowmelt !! Snow melt (liquid water) (m)
1385	REAL(r_std),DIMENSION (nsnow) :: zsnowtemp
1386	REAL(r_std),DIMENSION (nsnow) :: zmeltxs !! Excess melt
1387	REAL(r_std),DIMENSION (nsnow) :: zwholdmax !! Maximum liquid water holding (m)
1388	REAL(r_std),DIMENSION (nsnow) :: zcmprsfact !! Compression factor due to densification from melting
1389	REAL(r_std),DIMENSION (nsnow) :: zscap !! Snow heat capacity (J/m3 K)
1390	REAL(r_std),DIMENSION (nsnow) :: zsnowliq !! (m)
1391	REAL(r_std),DIMENSION (nsnow) :: snowtemp_old
1392	REAL(r_std),DIMENSION (0:nsnow) :: zflowliqt !!(m)
1393	REAL(r_std) :: zrainfall,zpcpxs
1394	REAL(r_std) :: ztotwcap
1395	REAL(r_std),DIMENSION(kjpindex,nsnow) :: snowdz_old,snowliq_old
1396	INTEGER(i_std) :: jj,ji, iv
1397	REAL(r_std),DIMENSION(nsnow) :: snowdz_old2
1398	REAL(r_std),DIMENSION(nsnow) :: zsnowrho
1399
1400	!initialize
1401	snowdz_old = snowdz
1402	snowliq_old = snowliq
1403
1404	DO ji = 1, kjpindex
1405
1406
1407	snowmass(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
1408	IF ((snowmass(ji) .GT. min_sechiba)) THEN
1409
1410	!! 1 snow melting due to positive snowpack snow temperature
1411
1412	!! 1.0 total liquid equivalent water content of each snow layer
1413
1414	zsnowlwe(:) = snowrho(ji,:) * snowdz(ji,:)/ ph2o
1415
1416	!! 1.1 phase change (J/m2)
1417
1418	pcapa_snow(ji,:) = snowrho(ji,:)*xci
1419
1420	zphase(:) = MIN(pcapa_snow(ji,:)MAX(0.0, snowtemp(ji,:)-tp_00) &
1421	snowdz(ji,:), &
1422	MAX(0.0,zsnowlwe(:)-snowliq(ji,:))chalfu0ph2o)
1423
1424	!! 1.2 update snow liq water and temperature if melting
1425
1426	zsnowmelt(:) = zphase(:)/(chalfu0*ph2o)
1427
1428	!! 1.3 cool off the snow layer temperature due to melt
1429
1430	zsnowtemp(:) = snowtemp(ji,:) - zphase(:)/(pcapa_snow(ji,:)* snowdz(ji,:))
1431
1432	snowtemp(ji,:) = MIN(tp_00, zsnowtemp(:))
1433
1434	zmeltxs(:) = (zsnowtemp(:)-snowtemp(ji,:))pcapa_snow(ji,:)snowdz(ji,:)
1435
1436	!! 1.4 loss of snowpack depth and liquid equivalent water
1437
1438	zwholdmax(:) = snow3lhold_1d(snowrho(ji,:),snowdz(ji,:)) ! 1 dimension
1439
1440	zcmprsfact(:) = (zsnowlwe(:)-MIN(snowliq(ji,:)+zsnowmelt(:),zwholdmax(:)))/ &
1441	(zsnowlwe(:)-MIN(snowliq(ji,:),zwholdmax(:)))
1442
1443	snowdz(ji,:) = snowdz(ji,:)*zcmprsfact(:)
1444
1445	snowrho(ji,:) = zsnowlwe(:)*ph2o/snowdz(ji,:)
1446
1447	snowliq(ji,:) = snowliq(ji,:) + zsnowmelt(:)
1448
1449
1450	!! 2 snow refreezing process
1451	!! 2.1 freeze liquid water in any layer
1452	zscap(:) = snowrho(ji,:)*xci !J K-1 m-3
1453	zphase2(:) = MIN(zscap(:)* &
1454	MAX(0.0, tp_00 - snowtemp(ji,:))*snowdz(ji,:), &
1455	snowliq(ji,:)chalfu0ph2o)
1456
1457	! warm layer and reduce liquid if freezing occurs
1458	zsnowdz(:) = MAX(xsnowdmin/nsnow, snowdz(ji,:))
1459	snowtemp_old(:) = snowtemp(ji,:)
1460	snowtemp(ji,:) = snowtemp(ji,:) + zphase2(:)/(zscap(:)*zsnowdz(:))
1461
1462	! Reduce liquid portion if freezing occurs:
1463	snowliq(ji,:) = snowliq(ji,:) - ( (snowtemp(ji,:)-snowtemp_old(:))* &
1464	zscap(:)zsnowdz(:)/(chalfu0ph2o) )
1465	snowliq(ji,:) = MAX(snowliq(ji,:), 0.0)
1466	!! 3. thickness change due to snowmelt in excess of holding capacity
1467	zwholdmax(:) = snow3lhold_1d(snowrho(ji,:),snowdz(ji,:)) ! 1 dimension
1468	flowliq(:) = MAX(0.,(snowliq(ji,:)-zwholdmax(:)))
1469	snowliq(ji,:) = snowliq(ji,:) - flowliq(:)
1470	snowdz(ji,:) = snowdz(ji,:) - flowliq(:)*ph2o/snowrho(ji,:)
1471	! to prevent possible very small negative values (machine
1472	! prescision as snow vanishes
1473	snowdz(ji,:) = MAX(0.0, snowdz(ji,:))
1474
1475	!! 4. liquid water flow
1476	ztotwcap = SUM(zwholdmax(:))
1477	! Rain entering snow (m):
1478	!ORCHIDEE has assumed that all precipitation entering snow has
1479	!left the snowpack and finally become runoff in hydrol_soil. Here we just put zrainfall as 0.0
1480	!!zrainfall = precip_rain(ji)/ph2o ! rainfall (m)
1481	zrainfall = 0.0
1482
1483	zflowliqt(0) = MIN(zrainfall,ztotwcap)
1484	! Rain assumed to directly pass through the pack to runoff (m):
1485	zpcpxs = zrainfall - zflowliqt(0)
1486	!
1487	DO jj=1,nsnow
1488	zflowliqt(jj) = flowliq(jj)
1489	ENDDO
1490
1491	! translated into a density increase:
1492	flowliq(:) = 0.0 ! clear this array for work
1493	zsnowliq(:) = snowliq(ji,:) ! reset liquid water content
1494	!
1495
1496	DO jj=1,nsnow
1497	snowliq(ji,jj) = snowliq(ji,jj) + zflowliqt(jj-1)
1498	flowliq(jj) = MAX(0.0, (snowliq(ji,jj)-zwholdmax(jj)))
1499	snowliq(ji,jj) = snowliq(ji,jj) - flowliq(jj)
1500
1501	!Modified by Tao Wang based on previous ISBA-ES scheme
1502	snowrho(ji,jj) = snowrho(ji,jj) + (snowliq(ji,jj) - zsnowliq(jj))* &
1503	& ph2o/MAX(xsnowdmin/nsnow,snowdz(ji,jj))
1504
1505	zflowliqt(jj) = zflowliqt(jj) + flowliq(jj)
1506	ENDDO
1507
1508	snowmelt(ji) = snowmelt(ji) + (zflowliqt(nsnow) + zpcpxs) * ph2o
1509
1510	! excess heat from melting, using it to warm underlying ground to conserve energy
1511	meltxs(ji) = SUM(zmeltxs(:))/dt_sechiba ! (W/m2)
1512
1513	! energy flux into the soil
1514	grndflux(ji) = grndflux(ji) + meltxs(ji)
1515
1516	ELSE
1517	snowdz(ji,:)=0.
1518	snowliq(ji,:)=0.
1519	snowmelt(ji)=snowmelt(ji)+SUM(snowrho(ji,:)*snowdz(ji,:))
1520	ENDIF
1521
1522	ENDDO
1523
1524	END SUBROUTINE explicitsnow_melt_refrz
1525
1526	!================================================================================================================================
1527	!! SUBROUTINE : explicitsnow_levels
1528	!!
1529	!>\BRIEF Computes snow discretization based on given total snow depth
1530	!!
1531	!! DESCRIPTION :
1532	!! RECENT CHANGE(S) : None
1533	!!
1534	!! MAIN OUTPUT VARIABLE(S): None
1535	!!
1536	!! REFERENCE(S) :
1537	!!
1538	!! FLOWCHART : None
1539	!! \n
1540	!_
1541	!================================================================================================================================
1542	SUBROUTINE explicitsnow_levels( kjpindex,snow_thick, snowdz)
1543
1544	!! 0.1 Input variables
1545	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1546	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow_thick !! Total snow depth
1547
1548	!! 0.2 Output variables
1549
1550	!! 0.3 Modified variables
1551
1552	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth [m]
1553
1554	!! 0.4 Local variables
1555
1556	INTEGER(i_std) :: il,it,ji
1557	INTEGER(i_std) :: i,j, ix
1558	REAL(r_std), DIMENSION(kjpindex) :: xi, xf
1559	REAL(r_std), PARAMETER, DIMENSION(3) :: ZSGCOEF1 = (/0.25, 0.50, 0.25/) !! Snow grid parameters
1560	REAL(r_std), PARAMETER, DIMENSION(2) :: ZSGCOEF2 = (/0.05, 0.34/) !! Snow grid parameters
1561	REAL(r_std), PARAMETER :: ZSNOWTRANS = 0.20 !! Minimum total snow depth at which surface layer thickness is constant (m)
1562	REAL(r_std), PARAMETER :: XSNOWCRITD = 0.03 !! (m)
1563
1564	DO ji=1,kjpindex
1565	IF ( snow_thick(ji) .LE. (XSNOWCRITD+0.01)) THEN
1566
1567	snowdz(ji,1) = MIN(0.01, snow_thick(ji)/nsnow)
1568	snowdz(ji,3) = MIN(0.01, snow_thick(ji)/nsnow)
1569	snowdz(ji,2) = snow_thick(ji) - snowdz(ji,1) - snowdz(ji,3)
1570
1571	ENDIF
1572	ENDDO
1573
1574	WHERE ( snow_thick(:) .LE. ZSNOWTRANS .AND. &
1575	snow_thick(:) .GT. (XSNOWCRITD+0.01) )
1576	!
1577	snowdz(:,1) = snow_thick(:)*ZSGCOEF1(1)
1578	snowdz(:,2) = snow_thick(:)*ZSGCOEF1(2)
1579	snowdz(:,3) = snow_thick(:)*ZSGCOEF1(3)
1580	!
1581	END WHERE
1582
1583	DO ji = 1,kjpindex
1584	IF (snow_thick(ji) .GT. ZSNOWTRANS) THEN
1585	snowdz(ji,1) = ZSGCOEF2(1)
1586	snowdz(ji,2) = (snow_thick(ji)-ZSGCOEF2(1))*ZSGCOEF2(2) + ZSGCOEF2(1)
1587	!When using simple finite differences, limit the thickness
1588	!factor between the top and 2nd layers to at most 10
1589	snowdz(ji,2) = MIN(10*ZSGCOEF2(1), snowdz(ji,2) )
1590	snowdz(ji,3) = snow_thick(ji) - snowdz(ji,2) - snowdz(ji,1)
1591
1592	ENDIF
1593	ENDDO
1594
1595
1596	END SUBROUTINE explicitsnow_levels
1597
1598	!!
1599	!================================================================================================================================
1600	!! SUBROUTINE : explicitsnow_profile
1601	!!
1602	!>\BRIEF
1603	!!
1604	!! DESCRIPTION : In this routine solves the numerical soil thermal scheme, ie calculates the new soil temperature profile.
1605	!!
1606	!! RECENT CHANGE(S) : None
1607	!!
1608	!! MAIN OUTPUT VARIABLE(S):
1609	!!
1610	!! REFERENCE(S) :
1611	!!
1612	!! FLOWCHART : None
1613	!! \n
1614	!_
1615	!================================================================================================================================
1616	SUBROUTINE explicitsnow_profile (kjpindex, cgrnd_snow,dgrnd_snow,lambda_snow,temp_sol_new, snowtemp,snowdz,temp_sol_add)
1617
1618	!! 0. Variables and parameter declaration
1619
1620	!! 0.1 Input variables
1621
1622	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1623	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! skin temperature
1624	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
1625	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
1626	REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
1627	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow layer thickness [m]
1628	REAL(r_std), DIMENSION (kjpindex),INTENT(inout) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
1629
1630	!! 0.3 Modified variable
1631
1632	!! 0.2 Output variables
1633
1634	!! 0.3 Modified variables
1635	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout) :: snowtemp
1636
1637	!! 0.4 Local variables
1638
1639	INTEGER(i_std) :: ji, jg
1640	!_
1641	!================================================================================================================================
1642	!! 1. Computes the snow temperatures ptn.
1643	DO ji = 1,kjpindex
1644	IF (SUM(snowdz(ji,:)) .GT. 0) THEN
1645	snowtemp(ji,1) = (lambda_snow(ji) * cgrnd_snow(ji,1) + (temp_sol_new(ji)+temp_sol_add(ji))) / &
1646	(lambda_snow(ji) * (un - dgrnd_snow(ji,1)) + un)
1647	temp_sol_add(ji) = zero
1648	DO jg = 1,nsnow-1
1649	snowtemp(ji,jg+1) = cgrnd_snow(ji,jg) + dgrnd_snow(ji,jg) * snowtemp(ji,jg)
1650	ENDDO
1651	ENDIF
1652	ENDDO
1653	IF (printlev>=3) WRITE (numout,*) ' explicitsnow_profile done '
1654
1655	END SUBROUTINE explicitsnow_profile
1656
1657	END MODULE explicitsnow

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

Context Navigation

source: branches/publications/ORCHIDEE_CAMEO_gmd_2022/src_sechiba/explicitsnow.f90 @ 7693

Download in other formats: