Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

explicitsnow.f90 @ 5143

Last change on this file since 5143 was 3965, checked in by jan.polcher, 8 years ago

Merge with trunk at revision3959.
This includes all the developments made for CMIP6 and passage to XIOS2.
All conflicts are resolved and the code compiles.

But it still does not link because of an "undefined reference to `_intel_fast_memmove'"

File size: 70.7 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_sechiba/explicitsnow.f90 @ 5143

Download in other formats: