Context Navigation

source: branches/publications/ORCHILEAK-Gommet/src_sechiba/explicitsnow.f90 @ 7746

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

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

Download in other formats: