Context Navigation

source: branches/publications/ORCHIDEE-ICE_SurfaceMassBalance/src_sechiba/explicitsnow.f90 @ 8398

Last change on this file since 8398 was 7642, checked in by christophe.dumas, 2 years ago
subroutine explicitsnow_maxmass : simplified code with test on nsnow (3 or 12) defining nsnowstart and nsnowstop (layers where snowmelt_from_maxmass can be applied) at the beginning.
Property svn:keywords set to `Date Revision HeadURL`
File size: 111.4 KB

Line
1	! ================================================================================================================================
2	! MODULE : explicitsnow
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Computes hydrologic Xsnow processes on continental points.
10	!!
11	!!\n DESCRIPTION: This module computes hydrologic snow processes on continental points.
12	!!
13	!! RECENT CHANGES:
14	!!
15	!! REFERENCES : None
16	!!
17	!! SVN :
18	!! $HeadURL$
19	!! $Date$
20	!! $Revision$
21	!! \n
22	!_ ================================================================================================================================
23
24	MODULE explicitsnow
25	USE ioipsl_para
26	USE constantes_soil
27	USE constantes
28	USE pft_parameters
29	USE pft_parameters_var !! -SC- Ajout pour ajout throughfall dans zrainfall
30	USE qsat_moisture
31	USE sechiba_io_p
32	USE xios_orchidee
33	USE grid !! cdc Ajout pour explicitsnow_read_reftempicefile
34
35	IMPLICIT NONE
36
37	! Public routines :
38	PRIVATE
39	PUBLIC :: explicitsnow_main, explicitsnow_initialize, explicitsnow_finalize
40
41	CONTAINS
42
43	!================================================================================================================================
44	!! SUBROUTINE : explicitsnow_initialize
45	!!
46	!>\BRIEF Read variables for explictsnow module from restart file
47	!!
48	!! DESCRIPTION : Read variables for explictsnow module from restart file
49	!!
50	!! \n
51	!_
52	!================================================================================================================================
53	SUBROUTINE explicitsnow_initialize( kjit, kjpindex, rest_id, snowrho, &
54	snowtemp, snowdz, snowheat, snowgrain, &
55	icetemp, icedz)
56
57	!! 0.1 Input variables
58	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
59	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
60	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
61
62	!! 0.2 Output variables
63	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowrho !! Snow density (Kg/m^3)
64	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowtemp !! Snow temperature
65	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowdz !! Snow layer thickness
66	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowheat !! Snow heat content (J/m^2)
67	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowgrain !! Snow grainsize (m)
68	REAL(r_std), DIMENSION (kjpindex,nice), INTENT(out) :: icetemp !! Ice temperature
69	REAL(r_std), DIMENSION (kjpindex,nice), INTENT(out) :: icedz !! Ice layer thickness [m]
70
71	!! local variables
72	integer(i_std) :: ji
73
74	!! 1. Read from restart file
75	CALL restget_p (rest_id, 'snowrho', nbp_glo, nsnow, 1, kjit,.TRUE.,snowrho, "gather", nbp_glo, index_g)
76	CALL setvar_p (snowrho, val_exp, 'Snow Density profile', xrhosmin)
77
78	CALL restget_p (rest_id, 'snowtemp', nbp_glo, nsnow, 1, kjit,.TRUE.,snowtemp, "gather", nbp_glo, index_g)
79	CALL setvar_p (snowtemp, val_exp, 'Snow Temperature profile', tp_00)
80
81	CALL restget_p (rest_id, 'snowdz', nbp_glo, nsnow, 1, kjit,.TRUE.,snowdz, "gather", nbp_glo, index_g)
82	CALL setvar_p (snowdz, val_exp, 'Snow depth profile', 0.0)
83
84	CALL restget_p (rest_id, 'snowheat', nbp_glo, nsnow, 1, kjit,.TRUE.,snowheat, "gather", nbp_glo, index_g)
85	CALL setvar_p (snowheat, val_exp, 'Snow Heat profile', 0.0)
86
87	CALL restget_p (rest_id, 'snowgrain', nbp_glo, nsnow, 1, kjit,.TRUE.,snowgrain, "gather", nbp_glo, index_g)
88	CALL setvar_p (snowgrain, val_exp, 'Snow grain profile', 0.0)
89
90	IF (ok_ice_sheet) THEN
91	CALL restget_p (rest_id, 'icetemp', nbp_glo, nice, 1, kjit,.TRUE.,icetemp, "gather", nbp_glo, index_g)
92
93	IF (ALL(icetemp(:,:)==val_exp)) THEN
94	! icetemp was not found in restart file
95	IF (read_reftempice) THEN
96	! Read variable ptn from file
97	CALL explicitsnow_read_reftempicefile(kjpindex,lalo,icetemp)
98	ELSE
99	! Initialize icetemp with a constant value which can be set in run.def
100
101	!Config Key = EXPLICITSNOW_TICEPRO
102	!Config Desc = Initial ice temperature profile if not found in restart
103	!Config Def = 273.
104	!Config If = OK_SECHIBA
105	!Config Help = The initial value of the temperature profile in the soil if
106	!Config its value is not found in the restart file. Here
107	!Config we only require one value as we will assume a constant
108	!Config throughout the column.
109	!Config Units = Kelvin [K]
110	CALL setvar_p (icetemp, val_exp,'EXPLICITSNOW_TICEPRO',273._r_std)
111	END IF
112	END IF
113	do ji=1,kjpindex
114	icedz(ji,:) = ZSICOEF1(:)
115	enddo
116	ENDIF
117
118
119
120	END SUBROUTINE explicitsnow_initialize
121
122
123	!================================================================================================================================
124	!! SUBROUTINE : explicitsnow_main
125	!!
126	!>\BRIEF Main call for snow calculations
127	!!
128	!! DESCRIPTION : Main routine for calculation of the snow processes with explictsnow module.
129	!!
130	!! RECENT CHANGE(S) : None
131	!!
132	!! MAIN OUTPUT VARIABLE(S): None
133	!!
134	!! REFERENCE(S) :
135	!!
136	!! FLOWCHART : None
137	!! \n
138	!_
139	!================================================================================================================================
140	SUBROUTINE explicitsnow_main(kjpindex, precip_rain, precip_snow, temp_air, pb, & ! in
141	u, v, temp_sol_new, soilcap, pgflux, & ! in
142	frac_nobio, totfrac_nobio,frac_snow_nobio,gtemp, & ! in
143	lambda_snow, cgrnd_snow, dgrnd_snow, & ! in
144	lambda_ice, cgrnd_ice, dgrnd_ice, njsc, & ! in
145	vevapsno, snow_age, snow_nobio_age,snow_nobio, snowrho, & ! inout
146	snowgrain, snowdz, snowtemp, snowheat, snow, & ! inout
147	temp_sol_add, icetemp, icedz, & ! inout
148	snowliq, subsnownobio, grndflux, snowmelt, tot_melt, & ! output
149	subsinksoil, zrainfall, frac_snow_veg, veget, veget_max) ! output
150
151
152	!! 0. Variable and parameter declaration
153
154	!! 0.1 Input variables
155	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
156	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall
157	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_snow !! Snowfall
158	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_air !! Air temperature
159	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: pb !! Surface pressure
160	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: u,v !! Horizontal wind speed
161	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_sol_new !! Surface temperature
162	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: soilcap !! Soil capacity
163	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: pgflux !! Net energy into snowpack
164	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
165	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+ ...
166	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
167	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegetation
168	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: gtemp !! First soil layer temperature
169	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
170	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
171	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
172	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: lambda_ice !! Coefficient of the linear extrapolation of ice surface temperature
173	REAL(r_std), DIMENSION (kjpindex,nice), INTENT (in) :: cgrnd_ice !! Integration coefficient for ice numerical scheme
174	REAL(r_std), DIMENSION (kjpindex,nice), INTENT (in) :: dgrnd_ice !! Integration coefficient for ice numerical scheme
175	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type
176	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type
177	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
178
179	!! 0.2 Output fields
180	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowliq !! Snow liquid content (m)
181	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(out) :: subsnownobio !! Sublimation of snow on other surface types (ice, lakes, ...)
182	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: grndflux !! Net flux into soil [W/m2]
183	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: snowmelt !! Snow melt [mm/dt_sechiba]
184	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_melt !! Total melt from ice and snow
185	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: subsinksoil !! Excess of sublimation as a sink for the soil
186	!cdc ajout zrainfall pour calcul smb
187	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: zrainfall !! Rain precipitation on snow
188	!! @tex $(kg m^{-2})$ @endtex
189
190	!! 0.3 Modified fields
191	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapsno !! Snow evaporation @tex ($kg m^{-2}$) @endtex
192	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow_age !! Snow age
193	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio !! Ice water balance
194	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio_age !! Snow age on ice, lakes, ...
195	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
196	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowgrain !! Snow grainsize (m)
197	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow layer thickness [m]
198	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature
199	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowheat !! Snow heat content (J/m^2)
200	REAL(r_std), DIMENSION (kjpindex,nice), INTENT(inout) :: icetemp !! Ice temperature
201	REAL(r_std), DIMENSION (kjpindex,nice), INTENT(inout) :: icedz !! Ice layer thickness [m]
202	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow !! Snow mass [Kg/m^2]
203	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
204
205
206	!! 0.4 Local declaration
207	INTEGER(i_std) :: ji, iv, jj,m,jv
208	REAL(r_std),DIMENSION (kjpindex) :: snow_depth_tmp
209	REAL(r_std),DIMENSION (kjpindex) :: snowmelt_from_maxmass
210	REAL(r_std) :: snowdzm1
211	REAL(r_std), DIMENSION (kjpindex) :: thrufal !! Water leaving snow pack [kg/m2/s]
212	REAL(r_std), DIMENSION (kjpindex) :: snowmelt_tmp,snowmelt_ice,temp_sol_new_old
213	REAL(r_std), DIMENSION (kjpindex,nsnow) :: snowdz_old
214	REAL(r_std), DIMENSION (kjpindex) :: ZLIQHEATXS
215	REAL(r_std) :: grndflux_tmp
216	REAL(r_std), DIMENSION (nsnow) :: snowtemp_tmp
217	REAL(r_std) :: s2flux_tmp,fromsoilflux
218	REAL(r_std), DIMENSION (kjpindex,nsnow) :: pcapa_snow
219	REAL(r_std), DIMENSION (kjpindex) :: psnowhmass
220	REAL(r_std), PARAMETER :: XP00 = 1.E5
221	REAL(r_std), DIMENSION (kjpindex) :: ice_sheet_melt !! Ice melt [mm/dt_sechiba]
222
223	!! 1. Initialization
224
225	temp_sol_new_old = temp_sol_new
226	snowmelt_ice(:) = zero
227	tot_melt(:) = zero
228	snowmelt(:) = zero
229	snowmelt_from_maxmass(:) = zero
230
231
232	!! 2. on Vegetation
233	! 2.1 Snow fall
234	CALL explicitsnow_fall(kjpindex,precip_snow,temp_air,u,v,totfrac_nobio,snowrho,snowdz,&
235	& snowheat,snowgrain,snowtemp,psnowhmass)
236
237	! 2.2 calculate the new snow discretization
238	snow_depth_tmp(:) = SUM(snowdz(:,:),2)
239
240	snowdz_old = snowdz
241
242	! update snowdz
243	CALL explicitsnow_levels(kjpindex,snow_depth_tmp, snowdz)
244
245	! 2.3 Snow heat redistribution
246	CALL explicitsnow_transf(kjpindex,snowdz_old,snowdz,snowrho,snowheat,snowgrain)
247
248	! 2.4 Diagonize water portion of the snow from snow heat content:
249	DO ji=1, kjpindex
250	IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
251	snowtemp(ji,:) = snow3ltemp_1d(snowheat(ji,:),snowrho(ji,:),snowdz(ji,:))
252	snowliq(ji,:) = snow3lliq_1d(snowheat(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:))
253	ELSE
254	snowliq(ji,:) = zero
255	snowtemp(ji,:) = tp_00
256	ENDIF
257	END DO
258
259	! 2.5 snow compaction
260
261	!cdc snow compaction : original scheme
262	CALL explicitsnow_compactn(kjpindex,snowtemp,snowrho,snowdz)
263
264	!cdc snow compaction resulting from changes in snow viscosity : compaction Decharme 2016
265	! CALL explicitsnow_compactn_up(kjpindex,u,v,snowtemp,snowrho,snowdz,snowliq)
266	!cdc snow compaction due to snowdrift (used with explicitsnow_compactn_up)
267	! CALL explicitsnow_drift(kjpindex,u,v,snowrho,snowdz)
268
269	! Update snow heat
270	DO ji = 1, kjpindex
271	snowheat(ji,:) = snow3lheat_1d(snowliq(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:))
272	ENDDO
273
274	!2.6 Calculate the snow temperature profile based on heat diffusion
275	CALL explicitsnow_profile (kjpindex,cgrnd_snow,dgrnd_snow,lambda_snow,temp_sol_new, snowtemp,snowdz,temp_sol_add)
276
277	!2.7 Test whether snow is existed on the ground or not
278	grndflux(:)=0.0
279	CALL explicitsnow_gone(kjpindex,pgflux,&
280	& snowheat,snowtemp,snowdz,snowrho,snowliq,grndflux,snowmelt)
281
282	!2.8 Calculate snow melt/refreezing processes
283	CALL explicitsnow_melt_refrz(kjpindex,precip_rain,pgflux,soilcap,totfrac_nobio,frac_snow_nobio, &
284	& snowtemp,snowdz,snowrho,snowliq,snowmelt,grndflux, temp_air,frac_snow_veg,veget,veget_max, zrainfall)
285
286	IF (ok_ice_sheet) THEN
287	!cdc Calculate ice temperature
288	CALL explicitsnow_iceprofile(kjpindex, cgrnd_ice,dgrnd_ice,lambda_ice,temp_sol_new,icedz, &
289	& njsc, snowdz, cgrnd_snow, dgrnd_snow, snowtemp, temp_sol_add, icetemp)
290
291	!cdc Calculate ice melt
292	CALL explicitsnow_icemelt(kjpindex,njsc,snowdz,precip_snow, zrainfall, snowmelt, vevapsno, icetemp,icedz, &
293	& grndflux,ice_sheet_melt)
294
295	!cdc Update icetemp and icedz
296	CALL explicitsnow_icelevels(kjpindex,njsc,ice_sheet_melt,icetemp,icedz)
297	ENDIF
298
299
300	! 2.9 Snow sublimation changing snow thickness
301	CALL explicitsnow_subli(kjpindex,snowrho,snowdz,vevapsno, &
302	snowliq,snowtemp,snow,subsnownobio,subsinksoil)
303
304	! 2.9.1 Snowmelt_from_maxmass
305	CALL explicitsnow_maxmass(kjpindex,snowrho,soilcap,snow,snowdz,snowmelt_from_maxmass)
306
307
308	!2.10 Calculate snow grain size using the updated thermal gradient
309	CALL explicitsnow_grain(kjpindex,snowliq,snowdz,gtemp,snowtemp,pb,snowgrain)
310
311	!2.11 Update snow heat
312	! Update the heat content (variable stored each time step)
313	! using current snow temperature and liquid water content:
314	!
315	! First, make check to make sure heat content not too large
316	! (this can result due to signifigant heating of thin snowpacks):
317	! add any excess heat to ground flux:
318	!
319	DO ji=1,kjpindex
320	DO jj=1,nsnow
321	ZLIQHEATXS(ji) = MAX(0.0, snowliq(ji,jj)ph2o - 0.10snowdz(ji,jj)snowrho(ji,jj))chalfu0/dt_sechiba
322	snowliq(ji,jj) = snowliq(ji,jj) - ZLIQHEATXS(ji)dt_sechiba/(ph2ochalfu0)
323	snowliq(ji,jj) = MAX(0.0, snowliq(ji,jj))
324	grndflux(ji) = grndflux(ji) + ZLIQHEATXS(ji)
325	ENDDO
326	ENDDO
327
328	snow(:) = 0.0
329	DO ji=1,kjpindex !domain size
330	snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
331	ENDDO
332
333	DO ji = 1, kjpindex
334	snowheat(ji,:) = snow3lheat_1d(snowliq(ji,:),snowrho(ji,:),snowdz(ji,:),snowtemp(ji,:))
335	ENDDO
336
337
338	!! 4. On other surface types - not done yet
339
340	IF ( nnobio .GT. 1 ) THEN
341	WRITE(numout,*) 'WE HAVE',nnobio-1,' SURFACE TYPES I DO NOT KNOW'
342	WRITE(numout,*) 'CANNOT TREAT SNOW ON THESE SURFACE TYPES'
343	CALL ipslerr_p(3,'explicitsnow_main','Cannot treat snow on these surface types.','','')
344	ENDIF
345
346	!! 5. Computes snow age on land and land ice (for albedo)
347	call explicitsnow_age(kjpindex,snow,precip_snow,precip_rain,frac_snow_nobio, &
348	temp_sol_new,snow_age,snow_nobio_age)
349
350
351	!! 6. Check the snow on land
352	DO ji=1,kjpindex
353	IF (snow(ji) .EQ. 0) THEN
354	snowrho(ji,:)=50.0
355	snowgrain(ji,:)=0.0
356	snowdz(ji,:)=0.0
357	snowliq(ji,:)=0.0
358	ENDIF
359	ENDDO
360
361
362	! Snow melt only if there is more than a given mass : maxmass_snow
363	! Here I suggest to remove the snow based on a certain threshold of snow
364	! depth instead of snow mass because it is quite difficult for
365	! explictsnow module to calculate other snow properties following the
366	! removal of snow mass
367	! to define the threshold of snow depth based on old snow density (330
368	! kg/m3)
369	! maxmass_snowdepth=maxmass_snow/sn_dens
370	! snowmelt_from_maxmass(:) = 0.0
371
372	!! 7. compute total melt
373
374	DO ji=1,kjpindex
375	!cdc tot_melt(ji) = icemelt(ji) + snowmelt(ji) + snowmelt_ice(ji) + snowmelt_from_maxmass(ji)
376	tot_melt(ji) = snowmelt(ji) + snowmelt_from_maxmass(ji)
377	ENDDO
378
379	!cdc modif Fabienne pour TWBR
380	snow_depth_tmp(:) = SUM(snowdz(:,:),2)
381	snowdz_old = snowdz
382	! update snowdz
383	CALL explicitsnow_levels(kjpindex,snow_depth_tmp, snowdz)
384	! snow heat redistribution
385	CALL explicitsnow_transf(kjpindex,snowdz_old,snowdz,snowrho,snowheat,snowgrain)
386
387	IF (printlev>=3) WRITE(numout,*) 'explicitsnow_main done'
388
389	! Write output variables
390	CALL xios_orchidee_send_field("snowliq",snowliq)
391	CALL xios_orchidee_send_field("snowrho",snowrho)
392	CALL xios_orchidee_send_field("snowheat",snowheat)
393	!cdc CALL xios_orchidee_send_field("snowgrain",snowgrain) ! plantage parfois avec sortie snowgrain
394	CALL xios_orchidee_send_field("snowmelt_from_maxmass",snowmelt_from_maxmass/dt_sechiba)
395	CALL xios_orchidee_send_field("soilcap",soilcap)
396
397	END SUBROUTINE explicitsnow_main
398
399
400	!================================================================================================================================
401	!! SUBROUTINE : explicitsnow_finalize
402	!!
403	!>\BRIEF Write variables for explictsnow module to restart file
404	!!
405	!! DESCRIPTION : Write variables for explictsnow module to restart file
406	!!
407	!! \n
408	!_
409	!================================================================================================================================
410	SUBROUTINE explicitsnow_finalize ( kjit, kjpindex, rest_id, snowrho, &
411	snowtemp, snowdz, snowheat, snowgrain, icetemp)
412
413	!! 0.1 Input variables
414	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
415	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
416	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
417	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density (Kg/m^3)
418	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature
419	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow layer thickness
420	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowheat !! Snow heat content (J/m^2)
421	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowgrain !! Snow grainsize (m)
422	REAL(r_std), DIMENSION (kjpindex,nice), INTENT(in) :: icetemp !! Ice temperature (K)
423
424	!! 1. Write to restart file
425	CALL restput_p(rest_id, 'snowrho', nbp_glo, nsnow, 1, kjit, snowrho, 'scatter', nbp_glo, index_g)
426	CALL restput_p(rest_id, 'snowtemp', nbp_glo, nsnow, 1, kjit, snowtemp, 'scatter', nbp_glo, index_g)
427	CALL restput_p(rest_id, 'snowdz', nbp_glo, nsnow, 1, kjit, snowdz, 'scatter', nbp_glo, index_g)
428	CALL restput_p(rest_id, 'snowheat', nbp_glo, nsnow, 1, kjit, snowheat, 'scatter', nbp_glo, index_g)
429	CALL restput_p(rest_id, 'snowgrain', nbp_glo, nsnow, 1, kjit, snowgrain, 'scatter', nbp_glo, index_g)
430	IF ( ok_ice_sheet) CALL restput_p(rest_id, 'icetemp', nbp_glo, nice, 1, kjit, icetemp, 'scatter', nbp_glo, index_g)
431
432	END SUBROUTINE explicitsnow_finalize
433
434
435	!================================================================================================================================
436	!! SUBROUTINE : explicitsnow_grain
437	!!
438	!>\BRIEF Compute evolution of snow grain size
439	!!
440	!! DESCRIPTION :
441	!!
442	!! RECENT CHANGE(S) : None
443	!!
444	!! MAIN OUTPUT VARIABLE(S): None
445	!!
446	!! REFERENCE(S) : R. Jordan (1991)
447	!!
448	!! FLOWCHART : None
449	!! \n
450	!_
451	!================================================================================================================================
452
453
454	SUBROUTINE explicitsnow_grain(kjpindex,snowliq,snowdz,gtemp,snowtemp,pb,snowgrain)
455
456	!! 0.1 Input variables
457	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
458	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowliq !! Liquid water content
459	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth (m)
460	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: gtemp !! First soil layer temperature
461	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature
462	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: pb !! Surface pressure (hpa)
463
464	!! 0.2 Output variables
465
466	!! 0.3 Modified variables
467	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m)
468
469	!! 0.4 Local variables
470	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowdz,zdz,ztheta
471	REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: ztemp,zdiff,ztgrad,zthetaa,zfrac,&
472	zexpo,zckt_liq,zckt_ice,zckt
473	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zrhomin,zgrainmin
474	INTEGER(i_std) :: ji,jj
475
476	!! 0.5 Local parameters
477	REAL(r_std), PARAMETER :: ztheta_crit = 0.02 !! m3 m-3
478	REAL(r_std), PARAMETER :: zc1_ice = 8.047E+9 !! kg m-3 K
479	REAL(r_std), PARAMETER :: zc1_liq = 5.726E+8 !! kg m-3 K
480	REAL(r_std), PARAMETER :: zdeos = 0.92E-4 !! effective diffusion
481	!! coef for water vapor in snow
482	!! at 0C and 1000 mb (m2 s-1)
483	REAL(r_std), PARAMETER :: zg1 = 5.0E-7 !! m4 kg-1
484	REAL(r_std), PARAMETER :: zg2 = 4.0E-12 !! m2 s-1
485	REAL(r_std), PARAMETER :: ztheta_w = 0.05 !! m3 m-3
486	REAL(r_std), PARAMETER :: ztheta_crit_w = 0.14 !! m3 m-3
487	REAL(r_std), PARAMETER :: zdzmin = 0.01 !! m : minimum thickness
488	!! for thermal gradient evaluation:
489	!! to prevent excessive gradients
490	!! for vanishingly thin snowpacks.
491	REAL(r_std), PARAMETER :: xp00=1.E5
492
493	!! 1. initialize
494
495	DO ji=1,kjpindex
496
497
498	zsnowdz(ji,:) = MAX(xsnowdmin/nsnow, snowdz(ji,:))
499
500	DO jj=1,nsnow-1
501	zdz(ji,jj) = zsnowdz(ji,jj) + zsnowdz(ji,jj+1)
502	ENDDO
503	zdz(ji,nsnow) = zsnowdz(ji,nsnow)
504
505	! compute interface average volumetric water content (m3 m-3):
506	! first, layer avg VWC:
507	!
508	ztheta(ji,:) = snowliq(ji,:)/MAX(xsnowdmin, zsnowdz(ji,:))
509
510	! at interfaces:
511	zthetaa(ji,0) = ztheta(ji,1)
512	DO jj=1,nsnow-1
513	zthetaa(ji,jj) = (zsnowdz(ji,jj) *ztheta(ji,jj) + &
514	zsnowdz(ji,jj+1)*ztheta(ji,jj+1))/zdz(ji,jj)
515	ENDDO
516	zthetaa(ji,nsnow) = ztheta(ji,nsnow)
517	! compute interface average temperatures (K):
518	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
519	!
520	ztemp(ji,0) = snowtemp(ji,1)
521	DO jj=1,nsnow-1
522	ztemp(ji,jj) = (zsnowdz(ji,jj) *snowtemp(ji,jj) + &
523	zsnowdz(ji,jj+1)*snowtemp(ji,jj+1))/zdz(ji,jj)
524	ENDDO
525	ztemp(ji,nsnow) = snowtemp(ji,nsnow)
526	!
527	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
528	! compute variation of saturation vapor pressure with temperature
529	! for solid and liquid phases:
530	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
531	zexpo(ji,:) = chalsu0/(xrv*ztemp(ji,:))
532	zckt_ice(ji,:) = (zc1_ice/ztemp(ji,:)*2)(zexpo(ji,:) - 1.0)*EXP(-zexpo(ji,:))
533	!
534	zexpo(ji,:) = chalev0/(xrv*ztemp(ji,:))
535	zckt_liq(ji,:) = (zc1_liq/ztemp(ji,:)*2)(zexpo(ji,:) - 1.0)*EXP(-zexpo(ji,:))
536	!
537	! compute the weighted ice/liquid total variation (N m-2 K):
538	!
539	zfrac(ji,:) = MIN(1.0, zthetaa(ji,:)/ztheta_crit)
540	zckt(ji,:) = zfrac(ji,:)zckt_liq(ji,:) + (1.0 - zfrac(ji,:))zckt_ice(ji,:)
541
542	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
543	! Compute effective diffusion coefficient (m2 s-1):
544	! -diffudivity relative to a reference diffusion at 1000 mb and freezing point
545	! multiplied by phase energy coefficient
546	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
547	!
548	DO jj=0,nsnow
549	zdiff(ji,jj) = zdeos(xp00/(pb(ji)100.))((ztemp(ji,jj)/tp_00)6)zckt(ji,jj)
550	ENDDO
551
552	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
553	! Temperature gradient (K m-1):
554
555	ztgrad(ji,0) = 0.0 ! uppermost layer-mean and surface T's are assumed to be equal
556	DO jj=1,nsnow-1
557	ztgrad(ji,jj) = 2*(snowtemp(ji,jj)-snowtemp(ji,jj+1))/MAX(zdzmin, zdz(ji,jj))
558	ENDDO
559	!
560	! assume at base of snow, temperature is in equilibrium with soil
561	! (but obviously must be at or below freezing point!)
562	!
563	ztgrad(ji,nsnow) = 2*(snowtemp(ji,nsnow) - MIN(tp_00, gtemp(ji)))/MAX(zdzmin, zdz(ji,nsnow))
564	! prognostic grain size (m) equation:
565	!-------------------------------------------------------------------
566	! first compute the minimum grain size (m):
567	!
568	zrhomin(ji,:) = xrhosmin
569	zgrainmin(ji,:) = snow3lgrain_1d(zrhomin(ji,:))
570
571	! dry snow:
572	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
573	!
574	DO jj=1,nsnow
575
576	IF(ztheta(ji,jj) == 0.0) THEN
577
578	! only allow growth due to vapor flux INTO layer:
579	! aab add sublimation(only condensation) as upper BC...?
580
581	snowgrain(ji,jj) = snowgrain(ji,jj) + &
582	(dt_sechibazg1/MAX(zgrainmin(ji,jj),snowgrain(ji,jj))) &
583	( zdiff(ji,jj-1)*MAX(0.0,ztgrad(ji,jj-1)) - &
584	zdiff(ji,jj) *MIN(0.0,ztgrad(ji,jj)) )
585	ELSE
586
587	! wet snow
588	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
589	!
590	snowgrain(ji,jj) = snowgrain(ji,jj) + &
591	(dt_sechibazg2/MAX(zgrainmin(ji,jj),snowgrain(ji,jj))) &
592	MIN(ztheta_crit_w, ztheta(ji,jj) + ztheta_w)
593	END IF
594
595	ENDDO
596
597
598	ENDDO
599
600
601	END SUBROUTINE explicitsnow_grain
602
603	!================================================================================================================================
604	!! SUBROUTINE : explicitsnow_compactn
605	!!
606	!>\BRIEF Compute Compaction/Settling
607	!!
608	!! DESCRIPTION :
609	!! Snow compaction due to overburden and settling.
610	!! Mass is unchanged: layer thickness is reduced
611	!! in proportion to density increases. Method
612	!! of Anderson (1976): see Loth and Graf, 1993,
613	!! J. of Geophys. Res., 98, 10,451-10,464.
614	!!
615	!! RECENT CHANGE(S) : None
616	!!
617	!! MAIN OUTPUT VARIABLE(S): snowrho, snowdz
618	!!
619	!! REFERENCE(S) : Loth and Graf (1993), Mellor (1964) and Anderson (1976)
620	!!
621	!! FLOWCHART : None
622	!! \n
623	!_
624	!================================================================================================================================
625
626
627	SUBROUTINE explicitsnow_compactn(kjpindex,snowtemp,snowrho,snowdz)
628
629	!! 0.1 Input variables
630	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
631	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature
632
633	!! 0.2 Output variables
634
635	!! 0.3 Modified variables
636	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
637	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth
638
639	!! 0.4 Local variables
640	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zwsnowdz,zsmass,zsnowrho2,zviscocity,zsettle
641	REAL(r_std),DIMENSION(kjpindex) :: zsmassc !! cummulative snow mass (kg/m2)
642	REAL(r_std),DIMENSION(kjpindex) :: snowdepth_crit
643	INTEGER(i_std) :: ji,jj
644
645	!! 1. initialize
646
647	zsnowrho2 = snowrho
648	zsettle(:,:) = ZSNOWCMPCT_ACM
649	zviscocity(:,:) = ZSNOWCMPCT_V0
650
651	!! 2. Calculating Cumulative snow mass (kg/m2):
652
653	DO ji=1, kjpindex
654
655
656	IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
657
658	zwsnowdz(ji,:)= snowdz(ji,:)*snowrho(ji,:)
659
660	zsmassc (ji)= 0.0
661
662	DO jj=1,nsnow
663	zsmass(ji,jj) = zsmassc(ji) + zwsnowdz(ji,jj)
664	zsmassc(ji) = zsmassc(ji) + zwsnowdz(ji,jj)
665	ENDDO
666
667
668	!! 3. Computing compaction/Settling
669	! ----------------------
670	! Compaction/settling if density below upper limit
671	! (compaction is generally quite small above ~ 500 kg m-3):
672	!
673	DO jj=1,nsnow
674	IF (snowrho(ji,jj) .LT. xrhosmax) THEN
675
676	!
677	! First calculate settling due to freshly fallen snow: (NOTE:bug here for the snow temperature profile)
678	!
679
680	zsettle(ji,jj) = ZSNOWCMPCT_ACM*EXP( &
681	-ZSNOWCMPCT_BCM*(tp_00-MIN(tp_00,snowtemp(ji,jj))) &
682	-ZSNOWCMPCT_CCM*MAX(0.0, &
683	snowrho(ji,jj)-ZSNOWCMPCT_RHOD))
684	!
685	! Snow viscocity:
686	!
687
688	!cdc test MAR-ice60 viscosite plus faible si snowtemp proche de 0Â°
689	!cdc Si T<264K on reprend l'ancienne fonction
690
691	IF (snowtemp(ji,jj).GT.264.) THEN
692	zviscocity(ji,jj) = ZSNOWCMPCT_W0EXP( ZSNOWCMPCT_WT(tp_00-MIN(tp_00,snowtemp(ji,jj))) + &
693	ZSNOWCMPCT_WR*snowrho(ji,jj) )
694	ELSE
695	zviscocity(ji,jj) = ZSNOWCMPCT_V0EXP( ZSNOWCMPCT_VT(tp_00-MIN(tp_00,snowtemp(ji,jj))) + &
696	ZSNOWCMPCT_VR*snowrho(ji,jj) )
697	ENDIF
698
699	! Calculate snow density: compaction from weight/over-burden
700	! Anderson 1976 method:
701	zsnowrho2(ji,jj) = snowrho(ji,jj) + snowrho(ji,jj)dt_sechiba( &
702	(cte_grav*zsmass(ji,jj)/zviscocity(ji,jj)) &
703	+ zsettle(ji,jj) )
704
705	! Conserve mass by decreasing grid thicknesses in response
706	! to density increases
707	!
708	snowdz(ji,jj) = snowdz(ji,jj)*(snowrho(ji,jj)/zsnowrho2(ji,jj))
709
710	ENDIF
711	ENDDO
712
713	! Update density (kg m-3):
714	snowrho(ji,:) = zsnowrho2(ji,:)
715
716	ENDIF
717
718	ENDDO
719
720	END SUBROUTINE explicitsnow_compactn
721	!!
722	!================================================================================================================================
723	!! SUBROUTINE : explicitsnow_compactn_up
724	!!
725	!>\BRIEF Compute Compaction/Settling
726	!!
727	!! DESCRIPTION :
728	!! Snow compaction due to overburden and settling.
729	!! Mass is unchanged: see Decharme et al. (2016)
730	!!
731	!! RECENT CHANGE(S) : 12/07/2021
732	!!
733	!! MAIN OUTPUT VARIABLE(S): snowrho, snowdz
734	!!
735	!! REFERENCE(S) : Decharme et al. (2016)vi
736	!!
737	!! FLOWCHART : None
738	!! \n
739	!_
740	!================================================================================================================================
741
742
743	SUBROUTINE explicitsnow_compactn_up(kjpindex,u,v,snowtemp,snowrho,snowdz,snowliq)
744
745	!! 0.1 Input variables
746	INTEGER(i_std),INTENT(in) :: kjpindex
747	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u,v
748	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature
749	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowliq !! Liquid water content
750
751	!! 0.2 Output variables
752
753	!! 0.3 Modified variables
754	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
755	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth
756
757	!! 0.4 Local variables
758	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsmass, zsnowrho2, zviscocity
759	REAL(r_std),DIMENSION(kjpindex,nsnow) :: ztemp !! temperature dependance for snow viscosity
760	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zwholdmax !! max water liquid content
761	REAL(r_std),DIMENSION(kjpindex,nsnow) :: Fwliq !! describes the decrease of viscosity in presence of liquid water
762	INTEGER(i_std) :: ji,jj
763
764	!! 1. initialize
765	zsnowrho2(:,:) = snowrho(:,:)
766	zsmass(:,:) = 0.0
767
768	!! 2. Calculating Cumulative snow mass (kg/m2):
769
770	DO ji=1, kjpindex
771	IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
772
773	! 1. Cumulative snow mass (kg/m2):
774	DO jj=2,nsnow
775	zsmass(ji,jj) = zsmass(ji,jj-1) + snowdz(ji,jj) * snowrho(ji,jj-1)
776	ENDDO
777	! 1st layer : half the mass of the uppermost layer applied to itself
778	zsmass(ji,1)= 0.5 * snowdz(ji,1) * snowrho(ji,1)
779
780	! 2. Compaction
781	! Liquid water effect
782	zwholdmax(ji,:) = snow3lhold_1d(snowrho(ji,:),snowdz(ji,:))
783	zwholdmax(ji,:) = max(1.e-10, zwholdmax(ji,:))
784	Fwliq(ji,:) = 1./ (1. + 10.*MIN(1.0,snowliq(ji,:)/zwholdmax(ji,:)))
785
786	! Snow viscocity, density and grid thicknesses
787	DO jj=1,nsnow
788	IF (snowrho(ji,jj) .LT. xrhosmax) THEN
789	! Temperature dependence limited to 5K: Schleef et al. (2014)
790	ztemp(ji,jj) = ZSNOWCMPCT_XT * MIN(DTref,tp_00-MIN(tp_00,snowtemp(ji,jj)))
791
792	! Calculate snow viscocity: Brun et al. (1989), Vionnet et al. (2012)
793	zviscocity(ji,jj) = ZSNOWCMPCT_X0 * Fwliq(ji,jj) * EXP(ztemp(ji,jj) &
794	+ ZSNOWCMPCT_XR * snowrho(ji,jj)) * snowrho(ji,jj)/xrhosref
795	! Calculate new snow density:
796	zsnowrho2(ji,jj) = snowrho(ji,jj) + dt_sechiba &
797	(snowrho(ji,jj)cte_grav*zsmass(ji,jj)/zviscocity(ji,jj))
798	! Conserve mass by decreasing grid thicknesses in response to density increases
799	snowdz(ji,jj) = snowdz(ji,jj)*(snowrho(ji,jj)/zsnowrho2(ji,jj))
800	ENDIF
801	ENDDO
802	! Update density (kg m-3):
803	snowrho(ji,:) = zsnowrho2(ji,:)
804	ENDIF
805	ENDDO
806
807	END SUBROUTINE explicitsnow_compactn_up
808
809	!!
810	!================================================================================================================================
811	!! SUBROUTINE : explicitsnow_drift
812	!!
813	!>\BRIEF Compute Compaction due to snow drift : wind induced densification of near-surface snow layers
814	!!
815	!! DESCRIPTION :
816	!! Snow compaction due to snowdrift
817	!! Mass is unchanged: see Decharme et al. (2016)
818	!!
819	!! RECENT CHANGE(S) : 04/10/2021
820	!!
821	!! MAIN OUTPUT VARIABLE(S): snowrho, snowdz
822	!!
823	!! REFERENCE(S) : Decharme et al. (2016)
824	!!
825	!! FLOWCHART : None
826	!! \n
827	!_
828	!================================================================================================================================
829
830
831	SUBROUTINE explicitsnow_drift(kjpindex,u,v,snowrho,snowdz)
832
833	!! 0.1 Input variables
834	INTEGER(i_std),INTENT(in) :: kjpindex
835	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u,v
836
837	!! 0.2 Output variables
838
839	!! 0.3 Modified variables
840	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
841	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth
842
843	!! 0.4 Local variables
844	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowrho2, zsettle
845	REAL(r_std),DIMENSION(kjpindex) :: zsmobc !! cumulative mobility index
846	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsmob !! mobility index
847	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zwmob !! for mobility index
848	REAL(r_std),DIMENSION(kjpindex,nsnow) :: Zmob !! Mobility index
849	REAL(r_std),DIMENSION(kjpindex,nsnow) :: Zwind !! Wind-driven compaction index
850	REAL(r_std),DIMENSION(kjpindex,nsnow) :: Ztau !! Function used for the calculation of teh compaction rate
851	REAL(r_std),DIMENSION(kjpindex,nsnow) :: tau_cmpct !! Compaction rate (s)
852	REAL(r_std) :: speed !! Wind speed
853	INTEGER(i_std) :: ji,jj
854	LOGICAL :: Zwindup !! True if Zwind(jj-1) > 0
855
856	!! 1. initialize
857	zsnowrho2(:,:) = snowrho(:,:)
858	zsettle(:,:) = 0.
859	zsmobc(:) = 0.0
860
861	!! 2. Calculating Cumulative snow mass (kg/m2):
862
863	DO ji=1, kjpindex
864	Zwindup = .true.
865	IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN
866
867	! Computation of zsmob : Analogous to zsmass (used in Equation B3 Decharme 2016)
868	DO jj=1,nsnow
869	! Mobility index : Equation B1 dans Appendix B Decharme 2016
870	Zmob(ji,jj) = Amob*(1-MAX(0.,(snowrho(ji,jj)-xrhosmin)/xrhosmob)) ! Amob= 1.25, xrhosmin= 50 kg.m-3, xrhosmob= 295 kg.m-3
871	zwmob(ji,jj) = snowdz(ji,jj)*(b_tau - Zmob(ji,jj))
872	zsmob(ji,jj) = zsmobc(ji) + zwmob(ji,jj)
873	zsmobc(ji) = zsmobc(ji) + zwmob(ji,jj)
874	ENDDO
875
876	! Snow viscocity, density and grid thicknesses
877	DO jj=1,nsnow
878	IF (snowrho(ji,jj) .LT. xrhosmax) THEN
879	! Calculate wind densification : Brun 1997
880	! Mobility index : Equation B1 dans Appendix B Decharme 2016
881	Zmob(ji,jj) = Amob*(1-MAX(0.,(snowrho(ji,jj)-xrhosmin)/xrhosmob)) ! Amob= 1.25, xrhosmin= 50 kg.m-3, xrhosmob= 295 kg.m-3
882	! Computation of zsmob : Analogous to zsmass (used in Equation B3 Decharme 2016)
883	zwmob(ji,jj) = snowdz(ji,jj)*(b_tau - Zmob(ji,jj))
884	zsmob(ji,jj) = zsmobc(ji) + zwmob(ji,jj)
885	zsmobc(ji) = zsmobc(ji) + zwmob(ji,jj)
886
887	speed=SQRT(u(ji)u(ji)+v(ji)v(ji)) ! wind speed
888
889	! Wind driven compaction index : Equation B2 in Appendix B Decharme 2016
890	Zwind(ji,jj) = 1 - ACMPCTEXP(-BCMPCTKCMPCT*speed) + Zmob(ji,jj) ! ACMPCT=2.868, BCMPCT=0.085 s.m-1, KCMPCT=1.25
891
892	IF ((Zwind(ji,jj).GT.0.).AND.Zwindup) THEN ! snowdrift occurs only if Zwind>0
893	Ztau(ji,jj) = MAX(0., Zwind(ji,jj))EXP(a_tauzsmob(ji,jj))
894
895	! Wind compaction rate (Equation B3 in Appendix B)
896	tau_cmpct(ji,jj) = 2KCMPCTone_day/Ztau(ji,jj)
897
898	! zsettle : Second term in the right-hand side of Equation 13 in Decharme et al. 2016
899	zsettle(ji,jj) = MAX(0.,(xrhoswind-snowrho(ji,jj))/tau_cmpct(ji,jj)) !xrhoswind=350
900	ELSE
901	Zwindup = .false.
902	Ztau(ji,jj) = 0.
903	tau_cmpct(ji,jj) = 0.
904	zsettle(ji,jj) = 0.
905	ENDIF
906
907	! Calculate new snow density:
908	! settling by wind transport only in case of not too dense snow
909	IF (snowrho(ji,jj).LT.xrhoswind) THEN
910	! settling by wind cannot lead to densities above xrhoswind
911	zsnowrho2(ji,jj) = min(xrhoswind, snowrho(ji,jj) + dt_sechiba * zsettle(ji,jj))
912	! Conserve mass by decreasing grid thicknesses in response to density increases
913	snowdz(ji,jj) = snowdz(ji,jj)*(snowrho(ji,jj)/zsnowrho2(ji,jj))
914	ENDIF
915	ENDIF
916	ENDDO
917	! Update density (kg m-3):
918	snowrho(ji,:) = zsnowrho2(ji,:)
919	ENDIF
920	ENDDO
921
922	END SUBROUTINE explicitsnow_drift
923
924
925	!!
926	!================================================================================================================================
927	!! SUBROUTINE : explicitsnow_transf
928	!!
929	!>\BRIEF Computing snow mass and heat redistribution due to grid thickness configuration resetting
930	!!
931	!! DESCRIPTION : Snow mass and heat redistibution due to grid thickness
932	!! configuration resetting. Total mass and heat content
933	!! of the overall snowpack unchanged/conserved within this routine.
934	!! RECENT CHANGE(S) : None
935	!!
936	!! MAIN OUTPUT VARIABLE(S): None
937	!!
938	!! REFERENCE(S) :
939	!!
940	!! FLOWCHART : None
941	!! \n
942	!_
943	!================================================================================================================================
944
945	SUBROUTINE explicitsnow_transf(kjpindex,snowdz_old,snowdz,snowrho,snowheat,snowgrain)
946
947	!! 0.1 Input variables
948	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
949	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowdz_old !! Snow depth at the previous time step
950
951	!! 0.2 Output variables
952
953	!! 0.3 Modified variables
954	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
955	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m)
956	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow depth (m)
957	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowheat !! Snow heat content/enthalpy (J/m2)
958
959	!! 0.4 Local varibles
960	REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: zsnowzo
961	REAL(r_std),DIMENSION(kjpindex,0:nsnow) :: zsnowzn
962	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowddz
963	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zdelta
964	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zsnowrhon,zsnowheatn,zsnowgrainn
965	REAL(r_std),DIMENSION(kjpindex) :: zsnowmix_delta
966	REAL(r_std),DIMENSION(kjpindex) :: zsumheat, zsumswe, zsumgrain
967	INTEGER(i_std),DIMENSION(nsnow,2) :: locflag
968	REAL(r_std) :: psnow
969	INTEGER(i_std) :: ji,jj,jjj
970
971	! Initialization
972	zsumheat(:) = 0.0
973	zsumswe(:) = 0.0
974	zsumgrain(:) = 0.0
975	zsnowmix_delta(:) = 0.0
976	locflag(:,:) = 0
977
978	DO ji=1, kjpindex
979	psnow = SUM(snowdz(ji,:))
980
981	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
982	zsnowzo(ji,0) = 0.
983	zsnowzn(ji,0) = 0.
984	zsnowzo(ji,1) = snowdz_old(ji,1)
985	zsnowzn(ji,1) = snowdz(ji,1)
986
987	DO jj=2,nsnow
988	zsnowzo(ji,jj) = zsnowzo(ji,jj-1) + snowdz_old(ji,jj)
989	zsnowzn(ji,jj) = zsnowzn(ji,jj-1) + snowdz(ji,jj)
990	ENDDO
991
992	DO jj=1,nsnow
993	IF (jj .EQ. 1) THEN
994	locflag(jj,1)=1
995	ELSE
996	DO jjj=nsnow,1,-1
997	!upper bound of the snow layer
998	IF (zsnowzn(ji,jj-1) .LE. zsnowzo(ji,jjj)) THEN
999	locflag(jj,1) = jjj
1000	ENDIF
1001	ENDDO
1002	ENDIF
1003
1004	IF (jj .EQ. nsnow) THEN
1005	locflag(jj,2)=nsnow
1006	ELSE
1007	DO jjj=nsnow,1,-1
1008	!lower bound of the snow layer
1009	IF (zsnowzn(ji,jj) .LE. zsnowzo(ji,jjj)) THEN
1010	locflag(jj,2) = jjj
1011	ENDIF
1012	ENDDO
1013	ENDIF
1014
1015	!to interpolate
1016	! when heavy snow occurred
1017	IF (locflag(jj,1) .EQ. locflag(jj,2)) THEN
1018	zsnowrhon(ji,jj) = snowrho(ji,locflag(jj,1))
1019
1020	zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,1))*snowdz(ji,jj)/snowdz_old(ji,locflag(jj,1))
1021
1022	zsnowgrainn(ji,jj) = snowgrain(ji,locflag(jj,1))
1023	ELSE
1024	!snow density
1025	zsnowrhon(ji,jj) = snowrho(ji,locflag(jj,1)) * &
1026	(zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1)) + &
1027	snowrho(ji,locflag(jj,2)) * &
1028	(zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))
1029
1030	DO jjj=locflag(jj,1),locflag(jj,2)-1
1031	zsnowrhon(ji,jj) = zsnowrhon(ji,jj) + &
1032	(jjj-locflag(jj,1))snowrho(ji,jjj)snowdz_old(ji,jjj)
1033	ENDDO
1034
1035	zsnowrhon(ji,jj) = zsnowrhon(ji,jj) / snowdz(ji,jj)
1036
1037	!snow heat
1038	IF (snowdz_old(ji,locflag(jj,1)) .GT. 0.0) THEN
1039	zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,1)) * &
1040	(zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1))/snowdz_old(ji,locflag(jj,1)) + &
1041	snowheat(ji,locflag(jj,2)) * &
1042	(zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))/snowdz_old(ji,locflag(jj,2))
1043	ELSE
1044	zsnowheatn(ji,jj) = snowheat(ji,locflag(jj,2))/snowdz_old(ji,locflag(jj,2))*(zsnowzn(ji,locflag(jj,1))-zsnowzo(ji,locflag(jj,1)))
1045	ENDIF
1046
1047	DO jjj=locflag(jj,1),locflag(jj,2)-1
1048	zsnowheatn(ji,jj) = zsnowheatn(ji,jj) + &
1049	(jjj-locflag(jj,1))*snowheat(ji,jjj)
1050	ENDDO
1051
1052	!snow grain
1053	zsnowgrainn(ji,jj) = snowgrain(ji,locflag(jj,1)) * &
1054	(zsnowzo(ji,locflag(jj,1))-zsnowzn(ji,jj-1)) + &
1055	snowgrain(ji,locflag(jj,2)) * &
1056	(zsnowzn(ji,jj)-zsnowzo(ji,locflag(jj,2)-1))
1057
1058	DO jjj=locflag(jj,1),locflag(jj,2)-1
1059	zsnowgrainn(ji,jj) = zsnowgrainn(ji,jj) + &
1060	(jjj-locflag(jj,1))snowgrain(ji,jjj)snowdz_old(ji,jjj)
1061	ENDDO
1062	zsnowgrainn(ji,jj) = zsnowgrainn(ji,jj) / snowdz(ji,jj)
1063	ENDIF
1064	ENDDO
1065	snowrho(ji,:) = zsnowrhon(ji,:)
1066	snowheat(ji,:) = zsnowheatn(ji,:)
1067	snowgrain(ji,:) = zsnowgrainn(ji,:)
1068	ENDIF
1069
1070	! Vanishing or very thin snowpack check:
1071	! -----------------------------------------
1072	!
1073	! NOTE: ONLY for very shallow snowpacks, mix properties (homogeneous):
1074	! this avoids problems related to heat and mass exchange for
1075	! thin layers during heavy snowfall or signifigant melt: one
1076	! new/old layer can exceed the thickness of several old/new layers.
1077	! Therefore, mix (conservative):
1078	!
1079	! modified by Tao Wang
1080	IF (psnow > 0 .AND. (psnow < xsnowcritd .OR. snowdz_old(ji,1) &
1081	& .eq. 0 .OR. snowdz_old(ji,2) .eq. 0 .OR. snowdz_old(ji,3) .eq. 0)) THEN
1082	zsumheat(ji) = SUM(snowheat(ji,:))
1083	zsumswe(ji) = SUM(snowrho(ji,:)*snowdz_old(ji,:))
1084	zsumgrain(ji)= SUM(snowgrain(ji,:)*snowdz_old(ji,:))
1085	zsnowmix_delta(ji) = 1.0
1086	DO jj=1,nsnow
1087	zsnowheatn(ji,jj) = zsnowmix_delta(ji)*(zsumheat(ji)/nsnow)
1088	snowdz(ji,jj) = zsnowmix_delta(ji)*(psnow/nsnow)
1089	zsnowrhon(ji,jj) = zsnowmix_delta(ji)*(zsumswe(ji)/psnow)
1090	zsnowgrainn(ji,jj) = zsnowmix_delta(ji)*(zsumgrain(ji)/psnow)
1091	ENDDO
1092	! Update mass (density and thickness), heat and grain size:
1093	! ------------------------------------------------------------
1094	snowrho(ji,:) = zsnowrhon(ji,:)
1095	snowheat(ji,:) = zsnowheatn(ji,:)
1096	snowgrain(ji,:) = zsnowgrainn(ji,:)
1097	ENDIF
1098	ENDDO
1099
1100	END SUBROUTINE explicitsnow_transf
1101
1102
1103	!!
1104	!================================================================================================================================
1105	!! SUBROUTINE : explicitsnow_fall
1106	!!
1107	!>\BRIEF Computes snowfall
1108	!!
1109	!! DESCRIPTION : Computes snowfall
1110	!routine.
1111	!! RECENT CHANGE(S) : None
1112	!!
1113	!! MAIN OUTPUT VARIABLE(S): None
1114	!!
1115	!! REFERENCE(S) :
1116	!!
1117	!! FLOWCHART : None
1118	!! \n
1119	!_
1120	!================================================================================================================================
1121
1122	SUBROUTINE explicitsnow_fall(kjpindex,precip_snow,temp_air,u,v,totfrac_nobio,snowrho,snowdz,&
1123	& snowheat,snowgrain,snowtemp,psnowhmass)
1124
1125	!! 0.1 Input variables
1126	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size
1127	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: precip_snow !! Snow rate (SWE) (kg/m2 per dt_sechiba)
1128	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_air !! Air temperature
1129	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u,v !! Horizontal wind speed
1130	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowtemp !! Snow temperature
1131	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: totfrac_nobio
1132
1133	!! 0.2 Output variables
1134	REAL(r_std), DIMENSION(kjpindex),INTENT(out) :: psnowhmass !! Heat content of snowfall (J/m2)
1135
1136	!! 0.3 Modified variables
1137	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowrho !! Snow density profile (kg/m3)
1138	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowdz !! Snow layer thickness profile (m)
1139	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowheat !! Snow heat content/enthalpy (J/m2)
1140	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(inout) :: snowgrain !! Snow grain size (m)
1141
1142	!! 0.4 Local variables
1143	REAL(r_std), DIMENSION(kjpindex) :: rhosnew !! Snowfall density
1144	REAL(r_std), DIMENSION(kjpindex) :: dsnowfall !! Snowfall thickness (m)
1145	REAL(r_std), DIMENSION(kjpindex,nsnow) :: snowdz_old
1146	REAL(r_std), DIMENSION(kjpindex) :: snow_depth,snow_depth_old,newgrain
1147	REAL(r_std) :: snowfall_delta,speed
1148	INTEGER(i_std) :: ji,jj
1149
1150	!! 1. initialize the variables
1151
1152	snowdz_old = snowdz
1153	DO ji=1,kjpindex
1154	snow_depth(ji) = SUM(snowdz(ji,:))
1155	ENDDO
1156
1157	snow_depth_old = snow_depth
1158	snowfall_delta = 0.0
1159
1160	!! 2. incorporate snowfall into snowpack
1161	DO ji = 1, kjpindex
1162
1163	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
1164
1165	! new snow fall on snowpack
1166	! NOTE: when the surface temperature is zero, it means that snowfall has no
1167	! heat content but it can be used for increasing the thickness and changing the density (maybe it is a bug)
1168	psnowhmass(ji) = 0.0
1169	IF ( (precip_snow(ji) .GT. 0.0) ) THEN
1170
1171	! calculate
1172	!cdc psnowhmass(ji) = precip_snow(ji)(xci(snowtemp(ji,1)-tp_00)-chalfu0)
1173	psnowhmass(ji) = precip_snow(ji)(xci(temp_air(ji)-tp_00)-chalfu0)
1174
1175	! Snowfall density: Following CROCUS (Pahaut 1976)
1176	rhosnew(ji) = MAX(xrhosmin, snowfall_a_sn + snowfall_b_sn*(temp_air(ji)-tp_00) + &
1177	snowfall_c_sn*SQRT(speed))
1178
1179	! Augment total pack depth:
1180	dsnowfall(ji) = precip_snow(ji)/rhosnew(ji) !snowfall thickness (m)
1181	snow_depth(ji) = snow_depth(ji) + dsnowfall(ji)
1182
1183	! Fresh snowfall changes the snowpack density and liquid content in uppermost layer
1184	IF (dsnowfall(ji) .GT. zero) THEN
1185	snowrho(ji,1) = (snowdz(ji,1)snowrho(ji,1) + dsnowfall(ji)rhosnew(ji))/ &
1186	(snowdz(ji,1)+dsnowfall(ji))
1187
1188	snowdz(ji,1) = snowdz(ji,1) + dsnowfall(ji)
1189
1190	! Add energy of snowfall to snowpack:
1191	! Update heat content (J/m2) (therefore the snow temperature
1192	! and liquid content):
1193	snowheat(ji,1) = snowheat(ji,1) + psnowhmass(ji)
1194
1195	! Incorporate snowfall grain size:
1196	newgrain(ji) = MIN(dgrain_new_max, snow3lgrain_0d(rhosnew(ji)))
1197	snowgrain(ji,1) = (snowdz_old(ji,1)snowgrain(ji,1) + dsnowfall(ji)newgrain(ji))/ &
1198	snowdz(ji,1)
1199	ENDIF
1200	ELSE
1201	dsnowfall(ji) = 0.
1202	ENDIF
1203
1204	! new snow fall on snow free surface.
1205	! we use the linearization for the new snow fall on snow-free ground
1206	IF ( (dsnowfall(ji) .GT. zero) .AND. (snow_depth_old(ji) .EQ. zero) ) THEN
1207	snowfall_delta = 1.0
1208	DO jj=1,nsnow
1209	snowdz(ji,jj) = snowfall_delta*(dsnowfall(ji)/nsnow) + &
1210	(1.0-snowfall_delta)*snowdz(ji,jj)
1211
1212	snowheat(ji,jj) = snowfall_delta*(psnowhmass(ji)/nsnow) + &
1213	(1.0-snowfall_delta)*snowheat(ji,jj)
1214
1215	snowrho(ji,jj) = snowfall_delta*rhosnew(ji) + &
1216	(1.0-snowfall_delta)*snowrho(ji,jj)
1217
1218	snowgrain(ji,jj) = snowfall_delta*newgrain(ji) + &
1219	(1.0-snowfall_delta)*snowgrain(ji,jj)
1220	ENDDO
1221	ENDIF
1222	ENDDO
1223
1224	END SUBROUTINE explicitsnow_fall
1225
1226	!!
1227	!================================================================================================================================
1228	!! SUBROUTINE : explicitsnow_gone
1229	!!
1230	!>\BRIEF Check whether snow is gone
1231	!!
1232	!! DESCRIPTION : If so, set thickness (and therefore mass and heat) and liquid
1233	!! content to zero, and adjust fluxes of water, evaporation and
1234	!! heat into underlying surface.
1235	!! RECENT CHANGE(S) : None
1236	!!
1237	!! MAIN OUTPUT VARIABLE(S): None
1238	!!
1239	!! REFERENCE(S) :
1240	!!
1241	!! FLOWCHART : None
1242	!! \n
1243	!_
1244	!================================================================================================================================
1245
1246	SUBROUTINE explicitsnow_gone(kjpindex,pgflux,&
1247	snowheat,snowtemp,snowdz,snowrho,snowliq,grndflux,snowmelt)
1248
1249	!! 0.1 Input variables
1250	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1251	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pgflux !! Net energy into snow pack(w/m2)
1252	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowheat !! Snow heat content (J/m^2)
1253
1254	!! 0.2 Output variables
1255
1256	!! 0.3 Modified variables
1257	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature
1258	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth
1259	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow density (Kg/m^3)
1260	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowliq !! Liquid water content
1261	REAL(r_std),DIMENSION(kjpindex), INTENT(inout) :: grndflux !! Soil/snow interface heat flux (W/m2)
1262	REAL(r_std),DIMENSION(kjpindex),INTENT(inout) :: snowmelt !! Snow melt
1263	REAL(r_std),DIMENSION(kjpindex) :: thrufal !! Water leaving snowpack(kg/m2/s)
1264
1265	!! 0.4 Local variables
1266	INTEGER(i_std) :: ji,jj
1267	REAL(r_std),DIMENSION(kjpindex) :: snowgone_delta
1268	REAL(r_std),DIMENSION (kjpindex) :: totsnowheat !!snow heat content at each layer
1269	REAL(r_std),DIMENSION(kjpindex) :: snowdepth_crit
1270
1271	! first caculate total snowpack snow heat content
1272	snowgone_delta(:) = un
1273	thrufal(:)=0.0
1274	snowmelt(:)=0
1275	totsnowheat(:) = SUM(snowheat(:,:),2)
1276
1277	DO ji = 1, kjpindex
1278
1279	IF ( (SUM(snowdz(ji,:)) .GT. min_sechiba)) THEN
1280
1281	IF ( pgflux(ji) >= (-totsnowheat(ji)/dt_sechiba) ) THEN
1282	grndflux(ji) = pgflux(ji) + (totsnowheat(ji)/dt_sechiba)
1283	thrufal(ji)=SUM(snowrho(ji,:)*snowdz(ji,:))
1284	snowgone_delta(ji) = 0.0
1285	snowmelt(ji) = snowmelt(ji)+thrufal(ji)
1286	ENDIF
1287
1288	! update of snow state (either still present or not)
1289	DO jj=1,nsnow
1290	snowdz(ji,jj) = snowdz(ji,jj) *snowgone_delta(ji)
1291	snowliq(ji,jj) = snowliq(ji,jj) *snowgone_delta(ji)
1292	snowtemp(ji,jj) = (1.0-snowgone_delta(ji))tp_00 + snowtemp(ji,jj)snowgone_delta(ji)
1293	ENDDO
1294
1295	ELSE
1296	snowdz(ji,:)=0.
1297	snowliq(ji,:)=0.
1298	snowmelt(ji)=snowmelt(ji)+SUM(snowrho(ji,:)*snowdz(ji,:))
1299	!cdc MODIF ICE
1300	grndflux(ji) = pgflux(ji)
1301	ENDIF
1302	ENDDO
1303
1304	END SUBROUTINE explicitsnow_gone
1305
1306	!================================================================================================================================
1307	!! SUBROUTINE : explicitsnow_melt_refrz
1308	!!
1309	!>\BRIEF Computes snow melt and refreezing processes within snowpack
1310	!!
1311	!! DESCRIPTION :
1312	!! RECENT CHANGE(S) : None
1313	!!
1314	!! MAIN OUTPUT VARIABLE(S): None
1315	!!
1316	!! REFERENCE(S) :
1317	!!
1318	!! FLOWCHART : None
1319	!! \n
1320	!_
1321	!================================================================================================================================
1322
1323	SUBROUTINE explicitsnow_melt_refrz(kjpindex,precip_rain,pgflux,soilcap,totfrac_nobio,frac_snow_nobio, &
1324	snowtemp,snowdz,snowrho,snowliq,snowmelt,grndflux,temp_air, frac_snow_veg, veget, veget_max, zrainfall)
1325
1326	!! 0.1 Input variables
1327	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
1328	REAL(r_std),DIMENSION (kjpindex,nsnow) :: pcapa_snow !! Heat capacity for snow
1329	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall
1330	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: temp_air !! Air temperature
1331	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: pgflux !! Net energy into snowpack(w/m2)
1332	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilcap !! Soil heat capacity
1333	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+ ...
1334	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
1335	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegetation
1336	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type
1337	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type
1338
1339	!! 0.2 Output variables
1340	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: zrainfall !! Rain precipitation on snow
1341	!! @tex $(kg m^{-2})$ @endtex
1342	!! 0.3 Modified variables
1343	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature
1344	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth
1345	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow layer density (Kg/m^3)
1346	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowliq !! Liquid water content
1347	REAL(r_std),DIMENSION(kjpindex),INTENT(inout) :: grndflux !! Net energy input to soil
1348	REAL(r_std),DIMENSION (kjpindex), INTENT(inout) :: snowmelt !! Snowmelt
1349
1350	!! 0.4 Local variables
1351	REAL(r_std),DIMENSION (kjpindex) :: meltxs !! Residual snowmelt energy applied to underlying soil
1352	REAL(r_std) :: enerin,melttot,hrain
1353	REAL(r_std),DIMENSION (nsnow) :: zsnowlwe
1354	REAL(r_std),DIMENSION (nsnow) :: flowliq
1355	REAL(r_std),DIMENSION (kjpindex) :: snowmass
1356	REAL(r_std),DIMENSION (nsnow) :: zphase !! Phase change (from ice to water) (J/m2)
1357	REAL(r_std),DIMENSION (nsnow) :: zphase2 !! Phase change (from water to ice)
1358	REAL(r_std),DIMENSION (nsnow) :: zphase3 !! Phase change related with net energy input to snowpack
1359	REAL(r_std),DIMENSION (nsnow) :: zsnowdz !! Snow layer depth
1360	REAL(r_std),DIMENSION (nsnow) :: zsnowmelt !! Snow melt (liquid water) (m)
1361	REAL(r_std),DIMENSION (nsnow) :: zsnowtemp
1362	REAL(r_std),DIMENSION (nsnow) :: zmeltxs !! Excess melt
1363	REAL(r_std),DIMENSION (nsnow) :: zwholdmax !! Maximum liquid water holding (m)
1364	REAL(r_std),DIMENSION (nsnow) :: zcmprsfact !! Compression factor due to densification from melting
1365	REAL(r_std),DIMENSION (nsnow) :: zscap !! Snow heat capacity (J/m3 K)
1366	REAL(r_std),DIMENSION (nsnow) :: zsnowliq !! (m)
1367	REAL(r_std),DIMENSION (nsnow) :: snowtemp_old
1368	REAL(r_std),DIMENSION (0:nsnow) :: zflowliqt !!(m)
1369	REAL(r_std),DIMENSION (kjpindex) :: frac_rain_veg
1370	REAL(r_std) :: zpcpxs
1371	REAL(r_std) :: ztotwcap
1372	REAL(r_std),DIMENSION(kjpindex,nsnow) :: snowdz_old,snowliq_old
1373	INTEGER(i_std) :: jj,ji, iv, jv
1374	REAL(r_std),DIMENSION(nsnow) :: snowdz_old2
1375	REAL(r_std),DIMENSION(nsnow) :: zsnowrho
1376	REAL(r_std),DIMENSION(kjpindex,nsnow) :: zrefrz !! (m)
1377
1378	!initialize
1379	snowdz_old = snowdz
1380	snowliq_old = snowliq
1381	zrainfall(:) = 0.
1382	frac_rain_veg(:) = 0.
1383	zrefrz(:,:) = 0.
1384
1385	DO ji = 1, kjpindex
1386
1387
1388	snowmass(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
1389	IF ((snowmass(ji) .GT. min_sechiba)) THEN
1390
1391	!! 1 snow melting due to positive snowpack snow temperature
1392
1393	!! 1.0 total liquid equivalent water content of each snow layer
1394
1395	zsnowlwe(:) = snowrho(ji,:) * snowdz(ji,:)/ ph2o
1396
1397	!! 1.1 phase change (J/m2)
1398
1399	pcapa_snow(ji,:) = snowrho(ji,:)*xci
1400
1401	zphase(:) = MIN(pcapa_snow(ji,:)MAX(0.0, snowtemp(ji,:)-tp_00) &
1402	snowdz(ji,:), &
1403	MAX(0.0,zsnowlwe(:)-snowliq(ji,:))chalfu0ph2o)
1404
1405	!! 1.2 update snow liq water and temperature if melting
1406
1407	zsnowmelt(:) = zphase(:)/(chalfu0*ph2o)
1408
1409	!! 1.3 cool off the snow layer temperature due to melt
1410
1411	zsnowtemp(:) = snowtemp(ji,:) - zphase(:)/(pcapa_snow(ji,:)* snowdz(ji,:))
1412
1413	snowtemp(ji,:) = MIN(tp_00, zsnowtemp(:))
1414
1415	zmeltxs(:) = (zsnowtemp(:)-snowtemp(ji,:))pcapa_snow(ji,:)snowdz(ji,:)
1416
1417	!! 1.4 loss of snowpack depth and liquid equivalent water
1418
1419	zwholdmax(:) = snow3lhold_1d(snowrho(ji,:),snowdz(ji,:)) ! 1 dimension
1420
1421	!!-SC- Modification de zwholdmax selon Vionnet_et_al_GMD_2012
1422	!!-SC- zwholdmax(:)= 0.05snowdz(ji,:)(1-(snowrho(ji,:)/920))
1423
1424
1425	zcmprsfact(:) = (zsnowlwe(:)-MIN(snowliq(ji,:)+zsnowmelt(:),zwholdmax(:)))/ &
1426	(zsnowlwe(:)-MIN(snowliq(ji,:),zwholdmax(:)))
1427
1428	snowdz(ji,:) = snowdz(ji,:)*zcmprsfact(:)
1429
1430	snowrho(ji,:) = zsnowlwe(:)*ph2o/snowdz(ji,:)
1431
1432	snowliq(ji,:) = snowliq(ji,:) + zsnowmelt(:)
1433
1434
1435	!! 2 snow refreezing process
1436	!! 2.1 freeze liquid water in any layer
1437	zscap(:) = snowrho(ji,:)*xci !J K-1 m-3
1438	zphase2(:) = MIN(zscap(:)* &
1439	MAX(0.0, tp_00 - snowtemp(ji,:))*snowdz(ji,:), &
1440	snowliq(ji,:)chalfu0ph2o)
1441
1442	! warm layer and reduce liquid if freezing occurs
1443	zsnowdz(:) = MAX(xsnowdmin/nsnow, snowdz(ji,:))
1444	snowtemp_old(:) = snowtemp(ji,:)
1445	snowtemp(ji,:) = snowtemp(ji,:) + zphase2(:)/(zscap(:)*zsnowdz(:))
1446
1447	! Reduce liquid portion if freezing occurs:
1448	snowliq(ji,:) = snowliq(ji,:) - ( (snowtemp(ji,:)-snowtemp_old(:))* &
1449	zscap(:)zsnowdz(:)/(chalfu0ph2o) )
1450	!!-SC- Refreezing :
1451	zrefrz(ji,:) = (snowtemp(ji,:)-snowtemp_old(:))zscap(:)zsnowdz(:)/(chalfu0*ph2o)
1452	!!-SC-
1453	snowliq(ji,:) = MAX(snowliq(ji,:), 0.0)
1454	!! 3. thickness change due to snowmelt in excess of holding capacity
1455	zwholdmax(:) = snow3lhold_1d(snowrho(ji,:),snowdz(ji,:)) ! 1 dimension
1456	flowliq(:) = MAX(0.,(snowliq(ji,:)-zwholdmax(:)))
1457	snowliq(ji,:) = snowliq(ji,:) - flowliq(:)
1458	snowdz(ji,:) = snowdz(ji,:) - flowliq(:)*ph2o/snowrho(ji,:)
1459	! to prevent possible very small negative values (machine
1460	! prescision as snow vanishes
1461	snowdz(ji,:) = MAX(0.0, snowdz(ji,:))
1462
1463	!! 4. liquid water flow
1464	ztotwcap = SUM(zwholdmax(:))
1465	! Rain entering snow (m):
1466	!cdc version svn zrainfall = precip_rain(ji)/ph2ofrac_snow_nobio(ji,iice)totfrac_nobio(ji) ! rainfall (m)
1467
1468	!!-SC-- zrainfall-06 (Le bon)
1469	!~ DO jv = 1,nvm
1470	!~ zrainfall = (frac_snow_nobio(ji,iice)*totfrac_nobio(ji) &
1471	!~ + frac_snow_veg(ji)throughfall_by_pft(jv)veget(ji,jv) &
1472	!~ + frac_snow_veg(ji)*(veget_max(ji,jv) - veget(ji,jv)) ) &
1473	!~ * precip_rain(ji)/ph2o
1474	!~ ENDDO
1475	!!-SC-- fin zrainfall-06
1476
1477	!cdc version Christophe zrainfall tableau et somme sur jv
1478	DO jv = 1,nvm
1479	frac_rain_veg(ji) = frac_snow_veg(ji)throughfall_by_pft(jv)veget(ji,jv) &
1480	+ frac_snow_veg(ji)*(veget_max(ji,jv) - veget(ji,jv)) &
1481	+ frac_rain_veg(ji)
1482	ENDDO
1483	zrainfall(ji) = (frac_snow_nobio(ji,iice)*totfrac_nobio(ji) &
1484	+ frac_rain_veg(ji) )* precip_rain(ji)/ph2o
1485
1486	zflowliqt(0) = MIN(zrainfall(ji),ztotwcap)
1487	! Rain assumed to directly pass through the pack to runoff (m):
1488	zpcpxs = zrainfall(ji) - zflowliqt(0)
1489	!
1490	DO jj=1,nsnow
1491	zflowliqt(jj) = flowliq(jj)
1492	ENDDO
1493
1494	! translated into a density increase:
1495	flowliq(:) = 0.0 ! clear this array for work
1496	zsnowliq(:) = snowliq(ji,:) ! reset liquid water content
1497	!
1498
1499	DO jj=1,nsnow
1500	snowliq(ji,jj) = snowliq(ji,jj) + zflowliqt(jj-1)
1501	flowliq(jj) = MAX(0.0, (snowliq(ji,jj)-zwholdmax(jj)))
1502	snowliq(ji,jj) = snowliq(ji,jj) - flowliq(jj)
1503
1504	!Modified by Tao Wang based on previous ISBA-ES scheme
1505	snowrho(ji,jj) = snowrho(ji,jj) + (snowliq(ji,jj) - zsnowliq(jj))* &
1506	& ph2o/MAX(xsnowdmin/nsnow,snowdz(ji,jj))
1507
1508	zflowliqt(jj) = zflowliqt(jj) + flowliq(jj)
1509	ENDDO
1510
1511	snowmelt(ji) = snowmelt(ji) + (zflowliqt(nsnow) + zpcpxs) * ph2o
1512
1513	! excess heat from melting, using it to warm underlying ground to conserve energy
1514	meltxs(ji) = SUM(zmeltxs(:))/dt_sechiba ! (W/m2)
1515
1516	! energy flux into the soil
1517	grndflux(ji) = grndflux(ji) + meltxs(ji)
1518
1519	ELSE
1520	snowdz(ji,:)=0.
1521	snowliq(ji,:)=0.
1522	snowmelt(ji)=snowmelt(ji)+SUM(snowrho(ji,:)*snowdz(ji,:))
1523
1524	!This addition is to get the precipitation that falls on the snow in case the snow has just totally disppeared in explicitsnow_gone.
1525	snowmelt(ji)=snowmelt(ji)+ precip_rain(ji)frac_snow_nobio(ji,iice)totfrac_nobio(ji)
1526	zrefrz(ji,:) = 0.
1527	ENDIF
1528
1529	ENDDO
1530
1531	!cdc sortie du zrainfall :
1532	CALL xios_orchidee_send_field("zrainfall",zrainfall)
1533	!!-SC-Refreezing:
1534	CALL xios_orchidee_send_field("zrefrz",zrefrz)
1535
1536	END SUBROUTINE explicitsnow_melt_refrz
1537
1538
1539	!================================================================================================================================
1540	!! SUBROUTINE : explicitsnow_icemelt
1541	!!
1542	!>\BRIEF Computes ice melt on ice sheet area (no refreezing process)
1543	!!
1544	!! DESCRIPTION :
1545	!! RECENT CHANGE(S) : None
1546	!!
1547	!! MAIN OUTPUT VARIABLE(S): ice_sheet_melt, icetemp, icedz, grndflux
1548	!!
1549	!! REFERENCE(S) :
1550	!!
1551	!! FLOWCHART : None
1552	!! \n
1553	!_
1554	!================================================================================================================================
1555
1556	SUBROUTINE explicitsnow_icemelt(kjpindex,njsc,snowdz, precip_snow, zrainfall, snowmelt, vevapsno, icetemp,icedz,&
1557	grndflux, ice_sheet_melt)
1558
1559	!! 0.1 Input variables
1560	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
1561	INTEGER(i_std),DIMENSION (kjpindex), INTENT(in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
1562	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth
1563	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: precip_snow !! Snow rate (SWE) (kg/m2 per dt_sechiba)
1564	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: zrainfall !! Rain precipitation on snow
1565	!! @tex $(kg m^{-2})$ @endtex
1566	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snowmelt !! Snowmelt
1567	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: vevapsno !! Snow evaporation @tex ($kg m^{-2}$) @endtex
1568
1569	!! 0.2 Output variables
1570	REAL(r_std),DIMENSION(kjpindex),INTENT(out) :: ice_sheet_melt !! Ice melt [mm/dt_sechiba]
1571
1572	!! 0.3 Modified variables
1573	REAL(r_std), DIMENSION (kjpindex,nice), INTENT(inout) :: icetemp !! Snow temperature
1574	REAL(r_std), DIMENSION (kjpindex,nice), INTENT(inout) :: icedz !! Snow depth
1575	REAL(r_std),DIMENSION(kjpindex),INTENT(inout) :: grndflux !! Net energy input to soil
1576
1577	!! 0.4 Local variables
1578	REAL(r_std),DIMENSION (kjpindex) :: meltxs !! Residual snowmelt energy applied to underlying soil
1579	REAL(r_std),DIMENSION (nice) :: zphase_ice !! Phase change (from ice to water) (J/m2)
1580	REAL(r_std),DIMENSION (nice) :: zicemelt !! Ice melt (liquid water) (m)
1581	REAL(r_std),DIMENSION (nice) :: zicetemp
1582	REAL(r_std),DIMENSION (nice) :: zmeltxs !! Excess melt
1583	INTEGER(i_std) :: ji
1584	REAL(r_std),DIMENSION (kjpindex,nice) :: pcapa_ice !! Heat capacity for ice
1585	REAL(r_std), DIMENSION (kjpindex) :: surf_massbal !! surface mass balance [mm/dt_sechiba]
1586
1587
1588	!initialize
1589	ice_sheet_melt(:) = 0.
1590
1591	DO ji = 1, kjpindex
1592
1593	IF (njsc(ji).EQ.13.) THEN ! ice grid point
1594
1595	!! 1 Ice melting due to positive ice temperature
1596
1597	!! 1.1 phase change (J/m2)
1598
1599	! pcapa_ice(ji,:) = 917. * xci
1600	pcapa_ice(ji,:) = rho_ice * (ZICETHRMHEAT1 + ZICETHRMHEAT2*(icetemp(ji,:)-tp_00)) ! formulation as in GRISLI prop_th_icetemp.f90
1601
1602	zphase_ice(:) = pcapa_ice(ji,:)MAX(0.0, icetemp(ji,:)-tp_00) icedz(ji,:)
1603
1604	!! 1.2 ice melt J/m2/(J/kg)=kg/m2 = mm
1605
1606	zicemelt(:) = zphase_ice(:)/chalfu0
1607
1608	!! 1.3 cool off the snow layer temperature due to melt
1609
1610	zicetemp(:) = icetemp(ji,:) - zphase_ice(:)/(pcapa_ice(ji,:)* icedz(ji,:))
1611
1612	icetemp(ji,:) = MIN(tp_00, zicetemp(:))
1613
1614	zmeltxs(:) = (zicetemp(:)-icetemp(ji,:))pcapa_ice(ji,:)icedz(ji,:)
1615
1616	icedz(ji,:) = max(0.,icedz(ji,:) - zicemelt(:)/ rho_ice)
1617
1618	ice_sheet_melt(ji) = sum(zicemelt(:))
1619
1620	! excess heat from melting, using it to warm underlying ground to conserve energy
1621	meltxs(ji) = SUM(zmeltxs(:))/dt_sechiba ! (W/m2)
1622
1623	! energy flux into the soil
1624	grndflux(ji) = grndflux(ji) + meltxs(ji)
1625
1626	ELSE
1627	ice_sheet_melt(ji) = 0.
1628	ENDIF
1629	ENDDO
1630	surf_massbal(:) = precip_snow(:) + zrainfall(:) - snowmelt(:) - ice_sheet_melt(:) - vevapsno(:)
1631	CALL xios_orchidee_send_field("surf_massbal",surf_massbal/dt_sechiba)
1632	CALL xios_orchidee_send_field("ice_sheet_melt",ice_sheet_melt/dt_sechiba)
1633
1634	END SUBROUTINE explicitsnow_icemelt
1635
1636
1637
1638	!================================================================================================================================
1639	!! SUBROUTINE : explicitsnow_icelevels
1640	!!
1641	!>\BRIEF Update icelevels (constant) and icetemp
1642	!!
1643	!! DESCRIPTION :
1644	!! RECENT CHANGE(S) : None
1645	!!
1646	!! MAIN OUTPUT VARIABLE(S): icetemp, icedz
1647	!!
1648	!! REFERENCE(S) :
1649	!!
1650	!! FLOWCHART : None
1651	!! \n
1652	!_
1653	!================================================================================================================================
1654
1655	SUBROUTINE explicitsnow_icelevels(kjpindex,njsc,ice_sheet_melt,icetemp,icedz)
1656
1657	!! 0.1 Input variables
1658
1659	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
1660	INTEGER(i_std),DIMENSION (kjpindex), INTENT(in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
1661	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: ice_sheet_melt !! Ice melt [mm/dt_sechiba]
1662
1663	!! 0.2 Output variables
1664
1665	!! 0.3 Modified variables
1666	REAL(r_std), DIMENSION (kjpindex,nice), INTENT(inout) :: icetemp !! Snow temperature
1667	REAL(r_std), DIMENSION (kjpindex,nice), INTENT(inout) :: icedz !! Snow depth
1668
1669	!! 0.4 Local variables
1670	INTEGER(i_std) :: ji,jj,xx,locxx
1671	REAL(r_std),DIMENSION(nice) :: znt_tmp !! Ice grid layer boundary after melt
1672	REAL(r_std),DIMENSION (kjpindex,nice) :: icetemp_new !! Update ice temperature
1673	REAL(r_std) :: totdz !! melt cumulated on all ice layers
1674
1675	DO ji = 1, kjpindex
1676	! on ne recalcule icetemp que pour le type Ice (classification usda nÂ°13)
1677	IF (njsc(ji).EQ.13) THEN
1678	IF (ice_sheet_melt(ji).GT.min_sechiba) THEN
1679	totdz = 0.
1680	DO jj=1,nice-1
1681	icetemp_new(ji,jj) = icetemp(ji,jj)icedz(ji,jj)/ZSICOEF1(jj) + icetemp(ji,jj+1)(ZSICOEF1(jj)-icedz(ji,jj))/ZSICOEF1(jj)
1682	totdz = totdz + ZSICOEF1(jj)-icedz(ji,jj) ! melt cumulated on all ice layers
1683	ENDDO
1684	icetemp_new(ji,nice) = icetemp(ji,nice)(ZSICOEF1(nice) - totdz)/ZSICOEF1(nice) + tp_00 totdz/ZSICOEF1(nice)
1685	icetemp(ji,:)=icetemp_new(ji,:)
1686	icedz(ji,:)=ZSICOEF1(:)
1687	!cdc temperature base glace = 0Â°C :
1688	! icetemp(ji,nice)= tp_00
1689	!cdc test en forcant la glace Ã 0Â°C sur toute l'epaisseur
1690	!cdc icetemp(ji,:)=tp_00
1691
1692	ENDIF
1693	ENDIF
1694	ENDDO
1695	END SUBROUTINE explicitsnow_icelevels
1696
1697	!================================================================================================================================
1698	!! SUBROUTINE : explicitsnow_levels
1699	!!
1700	!>\BRIEF Computes snow discretization based on given total snow depth
1701	!!
1702	!! DESCRIPTION :
1703	!! RECENT CHANGE(S) : compatible with 3 and 12 layers
1704	!!
1705	!! MAIN OUTPUT VARIABLE(S): snowdz
1706	!!
1707	!! REFERENCE(S) :
1708	!!
1709	!! FLOWCHART : None
1710	!! \n
1711	!_
1712	!================================================================================================================================
1713	SUBROUTINE explicitsnow_levels( kjpindex,snow_thick, snowdz)
1714
1715	!! 0.1 Input variables
1716	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1717	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow_thick !! Total snow depth
1718
1719	!! 0.2 Output variables
1720
1721	!! 0.3 Modified variables
1722	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth
1723
1724	!! 0.4 Local variables
1725	INTEGER(i_std) :: ji, jj
1726
1727	!! parameters for 3 layers snow model
1728	REAL(r_std), PARAMETER, DIMENSION(3) :: ZSGCOEF1 = (/0.25, 0.50, 0.25/) !! Snow grid parameters
1729	REAL(r_std), PARAMETER, DIMENSION(2) :: ZSGCOEF2 = (/0.05, 0.34/) !! Snow grid parameters
1730	REAL(r_std), PARAMETER :: ZSNOWTRANS = 0.20 !! Minimum total snow depth at which surface layer thickness is constant (m)
1731	REAL(r_std), PARAMETER :: XSNOWCRITD = 0.03 !! (m)
1732
1733	!! parameters for 12 layers snon model
1734	REAL(r_std), PARAMETER, DIMENSION(3) :: ZSGCOEF11 = (/0.3, 0.4, 0.3/) !! Snow grid parameters : distribution of snow on layers 6, 7 and 8
1735	REAL(r_std), PARAMETER, DIMENSION(12) :: Zmax = (/0.01, 0.05, 0.15, 0.5, 1., 0., 0., 0., 1., 0.5, 0.1, 0.02/) !! Snow grid parameters : max thickness of each layer
1736	REAL(r_std), PARAMETER, DIMENSION(3) :: ZSGCOEF22 = (/0.3, 0.3, 0.3/)
1737	LOGICAL, PARAMETER, DIMENSION(12) :: mask_d_r = (/.true.,.true.,.true.,.true.,.true.,.false.,.false.,.false.,.true.,.true.,.true.,.true./) ! True for nsnow=1,5 and nsnow=9,12
1738	REAL(r_std), DIMENSION(kjpindex) :: d_r ! sum of the snow depth of layers 6, 7 and 8.
1739
1740
1741	IF (nsnow .eq. 3) THEN
1742	DO ji=1,kjpindex
1743	IF ( snow_thick(ji) .LE. (XSNOWCRITD+0.01)) THEN
1744	snowdz(ji,1) = MIN(0.01, snow_thick(ji)/nsnow)
1745	snowdz(ji,3) = MIN(0.01, snow_thick(ji)/nsnow)
1746	snowdz(ji,2) = snow_thick(ji) - snowdz(ji,1) - snowdz(ji,3)
1747	ENDIF
1748	ENDDO
1749
1750	WHERE ( snow_thick(:) .LE. ZSNOWTRANS .AND. &
1751	snow_thick(:) .GT. (XSNOWCRITD+0.01) )
1752	snowdz(:,1) = snow_thick(:)*ZSGCOEF1(1)
1753	snowdz(:,2) = snow_thick(:)*ZSGCOEF1(2)
1754	snowdz(:,3) = snow_thick(:)*ZSGCOEF1(3)
1755	END WHERE
1756
1757	DO ji = 1,kjpindex
1758	IF (snow_thick(ji) .GT. ZSNOWTRANS) THEN
1759	snowdz(ji,1) = ZSGCOEF2(1)
1760	snowdz(ji,2) = (snow_thick(ji)-ZSGCOEF2(1))*ZSGCOEF2(2) + ZSGCOEF2(1)
1761	! When using simple finite differences, limit the thickness
1762	! factor between the top and 2nd layers to at most 10
1763	snowdz(ji,2) = MIN(10*ZSGCOEF2(1), snowdz(ji,2) )
1764	snowdz(ji,3) = snow_thick(ji) - snowdz(ji,2) - snowdz(ji,1)
1765	ENDIF
1766	ENDDO
1767
1768	ELSE IF (nsnow .eq. 12) THEN ! snow layering as in Decharme 2016 3.1.1
1769	DO ji=1,kjpindex
1770	!! test if snowdz must be updated
1771	IF ((snowdz(ji,1).LT.(0.5min(Zmax(1),snow_thick(ji)/12)).OR.(snowdz(ji,1).GT.(1.5min(Zmax(1),snow_thick(ji)/12)))) &
1772	.OR. (snowdz(ji,2).LT.(0.5min(Zmax(2),snow_thick(ji)/12)).OR.(snowdz(ji,2).GT.(1.5min(Zmax(2),snow_thick(ji)/12)))) &
1773	.OR. (snowdz(ji,12).LT.(0.5min(Zmax(12),snow_thick(ji)/12)).OR.(snowdz(ji,12).GT.(1.5min(Zmax(12),snow_thick(ji)/12))))) THEN
1774
1775	DO jj=1,5
1776	snowdz(ji,jj) = min(Zmax(jj),snow_thick(ji)/12)
1777	ENDDO
1778	DO jj=9,nsnow
1779	snowdz(ji,jj) = min(Zmax(jj),snow_thick(ji)/12)
1780	ENDDO
1781	! d_r(ji) = snow_thick(ji) - sum(snowdz(ji,:),mask=mask_d_r) ! version code std 12 couches
1782	!cdc version compatible avec 3 couches a la compilation
1783	d_r(ji) = snow_thick(ji)
1784	do jj=1,nsnow
1785	if (mask_d_r(jj)) then
1786	d_r(ji) = d_r(ji) - snowdz(ji,jj)
1787	endif
1788	enddo
1789	snowdz(ji,6) = ZSGCOEF11(1)d_r(ji) - min(0., ZSGCOEF22(1)d_r(ji) - snowdz(ji,5))
1790	snowdz(ji,7) = ZSGCOEF11(2)d_r(ji) + min(0., ZSGCOEF22(2)d_r(ji) - snowdz(ji,5)) + min(0.,ZSGCOEF22(2)*d_r(ji) - snowdz(ji,9))
1791	snowdz(ji,8) = ZSGCOEF11(3)d_r(ji) - min(0., ZSGCOEF22(3)d_r(ji) - snowdz(ji,9))
1792	ENDIF
1793	ENDDO
1794	ENDIF
1795
1796	END SUBROUTINE explicitsnow_levels
1797
1798
1799	!!
1800	!================================================================================================================================
1801	!! SUBROUTINE : explicitsnow_profile
1802	!!
1803	!>\BRIEF
1804	!!
1805	!! DESCRIPTION : In this routine solves the numerical soil thermal scheme, ie calculates the new soil temperature profile.
1806	!!
1807	!! RECENT CHANGE(S) : None
1808	!!
1809	!! MAIN OUTPUT VARIABLE(S): snowtemp, temp_sol_add
1810	!!
1811	!! REFERENCE(S) :
1812	!!
1813	!! FLOWCHART : None
1814	!! \n
1815	!_
1816	!================================================================================================================================
1817	SUBROUTINE explicitsnow_profile (kjpindex, cgrnd_snow,dgrnd_snow,lambda_snow,temp_sol_new, snowtemp,snowdz,temp_sol_add)
1818
1819	!! 0. Variables and parameter declaration
1820
1821	!! 0.1 Input variables
1822
1823	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1824	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! skin temperature
1825	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
1826	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
1827	REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
1828	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow layer thickness
1829
1830	!! 0.2 Output variables
1831
1832	!! 0.3 Modified variables
1833	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout) :: snowtemp
1834	REAL(r_std), DIMENSION (kjpindex),INTENT(inout) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K)
1835
1836	!! 0.4 Local variables
1837	INTEGER(i_std) :: ji, jg
1838	!_
1839	!================================================================================================================================
1840	!! 1. Computes the snow temperatures ptn.
1841	DO ji = 1,kjpindex
1842	IF (SUM(snowdz(ji,:)) .GT. 0) THEN
1843	snowtemp(ji,1) = (lambda_snow(ji) * cgrnd_snow(ji,1) + (temp_sol_new(ji)+temp_sol_add(ji))) / &
1844	(lambda_snow(ji) * (un - dgrnd_snow(ji,1)) + un)
1845	temp_sol_add(ji) = zero
1846	DO jg = 1,nsnow-1
1847	snowtemp(ji,jg+1) = cgrnd_snow(ji,jg) + dgrnd_snow(ji,jg) * snowtemp(ji,jg)
1848	ENDDO
1849	ENDIF
1850	ENDDO
1851	IF (printlev>=3) WRITE (numout,*) ' explicitsnow_profile done '
1852
1853	END SUBROUTINE explicitsnow_profile
1854
1855	!!
1856	!================================================================================================================================
1857	!! SUBROUTINE : explicitsnow_iceprofile
1858	!!
1859	!>\BRIEF
1860	!!
1861	!! DESCRIPTION : In this routine solves the numerical ice thermal scheme, ie calculates the new ice temperature profile.
1862	!!
1863	!! RECENT CHANGE(S) : None
1864	!!
1865	!! MAIN OUTPUT VARIABLE(S): icetemp, temp_sol_add
1866	!!
1867	!! REFERENCE(S) :
1868	!!
1869	!! FLOWCHART : None
1870	!! \n
1871	!_
1872	!================================================================================================================================
1873	SUBROUTINE explicitsnow_iceprofile (kjpindex, cgrnd_ice,dgrnd_ice,lambda_ice,temp_sol_new,icedz,&
1874	njsc, snowdz, cgrnd_snow, dgrnd_snow, snowtemp, temp_sol_add, icetemp)
1875
1876	!! 0. Variables and parameter declaration
1877
1878	!! 0.1 Input variables
1879
1880	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1881	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: temp_sol_new !! skin temperature
1882	REAL(r_std), DIMENSION (kjpindex,nice),INTENT(in) :: cgrnd_ice !! Integration coefficient for ice numerical scheme
1883	REAL(r_std), DIMENSION (kjpindex,nice),INTENT(in) :: dgrnd_ice !! Integration coefficient for ice numerical scheme
1884	REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: lambda_ice !! Coefficient of the linear extrapolation of surface temperature
1885	REAL(r_std), DIMENSION (kjpindex,nice),INTENT(in) :: icedz !! ice layer thickness
1886	INTEGER(i_std),DIMENSION (kjpindex), INTENT(in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
1887	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow layer thickness
1888	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT(in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
1889	REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT(in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
1890	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature
1891
1892	!! 0.2 Output variables
1893
1894	!! 0.3 Modified variables
1895	REAL(r_std), DIMENSION (kjpindex),INTENT(inout) :: temp_sol_add !! Additional energy to melt ice for ice ablation case (K)
1896	REAL(r_std), DIMENSION (kjpindex,nice), INTENT(inout) :: icetemp !! Ice temperature
1897
1898	!! 0.4 Local variables
1899	INTEGER(i_std) :: ji, jg
1900	!_
1901	!================================================================================================================================
1902	!! 1. Computes the snow temperatures ptn.
1903	DO ji = 1,kjpindex
1904	IF (njsc(ji).EQ.13.) THEN ! ice grid point
1905	IF (sum(snowdz(ji,:)).GT.0.) THEN ! grid point with snow
1906	icetemp(ji,1) = cgrnd_snow(ji,nsnow) + dgrnd_snow(ji,nsnow) * snowtemp(ji,nsnow)
1907	temp_sol_add(ji) = zero
1908	DO jg = 1,nice-1
1909	icetemp(ji,jg+1) = cgrnd_ice(ji,jg) + dgrnd_ice(ji,jg) * icetemp(ji,jg)
1910	ENDDO
1911	ELSE ! no snow
1912	icetemp(ji,1) = (lambda_ice(ji) * cgrnd_ice(ji,1) + (temp_sol_new(ji)+temp_sol_add(ji))) / &
1913	(lambda_ice(ji) * (un - dgrnd_ice(ji,1)) + un)
1914	temp_sol_add(ji) = zero
1915	DO jg = 1,nice-1
1916	icetemp(ji,jg+1) = cgrnd_ice(ji,jg) + dgrnd_ice(ji,jg) * icetemp(ji,jg)
1917	ENDDO
1918	ENDIF
1919	ELSE
1920	icetemp(ji,:) = tp_00 ! ice temperature on non ice-sheet area (could be an other value, to be tested)
1921	ENDIF
1922	ENDDO
1923	IF (printlev>=3) WRITE (numout,*) ' explicitsnow_iceprofile done '
1924
1925	END SUBROUTINE explicitsnow_iceprofile
1926
1927	!! ================================================================================================================================
1928	!! SUBROUTINE : explicitsnow_read_reftempicefile
1929	!!
1930	!>\BRIEF
1931	!!
1932	!! DESCRIPTION : Read file with longterm ice temperature
1933	!!
1934	!!
1935	!! RECENT CHANGE(S) : None
1936	!!
1937	!! MAIN OUTPUT VARIABLE(S): reftempice : Reference temperature for ice
1938	!!
1939	!! REFERENCE(S) :
1940	!!
1941	!! FLOWCHART : None
1942	!! \n
1943	!_ ================================================================================================================================
1944	SUBROUTINE explicitsnow_read_reftempicefile(kjpindex,lalo,reftempice)
1945
1946	USE interpweight
1947
1948	IMPLICIT NONE
1949
1950	!! 0. Variables and parameter declaration
1951
1952	!! 0.1 Input variables
1953	INTEGER(i_std), INTENT(in) :: kjpindex
1954	REAL(r_std), DIMENSION(kjpindex,2), INTENT(in) :: lalo
1955
1956	!! 0.2 Output variables
1957	REAL(r_std), DIMENSION(kjpindex, nice), INTENT(out) :: reftempice
1958
1959	!! 0.3 Local variables
1960	INTEGER(i_std) :: ib
1961	CHARACTER(LEN=80) :: filename
1962	REAL(r_std),DIMENSION(kjpindex) :: reftempice_file !! Horizontal temperature field interpolated from file [K]
1963	INTEGER(i_std),DIMENSION(kjpindex,8) :: neighbours
1964	REAL(r_std) :: vmin, vmax !! min/max values to use for the
1965	!! renormalization
1966	REAL(r_std), DIMENSION(kjpindex) :: areftemp !! Availability of data for the interpolation
1967	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
1968	!! the file
1969	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
1970	REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable
1971	!! (variabletypevals(1) = -un, not used)
1972	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
1973	!! 'XYKindTime': Input values are kinds
1974	!! of something with a temporal
1975	!! evolution on the dx*dy matrix'
1976	LOGICAL :: nonegative !! whether negative values should be removed
1977	CHARACTER(LEN=50) :: maskingtype !! Type of masking
1978	!! 'nomask': no-mask is applied
1979	!! 'mbelow': take values below maskvals(1)
1980	!! 'mabove': take values above maskvals(1)
1981	!! 'msumrange': take values within 2 ranges;
1982	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
1983	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
1984	!! (normalized by maskvals(3))
1985	!! 'var': mask values are taken from a
1986	!! variable inside the file (>0)
1987	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
1988	!! `maskingtype')
1989	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
1990	REAL(r_std) :: reftemp_norefinf
1991	REAL(r_std) :: reftemp_default !! Default value
1992
1993
1994	!Config Key = SOIL_REFTEMP_FILE
1995	!Config Desc = File with climatological soil temperature
1996	!Config If = READ_REFTEMPICE
1997	!Config Def = reftempice.nc
1998	!Config Help =
1999	!Config Units = [FILE]
2000	filename = 'reftempice.nc'
2001	CALL getin_p('REFTEMPICE_FILE',filename)
2002	variablename = 'icetemp'
2003
2004	IF (printlev >= 3) WRITE(numout,*) " in thermosoil_read_reftempicefile filename '" // TRIM(filename) // &
2005	"' variable name: '" //TRIM(variablename) // "'"
2006
2007	! For this case there are not types/categories. We have 'only' a continuos field
2008	! Assigning values to vmin, vmax
2009
2010	vmin = 0.
2011	vmax = 9999.
2012
2013	! For this file we do not need neightbours!
2014	neighbours = 0
2015
2016	!! Variables for interpweight
2017	! Type of calculation of cell fractions
2018	fractype = 'default'
2019	! Name of the longitude and latitude in the input file
2020	lonname = 'nav_lon'
2021	latname = 'nav_lat'
2022	! Default value when no value is get from input file
2023	reftemp_default = 1.
2024	! Reference value when no value is get from input file
2025	reftemp_norefinf = 1.
2026	! Should negative values be set to zero from input file?
2027	nonegative = .FALSE.
2028	! Type of mask to apply to the input data (see header for more details)
2029	maskingtype = 'nomask'
2030	! Values to use for the masking (here not used)
2031	maskvals = (/ undef_sechiba, undef_sechiba, undef_sechiba /)
2032	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
2033	namemaskvar = ''
2034
2035	CALL interpweight_2Dcont(kjpindex, 0, 0, lalo, resolution, neighbours, &
2036	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
2037	maskvals, namemaskvar, -1, fractype, reftemp_default, reftemp_norefinf, &
2038	reftempice_file, areftemp)
2039	IF (printlev >= 5) WRITE(numout,*)' thermosoil_read_reftempicefile after interpweight_2Dcont'
2040
2041	! Copy reftempice_file temperature to all ground levels
2042	DO ib=1, kjpindex
2043	reftempice(ib, :) = reftempice_file(ib)
2044	END DO
2045
2046	END SUBROUTINE explicitsnow_read_reftempicefile
2047
2048	!! ================================================================================================================================
2049	!! SUBROUTINE : explicitsnow_maxmass
2050	!!
2051	!>\BRIEF
2052	!!
2053	!! DESCRIPTION : limits the snow mass below a threshold (maxmass_snow)
2054	!!
2055	!!
2056	!! RECENT CHANGE(S) : adapted for nsnow=12
2057	!!
2058	!! MAIN OUTPUT VARIABLE(S): snow, snowdz & snowmelt_from_maxmass
2059	!!
2060	!! REFERENCE(S) :
2061	!!
2062	!! FLOWCHART : None
2063	!! \n
2064	!_ ================================================================================================================================
2065	SUBROUTINE explicitsnow_maxmass(kjpindex,snowrho,soilcap,snow,snowdz,snowmelt_from_maxmass)
2066
2067	IMPLICIT NONE
2068	! variables pas necessaires a transmettre ? : maxmass_snow, chalfu0
2069
2070	!! 0.1 Input variables
2071	INTEGER(i_std), INTENT(in) :: kjpindex
2072	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowrho !! Snow density (Kg/m^3)
2073	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilcap !! Soil heat capacity
2074
2075	!! 0.2 Output variables
2076	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: snow !! Snow mass [Kg/m^2]
2077	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth
2078	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: snowmelt_from_maxmass !! snowmelt to limit snow mass under maxmass_snow
2079
2080	!! 0.3 Local variables
2081	INTEGER(i_std) :: ji, jj
2082	INTEGER(i_std) :: nsnowstart, nsnowstop !! start and stop index layer
2083	INTEGER(i_std) :: locjj
2084	REAL(r_std) :: snow_d1k
2085	REAL(r_std),DIMENSION(kjpindex,nsnow) :: WSNOWDZ !! snow mass (kg/m2)
2086	REAL(r_std),DIMENSION(kjpindex) :: SMASSC !! snow mass cumulated
2087	REAL(r_std),DIMENSION(kjpindex,nsnow) :: SMASS !! snow mass cumulated starting from bottom layer
2088	REAL(r_std),DIMENSION(kjpindex) :: ZSNOWMELT_XS !! thickness of the layer that has been partially removed
2089	REAL(r_std),DIMENSION(kjpindex) :: ZSNOWDZ !! new thickness of the snow layer
2090	REAL(r_std),DIMENSION(kjpindex) :: snow_remove !! snow mass of layers that completely disapears
2091
2092
2093	! initialisation
2094	! snowmelt_from_maxmass(:) = zero
2095
2096	IF (nsnow.EQ.3) THEN ! snowmelt_from_maxmass can be applied to layers 3 to 2
2097	nsnowstart=3
2098	nsnowstop=2
2099	ELSEIF (nsnow.EQ.12) THEN ! snowmelt_from_maxmass can be applied to layers 9 to 6 because 10, 11 and 12 are thin
2100	nsnowstart=9
2101	nsnowstop=6
2102	ELSE
2103	WRITE(numout,*) 'explicitsnow_maxmass : nsnow has a non implemented value : ', nsnow
2104	STOP 'explicitsnow_maxmass'
2105	ENDIF
2106
2107	DO ji=1,kjpindex !domain size
2108	snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
2109	IF(snow(ji).GT.maxmass_snow)THEN
2110	IF (snow(ji).GT.(2.5*maxmass_snow)) THEN
2111	! melt much faster for points where accumulation is very very high
2112	snow_d1k = 10. * soilcap(ji) / chalfu0
2113	ELSE IF (snow(ji).GT.(1.5*maxmass_snow)) THEN
2114	! melt much faster for points where accumulation is very high
2115	snow_d1k = six * soilcap(ji) / chalfu0
2116	ELSE IF (snow(ji).GT.(1.2*maxmass_snow)) THEN
2117	! melt faster for points where accumulation is too high
2118	snow_d1k = trois * soilcap(ji) / chalfu0
2119	ELSE
2120	! standard melting
2121	snow_d1k = un * soilcap(ji) / chalfu0
2122	ENDIF
2123	snowmelt_from_maxmass(ji) = MIN((snow(ji) - maxmass_snow),snow_d1k)
2124	! Calculating the snow accumulation
2125	WSNOWDZ(ji,:)= snowdz(ji,:)*snowrho(ji,:) ! WSNOWDZ in kg/m2
2126	SMASSC(ji)= 0.0
2127	DO jj=nsnowstart,nsnowstop,-1
2128	SMASS(ji,jj) = SMASSC(ji) + WSNOWDZ(ji,jj)
2129	SMASSC(ji) = SMASSC(ji) + WSNOWDZ(ji,jj)
2130	ENDDO
2131
2132	! Finding the layer
2133	locjj = nsnowstart ! search the layer starting from nsnowstart
2134	DO jj=nsnowstart,nsnowstop,-1
2135	IF ((SMASS(ji,jj) .LE. snowmelt_from_maxmass(ji)) .AND. (SMASS(ji,jj-1) .GE. snowmelt_from_maxmass(ji)) ) THEN
2136	locjj=jj-1
2137	ENDIF
2138	ENDDO
2139
2140	! Calculating the removal of snow depth
2141	IF (locjj .EQ. nsnowstart) THEN
2142	! ZSNOWMELT_XS : Epaisseur de la couche qui a disparu partiellement
2143	ZSNOWMELT_XS(ji) = snowmelt_from_maxmass(ji)/snowrho(ji,nsnowstart)
2144	ZSNOWDZ(ji) = snowdz(ji,nsnowstart) - ZSNOWMELT_XS(ji)
2145	snowdz(ji,nsnowstart) = MAX(0.0, ZSNOWDZ(ji))
2146	ELSE IF (locjj .LT. nsnowstart) THEN
2147	snow_remove(ji)=0 ! masse de neige des couches qui disparaissent totalement
2148	DO jj=nsnowstart,locjj+1,-1
2149	snow_remove(ji)=snow_remove(ji)+snowdz(ji,jj)*snowrho(ji,jj)
2150	snowdz(ji,jj)=0
2151	ENDDO
2152	ZSNOWMELT_XS(ji) = (snowmelt_from_maxmass(ji) - snow_remove(ji))/snowrho(ji,locjj)
2153	ZSNOWDZ(ji) = snowdz(ji,locjj) - ZSNOWMELT_XS(ji)
2154	snowdz(ji,locjj) = MAX(0.0, ZSNOWDZ(ji))
2155	ENDIF
2156	ENDIF
2157	ENDDO
2158
2159	END SUBROUTINE explicitsnow_maxmass
2160
2161	!! ================================================================================================================================
2162	!! SUBROUTINE : explicitsnow_subli
2163	!!
2164	!>\BRIEF
2165	!!
2166	!! DESCRIPTION : sublimation on snow
2167	!!
2168	!!
2169	!! RECENT CHANGE(S) :
2170	!!
2171	!! MAIN OUTPUT VARIABLE(S): snow, snowdz, subsnownobio, subsinksoil, snowliq, snowtemp, vevapsno
2172	!!
2173	!! REFERENCE(S) :
2174	!!
2175	!! FLOWCHART : None
2176	!! \n
2177	!_ ================================================================================================================================
2178	SUBROUTINE explicitsnow_subli(kjpindex,snowrho,snowdz,vevapsno, &
2179	snowliq,snowtemp,snow,subsnownobio,subsinksoil)
2180
2181	IMPLICIT NONE
2182
2183	!! 0.1 Input variables
2184	INTEGER(i_std), INTENT(in) :: kjpindex
2185	REAL(r_std),DIMENSION(kjpindex,nsnow),INTENT(in) :: snowrho !! Snow density (Kg/m^3)
2186
2187	!! 0.2 Output variables
2188	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow depth
2189	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapsno !! Snow evaporation @tex ($kg m^{-2}$) @endtex
2190	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowliq !! Snow liquid content (m)
2191	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature
2192	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: snow !! Snow mass [Kg/m^2]
2193	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(out) :: subsnownobio !! Sublimation of snow on other surface types (ice, lakes, ...)
2194	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: subsinksoil !! Excess of sublimation as a sink for the soil
2195
2196	!! 0.3 Local variables
2197	INTEGER(i_std) :: ji, jj
2198	INTEGER(i_std) :: nsnowstart, nsnowstop !! start and stop index layer
2199	INTEGER(i_std) :: locjj
2200	REAL(r_std),DIMENSION(kjpindex,nsnow) :: WSNOWDZ !! snow mass (kg/m2)
2201	REAL(r_std),DIMENSION(kjpindex) :: SMASSC !! snow mass cumulated
2202	REAL(r_std),DIMENSION(kjpindex,nsnow) :: SMASS !! snow mass cumulated starting from bottom layer
2203	REAL(r_std),DIMENSION(kjpindex) :: ZSNOWEVAPS !! thickness of the layer that has been removed by sublimation
2204	REAL(r_std),DIMENSION(kjpindex) :: ZSNOWDZ !! new thickness of the snow layer
2205	REAL(r_std), DIMENSION (kjpindex) :: snowacc
2206
2207	snow(:) = 0.0
2208	DO ji=1,kjpindex !domain size
2209	snow(ji) = SUM(snowrho(ji,:) * snowdz(ji,:))
2210	ENDDO
2211
2212	subsnownobio(:,:) = zero
2213	subsinksoil(:) = zero
2214
2215	DO ji=1, kjpindex ! domain size
2216
2217	IF (vevapsno(ji) .GT. snow(ji)) THEN
2218	subsinksoil (ji) = vevapsno(ji) - snow(ji)
2219	vevapsno(ji) = snow(ji)
2220	snow(ji) = zero
2221	snowdz(ji,:) = 0
2222	snowliq(ji,:) = 0
2223	snowtemp(ji,:) = tp_00
2224	ELSE
2225	! Calculating the snow accumulation
2226	WSNOWDZ(ji,:)= snowdz(ji,:)*snowrho(ji,:)
2227	SMASSC (ji)= 0.0
2228	DO jj=1,nsnow
2229	SMASS(ji,jj) = SMASSC(ji) + WSNOWDZ(ji,jj)
2230	SMASSC(ji) = SMASSC(ji) + WSNOWDZ(ji,jj)
2231	ENDDO
2232	! Finding the layer
2233	locjj=1
2234	DO jj=1,nsnow-1
2235	IF ((SMASS(ji,jj) .LE. vevapsno(ji)) .AND. (SMASS(ji,jj+1) .GE. vevapsno(ji)) ) THEN
2236	locjj=jj+1
2237	ENDIF
2238	ENDDO
2239
2240	! Calculating the removal of snow depth
2241	IF (locjj .EQ. 1) THEN
2242	ZSNOWEVAPS(ji) = vevapsno(ji)/snowrho(ji,1)
2243	ZSNOWDZ(ji) = snowdz(ji,1) - ZSNOWEVAPS(ji)
2244	snowdz(ji,1) = MAX(0.0, ZSNOWDZ(ji))
2245	ELSE IF (locjj .GT. 1) THEN
2246	snowacc(ji)=0
2247	DO jj=1,locjj-1
2248	snowacc(ji)=snowacc(ji)+snowdz(ji,jj)*snowrho(ji,jj)
2249	snowdz(ji,jj)=0
2250	ENDDO
2251	ZSNOWEVAPS(ji) = (vevapsno(ji)-snowacc(ji))/snowrho(ji,locjj)
2252	ZSNOWDZ(ji) = snowdz(ji,locjj) - ZSNOWEVAPS(ji)
2253	snowdz(ji,locjj) = MAX(0.0, ZSNOWDZ(ji))
2254	! ELSE
2255	! ZSNOWEVAPS(ji) = subsnowveg(ji)/snowrho(ji,1)
2256	! ZSNOWDZ(ji) = snowdz(ji,1) - ZSNOWEVAPS(ji)
2257	! snowdz(ji,1) = MAX(0.0, ZSNOWDZ(ji))
2258	ENDIF
2259
2260	ENDIF
2261	ENDDO
2262	END SUBROUTINE explicitsnow_subli
2263
2264
2265	!! ================================================================================================================================
2266	!! SUBROUTINE : explicitsnow_age
2267	!!
2268	!>\BRIEF
2269	!!
2270	!! DESCRIPTION : compute snow age for albedo
2271	!!
2272	!!
2273	!! RECENT CHANGE(S) :
2274	!!
2275	!! MAIN OUTPUT VARIABLE(S): snowage, snownobioage
2276	!!
2277	!! REFERENCE(S) :
2278	!!
2279	!! FLOWCHART : None
2280	!! \n
2281	!_ ================================================================================================================================
2282	SUBROUTINE explicitsnow_age(kjpindex,snow,precip_snow,precip_rain,frac_snow_nobio, &
2283	temp_sol_new,snow_age,snow_nobio_age)
2284
2285	IMPLICIT NONE
2286
2287	!! 0.1 Input variables
2288	INTEGER(i_std), INTENT(in) :: kjpindex
2289	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass [Kg/m^2]
2290	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_snow !! Snowfall
2291	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall
2292	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area
2293	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_sol_new !! Surface temperature
2294
2295	!! 0.2 Output variables
2296	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow_age !! Snow age
2297	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio_age !! Snow age on ice, lakes, ...
2298
2299
2300	!! 0.3 Local variables
2301	INTEGER(i_std) :: ji
2302	REAL(r_std), DIMENSION (kjpindex) :: d_age !! Snow age change
2303	REAL(r_std), DIMENSION (kjpindex) :: xx !! Temporary
2304
2305
2306	DO ji = 1, kjpindex
2307
2308	!! 5.1. Snow age on land
2309
2310	IF (snow(ji) .LE. zero) THEN
2311	snow_age(ji) = zero
2312	ELSE
2313	snow_age(ji) =(snow_age(ji) + (un - snow_age(ji)/max_snow_age) * dt_sechiba/one_day) &
2314	& * EXP(-precip_snow(ji) / snow_trans)
2315	ENDIF
2316
2317	!! 5.2. Snow age on land ice
2318
2319	!! Age of snow on ice: a little bit different because in cold regions, we really
2320	!! cannot negect the effect of cold temperatures on snow metamorphism any more.
2321
2322	IF ( (frac_snow_nobio(ji,iice) .LE. zero) .OR. (snow(ji) .LE. zero) ) THEN
2323	snow_nobio_age(ji,iice) = zero
2324	ELSE
2325
2326	d_age(ji) = ( snow_nobio_age(ji,iice) + &
2327	& (un - snow_nobio_age(ji,iice)/max_snow_age) * dt_sechiba/one_day ) * &
2328	& EXP(-precip_snow(ji) / snow_trans_nobio) - snow_nobio_age(ji,iice)
2329
2330	IF (d_age(ji) .GT. 0. ) THEN
2331	xx(ji) = MAX( tp_00 - temp_sol_new(ji), zero )
2332	!!-SC xx(ji) = ( xx(ji) / 7._r_std ) ** 4._r_std ! paramÃ©trisation standard
2333	!!-SC Pour ralentir le vieillissemnt (surtout pour les basses tempÃ©ratures) :
2334	xx(ji) = ( xx(ji) / omg1 ) ** omg2 !! omg1,omg2 values in run.def
2335	d_age(ji) = d_age(ji) / (un+xx(ji))
2336
2337	!cdc vieillissement accelere par 2 si il pleut :
2338	IF (precip_rain(ji) .GT.0.) THEN
2339	d_age(ji) = d_age(ji)*2.
2340	ENDIF
2341	ENDIF
2342
2343	snow_nobio_age(ji,iice) = MAX( snow_nobio_age(ji,iice) + d_age(ji), zero )
2344	ENDIF
2345	ENDDO
2346
2347	END SUBROUTINE explicitsnow_age
2348
2349	END MODULE explicitsnow

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

Download in other formats: