Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

explicitsnow.f90 @ 5501

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

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

Context Navigation

source: CONFIG_DEVT/IPSLCM6.5_work_ENSEMBLES/modeles/ORCHIDEE/src_sechiba/explicitsnow.f90 @ 5501

Download in other formats: