Context Navigation

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

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

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

Download in other formats: