Context Navigation

source: branches/publications/ORCHIDEE_GLUC_r6545/src_sechiba/explicitsnow.f90 @ 7346

Last change on this file since 7346 was 4719, checked in by albert.jornet, 7 years ago

Merge: from revisions [4491:4695/trunk/ORCHIDEE]

Merge done in [4671:4718/perso/albert.jornet/MICT_MERGE]

Property svn:keywords set to Date Revision HeadURL

File size: 71.6 KB

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

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

Download in other formats: