SURFEX  V8_0
Surfex V8_0 release
 All Classes Files Functions Variables
preps_for_meb_ebud_rad.F90
Go to the documentation of this file.
1 !SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
3 !SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
4 !SFX_LIC for details. version 1.
5 ! ############################################################################
6 SUBROUTINE preps_for_meb_ebud_rad(PPS, &
7  plaicv,psnowrho,psnowswe,psnowheat, &
8  psnowtemp,psnowdz,pscond,pheatcaps,pemisnow,psigma_f,pchip, &
9  ptstep,psr,pta,pvmod,psnowage,ppermsnowfrac )
10 ! ############################################################################
11 !
12 !!**** *PREPS_FOR_MEB_EBUD_RAD*
13 !!
14 !! PURPOSE
15 !! -------
16 !
17 ! Get preliminary estimates of certain parameters needed for energy budget
18 ! solution of snowpack, and some other misc inputs needed by radiation
19 ! routines for MEB.
20 !
21 !!** METHOD
22 !! ------
23 !
24 !
25 !! EXTERNAL
26 !! --------
27 !!
28 !!
29 !! IMPLICIT ARGUMENTS
30 !! ------------------
31 !!
32 !!
33 !!
34 !! REFERENCE
35 !! ---------
36 !!
37 !!
38 !! AUTHOR
39 !! ------
40 !!
41 !! A. Boone * CNRM-GAME, Meteo-France *
42 !!
43 !! MODIFICATIONS
44 !! -------------
45 !! Original 02/2011
46 !-------------------------------------------------------------------------------
47 !
48 !* 0. DECLARATIONS
49 ! ------------
50 !
51 !
52 USE modd_csts, ONLY : xtt, xlmtt, xrholw
53 !
54 USE modd_snow_par, ONLY : xrhosmax_es, xrhosmin_es, xemissn
55 !
56 USE modd_snow_metamo, ONLY : xsnowdzmin
57 !
58 USE modd_surf_par, ONLY : xundef
59 !
61 !
62 USE mode_meb, ONLY : meb_shield_factor
63 !
64 USE yomhook ,ONLY : lhook, dr_hook
65 USE parkind1 ,ONLY : jprb
66 !
67 IMPLICIT NONE
68 !
69 !
70 !* 0.1 Declaration of Arguments
71 !
72 REAL :: ptstep ! time step (s)
73 REAL, DIMENSION(:), INTENT(IN) :: plaicv
74 REAL, DIMENSION(:), INTENT(IN) :: pps
75 REAL, DIMENSION(:), INTENT(IN) :: psr
76 REAL, DIMENSION(:), INTENT(IN) :: pta
77 REAL, DIMENSION(:), INTENT(IN) :: pvmod
78 REAL, DIMENSION(:), INTENT(IN) :: ppermsnowfrac
79 REAL, DIMENSION(:,:), INTENT(IN) :: psnowheat
80 
81 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowswe, psnowage, psnowrho
82 
83 REAL, DIMENSION(:), INTENT(OUT) :: psigma_f, pchip
84 REAL, DIMENSION(:), INTENT(OUT) :: pemisnow
85 REAL, DIMENSION(:,:), INTENT(OUT) :: psnowdz, pscond, pheatcaps, psnowtemp
86 !
87 !
88 !* 0.2 declarations of local variables
89 !
90 INTEGER :: ji, jk
91 REAL, DIMENSION(SIZE(PLAICV,1)) :: zpsna
92 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowliq, zsnowheat, zsnowdzn
93 REAL, DIMENSION(SIZE(PTA)) :: zsnow, zsnowhmass
94 !
95 REAL(KIND=JPRB) :: zhook_handle
96 !----------------------------------------------------
97 ! 0) Initialization
98 !
99 IF (lhook) CALL dr_hook('PREPS_FOR_MEB_EBUD_RAD',0,zhook_handle)
100 !
101 ! Here, as in snow3l (ISBA-ES), we account for several processes
102 ! on the snowpack before surface energy budget computations
103 ! (i.e. snowfall on albedo, density, thickness, and compaction etc...)
104 !
105 ! First, since snow might not exist, check using snow density as a
106 ! proxy (since several computations below depend on this value). Note
107 ! that if snow doesn't exist, the variables below depending on snow
108 ! density are either zero (e.g. PSNOWDZ) or unused (e.g. PSCOND)
109 !
110 WHERE(psnowrho(:,:)==xundef)
111  psnowrho(:,:) = xrhosmin_es
112 ENDWHERE
113 !
114 psnowdz(:,:) = psnowswe(:,:)/psnowrho(:,:) ! diagnose current layer thicknesses (m)
115 !
116 ! Local working:
117 !
118 zsnowheat(:,:) = psnowheat(:,:)*psnowswe(:,:)/psnowrho(:,:) ! J/m3 to J/m2
119 
120 zsnow(:) = 0.0
121 DO jk=1,SIZE(psnowdz,2)
122  DO ji=1,SIZE(psnowdz,1)
123  zsnow(ji) = zsnow(ji) + psnowdz(ji,jk)
124  ENDDO
125 ENDDO
126 !
127  CALL snow3lfall(ptstep,psr,pta,pvmod,zsnow,psnowrho,psnowdz, &
128  zsnowheat,zsnowhmass,psnowage,ppermsnowfrac )
129 
130  CALL snow3lgrid(zsnowdzn,zsnow,psnowdz_old=psnowdz)
131 
132  CALL snow3ltransf(zsnow,psnowdz,zsnowdzn,psnowrho,zsnowheat,psnowage)
133 
134 ! Snow heat capacity:
135 !
136 pheatcaps(:,:) = snow3lscap(psnowrho) ! J m-3 K-1
137 !
138 ! Snow temperature (K)
139 !
140 psnowtemp(:,:) = xtt + ( ((zsnowheat(:,:)/max(1.e-10,psnowdz(:,:))) &
141  + xlmtt*psnowrho(:,:))/pheatcaps(:,:) )
142 !
143 zsnowliq(:,:) = max(0.0,psnowtemp(:,:)-xtt)*pheatcaps(:,:)* &
144  psnowdz(:,:)/(xlmtt*xrholw)
145 
146 psnowtemp(:,:) = min(xtt,psnowtemp(:,:))
147 
148 ! SWE:
149 
150 psnowswe(:,:) = psnowdz(:,:)*psnowrho(:,:)
151 
152  CALL snow3lcompactn(ptstep,xsnowdzmin,psnowrho,psnowdz,psnowtemp,zsnow,zsnowliq)
153 
154 ! Snow thermal conductivity:
155 !
156  CALL snow3lthrm(psnowrho,pscond,psnowtemp,pps)
157 !
158 ! View factor: (1 - shielding factor)
159 !
160 zpsna(:) = 0.
161 pchip(:) = meb_shield_factor(plaicv,zpsna)
162 psigma_f(:) = 1.0 - pchip(:)
163 !
164 ! snow emissivity
165 !
166 pemisnow(:) = xemissn
167 !
168 IF (lhook) CALL dr_hook('PREPS_FOR_MEB_EBUD_RAD',1,zhook_handle)
169 !
170  CONTAINS
171 !####################################################################
172 SUBROUTINE snow3lfall(PTSTEP,PSR,PTA,PVMOD,PSNOW,PSNOWRHO,PSNOWDZ, &
173  psnowheat,psnowhmass,psnowage,ppermsnowfrac)
174 !
175 !! PURPOSE
176 !! -------
177 ! Calculate changes to snowpack resulting from snowfall.
178 ! Update mass and heat content of uppermost layer.
179 !
180 !
181 USE modd_csts, ONLY : xlmtt, xtt, xci
182 USE modd_snow_par, ONLY : xrhosmin_es, xsnowdmin, &
183  xsnowfall_a_sn, &
184  xsnowfall_b_sn, &
185  xsnowfall_c_sn
186 !
187 USE mode_snow3l
188 !
189 IMPLICIT NONE
190 !
191 !* 0.1 declarations of arguments
192 !
193 REAL, INTENT(IN) :: ptstep
194 !
195 REAL, DIMENSION(:), INTENT(IN) :: psr, pta, pvmod, ppermsnowfrac
196 !
197 REAL, DIMENSION(:), INTENT(INOUT) :: psnow
198 !
199 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowrho, psnowdz, psnowheat, psnowage
200 !
201 REAL, DIMENSION(:), INTENT(OUT) :: psnowhmass
202 !
203 !
204 !* 0.2 declarations of local variables
205 !
206 INTEGER :: jj, ji
207 !
208 INTEGER :: ini
209 INTEGER :: inlvls
210 !
211 REAL, DIMENSION(SIZE(PTA)) :: zsnowfall, zrhosnew, &
212  zsnow, zsnowtemp, &
213  zsnowfall_delta, zscap, &
214  zagenew
215 !
216 REAL(KIND=JPRB) :: zhook_handle
217 !
218 !-------------------------------------------------------------------------------
219 !
220 ! 0. Initialize:
221 ! ------------------
222 !
223 IF (lhook) CALL dr_hook('SNOW3LFALL',0,zhook_handle)
224 !
225 ini = SIZE(psnowdz(:,:),1)
226 inlvls = SIZE(psnowdz(:,:),2)
227 !
228 zrhosnew(:) = xrhosmin_es
229 zagenew(:) = 0.0
230 zsnowfall(:) = 0.0
231 zscap(:) = 0.0
232 zsnow(:) = psnow(:)
233 !
234 psnowhmass(:) = 0.0
235 !
236 ! 1. Incorporate snowfall into snowpack:
237 ! --------------------------------------
238 !
239 !
240 ! Heat content of newly fallen snow (J/m2):
241 ! NOTE for now we assume the snowfall has
242 ! the temperature of the snow surface upon reaching the snow.
243 ! This is done as opposed to using the air temperature since
244 ! this flux is quite small and has little to no impact
245 ! on the time scales of interest. If we use the above assumption
246 ! then, then the snowfall advective heat flux is zero.
247 !
248 zsnowtemp(:) = xtt
249 zscap(:) = snow3lscap(psnowrho(:,1))
250 !
251 WHERE (psr(:) > 0.0 .AND. psnowdz(:,1)>0.)
252  zsnowtemp(:) = xtt + (psnowheat(:,1) + &
253  xlmtt*psnowrho(:,1)*psnowdz(:,1))/ &
254  (zscap(:)*max(xsnowdmin/inlvls,psnowdz(:,1)))
255  zsnowtemp(:) = min(xtt, zsnowtemp(:))
256 END WHERE
257 !
258 WHERE (psr(:) > 0.0)
259 !
260  psnowhmass(:) = psr(:)*(xci*(zsnowtemp(:)-xtt)-xlmtt)*ptstep
261 !
262 ! Snowfall density: Following CROCUS (Pahaut 1976)
263 !
264  zrhosnew(:) = max(xrhosmin_es, xsnowfall_a_sn + xsnowfall_b_sn*(pta(:)-xtt)+ &
265  xsnowfall_c_sn*sqrt(pvmod(:)))
266 !
267 !
268 ! Fresh snowfall changes the snowpack age,
269 ! decreasing in uppermost snow layer (mass weighted average):
270 !
271  psnowage(:,1) = (psnowage(:,1)*psnowdz(:,1)*psnowrho(:,1)+zagenew(:)*psr(:)*ptstep) / &
272  (psnowdz(:,1)*psnowrho(:,1)+psr(:)*ptstep)
273 !
274 ! Augment total pack depth:
275 !
276  zsnowfall(:) = psr(:)*ptstep/zrhosnew(:) ! snowfall thickness (m)
277 !
278  psnow(:) = psnow(:) + zsnowfall(:)
279 !
280 ! Fresh snowfall changes the snowpack
281 ! density, increases the total liquid water
282 ! equivalent: in uppermost snow layer:
283 !
284  psnowrho(:,1) = (psnowdz(:,1)*psnowrho(:,1) + zsnowfall(:)*zrhosnew(:))/ &
285  (psnowdz(:,1)+zsnowfall(:))
286 !
287  psnowdz(:,1) = psnowdz(:,1) + zsnowfall(:)
288 !
289 ! Add energy of snowfall to snowpack:
290 ! Update heat content (J/m2) (therefore the snow temperature
291 ! and liquid content):
292 !
293  psnowheat(:,1) = psnowheat(:,1) + psnowhmass(:)
294 !
295 END WHERE
296 !
297 !
298 ! 2. Case of new snowfall on a previously snow-free surface:
299 ! ----------------------------------------------------------
300 !
301 ! When snow first falls on a surface devoid of snow,
302 ! redistribute the snow mass throughout the 3 layers:
303 ! (temperature already set in the calling routine
304 ! for this case)
305 !
306 zsnowfall_delta(:) = 0.0
307 WHERE(zsnow(:) == 0.0 .AND. psr(:) > 0.0)
308  zsnowfall_delta(:) = 1.0
309 END WHERE
310 !
311 DO jj=1,inlvls
312  DO ji=1,ini
313 !
314  psnowdz(ji,jj) = zsnowfall_delta(ji)*(zsnowfall(ji) /inlvls) + &
315  (1.0-zsnowfall_delta(ji))*psnowdz(ji,jj)
316 !
317  psnowheat(ji,jj) = zsnowfall_delta(ji)*(psnowhmass(ji)/inlvls) + &
318  (1.0-zsnowfall_delta(ji))*psnowheat(ji,jj)
319 !
320  psnowrho(ji,jj) = zsnowfall_delta(ji)*zrhosnew(ji) + &
321  (1.0-zsnowfall_delta(ji))*psnowrho(ji,jj)
322 !
323  psnowage(ji,jj) = zsnowfall_delta(ji)*(zagenew(ji)/inlvls) + &
324  (1.0-zsnowfall_delta(ji))*psnowage(ji,jj)
325 !
326  ENDDO
327 ENDDO
328 !
329 IF (lhook) CALL dr_hook('SNOW3LFALL',1,zhook_handle)
330 !
331 !
332 END SUBROUTINE snow3lfall
333 !####################################################################
334  SUBROUTINE snow3ltransf(PSNOW,PSNOWDZ,PSNOWDZN, &
335  psnowrho,psnowheat,psnowage)
336 !
337 !! PURPOSE
338 !! -------
339 ! Snow mass,heat and characteristics redistibution in case of
340 ! grid resizing. Total mass and heat content of the overall snowpack
341 ! unchanged/conserved within this routine.
342 ! Same method as in Crocus
343 !
344 USE modd_snow_par, ONLY : xsnowcritd
345 !
346 IMPLICIT NONE
347 !
348 !
349 !* 0.1 declarations of arguments
350 !
351 REAL, DIMENSION(: ), INTENT(IN) :: psnow
352 !
353 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowdzn
354 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowheat
355 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowrho
356 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowdz
357 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowage
358 !
359 !* 0.2 declarations of local variables
360 !
361 INTEGER :: ji, jl, jlo
362 !
363 INTEGER :: ini
364 INTEGER :: inlvls
365 !
366 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowrhon
367 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowheatn
368 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowagen
369 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowztop_new
370 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowzbot_new
371 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowrhoo
372 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowheato
373 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowageo
374 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowdzo
375 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowztop_old
376 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowzbot_old
377 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowhean
378 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowagn
379 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zmastotn
380 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zmassdzo
381 !
382 REAL, DIMENSION(SIZE(PSNOW)) :: zpsnow_old, zpsnow_new
383 REAL, DIMENSION(SIZE(PSNOW)) :: zsumheat, zsumswe, zsumage, zsnowmix_delta
384 !
385 REAL :: zpropor
386 !
387 REAL(KIND=JPRB) :: zhook_handle
388 !
389 !-------------------------------------------------------------------------------
390 !
391 ! 0. Initialization:
392 ! ------------------
393 !
394 !
395 IF (lhook) CALL dr_hook('SNOW3LTRANSF',0,zhook_handle)
396 !
397 ini = SIZE(psnowrho,1)
398 inlvls = SIZE(psnowrho,2)
399 !
400 zpsnow_new(:) = 0.0
401 zpsnow_old(:) = psnow(:)
402 !
403 DO jl=1,inlvls
404  DO ji=1,ini
405  zpsnow_new(ji)=zpsnow_new(ji)+psnowdzn(ji,jl)
406  ENDDO
407 ENDDO
408 !
409 ! initialization of variables describing the initial snowpack
410 !
411 zsnowdzo(:,:) = psnowdz(:,:)
412 zsnowrhoo(:,:) = psnowrho(:,:)
413 zsnowheato(:,:) = psnowheat(:,:)
414 zsnowageo(:,:) = psnowage(:,:)
415 zmassdzo(:,:) = xundef
416 !
417 ! 1. Calculate vertical grid limits (m):
418 ! --------------------------------------
419 !
420 zsnowztop_old(:,1) = zpsnow_old(:)
421 zsnowztop_new(:,1) = zpsnow_new(:)
422 zsnowzbot_old(:,1) = zsnowztop_old(:,1)-zsnowdzo(:,1)
423 zsnowzbot_new(:,1) = zsnowztop_new(:,1)-psnowdzn(:,1)
424 !
425 DO jl=2,inlvls
426  DO ji=1,ini
427  zsnowztop_old(ji,jl) = zsnowzbot_old(ji,jl-1)
428  zsnowztop_new(ji,jl) = zsnowzbot_new(ji,jl-1)
429  zsnowzbot_old(ji,jl) = zsnowztop_old(ji,jl )-zsnowdzo(ji,jl)
430  zsnowzbot_new(ji,jl) = zsnowztop_new(ji,jl )-psnowdzn(ji,jl)
431  ENDDO
432 ENDDO
433 zsnowzbot_old(:,inlvls)=0.0
434 zsnowzbot_new(:,inlvls)=0.0
435 
436 WHERE(psnowdzn(:,:)==0.)
437  zsnowztop_old(:,:) = 0.
438  zsnowztop_new(:,:) = 0.
439  zsnowzbot_old(:,:) = 0.
440  zsnowzbot_new(:,:) = 0.
441 END WHERE
442 !
443 ! 3. Calculate mass, heat, charcateristics mixing due to vertical grid resizing:
444 ! --------------------------------------------------------------------
445 !
446 ! loop over the new snow layers
447 ! Summ or avergage of the constituting quantities of the old snow layers
448 ! which are totally or partially inserted in the new snow layer
449 ! For snow age, mass weighted average is used.
450 !
451 zsnowhean(:,:)=0.0
452 zmastotn(:,:)=0.0
453 zsnowagn(:,:)=0.0
454 !
455 DO jl=1,inlvls
456  DO jlo=1, inlvls
457  DO ji=1,ini
458  IF((zsnowztop_old(ji,jlo)>zsnowzbot_new(ji,jl)).AND.(zsnowzbot_old(ji,jlo)<zsnowztop_new(ji,jl)))THEN
459 !
460  zpropor = (min(zsnowztop_old(ji,jlo), zsnowztop_new(ji,jl)) &
461  - max(zsnowzbot_old(ji,jlo), zsnowzbot_new(ji,jl)))&
462  / zsnowdzo(ji,jlo)
463 !
464  zmassdzo(ji,jlo)=zsnowrhoo(ji,jlo)*zsnowdzo(ji,jlo)*zpropor
465 !
466  zmastotn(ji,jl)=zmastotn(ji,jl)+zmassdzo(ji,jlo)
467  zsnowagn(ji,jl)=zsnowagn(ji,jl)+zsnowageo(ji,jlo)*zmassdzo(ji,jlo)
468 !
469  zsnowhean(ji,jl)=zsnowhean(ji,jl)+zsnowheato(ji,jlo)*zpropor
470 !
471  ENDIF
472  ENDDO
473  ENDDO
474 ENDDO
475 !
476 ! the new layer inherits from the weighted average properties of the old ones
477 ! heat and mass
478 !
479 zsnowheatn(:,:)= zsnowhean(:,:)
480 WHERE(psnowdzn(:,:)==0.)
481  zsnowagen(:,:)= psnowage(:,:)
482  zsnowrhon(:,:)= psnowrho(:,:)
483  zsnowheatn(:,:)= psnowheat(:,:)
484 ELSEWHERE
485  zsnowagen(:,:)= zsnowagn(:,:)/zmastotn(:,:)
486  zsnowrhon(:,:)= zmastotn(:,:)/psnowdzn(:,:)
487 END WHERE
488 !
489 ! 4. Vanishing or very thin snowpack check:
490 ! -----------------------------------------
491 !
492 ! NOTE: ONLY for very shallow snowpacks, mix properties (homogeneous):
493 ! this avoids problems related to heat and mass exchange for
494 ! thin layers during heavy snowfall or signifigant melt: one
495 ! new/old layer can exceed the thickness of several old/new layers.
496 ! Therefore, mix (conservative):
497 !
498 zsumheat(:) = 0.0
499 zsumswe(:) = 0.0
500 zsumage(:) = 0.0
501 zsnowmix_delta(:) = 0.0
502 !
503 DO jl=1,inlvls
504  DO ji=1,ini
505  IF(psnow(ji) < xsnowcritd)THEN
506  zsumheat(ji) = zsumheat(ji) + psnowheat(ji,jl)
507  zsumswe(ji) = zsumswe(ji) + psnowrho(ji,jl)*psnowdz(ji,jl)
508  zsumage(ji) = zsumage(ji) + psnowage(ji,jl)
509  zsnowmix_delta(ji) = 1.0
510  ENDIF
511  ENDDO
512 ENDDO
513 !
514 ! Heat and mass are evenly distributed vertically:
515 ! heat and mass (density and thickness) are constant
516 ! in profile:
517 !
518 DO jl=1,inlvls
519  DO ji=1,ini
520 
521  IF(psnowdzn(ji,jl)>0.)THEN
522 
523 !
524  zsnowheatn(ji,jl) = zsnowmix_delta(ji)*(zsumheat(ji)/inlvls) + &
525  (1.0-zsnowmix_delta(ji))*zsnowheatn(ji,jl)
526 !
527  psnowdzn(ji,jl) = zsnowmix_delta(ji)*(psnow(ji)/inlvls) + &
528  (1.0-zsnowmix_delta(ji))*psnowdzn(ji,jl)
529 !
530  zsnowrhon(ji,jl) = zsnowmix_delta(ji)*(zsumswe(ji)/psnow(ji)) + &
531  (1.0-zsnowmix_delta(ji))*zsnowrhon(ji,jl)
532 !
533  zsnowagen(ji,jl) = zsnowmix_delta(ji)*(zsumage(ji)/inlvls) + &
534  (1.0-zsnowmix_delta(ji))*zsnowagen(ji,jl)
535 !
536  ENDIF
537  ENDDO
538 ENDDO
539 !
540 ! 5. Update mass (density and thickness) and heat:
541 ! ------------------------------------------------
542 !
543 psnowdz(:,:) = psnowdzn(:,:)
544 psnowrho(:,:) = zsnowrhon(:,:)
545 psnowheat(:,:) = zsnowheatn(:,:)
546 psnowage(:,:) = zsnowagen(:,:)
547 !
548 IF (lhook) CALL dr_hook('SNOW3LTRANSF',1,zhook_handle)
549 !
550 !
551 !-------------------------------------------------------------------------------
552 !
553 END SUBROUTINE snow3ltransf
554 !####################################################################
555  SUBROUTINE snow3lgrid(PSNOWDZ,PSNOW,PSNOWDZ_OLD)
556 !
557 !! PURPOSE
558 !! -------
559 ! Once during each time step, update grid to maintain
560 ! grid proportions. Similar to approach of Lynch-Steiglitz,
561 ! 1994, J. Clim., 7, 1842-1855. Corresponding mass and
562 ! heat adjustments made directly after the call to this
563 ! routine. 3 grid configurations:
564 ! 1) for very thin snow, constant grid spacing
565 ! 2) for intermediate thicknesses, highest resolution at soil/snow
566 ! interface and at the snow/atmosphere interface
567 ! 3) for deep snow, vertical resoution finest at snow/atmosphere
568 ! interface (set to a constant value) and increases with snow depth.
569 ! Second layer can't be more than an order of magnitude thicker
570 ! than surface layer.
571 !
572 !
573 USE modd_surf_par, ONLY : xundef
574 USE modd_snow_par, ONLY : xsnowcritd
575 !
576 USE yomhook ,ONLY : lhook, dr_hook
577 USE parkind1 ,ONLY : jprb
578 !
579 IMPLICIT NONE
580 !
581 !* 0.1 declarations of arguments
582 !
583 REAL, DIMENSION(: ), INTENT(IN ) :: psnow
584 REAL, DIMENSION(:,:), INTENT(OUT) :: psnowdz
585 REAL, DIMENSION(:,:), INTENT(IN ), OPTIONAL :: psnowdz_old
586 !
587 !* 0.1 declarations of local variables
588 !
589 INTEGER :: jj, ji
590 !
591 INTEGER :: inlvls, ini
592 !
593 REAL, DIMENSION(SIZE(PSNOW)) :: zwork
594 !
595 LOGICAL , DIMENSION(SIZE(PSNOW)) :: gregrid
596 
597 ! ISBA-ES snow grid parameters
598 !
599 REAL, PARAMETER, DIMENSION(3) :: zsgcoef1 = (/0.25, 0.50, 0.25/)
600 REAL, PARAMETER, DIMENSION(2) :: zsgcoef2 = (/0.05, 0.34/)
601 !
602 REAL, PARAMETER, DIMENSION(3) :: zsgcoef = (/0.3, 0.4, 0.3/)
603 !
604 ! Minimum total snow depth at which surface layer thickness is constant:
605 !
606 REAL, PARAMETER :: zsnowtrans = 0.20 ! (m)
607 !
608 ! Minimum snow depth by layer for 6-L or 12-L configuration :
609 !
610 REAL, PARAMETER :: zdz1=0.01
611 REAL, PARAMETER :: zdz2=0.05
612 REAL, PARAMETER :: zdz3=0.15
613 REAL, PARAMETER :: zdz4=0.50
614 REAL, PARAMETER :: zdz5=1.00
615 REAL, PARAMETER :: zdzn0=0.02
616 REAL, PARAMETER :: zdzn1=0.1
617 REAL, PARAMETER :: zdzn2=0.5
618 REAL, PARAMETER :: zdzn3=1.0
619 !
620 REAL, PARAMETER :: zcoef1 = 0.5
621 REAL, PARAMETER :: zcoef2 = 1.5
622 !
623 REAL(KIND=JPRB) :: zhook_handle
624 !
625 !-------------------------------------------------------------------------------
626 !
627 ! 0. Initialization:
628 ! ------------------
629 !
630 IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LGRID_2D',0,zhook_handle)
631 !
632 inlvls = SIZE(psnowdz(:,:),2)
633 ini = SIZE(psnowdz(:,:),1)
634 !
635 zwork(:) = 0.0
636 gregrid(:) = .true.
637 !
638 ! 1. Calculate current grid for 3-layer (default) configuration):
639 ! ---------------------------------------------------------------
640 ! Based on formulation of Lynch-Stieglitz (1994)
641 ! except for 3 modifications:
642 ! i) smooth transition here at ZSNOWTRANS
643 ! ii) constant ratio for very thin snow:
644 ! iii) ratio of layer 2 to surface layer <= 10
645 !
646 IF(inlvls == 1)THEN
647 !
648  DO ji=1,ini
649  psnowdz(ji,1) = psnow(ji)
650  ENDDO
651 !
652 ELSEIF(inlvls == 3)THEN
653 !
654  WHERE(psnow <= xsnowcritd+0.01)
655  psnowdz(:,1) = min(0.01, psnow(:)/inlvls)
656  psnowdz(:,3) = min(0.01, psnow(:)/inlvls)
657  psnowdz(:,2) = psnow(:) - psnowdz(:,1) - psnowdz(:,3)
658  END WHERE
659 !
660  WHERE(psnow <= zsnowtrans .AND. psnow > xsnowcritd+0.01)
661  psnowdz(:,1) = psnow(:)*zsgcoef1(1)
662  psnowdz(:,2) = psnow(:)*zsgcoef1(2)
663  psnowdz(:,3) = psnow(:)*zsgcoef1(3)
664  END WHERE
665 !
666  WHERE(psnow > zsnowtrans)
667  psnowdz(:,1) = zsgcoef2(1)
668  psnowdz(:,2) = (psnow(:)-zsgcoef2(1))*zsgcoef2(2) + zsgcoef2(1)
669 !
670 ! When using simple finite differences, limit the thickness
671 ! factor between the top and 2nd layers to at most 10
672 !
673  psnowdz(:,2) = min(10*zsgcoef2(1), psnowdz(:,2))
674  psnowdz(:,3) = psnow(:) - psnowdz(:,2) - psnowdz(:,1)
675  END WHERE
676 !
677 !
678 ! 2. Calculate current grid for 6-layer :
679 ! ---------------------------------------------------------------
680 !
681 ELSEIF(inlvls == 6)THEN
682 !
683 ! critere a satisfaire pour remaillage
684 !
685  IF(present(psnowdz_old))THEN
686  gregrid(:) = psnowdz_old(:,1) < zcoef1 * min(zdz1 ,psnow(:)/inlvls) .OR. &
687  & psnowdz_old(:,1) > zcoef2 * min(zdz1 ,psnow(:)/inlvls) .OR. &
688  & psnowdz_old(:,2) < zcoef1 * min(zdz2 ,psnow(:)/inlvls) .OR. &
689  & psnowdz_old(:,2) > zcoef2 * min(zdz2 ,psnow(:)/inlvls) .OR. &
690  & psnowdz_old(:,6) < zcoef1 * min(zdzn1,psnow(:)/inlvls) .OR. &
691  & psnowdz_old(:,6) > zcoef2 * min(zdzn1,psnow(:)/inlvls)
692  ENDIF
693 !
694  WHERE(gregrid(:))
695 ! top layers
696  psnowdz(:,1) = min(zdz1,psnow(:)/inlvls)
697  psnowdz(:,2) = min(zdz2,psnow(:)/inlvls)
698 ! last layers
699  psnowdz(:,6) = min(zdzn1,psnow(:)/inlvls)
700 ! remaining snow for remaining layers
701  zwork(:) = psnow(:) - psnowdz(:,1) - psnowdz(:,2) - psnowdz(:,6)
702  psnowdz(:,3) = zwork(:)*zsgcoef(1)
703  psnowdz(:,4) = zwork(:)*zsgcoef(2)
704  psnowdz(:,5) = zwork(:)*zsgcoef(3)
705 ! layer 3 tickness >= layer 2 tickness
706  zwork(:)=min(0.0,psnowdz(:,3)-psnowdz(:,2))
707  psnowdz(:,3)=psnowdz(:,3)-zwork(:)
708  psnowdz(:,4)=psnowdz(:,4)+zwork(:)
709 ! layer 5 tickness >= layer 6 tickness
710  zwork(:)=min(0.0,psnowdz(:,5)-psnowdz(:,6))
711  psnowdz(:,5)=psnowdz(:,5)-zwork(:)
712  psnowdz(:,4)=psnowdz(:,4)+zwork(:)
713  ENDWHERE
714 !
715 ! 3. Calculate current grid for 9-layer :
716 ! ---------------------------------------------------------------
717 !
718 ELSEIF(inlvls == 9)THEN
719 !
720 ! critere a satisfaire pour remaillage
721 !
722  IF(present(psnowdz_old))THEN
723  gregrid(:) = psnowdz_old(:,1) < zcoef1 * min(zdz1 ,psnow(:)/inlvls) .OR. &
724  & psnowdz_old(:,1) > zcoef2 * min(zdz1 ,psnow(:)/inlvls) .OR. &
725  & psnowdz_old(:,2) < zcoef1 * min(zdz2 ,psnow(:)/inlvls) .OR. &
726  & psnowdz_old(:,2) > zcoef2 * min(zdz2 ,psnow(:)/inlvls) .OR. &
727  & psnowdz_old(:,9) < zcoef1 * min(zdzn0,psnow(:)/inlvls) .OR. &
728  & psnowdz_old(:,9) > zcoef2 * min(zdzn0,psnow(:)/inlvls)
729  ENDIF
730 !
731  WHERE(gregrid(:))
732 ! top layers
733  psnowdz(:,1) = min(zdz1,psnow(:)/inlvls)
734  psnowdz(:,2) = min(zdz2,psnow(:)/inlvls)
735  psnowdz(:,3) = min(zdz3,psnow(:)/inlvls)
736 ! last layers
737  psnowdz(:,9)= min(zdzn0,psnow(:)/inlvls)
738  psnowdz(:,8)= min(zdzn1,psnow(:)/inlvls)
739  psnowdz(:,7)= min(zdzn2,psnow(:)/inlvls)
740 ! remaining snow for remaining layers
741  zwork(:) = psnow(:) - psnowdz(:, 1) - psnowdz(:, 2) - psnowdz(:, 3) &
742  - psnowdz(:, 7) - psnowdz(:, 8) - psnowdz(:, 9)
743  psnowdz(:,4) = zwork(:)*zsgcoef(1)
744  psnowdz(:,5) = zwork(:)*zsgcoef(2)
745  psnowdz(:,6) = zwork(:)*zsgcoef(3)
746 ! layer 4 tickness >= layer 3 tickness
747  zwork(:)=min(0.0,psnowdz(:,4)-psnowdz(:,3))
748  psnowdz(:,4)=psnowdz(:,4)-zwork(:)
749  psnowdz(:,5)=psnowdz(:,5)+zwork(:)
750 ! layer 6 tickness >= layer 7 tickness
751  zwork(:)=min(0.0,psnowdz(:,6)-psnowdz(:,7))
752  psnowdz(:,6)=psnowdz(:,6)-zwork(:)
753  psnowdz(:,5)=psnowdz(:,5)+zwork(:)
754  ENDWHERE
755 !
756 ! 4. Calculate current grid for 12-layer :
757 ! ---------------------------------------------------------------
758 !
759 ELSEIF(inlvls == 12)THEN
760 !
761 ! critere a satisfaire pour remaillage
762 !
763  IF(present(psnowdz_old))THEN
764  gregrid(:) = psnowdz_old(:, 1) < zcoef1 * min(zdz1 ,psnow(:)/inlvls) .OR. &
765  & psnowdz_old(:, 1) > zcoef2 * min(zdz1 ,psnow(:)/inlvls) .OR. &
766  & psnowdz_old(:, 2) < zcoef1 * min(zdz2 ,psnow(:)/inlvls) .OR. &
767  & psnowdz_old(:, 2) > zcoef2 * min(zdz2 ,psnow(:)/inlvls) .OR. &
768  & psnowdz_old(:,12) < zcoef1 * min(zdzn0,psnow(:)/inlvls) .OR. &
769  & psnowdz_old(:,12) > zcoef2 * min(zdzn0,psnow(:)/inlvls)
770  ENDIF
771 !
772  WHERE(gregrid(:))
773 ! top layers
774  psnowdz(:,1) = min(zdz1,psnow(:)/inlvls)
775  psnowdz(:,2) = min(zdz2,psnow(:)/inlvls)
776  psnowdz(:,3) = min(zdz3,psnow(:)/inlvls)
777  psnowdz(:,4) = min(zdz4,psnow(:)/inlvls)
778  psnowdz(:,5) = min(zdz5,psnow(:)/inlvls)
779 ! last layers
780  psnowdz(:,12)= min(zdzn0,psnow(:)/inlvls)
781  psnowdz(:,11)= min(zdzn1,psnow(:)/inlvls)
782  psnowdz(:,10)= min(zdzn2,psnow(:)/inlvls)
783  psnowdz(:, 9)= min(zdzn3,psnow(:)/inlvls)
784 ! remaining snow for remaining layers
785  zwork(:) = psnow(:) - psnowdz(:, 1) - psnowdz(:, 2) - psnowdz(:, 3) &
786  - psnowdz(:, 4) - psnowdz(:, 5) - psnowdz(:, 9) &
787  - psnowdz(:,10) - psnowdz(:,11) - psnowdz(:,12)
788  psnowdz(:,6) = zwork(:)*zsgcoef(1)
789  psnowdz(:,7) = zwork(:)*zsgcoef(2)
790  psnowdz(:,8) = zwork(:)*zsgcoef(3)
791 ! layer 6 tickness >= layer 5 tickness
792  zwork(:)=min(0.0,psnowdz(:,6)-psnowdz(:,5))
793  psnowdz(:,6)=psnowdz(:,6)-zwork(:)
794  psnowdz(:,7)=psnowdz(:,7)+zwork(:)
795 ! layer 8 tickness >= layer 9 tickness
796  zwork(:)=min(0.0,psnowdz(:,8)-psnowdz(:,9))
797  psnowdz(:,8)=psnowdz(:,8)-zwork(:)
798  psnowdz(:,7)=psnowdz(:,7)+zwork(:)
799  ENDWHERE
800 !
801 ! 4. Calculate other non-optimized grid :
802 ! ---------------------------------------
803 !
804 ELSEIF(inlvls<10.AND.inlvls/=3.AND.inlvls/=6.AND.inlvls/=9) THEN
805 !
806  DO jj=1,inlvls
807  DO ji=1,ini
808  psnowdz(ji,jj) = psnow(ji)/inlvls
809  ENDDO
810  ENDDO
811 !
812  psnowdz(:,inlvls) = psnowdz(:,inlvls) + (psnowdz(:,1) - min(0.05, psnowdz(:,1)))
813  psnowdz(:,1) = min(0.05, psnowdz(:,1))
814 !
815 ELSE !(INLVLS>=10 and /=12)
816 !
817 ! critere a satisfaire pour remaillage
818 !
819  IF(present(psnowdz_old))THEN
820  gregrid(:) = psnowdz_old(:, 1) < zcoef1 * min(zdz1 ,psnow(:)/inlvls) .OR. &
821  & psnowdz_old(:, 1) > zcoef2 * min(zdz1 ,psnow(:)/inlvls) .OR. &
822  & psnowdz_old(:, 2) < zcoef1 * min(zdz2 ,psnow(:)/inlvls) .OR. &
823  & psnowdz_old(:, 2) > zcoef2 * min(zdz2 ,psnow(:)/inlvls) .OR. &
824  & psnowdz_old(:,inlvls) < zcoef1 * min(0.05*psnow(:),psnow(:)/inlvls) .OR. &
825  & psnowdz_old(:,inlvls) > zcoef2 * min(0.05*psnow(:),psnow(:)/inlvls)
826  ENDIF
827 !
828  WHERE(gregrid(:))
829  psnowdz(:,1 ) = min(zdz1 ,psnow(:)/inlvls)
830  psnowdz(:,2 ) = min(zdz2 ,psnow(:)/inlvls)
831  psnowdz(:,3 ) = min(zdz3 ,psnow(:)/inlvls)
832  psnowdz(:,4 ) = min(zdz4 ,psnow(:)/inlvls)
833  psnowdz(:,5 ) = min(zdz5 ,psnow(:)/inlvls)
834  psnowdz(:,inlvls) = min(0.05*psnow(:),psnow(:)/inlvls)
835  ENDWHERE
836 !
837  DO jj=6,inlvls-1,1
838  DO ji=1,ini
839  IF(gregrid(ji))THEN
840  zwork(ji) = psnowdz(ji,1)+psnowdz(ji,2)+psnowdz(ji,3)+psnowdz(ji,4)+psnowdz(ji,5)
841  psnowdz(ji,jj) = (psnow(ji)-zwork(ji)-psnowdz(ji,inlvls))/(inlvls-6)
842  ENDIF
843  ENDDO
844  ENDDO
845 !
846 ENDIF
847 !
848 DO jj=1,inlvls
849  DO ji=1,ini
850  IF(psnow(ji)==xundef)THEN
851  psnowdz(ji,jj) = xundef
852  ELSEIF(.NOT.gregrid(ji))THEN
853  psnowdz(ji,jj)=psnowdz_old(ji,jj)
854  ENDIF
855  ENDDO
856 ENDDO
857 !
858 IF (lhook) CALL dr_hook('MODE_SNOW3L:SNOW3LGRID_2D',1,zhook_handle)
859 !
860 END SUBROUTINE snow3lgrid
861 !####################################################################
862 SUBROUTINE snow3lcompactn(PTSTEP,PSNOWDZMIN,PSNOWRHO,PSNOWDZ,PSNOWTEMP,PSNOW,PSNOWLIQ)
863 !
864 !! PURPOSE
865 !! -------
866 ! Snow compaction due to overburden and settling.
867 ! Mass is unchanged: layer thickness is reduced
868 ! in proportion to density increases. Method
869 ! of Brun et al (1989) and Vionnet et al. (2012)
870 !
871 !
872 USE modd_surf_par, ONLY : xundef
873 !
874 USE modd_csts, ONLY : xtt, xg
875 USE modd_snow_par, ONLY : xrhosmax_es
876 !
877 USE modd_snow_metamo, ONLY : xvvisc1,xvvisc3,xvvisc4, &
878  xvvisc5,xvvisc6,xvro11
879 !
880 USE mode_snow3l, ONLY : snow3lhold
881 !
882 IMPLICIT NONE
883 !
884 !* 0.1 declarations of arguments
885 !
886 REAL, INTENT(IN) :: ptstep
887 REAL, INTENT(IN) :: psnowdzmin
888 !
889 REAL, DIMENSION(:,:), INTENT(IN) :: psnowtemp, psnowliq
890 !
891 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowrho, psnowdz
892 !
893 REAL, DIMENSION(:), INTENT(OUT) :: psnow
894 !
895 !
896 !* 0.2 declarations of local variables
897 !
898 INTEGER :: jj, ji
899 !
900 INTEGER :: ini
901 INTEGER :: inlvls
902 !
903 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowrho2, zviscocity, zf1, &
904  ztemp, zsmass, zsnowdz, &
905  zwsnowdz, zwholdmax
906 !
907 !
908 REAL(KIND=JPRB) :: zhook_handle
909 !
910 !-------------------------------------------------------------------------------
911 !
912 ! 0. Initialization:
913 ! ------------------
914 !
915 IF (lhook) CALL dr_hook('SNOW3LCOMPACTN',0,zhook_handle)
916 !
917 ini = SIZE(psnowdz(:,:),1)
918 inlvls = SIZE(psnowdz(:,:),2)
919 !
920 zsnowrho2(:,:) = psnowrho(:,:)
921 zsnowdz(:,:) = max(psnowdzmin,psnowdz(:,:))
922 zviscocity(:,:) = 0.0
923 ztemp(:,:) = 0.0
924 !
925 ! 1. Cumulative snow mass (kg/m2):
926 ! --------------------------------
927 !
928 zsmass(:,:) = 0.0
929 DO jj=2,inlvls
930  DO ji=1,ini
931  zsmass(ji,jj) = zsmass(ji,jj-1) + psnowdz(ji,jj-1)*psnowrho(ji,jj-1)
932  ENDDO
933 ENDDO
934 ! overburden of half the mass of the uppermost layer applied to itself
935 zsmass(:,1) = 0.5 * psnowdz(:,1) * psnowrho(:,1)
936 !
937 ! 2. Compaction
938 ! -------------
939 !
940 !Liquid water effect
941 !
942 zwholdmax(:,:) = snow3lhold(psnowrho,psnowdz)
943 zwholdmax(:,:) = max(1.e-10, zwholdmax(:,:))
944 zf1(:,:) = 1.0/(xvvisc5+10.*min(1.0,psnowliq(:,:)/zwholdmax(:,:)))
945 !
946 !Snow viscocity, density and grid thicknesses
947 !
948 DO jj=1,inlvls
949  DO ji=1,ini
950 !
951  IF(psnowrho(ji,jj) < xrhosmax_es)THEN
952 !
953 ! temperature dependence limited to 10K: Schleef et al. (2014)
954  ztemp(ji,jj) = xvvisc4*min(10.,abs(xtt-psnowtemp(ji,jj)))
955 !
956 ! Calculate snow viscocity: Brun et al. (1989), Vionnet et al. (2012)
957  zviscocity(ji,jj) = xvvisc1*zf1(ji,jj)*exp(xvvisc3*psnowrho(ji,jj)+ztemp(ji,jj))*psnowrho(ji,jj)/xvro11
958 !
959 ! Calculate snow density:
960  zsnowrho2(ji,jj) = psnowrho(ji,jj) + psnowrho(ji,jj)*ptstep &
961  * ( (xg*zsmass(ji,jj)/zviscocity(ji,jj)) )
962 !
963 ! Conserve mass by decreasing grid thicknesses in response to density increases
964  psnowdz(ji,jj) = psnowdz(ji,jj)*(psnowrho(ji,jj)/zsnowrho2(ji,jj))
965 !
966  ENDIF
967 !
968  ENDDO
969 ENDDO
970 !
971 ! 3. Update total snow depth and density profile:
972 ! -----------------------------------------------
973 !
974 ! Compaction/augmentation of total snowpack depth
975 !
976 psnow(:) = 0.
977 DO jj=1,inlvls
978  DO ji=1,ini
979  psnow(ji) = psnow(ji) + psnowdz(ji,jj)
980  ENDDO
981 ENDDO
982 !
983 ! Update density (kg m-3):
984 !
985 psnowrho(:,:) = zsnowrho2(:,:)
986 !
987 IF (lhook) CALL dr_hook('SNOW3LCOMPACTN',1,zhook_handle)
988 !
989 !-------------------------------------------------------------------------------
990 !
991 END SUBROUTINE snow3lcompactn
992 !####################################################################
993 
994 
995 END SUBROUTINE preps_for_meb_ebud_rad
subroutine snow3lcompactn(PTSTEP, PSNOWDZMIN, PSNOWRHO, PSNOWDZ, PSNOWTEMP, PSNOW, PSNOWLIQ)
subroutine preps_for_meb_ebud_rad(PPS, PLAICV, PSNOWRHO, PSNOWSWE, PSNOWHEAT, PSNOWTEMP, PSNOWDZ, PSCOND, PHEATCAPS, PEMISNOW, PSIGMA_F, PCHIP, PTSTEP, PSR, PTA, PVMOD, PSNOWAGE, PPERMSNOWFRAC)
subroutine snow3lfall(PTSTEP, PSR, PTA, PVMOD, PSNOW, PSNOWRHO, PSNOWDZ, PSNOWHEAT, PSNOWHMASS, PSNOWAGE, PPERMSNOWFRAC)
subroutine snow3ltransf(PSNOW, PSNOWDZ, PSNOWDZN, PSNOWRHO, PSNOWHEAT, PSNOWAGE)