Context Navigation

source: branches/publications/ORCHIDEE-PEAT_r5488/src_sechiba/explicitsnow.f90 @ 8787

Last change on this file since 8787 was 3924, checked in by albert.jornet, 8 years ago

Merge: from trunk revisions [3643:3883/trunk/ORCHIDEE]

Done in [3887:3922/perso/albert.jornet/MICT-MERGING]

File size: 72.2 KB

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

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

Download in other formats: