Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

explicitsnow.f90 @ 7541

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

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

Context Navigation

source: branches/publications/ORCHIDEE_2.2_r7266/ORCHIDEE/src_sechiba/explicitsnow.f90 @ 7541

Download in other formats: