SURFEX  V8_0
Surfex V8_0 release
 All Classes Files Functions Variables
snow3l.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 snow3l(HSNOWRES, TPTIME, OMEB, HIMPLICIT_WIND, &
7  ppew_a_coef, ppew_b_coef, &
8  ppet_a_coef, ppeq_a_coef, ppet_b_coef, ppeq_b_coef, &
9  psnowswe,psnowrho,psnowheat,psnowalb, &
10  psnowgran1,psnowgran2,psnowhist,psnowage, &
11  ptstep,pps,psr,prr,ppsn3l, &
12  pta,ptg,psw_rad,pqa,pvmod,plw_rad, prhoa, &
13  puref,pexns,pexna,pdircoszw, &
14  pzref,pz0,pz0eff,pz0h,palb, &
15  psoilcond,pd_g,plvtt,plstt, &
16  psnowliq,psnowtemp,psnowdz, &
17  pthrufal,pgrndflux,pevapcor,psoilcor, &
18  pgflxcor,psnowsfch, pdelheatn, pdelheatn_sfc, &
19  pswnetsnow,pswnetsnows,plwnetsnow,psnowflux, &
20  prnsnow,phsnow,pgfluxsnow, &
21  phpsnow,ples3l,plel3l,pevap,psndrift,pri, &
22  pemisnow,pcdsnow,pustar,pchsnow,psnowhmass,pqs, &
23  ppermsnowfrac,pforestfrac,pzenith,pxlat,pxlon, &
24  osnowdrift,osnowdrift_sublim )
25 ! ##########################################################################
26 !
27 !!**** *SNOW3L*
28 !!
29 !! PURPOSE
30 !! -------
31 !
32 ! 3(N)-Layer snow scheme option (Boone and Etchevers, J Hydrometeor., 2000)
33 !
34 !!** METHOD
35 !! ------
36 !
37 ! Direct calculation
38 !
39 !! EXTERNAL
40 !! --------
41 !
42 ! None
43 !!
44 !! IMPLICIT ARGUMENTS
45 !! ------------------
46 !!
47 !!
48 !!
49 !! REFERENCE
50 !! ---------
51 !!
52 !! ISBA-ES: Boone and Etchevers (2001)
53 !! ISBA: Belair (1995)
54 !! ISBA: Noilhan and Planton (1989)
55 !! ISBA: Noilhan and Mahfouf (1996)
56 !!
57 !! AUTHOR
58 !! ------
59 !! A. Boone * Meteo-France *
60 !!
61 !! MODIFICATIONS
62 !! -------------
63 !! Original 7/99
64 !! Modified by A.Boone 05/02 (code, not physics)
65 !! Modified by A.Boone 11/04 i) maximum density limit imposed (although
66 !! rarely if ever reached), ii) check to
67 !! see if upermost layer completely sublimates
68 !! during a timestep (as snowpack becomes vanishly
69 !! thin), iii) impose maximum grain size limit
70 !! in radiation transmission computation.
71 !!
72 !! Modified by B. Decharme (03/2009): Consistency with Arpege permanent
73 !! snow/ice treatment (LGLACIER for alb)
74 !! Modified by A. Boone (04/2010): Implicit coupling and replace Qsat and DQsat
75 !! by Qsati and DQsati, respectively.
76 !! Modified by E. Brun (08/2012): Mass conservation in SNOW3LEVAPGONE
77 !! Modified by B. Decharme (08/2012): Loop optimization
78 !! Modified by B. Decharme (09/2012): New wind implicitation
79 !! Modified by E. Brun (10/2012): Bug in vertical snow redistribution
80 !! Modified by B. Decharme (08/2013): Qsat as argument (needed for coupling with atm)
81 !! Add snow drift effect like in CROCUS
82 !! Change snow albedo (CROCUS spectral albedo)
83 !! Change snow viscosity like in CROCUS
84 !! Modified by P. Samuelsson(10/2014): MEB added an optional snow albedo calculation
85 !! for litter
86 !! Modified by A. Boone (10/2014): MEB modifs: removed above since spectral albedo
87 !! method from CROCUS added...can eventually add above
88 !! back in a manner consistent with spectral bands...
89 !! Modified by A. Boone (10/2014): MEB modifs: permit option to impose fluxes at sfc
90 !! Modified by A. Boone (10/2014): SNOW3LREFRZ and SNOW3LEVAPN edited to give
91 !! better enthalpy conservation.
92 !! Modified by B. Decharme (03/2016): No snowdrift under forest
93 !!
94 !!
95 !-------------------------------------------------------------------------------
96 !
97 !* 0. DECLARATIONS
98 ! ------------
99 !
100 USE modd_csts, ONLY : xtt, xrholw, xlmtt, xcl, xday
101 USE modd_surf_par, ONLY : xundef
102 !
103 USE modd_snow_metamo, ONLY : xsnowdzmin
104 USE modd_snow_par, ONLY : xsnowdmin, nspec_band_snow
105 !
106 USE mode_snow3l
107 USE mode_thermos, ONLY : qsati
108 !
109 USE modd_type_date_surf, ONLY: date_time
110 !
111 USE yomhook ,ONLY : lhook, dr_hook
112 USE parkind1 ,ONLY : jprb
113 !
114 IMPLICIT NONE
115 !
116 !
117 !* 0.1 declarations of arguments
118 !
119 REAL, INTENT(IN) :: ptstep
120 ! PTSTEP = time step of the integration
121 TYPE(date_time), INTENT(IN) :: tptime ! current date and time
122 !
123  CHARACTER(LEN=*), INTENT(IN) :: hsnowres
124 ! HSNOWRES = ISBA-SNOW3L turbulant exchange option
125 ! 'DEF' = Default: Louis (ISBA: Noilhan and Mahfouf 1996)
126 ! 'RIL' = Limit Richarson number under very stable
127 ! conditions (currently testing)
128 !
129 LOGICAL, INTENT(IN) :: omeb ! True = coupled to MEB. This means surface fluxes ae IMPOSED
130 ! ! as an upper boundary condition to the explicit snow schemes.
131 ! ! If = False, then energy
132 ! ! budget and fluxes are computed herein.
133 !
134  CHARACTER(LEN=*), INTENT(IN) :: himplicit_wind ! wind implicitation option
135 ! ! 'OLD' = direct
136 ! ! 'NEW' = Taylor serie, order 1
137 !
138 REAL, DIMENSION(:), INTENT(IN) :: pps, pta, psw_rad, pqa, &
139  pvmod, plw_rad, psr, prr
140 ! PSW_RAD = incoming solar radiation (W/m2)
141 ! PLW_RAD = atmospheric infrared radiation (W/m2)
142 ! PRR = rain rate [kg/(m2 s)]
143 ! PSR = snow rate (SWE) [kg/(m2 s)]
144 ! PTA = atmospheric temperature at level za (K)
145 ! PVMOD = modulus of the wind parallel to the orography (m/s)
146 ! PPS = surface pressure
147 ! PQA = atmospheric specific humidity
148 ! at level za
149 !
150 REAL, DIMENSION(:), INTENT(IN) :: psoilcond, pd_g, ppsn3l
151 ! PSOILCOND = soil thermal conductivity [W/(m K)]
152 ! PD_G = Assumed first soil layer thickness (m)
153 ! Used to calculate ground/snow heat flux
154 ! PPSN3L = snow fraction
155 !
156 REAL, DIMENSION(:), INTENT(IN) :: pzref, puref, pexns, pexna, pdircoszw, prhoa, pz0, pz0eff, &
157  palb, pz0h, ppermsnowfrac, pforestfrac
158 ! PZ0EFF = roughness length for momentum
159 ! PZ0 = grid box average roughness length
160 ! PZ0H = grid box average roughness length for heat
161 ! PZREF = reference height of the first
162 ! atmospheric level
163 ! PUREF = reference height of the wind
164 ! PRHOA = air density
165 ! PEXNS = Exner function at surface
166 ! PEXNA = Exner function at lowest atmos level
167 ! PDIRCOSZW = Cosinus of the angle between the
168 ! normal to the surface and the vertical
169 ! PALB = soil/vegetation albedo
170 ! PPERMSNOWFRAC = fraction of permanet snow/ice
171 ! PFORESTFRAC = fraction of forest
172 !
173 REAL, DIMENSION(:), INTENT(IN) :: ppew_a_coef, ppew_b_coef, &
174  ppet_a_coef, ppeq_a_coef, ppet_b_coef, &
175  ppeq_b_coef
176 ! PPEW_A_COEF = wind coefficient (m2s/kg)
177 ! PPEW_B_COEF = wind coefficient (m/s)
178 ! PPET_A_COEF = A-air temperature coefficient
179 ! PPET_B_COEF = B-air temperature coefficient
180 ! PPEQ_A_COEF = A-air specific humidity coefficient
181 ! PPEQ_B_COEF = B-air specific humidity coefficient
182 !
183 REAL, DIMENSION(:), INTENT(IN) :: ptg
184 ! PTG = Surface soil temperature (effective
185 ! temperature the of layer lying below snow)
186 REAL, DIMENSION(:), INTENT(IN) :: plvtt, plstt ! = latent heats for hydrology
187 !
188 REAL, DIMENSION(:), INTENT(INOUT) :: psnowalb
189 ! PSNOWALB = Prognostic surface snow albedo
190 ! (does not include anything but
191 ! the actual snow cover)
192 !
193 REAL, DIMENSION(:,:), INTENT(INOUT):: psnowheat, psnowrho, psnowswe
194 ! PSNOWHEAT = Snow layer(s) heat content (J/m2)
195 ! PSNOWRHO = Snow layer(s) averaged density (kg/m3)
196 ! PSNOWSWE = Snow layer(s) liquid Water Equivalent (SWE:kg m-2)
197 !
198 REAL, DIMENSION(:,:), INTENT(INOUT):: psnowgran1, psnowgran2, psnowhist
199 ! PSNOWGRAN1 = Snow layers grain feature 1
200 ! PSNOWGRAN2 = Snow layer grain feature 2
201 ! PSNOWHIST = Snow layer grain historical
202 ! parameter (only for non
203 ! dendritic snow)
204 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowage ! Snow grain age
205 !
206 REAL, DIMENSION(:), INTENT(INOUT) :: prnsnow, phsnow, ples3l, plel3l, &
207  phpsnow, pevap, pgrndflux, pemisnow
208 ! PLES3L = evaporation heat flux from snow (W/m2)
209 ! PLEL3L = sublimation (W/m2)
210 ! PHPSNOW = heat release from rainfall (W/m2)
211 ! PRNSNOW = net radiative flux from snow (W/m2)
212 ! PHSNOW = sensible heat flux from snow (W/m2)
213 ! PEVAP = total evaporative flux (kg/m2/s)
214 ! PGRNDFLUX = soil/snow interface heat flux (W/m2)
215 ! PEMISNOW = snow surface emissivity
216 !
217 REAL, DIMENSION(:), INTENT(OUT) :: pgfluxsnow
218 ! PGFLUXSNOW = net heat flux from snow (W/m2)
219 !
220 REAL, DIMENSION(:), INTENT(INOUT) :: pswnetsnow, plwnetsnow, pswnetsnows
221 ! PSWNETSNOW = net shortwave radiation entering top of snowpack
222 ! (W m-2) Imposed if MEB=T, diagnosed herein if MEB=F
223 ! PSWNETSNOWS= net shortwave radiation in uppermost layer of snowpack
224 ! (W m-2) Imposed if MEB=T, diagnosed herein if MEB=F
225 ! Used for surface energy budget diagnostics
226 ! PLWNETSNOW = net longwave radiation entering top of snowpack
227 ! (W m-2) Imposed if MEB=T, diagnosed herein if MEB=F
228 !
229 REAL, DIMENSION(:), INTENT(INOUT) :: pustar, pcdsnow, pchsnow, pri
230 ! PCDSNOW = drag coefficient for momentum over snow (-)
231 ! PUSTAR = friction velocity over snow (m/s)
232 ! PCHSNOW = drag coefficient for heat over snow (-)
233 ! PRI = Richardson number (-)
234 !
235 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowtemp
236 REAL, DIMENSION(:,:), INTENT(OUT) :: psnowliq, psnowdz
237 ! PSNOWLIQ = Snow layer(s) liquid water content (m)
238 ! PSNOWTEMP = Snow layer(s) temperature (m)
239 ! PSNOWDZ = Snow layer(s) thickness (m)
240 !
241 REAL, DIMENSION(:), INTENT(OUT) :: pthrufal, pevapcor, psoilcor, pgflxcor, &
242  psnowflux, psnowsfch, pdelheatn, pdelheatn_sfc
243 ! PTHRUFAL = rate that liquid water leaves snow pack:
244 ! paritioned into soil infiltration/runoff
245 ! by ISBA [kg/(m2 s)]
246 ! PEVAPCOR = evaporation/sublimation correction term:
247 ! extract any evaporation exceeding the
248 ! actual snow cover (as snow vanishes)
249 ! and apply it as a surface soil water
250 ! sink. [kg/(m2 s)]
251 ! PSOILCOR = for vanishingy thin snow cover,
252 ! allow any excess evaporation
253 ! to be extracted from the soil
254 ! to maintain an accurate water
255 ! balance [kg/(m2 s)]
256 ! PGFLXCOR = flux correction to underlying soil for vanishing snowpack
257 ! (to put any energy excess from snow to soil) (W/m2)
258 ! PSNOWFLUX = heat flux between the surface and sub-surface
259 ! snow layers (W/m2)
260 ! PSNOWSFCH = snow surface layer pseudo-heating term owing to
261 ! changes in grid thickness (W m-2)
262 ! PDELHEATN = total snow heat content change in the surface layer (W m-2)
263 ! PDELHEATN_SFC = total snow heat content change during the timestep (W m-2)
264 !
265 REAL, DIMENSION(:), INTENT(OUT) :: psndrift
266 ! PSNDRIFT = blowing snow sublimation (kg/m2/s)
267 !
268 REAL, DIMENSION(:), INTENT(OUT) :: psnowhmass
269 ! PSNOWHMASS = heat content change due to mass
270 ! changes in snowpack (J/m2): for budget
271 ! calculations only.
272 !
273 REAL, DIMENSION(:), INTENT(OUT) :: pqs
274 ! PQS = surface humidity
275 !
276 REAL, DIMENSION(:), INTENT(IN) :: pzenith ! solar zenith angle
277 REAL, DIMENSION(:), INTENT(IN) :: pxlat,pxlon ! LAT/LON after packing
278 !
279 LOGICAL, INTENT(IN) :: osnowdrift, osnowdrift_sublim ! activate snowdrift, sublimation during drift
280 !
281 !* 0.2 declarations of local variables
282 !
283 INTEGER :: jj, ji ! Loop control
284 !
285 INTEGER :: ini ! number of point
286 INTEGER :: inlvls ! number of snow layers
287 !
288 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowtemp, zscap, zsnowdzn, zscond, &
289  zradsink, zwork2d, zsnowtempo
290 ! ZSNOWTEMP = Snow layer(s) averaged temperature (K)
291 ! ZSCAP = Snow layer(s) heat capacity [J/(K m3)]
292 ! ZSNOWDZN = Updated snow layer thicknesses (m)
293 ! ZSCOND = Snow layer(s) thermal conducivity [W/(m K)]
294 ! ZRADSINK = Snow solar Radiation source terms (W/m2)
295 ! ZWORK2D = working variable (*)
296 !
297 REAL, DIMENSION(SIZE(PTA)) :: zsnow, zsfcfrz, ztsterm1, ztsterm2, &
298  zct, zra, zsnowtempo1
299 ! ZSNOW = Total snow depth (m)
300 ! ZCT = inverse of the product of snow heat capacity
301 ! and layer thickness [(m2 K)/J]
302 ! ZRA = Surface aerodynamic resistance
303 ! ZTSTERM1,ZTSTERM2 = Surface energy budget coefficients
304 ! ZSNOWTEMPO1= value of uppermost snow temperature
305 ! before time integration (K)
306 !
307 LOGICAL, DIMENSION(SIZE(PTA)) :: gsfcmelt
308 ! GSFCMELT = FLAG if surface melt is occurring, used
309 ! for surface albedo calculation.
310 !
311 REAL, DIMENSION(SIZE(PTA)) :: zrsra, zdqsat, zqsat, zradxs, zmeltxs, zliqheatxs, &
312  zlwupsnow, zgrndflux, zgrndfluxo, zgrndfluxi, zpsn3l
313 ! ZRSRA = air density over aerodynamic resistance
314 ! ZDQSAT = derrivative of saturation specific humidity
315 ! ZQSAT = saturation specific humidity
316 ! ZRADXS = shortwave radiation absorbed by soil surface
317 ! (for thin snow sover) (W m-2)
318 ! ZMELTXS = excess energy for snowmelt for vanishingly
319 ! thin snowpacks: used to heat underlying surface
320 ! in order to maintain energy conservation (W m-2)
321 ! ZLIQHEATXS = excess snowpack heating for vanishingly thin
322 ! snow cover: add energy to snow/ground heat
323 ! flux (W m-2)
324 ! ZLWUPSNOW = upwelling longwave radiative flux (W m-2)
325 ! ZGRNDFLUXO= snow-ground flux before correction (W m-2)
326 ! ZGRNDFLUX = snow-ground flux after correction (W m-2)
327 ! ZGRNDFLUXI= for the case where the ground flux is imposed,
328 ! this is the actual imposed value.
329 ! ZPSN3L = snow fraction: different use if MEB "on".
330 ! In this case, it is only used for Tg update
331 ! since only this variable has a sub-grid relevance.
332 !
333 REAL, DIMENSION(SIZE(PTA)) :: zustar2_ic, zta_ic, zqa_ic, zwork, zwork2, zwork3, &
334  zpet_a_coef_t, zpeq_a_coef_t, zpet_b_coef_t, zpeq_b_coef_t
335 ! ZUSTAR2_IC = implicit lowest atmospheric level friction (m2/s2)
336 ! ZTA_IC = implicit lowest atmospheric level air temperature
337 ! ZQA_IC = implicit lowest atmospheric level specific humidity
338 ! ZPET_A_COEF_T = transformed A-air temperature coefficient
339 ! ZPET_B_COEF_T = transformed B-air temperature coefficient
340 ! ZPEQ_A_COEF_T = transformed A-air specific humidity coefficient
341 ! ZPEQ_B_COEF_T = transformed B-air specific humidity coefficient
342 !
343 REAL, DIMENSION(SIZE(PSNOWRHO,1),NSPEC_BAND_SNOW) :: zspectralalbedo, zspectralwork
344 ! ZSPECTRALALBEDO=spectral albedo
345 !
346 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowheat0
347 !
348 REAL(KIND=JPRB) :: zhook_handle
349 !
350 ! - - ---------------------------------------------------
351 !
352 ! 0. Initialization
353 ! --------------
354 ! NOTE that snow layer thickness is used throughout this code: SWE
355 ! is only used to diagnose the thickness at the beginning of this routine
356 ! and it is updated at the end of this routine.
357 !
358 IF (lhook) CALL dr_hook('SNOW3L',0,zhook_handle)
359 !
360 psnowdz(:,:) = psnowswe(:,:)/psnowrho(:,:)
361 !
362 ini = SIZE(psnowswe(:,:),1)
363 inlvls = SIZE(psnowswe(:,:),2) ! total snow layers
364 !
365 zustar2_ic = 0.0
366 zta_ic = 0.0
367 zqa_ic = 0.0
368 !
369 zgrndfluxo = 0.0
370 zgrndflux = pgrndflux
371 !
372 zsnowheat0(:,:) = psnowheat(:,:) ! save initial heat content
373 !
374 zsnowtempo(:,:) = psnowtemp(:,:) ! for MEB, this is the updated T profile
375 !
376 !* 1. Snow total depth
377 ! ----------------
378 !
379 zsnow(:) = 0.
380 DO jj=1,inlvls
381  DO ji=1,ini
382  zsnow(ji) = zsnow(ji) + psnowdz(ji,jj) ! m
383  ENDDO
384 END DO
385 !
386 zwork(:)=zsnow(:)
387 !
388 ! Caluclate new snow albedo at time t
389 !
390 zwork2(:)=psnowalb(:)
391 !
392  CALL snow3lalb(zwork2,zspectralalbedo,psnowrho(:,1),psnowage(:,1),ppermsnowfrac,pps)
393 zwork3(:) = psnowalb(:)/zwork2(:)
394 DO jj=1,SIZE(zspectralalbedo,2)
395  DO ji=1,ini
396  zspectralalbedo(ji,jj)=zspectralalbedo(ji,jj)*zwork3(ji)
397  ENDDO
398 ENDDO
399 !
400 !* 2. Snowfall
401 ! --------
402 !
403 ! Caluclate uppermost density and thickness changes due to snowfall,
404 ! and add heat content of falling snow
405 !
406  CALL snow3lfall(ptstep,psr,pta,pvmod,zsnow,psnowrho,psnowdz, &
407  psnowheat,psnowhmass,psnowage,ppermsnowfrac )
408 !
409 ! Caluclate new snow albedo at time t if snowfall
410 !
411  CALL snow3lalb(zwork2,zspectralwork,psnowrho(:,1),psnowage(:,1),ppermsnowfrac,pps)
412 !
413 DO jj=1,SIZE(zspectralalbedo,2)
414  DO ji=1,ini
415  IF(zwork(ji)==0.0.AND.psr(ji)>0.0)THEN
416  zspectralalbedo(ji,jj) = zspectralwork(ji,jj)
417  ENDIF
418  ENDDO
419 ENDDO
420 WHERE(zwork(:)==0.0.AND.psr(:)>0.0)
421  psnowalb(:) = zwork2(:)
422 ENDWHERE
423 !
424 !* 3. Update grid/discretization
425 ! --------------------------
426 ! Reset grid to conform to model specifications:
427 !
428  CALL snow3lgrid(zsnowdzn,zsnow,psnowdz_old=psnowdz)
429 !
430 ! Mass/Heat redistribution:
431 !
432  CALL snow3ltransf(zsnow,psnowdz,zsnowdzn,psnowrho,psnowheat,psnowage)
433 !
434 !
435 !* 4. Liquid water content and snow temperature
436 ! -----------------------------------------
437 !
438 ! First diagnose snow temperatures and liquid
439 ! water portion of the snow from snow heat content:
440 !
441 zscap(:,:) = snow3lscap(psnowrho)
442 !
443 zsnowtemp(:,:) = xtt + ( ((psnowheat(:,:)/psnowdz(:,:)) &
444  + xlmtt*psnowrho(:,:))/zscap(:,:) )
445 !
446 psnowliq(:,:) = max(0.0,zsnowtemp(:,:)-xtt)*zscap(:,:)* &
447  psnowdz(:,:)/(xlmtt*xrholw)
448 !
449 zsnowtemp(:,:) = min(xtt,zsnowtemp(:,:))
450 !
451 !
452 !* 5. Snow Compaction
453 ! ---------------
454 !
455 ! Calculate snow density: compaction/aging: density increases
456 !
457  CALL snow3lcompactn(ptstep,xsnowdzmin,psnowrho,psnowdz,zsnowtemp,zsnow,psnowliq)
458 !
459 ! Snow compaction and metamorphism due to drift
460 !
461 psndrift(:) = 0.0
462 IF (osnowdrift) THEN
463  CALL snow3ldrift(ptstep,pforestfrac,pvmod,pta,pqa,pps,prhoa,&
464  psnowrho,psnowdz,zsnow,osnowdrift_sublim,psndrift)
465 ENDIF
466 !
467 ! Update snow heat content (J/m2):
468 !
469 zscap(:,:) = snow3lscap(psnowrho)
470 psnowheat(:,:) = psnowdz(:,:)*( zscap(:,:)*(zsnowtemp(:,:)-xtt) &
471  - xlmtt*psnowrho(:,:) ) + xlmtt*xrholw*psnowliq(:,:)
472 !
473 !
474 !* 6. Solar radiation transmission
475 ! -----------------------------
476 !
477 ! Heat source (-sink) term due to shortwave
478 ! radiation transmission within the snowpack:
479 !
480  CALL snow3lrad(omeb,xsnowdzmin,psw_rad,psnowalb, &
481  zspectralalbedo,psnowdz,psnowrho,palb, &
482  ppermsnowfrac,pzenith,pswnetsnow, &
483  pswnetsnows,zradsink,zradxs,psnowage)
484 !
485 !
486 !* 7. Heat transfer and surface energy budget
487 ! ---------------------------------------
488 ! Snow thermal conductivity:
489 !
490  CALL snow3lthrm(psnowrho,zscond,zsnowtemp,pps)
491 !
492 ! Precipitation heating term:
493 ! Rainfall renders it's heat to the snow when it enters
494 ! the snowpack:
495 ! NOTE: for now set to zero because we'd need to remove this heat from
496 ! the atmosphere to conserve energy. This can be done, but should
497 ! be tested within an atmos model first probably...NOTE, this
498 ! could be actiavted for use in OFFLINE mode however.
499 !
500 !PHPSNOW(:) = PRR(:)*XCL*(MAX(XTT,PTA(:))-XTT) ! (W/m2)
501 phpsnow(:) = 0.0
502 !
503 ! Surface Energy Budget calculations using ISBA linearized form
504 ! and standard ISBA turbulent transfer formulation
505 !
506 IF(omeb)THEN
507 
508 ! - surface fluxes prescribed:
509 
510  CALL snow3lebudmeb(ptstep,xsnowdzmin, &
511  zsnowtemp(:,1),psnowdz(:,1),psnowdz(:,2), &
512  zscond(:,1),zscond(:,2),zscap(:,1), &
513  pswnetsnows,plwnetsnow, &
514  phsnow,ples3l,plel3l,phpsnow, &
515  zct,ztsterm1,ztsterm2,pgfluxsnow)
516 
517 ! - for mass fluxes, snow fraction notion removed
518 
519  zpsn3l(:) = 1.0
520 
521 ELSE
522 
523  zpsn3l(:) = ppsn3l(:)
524 
525 ! - compute surface energy budget and fluxes:
526 
527  CALL snow3lebud(hsnowres, himplicit_wind, &
528  ppew_a_coef, ppew_b_coef, &
529  ppet_a_coef, ppeq_a_coef, ppet_b_coef, ppeq_b_coef, &
530  xsnowdzmin, &
531  pzref,zsnowtemp(:,1),psnowrho(:,1),psnowliq(:,1),zscap(:,1), &
532  zscond(:,1),zscond(:,2), &
533  puref,pexns,pexna,pdircoszw,pvmod, &
534  plw_rad,psw_rad,pta,pqa,pps,ptstep, &
535  psnowdz(:,1),psnowdz(:,2),psnowalb,pz0,pz0eff,pz0h, &
536  zsfcfrz,zradsink(:,1),phpsnow, &
537  zct,pemisnow,prhoa,ztsterm1,ztsterm2,zra,pcdsnow,pchsnow, &
538  zqsat, zdqsat, zrsra, zustar2_ic, pri, &
539  zpet_a_coef_t,zpeq_a_coef_t,zpet_b_coef_t,zpeq_b_coef_t )
540 !
541 ENDIF
542 !
543 ! Heat transfer: simple diffusion along the thermal gradient
544 !
545 zsnowtempo1(:) = zsnowtemp(:,1) ! save surface snow temperature before update
546 !
547 zgrndfluxi(:) = zgrndflux(:)
548 !
549  CALL snow3lsolvt(omeb,ptstep,xsnowdzmin,psnowdz,zscond,zscap,ptg, &
550  psoilcond,pd_g,zradsink,zct,ztsterm1,ztsterm2, &
551  zpet_a_coef_t,zpeq_a_coef_t,zpet_b_coef_t,zpeq_b_coef_t, &
552  zta_ic,zqa_ic,zgrndflux,zgrndfluxo,zsnowtemp,psnowflux )
553 !
554 !
555 !* 8. Surface fluxes
556 ! --------------
557 !
558 ! Since surface fluxes already computed under MEB option (and already
559 ! recomputed for the case when T>Tf), Only call if MEB not in use:
560 !
561 IF(.NOT.omeb)THEN
562 
563  CALL snow3lflux(zsnowtemp(:,1),psnowdz(:,1),pexns,pexna, &
564  zustar2_ic, &
565  ptstep,psnowalb,psw_rad, &
566  pemisnow,zlwupsnow,plw_rad,plwnetsnow, &
567  zta_ic,zsfcfrz,zqa_ic,phpsnow, &
568  zsnowtempo1,psnowflux,zct,zradsink(:,1), &
569  zqsat,zdqsat,zrsra, &
570  prnsnow,phsnow,pgfluxsnow,ples3l,plel3l,pevap, &
571  pustar,gsfcmelt )
572 !
573 ENDIF
574 !
575 !
576 !* 9. Snow melt
577 ! ---------
578 !
579 ! First Test to see if snow pack vanishes during this time step:
580 !
581  CALL snow3lgone(ptstep,plel3l,ples3l,psnowrho, &
582  psnowheat,zradsink(:,inlvls),pevapcor,pthrufal,zgrndflux, &
583  pgfluxsnow,zgrndfluxo,psnowdz,psnowliq,zsnowtemp, &
584  plvtt,plstt,zradxs )
585 !
586 ! For "normal" melt: transform excess heat content into snow liquid:
587 !
588  CALL snow3lmelt(ptstep,zscap,zsnowtemp,psnowdz,psnowrho,psnowliq,zmeltxs)
589 !
590 !
591 !* 10. Snow water flow and refreezing
592 ! ------------------------------
593 ! Liquid water vertical transfer and possible snowpack runoff
594 ! And refreezing/freezing of meltwater/rainfall (ripening of the snow)
595 !
596  CALL snow3lrefrz(ptstep,prr,psnowrho,zsnowtemp,psnowdz,psnowliq,pthrufal)
597 !
598 zscap(:,:) = snow3lscap(psnowrho)
599 psnowheat(:,:) = psnowdz(:,:)*( zscap(:,:)*(zsnowtemp(:,:)-xtt) &
600  - xlmtt*psnowrho(:,:) ) + xlmtt*xrholw*psnowliq(:,:)
601 !
602 !
603 !* 11. Snow Evaporation/Sublimation mass updates:
604 ! ------------------------------------------
605 !
606  CALL snow3levapn(zpsn3l,ples3l,plel3l,ptstep,zsnowtemp(:,1),psnowrho(:,1), &
607  psnowdz,psnowliq(:,1),pta,plvtt,plstt,psnowheat,psoilcor )
608 !
609 ! Update snow temperatures and liquid
610 ! water portion of the snow from snow heat content
611 ! since vertical distribution of enthalpy might have been
612 ! modified in the previous routine:
613 !
614 zscap(:,:) = snow3lscap(psnowrho)
615 zwork2d(:,:) = min(1.0, psnowdz(:,:)/xsnowdmin)
616 zsnowtemp(:,:) = xtt + zwork2d(:,:)*( ((psnowheat(:,:)/max(xsnowdmin,psnowdz(:,:))) &
617  + xlmtt*psnowrho(:,:))/zscap(:,:) )
618 psnowliq(:,:) = max(0.0,zsnowtemp(:,:)-xtt)*zscap(:,:)*psnowdz(:,:)/(xlmtt*xrholw)
619 zsnowtemp(:,:) = min(xtt,zsnowtemp(:,:))
620 !
621 !
622 ! If all snow in uppermost layer evaporates/sublimates, re-distribute
623 ! grid (below could be evoked for vanishingly thin snowpacks):
624 !
625  CALL snow3levapgone(psnowheat,psnowdz,psnowrho,zsnowtemp,psnowliq)
626 !
627 !
628 !* 12. Update surface albedo:
629 ! ----------------------
630 ! Snow clear sky albedo:
631 !
632  CALL snow3lalb(psnowalb,zspectralalbedo,psnowrho(:,1),psnowage(:,1),ppermsnowfrac,pps)
633 !
634 !
635 !* 13. Update snow heat content:
636 ! -------------------------
637 ! Update the heat content (variable stored each time step)
638 ! using current snow temperature and liquid water content:
639 !
640 ! First, make check to make sure heat content not too large
641 ! (this can result due to signifigant heating of thin snowpacks):
642 !
643 DO jj=1,inlvls
644  DO ji=1,ini
645  zliqheatxs(ji) = max(0.0,psnowliq(ji,jj)*xrholw-psnowdz(ji,jj)*psnowrho(ji,jj))*xlmtt/ptstep
646  psnowliq(ji,jj)= psnowliq(ji,jj) - zliqheatxs(ji)*ptstep/(xrholw*xlmtt)
647  psnowliq(ji,jj)= max(0.0, psnowliq(ji,jj))
648  ENDDO
649 ENDDO
650 !
651 psnowtemp(:,:) = zsnowtemp(:,:)
652 !
653 zscap(:,:) = snow3lscap(psnowrho)
654 !
655 psnowheat(:,:) = psnowdz(:,:)*( zscap(:,:)*(psnowtemp(:,:)-xtt) &
656  - xlmtt*psnowrho(:,:) ) + xlmtt*xrholw*psnowliq(:,:)
657 !
658 !
659 !* 14. Snow/Ground heat flux:
660 ! ----------------------
661 !
662 !Add radiation not absorbed by snow to soil/vegetation interface flux
663 !
664 pgrndflux(:) = zgrndfluxo(:)+zradxs(:)
665 !
666 !Comput soil/snow correction heat flux to prevent TG1 numerical oscillations
667 !Add any excess heat for melting and liquid water :
668 !
669 pgflxcor(:) = (zgrndflux(:)-zgrndfluxo(:))+zmeltxs(:)+zliqheatxs(:)
670 !
671 !
672 !* 15. Update snow mass and age:
673 ! -------------------------
674 !
675 psnowswe(:,:) = psnowdz(:,:)*psnowrho(:,:)
676 !
677 WHERE (psnowswe(:,:)>0)
678  psnowage(:,:)=psnowage(:,:)+(ptstep/xday)
679 ELSEWHERE
680  psnowage(:,:)= xundef
681 ENDWHERE
682 !
683 !
684 !* 16. Update surface specific humidity:
685 ! ---------------------------------
686 !
687 IF(omeb)THEN
688  pqs(:) = qsati(psnowtemp(:,1),pps) ! purely diagnostic
689 ELSE
690  pqs(:) = zqsat(:)
691 ENDIF
692 !
693 !
694 !* 17. MEB related checks
695 ! ------------------
696 !
697 IF(omeb)THEN
698 
699 ! Final adjustment: If using MEB, then all excess heat (above initial guess ground flux)
700 ! is used herein to adjust the ground temperature since for the MEB case, boundary
701 ! fluxes are imposed (the energy budget of the ground has already been computed).
702 ! This correction is needed since snow-soil flux here is implicit, while guess
703 ! used in surface energy soil budget was semi-implicit. This forces the flux
704 ! seen by the ground to be *the same* as that leaving the snow (energy conservation).
705 ! Also, add any excessing cooling from sublimation as snowpack becomes vanishingly thin.
706 ! This is added back to total heat content of the snowpack (and distributed among
707 ! all snow layers while conserving total heat content plus correction)
708 
709  zwork(:) = 0.
710  DO jj=1,inlvls
711  DO ji=1,ini
712  zwork(ji) = zwork(ji) + psnowheat(ji,jj)
713  ENDDO
714  ENDDO
715  zwork2(:) = min(0.0, zwork(:) + pgrndflux(:) - zgrndfluxi(:))
716  pgflxcor(:) = pgflxcor(:) + max(0., zwork2(:)) ! add any possible (rare!) excess to soil
717 
718  WHERE(zwork(:) > -1.e-10)
719  zwork(:) = 1. ! i.e. no modifs to H profile
720  ELSEWHERE
721  zwork(:) = zwork2(:)/zwork(:)
722  END WHERE
723 
724  DO jj=1,inlvls
725  DO ji=1,ini
726  psnowheat(ji,jj) = psnowheat(ji,jj)*zwork(ji)
727  ENDDO
728  ENDDO
729 
730 ! these are just updated diagnostics at this point:
731 
732  zwork2d(:,:) = min(1.0, psnowdz(:,:)/xsnowdmin)/max(xsnowdmin,psnowdz(:,:))
733  psnowtemp(:,:) = xtt + zwork2d(:,:)*( (psnowheat(:,:) + xlmtt*psnowswe(:,:))/zscap(:,:) )
734  psnowliq(:,:) = max(0.0,psnowtemp(:,:)-xtt)*zscap(:,:)*psnowdz(:,:)/(xlmtt*xrholw)
735  psnowtemp(:,:) = min(xtt,psnowtemp(:,:))
736 
737 ENDIF
738 !
739 !
740 !* 18. Energy Budget Diagnostics:
741 ! --------------------------
742 !
743 ! NOTE: To check enthalpy conservation, the error in W m-2 is defined HERE as
744 !
745 ! error = (SUM(PSNOWHEAT(:,:),2)-SUM(ZSNOWHEAT0(:,:),2))/PTSTEP &
746 ! - PSNOWHMASS(:)/PTSTEP - (PGFLUXSNOW(:) - PGRNDFLUX(:))
747 !
748 ! where ZSNOWHEAT0 is the value of PSNOWHEAT that enters this routine before any adjustments/processes.
749 ! Errors should be VERY small ( ~0) at each time step.
750 ! The above is "snow relative"...the actual tile or grid box error is multiplied by PPSN3L.
751 !
752 ! Snow heat storage change over the time step (W m-2):
753 
754 pdelheatn(:) = 0.0
755 DO jj=1,inlvls
756  DO ji=1,ini
757  pdelheatn(ji) = pdelheatn(ji) + (psnowheat(ji,jj)-zsnowheat0(ji,jj))
758  ENDDO
759 END DO
760 pdelheatn(:) = pdelheatn(:) /ptstep
761 pdelheatn_sfc(:) = (psnowheat(:,1)-zsnowheat0(:,1))/ptstep
762 
763 !
764 ! The pseudo-heating (W m-2) is derrived as a diagnostic from the surface energy budget
765 ! here...NOTE that the total snow energy budget is explicitly computed, and since
766 ! it balances, then simply compute this term as a residue (which ensures sfc energy balance):
767 !
768 psnowsfch(:) = pdelheatn_sfc(:) - (pswnetsnows(:) +plwnetsnow(:) - phsnow(:) -ples3l(:)-plel3l(:)) &
769  + psnowflux(:) - psnowhmass(:)/ptstep
770 !
771 IF (lhook) CALL dr_hook('SNOW3L',1,zhook_handle)
772 !
773 !
774  CONTAINS
775 !
776 !
777 !
778 !####################################################################
779 !####################################################################
780 !####################################################################
781 SUBROUTINE snow3lfall(PTSTEP,PSR,PTA,PVMOD,PSNOW,PSNOWRHO,PSNOWDZ, &
782  psnowheat,psnowhmass,psnowage,ppermsnowfrac)
783 !
784 !! PURPOSE
785 !! -------
786 ! Calculate changes to snowpack resulting from snowfall.
787 ! Update mass and heat content of uppermost layer.
788 !
789 !
790 USE modd_csts, ONLY : xlmtt, xtt, xci
791 USE modd_snow_par, ONLY : xrhosmin_es, xsnowdmin, &
792  xsnowfall_a_sn, &
793  xsnowfall_b_sn, &
794  xsnowfall_c_sn
795 !
796 USE mode_snow3l
797 !
798 IMPLICIT NONE
799 !
800 !* 0.1 declarations of arguments
801 !
802 REAL, INTENT(IN) :: ptstep
803 !
804 REAL, DIMENSION(:), INTENT(IN) :: psr, pta, pvmod, ppermsnowfrac
805 !
806 REAL, DIMENSION(:), INTENT(INOUT) :: psnow
807 !
808 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowrho, psnowdz, psnowheat, psnowage
809 !
810 REAL, DIMENSION(:), INTENT(OUT) :: psnowhmass
811 !
812 !
813 !* 0.2 declarations of local variables
814 !
815 INTEGER :: jj, ji
816 !
817 INTEGER :: ini
818 INTEGER :: inlvls
819 !
820 REAL, DIMENSION(SIZE(PTA)) :: zsnowfall, zrhosnew, &
821  zsnow, zsnowtemp, &
822  zsnowfall_delta, zscap, &
823  zagenew
824 !
825 REAL(KIND=JPRB) :: zhook_handle
826 !
827 !-------------------------------------------------------------------------------
828 !
829 ! 0. Initialize:
830 ! ------------------
831 !
832 IF (lhook) CALL dr_hook('SNOW3LFALL',0,zhook_handle)
833 !
834 ini = SIZE(psnowdz(:,:),1)
835 inlvls = SIZE(psnowdz(:,:),2)
836 !
837 zrhosnew(:) = xrhosmin_es
838 zagenew(:) = 0.0
839 zsnowfall(:) = 0.0
840 zscap(:) = 0.0
841 zsnow(:) = psnow(:)
842 !
843 psnowhmass(:) = 0.0
844 !
845 ! 1. Incorporate snowfall into snowpack:
846 ! --------------------------------------
847 !
848 !
849 ! Heat content of newly fallen snow (J/m2):
850 ! NOTE for now we assume the snowfall has
851 ! the temperature of the snow surface upon reaching the snow.
852 ! This is done as opposed to using the air temperature since
853 ! this flux is quite small and has little to no impact
854 ! on the time scales of interest. If we use the above assumption
855 ! then, then the snowfall advective heat flux is zero.
856 !
857 zsnowtemp(:) = xtt
858 zscap(:) = snow3lscap(psnowrho(:,1))
859 !
860 WHERE (psr(:) > 0.0 .AND. psnowdz(:,1)>0.)
861  zsnowtemp(:) = xtt + (psnowheat(:,1) + &
862  xlmtt*psnowrho(:,1)*psnowdz(:,1))/ &
863  (zscap(:)*max(xsnowdmin/inlvls,psnowdz(:,1)))
864  zsnowtemp(:) = min(xtt, zsnowtemp(:))
865 END WHERE
866 !
867 WHERE (psr(:) > 0.0)
868 !
869  psnowhmass(:) = psr(:)*(xci*(zsnowtemp(:)-xtt)-xlmtt)*ptstep
870 !
871 ! Snowfall density: Following CROCUS (Pahaut 1976)
872 !
873  zrhosnew(:) = max(xrhosmin_es, xsnowfall_a_sn + xsnowfall_b_sn*(pta(:)-xtt)+ &
874  xsnowfall_c_sn*sqrt(pvmod(:)))
875 !
876 !
877 ! Fresh snowfall changes the snowpack age,
878 ! decreasing in uppermost snow layer (mass weighted average):
879 !
880  psnowage(:,1) = (psnowage(:,1)*psnowdz(:,1)*psnowrho(:,1)+zagenew(:)*psr(:)*ptstep) / &
881  (psnowdz(:,1)*psnowrho(:,1)+psr(:)*ptstep)
882 !
883 ! Augment total pack depth:
884 !
885  zsnowfall(:) = psr(:)*ptstep/zrhosnew(:) ! snowfall thickness (m)
886 !
887  psnow(:) = psnow(:) + zsnowfall(:)
888 !
889 ! Fresh snowfall changes the snowpack
890 ! density, increases the total liquid water
891 ! equivalent: in uppermost snow layer:
892 !
893  psnowrho(:,1) = (psnowdz(:,1)*psnowrho(:,1) + zsnowfall(:)*zrhosnew(:))/ &
894  (psnowdz(:,1)+zsnowfall(:))
895 !
896  psnowdz(:,1) = psnowdz(:,1) + zsnowfall(:)
897 !
898 ! Add energy of snowfall to snowpack:
899 ! Update heat content (J/m2) (therefore the snow temperature
900 ! and liquid content):
901 !
902  psnowheat(:,1) = psnowheat(:,1) + psnowhmass(:)
903 !
904 END WHERE
905 !
906 !
907 ! 2. Case of new snowfall on a previously snow-free surface:
908 ! ----------------------------------------------------------
909 !
910 ! When snow first falls on a surface devoid of snow,
911 ! redistribute the snow mass throughout the 3 layers:
912 ! (temperature already set in the calling routine
913 ! for this case)
914 !
915 zsnowfall_delta(:) = 0.0
916 WHERE(zsnow(:) == 0.0 .AND. psr(:) > 0.0)
917  zsnowfall_delta(:) = 1.0
918 END WHERE
919 !
920 DO jj=1,inlvls
921  DO ji=1,ini
922 !
923  psnowdz(ji,jj) = zsnowfall_delta(ji)*(zsnowfall(ji) /inlvls) + &
924  (1.0-zsnowfall_delta(ji))*psnowdz(ji,jj)
925 !
926  psnowheat(ji,jj) = zsnowfall_delta(ji)*(psnowhmass(ji)/inlvls) + &
927  (1.0-zsnowfall_delta(ji))*psnowheat(ji,jj)
928 !
929  psnowrho(ji,jj) = zsnowfall_delta(ji)*zrhosnew(ji) + &
930  (1.0-zsnowfall_delta(ji))*psnowrho(ji,jj)
931 !
932  psnowage(ji,jj) = zsnowfall_delta(ji)*(zagenew(ji)/inlvls) + &
933  (1.0-zsnowfall_delta(ji))*psnowage(ji,jj)
934 !
935  ENDDO
936 ENDDO
937 !
938 IF (lhook) CALL dr_hook('SNOW3LFALL',1,zhook_handle)
939 !
940 !
941 END SUBROUTINE snow3lfall
942 !####################################################################
943 !####################################################################
944 !####################################################################
945 SUBROUTINE snow3lcompactn(PTSTEP,PSNOWDZMIN,PSNOWRHO,PSNOWDZ,PSNOWTEMP,PSNOW,PSNOWLIQ)
946 !
947 !! PURPOSE
948 !! -------
949 ! Snow compaction due to overburden and settling.
950 ! Mass is unchanged: layer thickness is reduced
951 ! in proportion to density increases. Method
952 ! of Brun et al (1989) and Vionnet et al. (2012)
953 !
954 !
955 USE modd_surf_par, ONLY : xundef
956 !
957 USE modd_csts, ONLY : xtt, xg
958 USE modd_snow_par, ONLY : xrhosmax_es
959 !
960 USE modd_snow_metamo, ONLY : xvvisc1,xvvisc3,xvvisc4, &
961  xvvisc5,xvvisc6,xvro11
962 !
963 IMPLICIT NONE
964 !
965 !* 0.1 declarations of arguments
966 !
967 REAL, INTENT(IN) :: ptstep
968 REAL, INTENT(IN) :: psnowdzmin
969 !
970 REAL, DIMENSION(:,:), INTENT(IN) :: psnowtemp, psnowliq
971 !
972 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowrho, psnowdz
973 !
974 REAL, DIMENSION(:), INTENT(OUT) :: psnow
975 !
976 !
977 !* 0.2 declarations of local variables
978 !
979 INTEGER :: jj, ji
980 !
981 INTEGER :: ini
982 INTEGER :: inlvls
983 !
984 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowrho2, zviscocity, zf1, &
985  ztemp, zsmass, zsnowdz, &
986  zwsnowdz, zwholdmax
987 !
988 !
989 REAL(KIND=JPRB) :: zhook_handle
990 !
991 !-------------------------------------------------------------------------------
992 !
993 ! 0. Initialization:
994 ! ------------------
995 !
996 IF (lhook) CALL dr_hook('SNOW3LCOMPACTN',0,zhook_handle)
997 !
998 ini = SIZE(psnowdz(:,:),1)
999 inlvls = SIZE(psnowdz(:,:),2)
1000 !
1001 zsnowrho2(:,:) = psnowrho(:,:)
1002 zsnowdz(:,:) = max(psnowdzmin,psnowdz(:,:))
1003 zviscocity(:,:) = 0.0
1004 ztemp(:,:) = 0.0
1005 !
1006 ! 1. Cumulative snow mass (kg/m2):
1007 ! --------------------------------
1008 !
1009 zsmass(:,:) = 0.0
1010 DO jj=2,inlvls
1011  DO ji=1,ini
1012  zsmass(ji,jj) = zsmass(ji,jj-1) + psnowdz(ji,jj-1)*psnowrho(ji,jj-1)
1013  ENDDO
1014 ENDDO
1015 ! overburden of half the mass of the uppermost layer applied to itself
1016 zsmass(:,1) = 0.5 * psnowdz(:,1) * psnowrho(:,1)
1017 !
1018 ! 2. Compaction
1019 ! -------------
1020 !
1021 !Liquid water effect
1022 !
1023 zwholdmax(:,:) = snow3lhold(psnowrho,psnowdz)
1024 zf1(:,:) = 1.0/(xvvisc5+10.*min(1.0,psnowliq(:,:)/zwholdmax(:,:)))
1025 !
1026 !Snow viscocity, density and grid thicknesses
1027 !
1028 DO jj=1,inlvls
1029  DO ji=1,ini
1030 !
1031  IF(psnowrho(ji,jj) < xrhosmax_es)THEN
1032 !
1033 ! temperature dependence limited to 5K: Schleef et al. (2014)
1034  ztemp(ji,jj) = xvvisc4*min(5.0,abs(xtt-psnowtemp(ji,jj)))
1035 !
1036 ! Calculate snow viscocity: Brun et al. (1989), Vionnet et al. (2012)
1037  zviscocity(ji,jj) = xvvisc1*zf1(ji,jj)*exp(xvvisc3*psnowrho(ji,jj)+ztemp(ji,jj))*psnowrho(ji,jj)/xvro11
1038 !
1039 ! Calculate snow density:
1040  zsnowrho2(ji,jj) = psnowrho(ji,jj) + psnowrho(ji,jj)*ptstep &
1041  * ( (xg*zsmass(ji,jj)/zviscocity(ji,jj)) )
1042 !
1043 ! Conserve mass by decreasing grid thicknesses in response to density increases
1044  psnowdz(ji,jj) = psnowdz(ji,jj)*(psnowrho(ji,jj)/zsnowrho2(ji,jj))
1045 !
1046  ENDIF
1047 !
1048  ENDDO
1049 ENDDO
1050 !
1051 ! 3. Update total snow depth and density profile:
1052 ! -----------------------------------------------
1053 !
1054 ! Compaction/augmentation of total snowpack depth
1055 !
1056 psnow(:) = 0.
1057 DO jj=1,inlvls
1058  DO ji=1,ini
1059  psnow(ji) = psnow(ji) + psnowdz(ji,jj)
1060  ENDDO
1061 ENDDO
1062 !
1063 ! Update density (kg m-3):
1064 !
1065 psnowrho(:,:) = zsnowrho2(:,:)
1066 !
1067 IF (lhook) CALL dr_hook('SNOW3LCOMPACTN',1,zhook_handle)
1068 !
1069 !-------------------------------------------------------------------------------
1070 !
1071 END SUBROUTINE snow3lcompactn
1072 !####################################################################
1073 !####################################################################
1074 !####################################################################
1075 SUBROUTINE snow3ldrift(PTSTEP,PFORESTFRAC,PVMOD,PTA,PQA,PPS,PRHOA,&
1076  psnowrho,psnowdz,psnow,osnowdrift_sublim,psndrift)
1077 !
1078 USE modd_surf_par, ONLY : xundef
1079 !
1080 USE modd_snow_par, ONLY : xvtime, xvromax, xvromin, xvmob1, &
1081  xvdrift1, xvdrift2, xvdrift3, &
1082  xcoef_ff, xcoef_effect, xqs_ref
1083 !
1084 USE mode_thermos
1085 !
1086 !
1087 !! PURPOSE
1088 !! -------
1089 ! Snow compaction and metamorphism due to drift
1090 ! Mass is unchanged: layer thickness is reduced
1091 ! in proportion to density increases. Method inspired from
1092 ! Brun et al. (1997) and Guyomarch
1093 !
1094 ! - computes a mobility index of each snow layer from its density
1095 ! - computes a drift index of each layer from its mobility and wind speed
1096 ! - computes a transport index
1097 ! - increases density in case of transport
1098 !
1099 ! HISTORY:
1100 ! Basic parameterization from Crocus/ARPEGE Coupling (1997)
1101 ! Adaptation from Crocus SNOWDRIFT subroutine
1102 !
1103 !! PURPOSE
1104 !! -------
1105 ! Snow compaction due to overburden and settling.
1106 ! Mass is unchanged: layer thickness is reduced
1107 ! in proportion to density increases.
1108 !
1109 IMPLICIT NONE
1110 !
1111 !* 0.1 declarations of arguments
1112 !
1113 REAL, INTENT(IN) :: ptstep
1114 !
1115 REAL, DIMENSION(:), INTENT(IN) :: pforestfrac
1116 REAL, DIMENSION(:), INTENT(IN) :: pvmod
1117 REAL, DIMENSION(:), INTENT(IN) :: pta
1118 REAL, DIMENSION(:), INTENT(IN) :: pqa
1119 REAL, DIMENSION(:), INTENT(IN) :: pps
1120 REAL, DIMENSION(:), INTENT(IN) :: prhoa
1121 !
1122 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowrho, psnowdz
1123 !
1124 REAL, DIMENSION(:), INTENT(OUT) :: psnow
1125 !
1126 LOGICAL, INTENT(IN) :: osnowdrift_sublim
1127 REAL, DIMENSION(:), INTENT(OUT) :: psndrift !blowing snow sublimation (kg/m2/s)
1128 !
1129 !* 0.2 declarations of local variables
1130 !
1131 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowrho2
1132 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zrmob
1133 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zrdrift
1134 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zrt
1135 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zdro
1136 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zqs_effect
1137 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zdrift_effect
1138 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: zprofequ
1139 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: zwind
1140 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: zqsati
1141 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: zvt
1142 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: zqs
1143 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: zw
1144 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: zt
1145 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: zsnowdz1
1146 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: zforest_effect
1147 !
1148 LOGICAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: gdrift
1149 !
1150 INTEGER :: jj, ji
1151 !
1152 INTEGER :: ini
1153 INTEGER :: inlvls
1154 !
1155 REAL(KIND=JPRB) :: zhook_handle
1156 !
1157 !-------------------------------------------------------------------------------
1158 !
1159 IF (lhook) CALL dr_hook('SNOW3LDRIFT',0,zhook_handle)
1160 !
1161 ! 0. Initialization:
1162 ! ------------------
1163 !
1164 ini = SIZE(psnowdz(:,:),1)
1165 inlvls = SIZE(psnowdz(:,:),2)
1166 !
1167 zsnowrho2(:,:) = psnowrho(:,:)
1168 zrdrift(:,:) = xundef
1169 zrt(:,:) = xundef
1170 zdro(:,:) = xundef
1171 zqs_effect(:,:) = 0.0
1172 zdrift_effect(:,:) = 0.0
1173 gdrift(:,:) = .false.
1174 !
1175 zprofequ(:) = 0.0
1176 !
1177 zsnowdz1(:) = psnowdz(:,1)
1178 !
1179 ! 1. Initialazation of drift and induced settling
1180 ! -----------------------------------------------
1181 !
1182 !Limitation of wind speed in Forest : ~ 15% of wind in open area
1183 zforest_effect(:) = 1.0 - 0.85 * pforestfrac(:)
1184 !
1185 zwind(:) = xcoef_ff * zforest_effect(:) * pvmod(:)
1186 !
1187 DO jj=1,inlvls
1188  DO ji=1,ini
1189 !
1190 ! mobility index computation of a layer as a function of its density
1191  zrmob(ji,jj)= 1.25-1.25e-3*(max(psnowrho(ji,jj),xvromin)-xvromin) / xvmob1
1192 !
1193 ! computation of the drift index inclunding the decay by overburden snow
1194  zrdrift(ji,jj) = zrmob(ji,jj)-(xvdrift1*exp(-xvdrift2*zwind(ji))-1.0)
1195 !
1196  ENDDO
1197 ENDDO
1198 !
1199 !Compaction from the surface to layers beneath until
1200 !a snow layer having negative drift index
1201 gdrift(:,:) = (zrdrift(:,:)>0.0)
1202 DO ji=1,ini
1203  DO jj=1,inlvls
1204  IF(.NOT.gdrift(ji,jj))THEN
1205  gdrift(ji,jj:inlvls)=.false.
1206  EXIT
1207  ENDIF
1208  ENDDO
1209 ENDDO
1210 
1211 ! 2. Specific case for blowing snow sublimation
1212 ! ---------------------------------------------
1213 !
1214 IF(osnowdrift_sublim)THEN
1215 !
1216  zqsati(:)= qsati(pta(:),pps(:))
1217 !
1218 ! computation of wind speed threshold QSATI and RH withe respect to ice
1219 !
1220  zw(:)=max(-0.99,zrmob(:,1))
1221  zvt(:)=log(xvdrift1/(1.0+zw(:)))/xvdrift2
1222 !
1223  WHERE(zwind(:)>0.0)
1224  zw(:)=log(zwind(:)/zvt(:))
1225  zw(:)=exp(3.6*zw(:))
1226  ELSEWHERE
1227  zw(:)=0.0
1228  ENDWHERE
1229 !
1230 ! computation of sublimation rate according to Gordon's PhD
1231 !
1232  zt(:)=xtt/pta(:)
1233  zt(:)=0.0018*(zt(:)**4)
1234 !
1235  zqs(:)=zt(:)*zvt(:)*prhoa(:)*zqsati(:)*(1.-pqa(:)/zqsati(:))*zw(:)
1236 !
1237 ! surface depth decrease in case of blowing snow sublimation
1238 !
1239  psnowdz(:,1)=max(0.5*psnowdz(:,1),psnowdz(:,1)-max(0.,zqs(:))*ptstep/(xcoef_ff*psnowrho(:,1)))
1240 !
1241  psndrift(:) = (zsnowdz1(:)-psnowdz(:,1))*psnowrho(:,1)/ptstep
1242 !
1243  zqs_effect(:,1)=min(3.,max(0.,zqs(:))/xqs_ref)
1244 !
1245 ENDIF
1246 !
1247 ! 3. Computation of drift and induced settling
1248 ! --------------------------------------------
1249 !
1250 DO jj=1,inlvls
1251  DO ji=1,ini
1252 
1253 ! update the decay coeff by half the current layer
1254  zprofequ(ji) = zprofequ(ji) + 0.5 * psnowdz(ji,jj) * 0.1 * (xvdrift3-zrdrift(ji,jj))
1255 !
1256  IF(gdrift(ji,jj).AND.psnowrho(ji,jj)<xvromax)THEN
1257 !
1258 ! computation of the drift index inclunding the decay by overburden snow
1259  zrt(ji,jj) = max(0.0,zrdrift(ji,jj)*exp(-zprofequ(ji)*100.0))
1260 !
1261  zdrift_effect(ji,jj) = (zqs_effect(ji,jj)+xcoef_effect)*zrt(ji,jj)/(xvtime*xcoef_ff)
1262 !
1263 ! settling by wind transport only in case of not too dense snow
1264  zdro(ji,jj) = (xvromax - psnowrho(ji,jj)) * zdrift_effect(ji,jj) * ptstep
1265 !
1266 ! Calculate new snow density:
1267  zsnowrho2(ji,jj) = min(xvromax,psnowrho(ji,jj)+zdro(ji,jj))
1268 
1269 ! Conserve mass by decreasing grid thicknesses in response to density increases
1270  psnowdz(ji,jj) = psnowdz(ji,jj)*(psnowrho(ji,jj)/zsnowrho2(ji,jj))
1271 !
1272  ENDIF
1273 !
1274 ! update the decay coeff by half the current layer
1275  zprofequ(ji) = zprofequ(ji) + 0.5 * psnowdz(ji,jj) * 0.1 * (xvdrift3-zrdrift(ji,jj))
1276 !
1277  ENDDO
1278 ENDDO
1279 !
1280 ! 4. Update total snow depth and density profile:
1281 ! -----------------------------------------------
1282 !
1283 ! Compaction/augmentation of total snowpack depth
1284 !
1285 psnow(:) = 0.
1286 DO jj=1,inlvls
1287  DO ji=1,ini
1288  psnow(ji) = psnow(ji) + psnowdz(ji,jj)
1289  ENDDO
1290 ENDDO
1291 !
1292 ! Update density (kg m-3):
1293 !
1294 psnowrho(:,:) = zsnowrho2(:,:)
1295 !
1296 IF (lhook) CALL dr_hook('SNOW3LDRIFT',1,zhook_handle)
1297 !
1298 !-------------------------------------------------------------------------------
1299 !
1300 END SUBROUTINE snow3ldrift
1301 !####################################################################
1302 !####################################################################
1303 !####################################################################
1304  SUBROUTINE snow3ltransf(PSNOW,PSNOWDZ,PSNOWDZN, &
1305  psnowrho,psnowheat,psnowage)
1306 !
1307 !! PURPOSE
1308 !! -------
1309 ! Snow mass,heat and characteristics redistibution in case of
1310 ! grid resizing. Total mass and heat content of the overall snowpack
1311 ! unchanged/conserved within this routine.
1312 ! Same method as in Crocus
1313 !
1314 USE modd_snow_par, ONLY : xsnowcritd
1315 !
1316 IMPLICIT NONE
1317 !
1318 !
1319 !* 0.1 declarations of arguments
1320 !
1321 REAL, DIMENSION(: ), INTENT(IN) :: psnow
1322 !
1323 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowdzn
1324 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowheat
1325 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowrho
1326 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowdz
1327 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowage
1328 !
1329 !* 0.2 declarations of local variables
1330 !
1331 INTEGER :: ji, jl, jlo
1332 !
1333 INTEGER :: ini
1334 INTEGER :: inlvls
1335 !
1336 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowrhon
1337 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowheatn
1338 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowagen
1339 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowztop_new
1340 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowzbot_new
1341 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowrhoo
1342 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowheato
1343 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowageo
1344 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowdzo
1345 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowztop_old
1346 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowzbot_old
1347 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowhean
1348 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowagn
1349 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zmastotn
1350 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zmassdzo
1351 !
1352 REAL, DIMENSION(SIZE(PSNOW)) :: zpsnow_old, zpsnow_new
1353 REAL, DIMENSION(SIZE(PSNOW)) :: zsumheat, zsumswe, zsumage, zsnowmix_delta
1354 !
1355 REAL :: zpropor
1356 !
1357 REAL(KIND=JPRB) :: zhook_handle
1358 !
1359 !-------------------------------------------------------------------------------
1360 !
1361 ! 0. Initialization:
1362 ! ------------------
1363 !
1364 !
1365 IF (lhook) CALL dr_hook('SNOW3LTRANSF',0,zhook_handle)
1366 !
1367 ini = SIZE(psnowrho,1)
1368 inlvls = SIZE(psnowrho,2)
1369 !
1370 zpsnow_new(:) = 0.0
1371 zpsnow_old(:) = psnow(:)
1372 !
1373 DO jl=1,inlvls
1374  DO ji=1,ini
1375  zpsnow_new(ji)=zpsnow_new(ji)+psnowdzn(ji,jl)
1376  ENDDO
1377 ENDDO
1378 !
1379 ! initialization of variables describing the initial snowpack
1380 !
1381 zsnowdzo(:,:) = psnowdz(:,:)
1382 zsnowrhoo(:,:) = psnowrho(:,:)
1383 zsnowheato(:,:) = psnowheat(:,:)
1384 zsnowageo(:,:) = psnowage(:,:)
1385 zmassdzo(:,:) = xundef
1386 !
1387 ! 1. Calculate vertical grid limits (m):
1388 ! --------------------------------------
1389 !
1390 zsnowztop_old(:,1) = zpsnow_old(:)
1391 zsnowztop_new(:,1) = zpsnow_new(:)
1392 zsnowzbot_old(:,1) = zsnowztop_old(:,1)-zsnowdzo(:,1)
1393 zsnowzbot_new(:,1) = zsnowztop_new(:,1)-psnowdzn(:,1)
1394 !
1395 DO jl=2,inlvls
1396  DO ji=1,ini
1397  zsnowztop_old(ji,jl) = zsnowzbot_old(ji,jl-1)
1398  zsnowztop_new(ji,jl) = zsnowzbot_new(ji,jl-1)
1399  zsnowzbot_old(ji,jl) = zsnowztop_old(ji,jl )-zsnowdzo(ji,jl)
1400  zsnowzbot_new(ji,jl) = zsnowztop_new(ji,jl )-psnowdzn(ji,jl)
1401  ENDDO
1402 ENDDO
1403 zsnowzbot_old(:,inlvls)=0.0
1404 zsnowzbot_new(:,inlvls)=0.0
1405 !
1406 ! 3. Calculate mass, heat, charcateristics mixing due to vertical grid resizing:
1407 ! --------------------------------------------------------------------
1408 !
1409 ! loop over the new snow layers
1410 ! Summ or avergage of the constituting quantities of the old snow layers
1411 ! which are totally or partially inserted in the new snow layer
1412 ! For snow age, mass weighted average is used.
1413 !
1414 zsnowhean(:,:)=0.0
1415 zmastotn(:,:)=0.0
1416 zsnowagn(:,:)=0.0
1417 !
1418 DO jl=1,inlvls
1419  DO jlo=1, inlvls
1420  DO ji=1,ini
1421  IF((zsnowztop_old(ji,jlo)>zsnowzbot_new(ji,jl)).AND.(zsnowzbot_old(ji,jlo)<zsnowztop_new(ji,jl)))THEN
1422 !
1423  zpropor = (min(zsnowztop_old(ji,jlo), zsnowztop_new(ji,jl)) &
1424  - max(zsnowzbot_old(ji,jlo), zsnowzbot_new(ji,jl)))&
1425  / zsnowdzo(ji,jlo)
1426 !
1427  zmassdzo(ji,jlo)=zsnowrhoo(ji,jlo)*zsnowdzo(ji,jlo)*zpropor
1428 !
1429  zmastotn(ji,jl)=zmastotn(ji,jl)+zmassdzo(ji,jlo)
1430  zsnowagn(ji,jl)=zsnowagn(ji,jl)+zsnowageo(ji,jlo)*zmassdzo(ji,jlo)
1431 !
1432  zsnowhean(ji,jl)=zsnowhean(ji,jl)+zsnowheato(ji,jlo)*zpropor
1433 !
1434  ENDIF
1435  ENDDO
1436  ENDDO
1437 ENDDO
1438 !
1439 ! the new layer inherits from the weighted average properties of the old ones
1440 ! heat and mass
1441 !
1442 zsnowheatn(:,:)= zsnowhean(:,:)
1443 zsnowagen(:,:)= zsnowagn(:,:)/zmastotn(:,:)
1444 zsnowrhon(:,:)= zmastotn(:,:)/psnowdzn(:,:)
1445 !
1446 ! 4. Vanishing or very thin snowpack check:
1447 ! -----------------------------------------
1448 !
1449 ! NOTE: ONLY for very shallow snowpacks, mix properties (homogeneous):
1450 ! this avoids problems related to heat and mass exchange for
1451 ! thin layers during heavy snowfall or signifigant melt: one
1452 ! new/old layer can exceed the thickness of several old/new layers.
1453 ! Therefore, mix (conservative):
1454 !
1455 zsumheat(:) = 0.0
1456 zsumswe(:) = 0.0
1457 zsumage(:) = 0.0
1458 zsnowmix_delta(:) = 0.0
1459 !
1460 DO jl=1,inlvls
1461  DO ji=1,ini
1462  IF(psnow(ji) < xsnowcritd)THEN
1463  zsumheat(ji) = zsumheat(ji) + psnowheat(ji,jl)
1464  zsumswe(ji) = zsumswe(ji) + psnowrho(ji,jl)*psnowdz(ji,jl)
1465  zsumage(ji) = zsumage(ji) + psnowage(ji,jl)
1466  zsnowmix_delta(ji) = 1.0
1467  ENDIF
1468  ENDDO
1469 ENDDO
1470 !
1471 ! Heat and mass are evenly distributed vertically:
1472 ! heat and mass (density and thickness) are constant
1473 ! in profile:
1474 !
1475 DO jl=1,inlvls
1476  DO ji=1,ini
1477 !
1478  zsnowheatn(ji,jl) = zsnowmix_delta(ji)*(zsumheat(ji)/inlvls) + &
1479  (1.0-zsnowmix_delta(ji))*zsnowheatn(ji,jl)
1480 !
1481  psnowdzn(ji,jl) = zsnowmix_delta(ji)*(psnow(ji)/inlvls) + &
1482  (1.0-zsnowmix_delta(ji))*psnowdzn(ji,jl)
1483 !
1484  zsnowrhon(ji,jl) = zsnowmix_delta(ji)*(zsumswe(ji)/psnow(ji)) + &
1485  (1.0-zsnowmix_delta(ji))*zsnowrhon(ji,jl)
1486 !
1487  zsnowagen(ji,jl) = zsnowmix_delta(ji)*(zsumage(ji)/inlvls) + &
1488  (1.0-zsnowmix_delta(ji))*zsnowagen(ji,jl)
1489 !
1490  ENDDO
1491 ENDDO
1492 !
1493 ! 5. Update mass (density and thickness) and heat:
1494 ! ------------------------------------------------
1495 !
1496 psnowdz(:,:) = psnowdzn(:,:)
1497 psnowrho(:,:) = zsnowrhon(:,:)
1498 psnowheat(:,:) = zsnowheatn(:,:)
1499 psnowage(:,:) = zsnowagen(:,:)
1500 !
1501 IF (lhook) CALL dr_hook('SNOW3LTRANSF',1,zhook_handle)
1502 !
1503 !
1504 !-------------------------------------------------------------------------------
1505 !
1506 END SUBROUTINE snow3ltransf
1507 !####################################################################
1508 !####################################################################
1509 !####################################################################
1510  SUBROUTINE snow3lrad(OMEB, PSNOWDZMIN, PSW_RAD, PSNOWALB, &
1511  pspectralalbedo, psnowdz, psnowrho, palb, &
1512  ppermsnowfrac, pzenith, pswnetsnow, &
1513  pswnetsnows, pradsink, pradxs, psnowage )
1514 !
1515 !! PURPOSE
1516 !! -------
1517 ! Calculate the transmission of shortwave (solar) radiation
1518 ! through the snowpack (using a form of Beer's Law: exponential
1519 ! decay of radiation with increasing snow depth).
1520 !
1521 USE modd_snow_par, ONLY : xvspec1,xvspec2,xvspec3,xvbeta1,xvbeta2, &
1522  xvbeta4,xvbeta3,xvbeta5, xmincoszen
1523 USE modd_meb_par, ONLY : xsw_wght_vis, xsw_wght_nir
1524 !
1525 USE mode_snow3l, ONLY : snow3ldopt
1526 !
1527 IMPLICIT NONE
1528 !
1529 !* 0.1 declarations of arguments
1530 !
1531 LOGICAL, INTENT(IN) :: omeb ! if=T, then uppermost abs is diagnosed
1532 ! ! since fluxes known
1533 !
1534 REAL, INTENT(IN) :: psnowdzmin
1535 !
1536 REAL, DIMENSION(:), INTENT(IN) :: psw_rad
1537 REAL, DIMENSION(:), INTENT(IN) :: psnowalb
1538 REAL, DIMENSION(:), INTENT(IN) :: palb
1539 REAL, DIMENSION(:), INTENT(IN) :: ppermsnowfrac
1540 REAL, DIMENSION(:), INTENT(IN) :: pzenith
1541 !
1542 REAL, DIMENSION(:,:), INTENT(IN) :: psnowrho, psnowdz, psnowage
1543 REAL, DIMENSION(:,:), INTENT(IN) :: pspectralalbedo
1544 !
1545 REAL, DIMENSION(:), INTENT(INOUT) :: pswnetsnow, pswnetsnows
1546 !
1547 REAL, DIMENSION(:), INTENT(OUT) :: pradxs
1548 !
1549 REAL, DIMENSION(:,:), INTENT(OUT) :: pradsink
1550 !
1551 !
1552 !* 0.2 declarations of local variables
1553 !
1554 INTEGER :: jj, ji
1555 !
1556 INTEGER :: ini
1557 INTEGER :: inlvls
1558 !
1559 REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: zradtot, zprojlat, zcoszen
1560 REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: zopticalpath1, zopticalpath2, zopticalpath3
1561 !
1562 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zdsgrain, zcoef, zsnowdz, zage
1563 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zbeta1, zbeta2, zbeta3, zwork
1564 REAL, DIMENSION(SIZE(PSPECTRALALBEDO,1),SIZE(PSPECTRALALBEDO,2)) :: zspectralalbedo
1565 !
1566 REAL(KIND=JPRB) :: zhook_handle
1567 !-------------------------------------------------------------------------------
1568 !
1569 ! 0. Initialization:
1570 ! ------------------
1571 !
1572 IF (lhook) CALL dr_hook('SNOW3LRAD',0,zhook_handle)
1573 !
1574 ini = SIZE(psnowdz(:,:),1)
1575 inlvls = SIZE(psnowdz(:,:),2)
1576 !
1577 zspectralalbedo(:,:) = 0. ! Init
1578 !
1579 ! 1. Vanishingly thin snowpack check:
1580 ! -----------------------------------
1581 ! For vanishingly thin snowpacks, much of the radiation
1582 ! can pass through snowpack into underlying soil, making
1583 ! a large (albeit temporary) thermal gradient: by imposing
1584 ! a minimum thickness, this increases the radiation absorbtion
1585 ! for vanishingly thin snowpacks.
1586 !
1587 zsnowdz(:,:) = max(psnowdzmin, psnowdz(:,:))
1588 !
1589 !
1590 ! 2. Extinction of net shortwave radiation
1591 ! ----------------------------------------
1592 ! Fn of snow depth and density (Loth and Graf 1993:
1593 ! SNOWCVEXT => from Bohren and Barkstrom 1974
1594 ! SNOWAGRAIN and SNOWBGRAIN=> from Jordan 1976)
1595 !
1596 ! Coefficient for taking into account the increase of path length of rays
1597 ! in snow due to zenithal angle
1598 !
1599 zcoszen(:)=max(xmincoszen,cos(pzenith(:)))
1600 !
1601 ! This formulation is incorrect but it compensate partly the fact that
1602 ! the albedo formulation does not account for zenithal angle.
1603 ! Only for polar or glacier regions
1604 !
1605 zprojlat(:)=(1.0-ppermsnowfrac(:))+ppermsnowfrac(:)/zcoszen(:)
1606 !
1607 ! Snow optical grain diameter (no age dependency over polar regions):
1608 !
1609 DO jj=1,inlvls
1610  DO ji=1,ini
1611  zage(ji,jj) = (1.0-ppermsnowfrac(ji))*psnowage(ji,jj)
1612  ENDDO
1613 ENDDO
1614 !
1615 zdsgrain(:,:) = snow3ldopt(psnowrho,zage)
1616 !
1617 ! Extinction coefficient from Brun et al. (1989):
1618 !
1619 zwork(:,:)=sqrt(zdsgrain(:,:))
1620 !
1621 zbeta1(:,:)=max(xvbeta1*psnowrho(:,:)/zwork(:,:),xvbeta2)
1622 zbeta2(:,:)=max(xvbeta3*psnowrho(:,:)/zwork(:,:),xvbeta4)
1623 zbeta3(:,:)=xvbeta5
1624 !
1625 zopticalpath1(:) = 0.0
1626 zopticalpath2(:) = 0.0
1627 zopticalpath3(:) = 0.0
1628 !
1629 IF(omeb)THEN
1630 
1631 ! Only 2 bands currently considered (for vegetation and soil)
1632 ! thus eliminate a band and renormalize (as was done for the surface snow albedo
1633 ! for the MEB snow surface energy budget computations)
1634 
1635  zspectralalbedo(:,1) = pspectralalbedo(:,1)
1636  zspectralalbedo(:,2) = (psnowalb(:) - xsw_wght_vis*zspectralalbedo(:,1))/xsw_wght_nir
1637 
1638  DO jj=1,inlvls
1639  DO ji=1,ini
1640  !
1641  zopticalpath1(ji) = zopticalpath1(ji) + zbeta1(ji,jj)*zsnowdz(ji,jj)
1642  zopticalpath2(ji) = zopticalpath2(ji) + zbeta2(ji,jj)*zsnowdz(ji,jj)
1643  !
1644  zcoef(ji,jj) = xsw_wght_vis*(1.0-zspectralalbedo(ji,1))*exp(-zopticalpath1(ji)*zprojlat(ji)) &
1645  + xsw_wght_nir*(1.0-zspectralalbedo(ji,2))*exp(-zopticalpath2(ji)*zprojlat(ji))
1646  !
1647  ENDDO
1648  ENDDO
1649 
1650 ! diagnose surface layer coef (should be very close/identical to ZCOEF(:,1) computed above)
1651 
1652  zcoef(:,1) = 1.0 - pswnetsnows(:)/max(1.e-4,pswnetsnow(:))
1653 
1654 ELSE
1655 !
1656  DO jj=1,inlvls
1657  DO ji=1,ini
1658  !
1659  zopticalpath1(ji) = zopticalpath1(ji) + zbeta1(ji,jj)*zsnowdz(ji,jj)
1660  zopticalpath2(ji) = zopticalpath2(ji) + zbeta2(ji,jj)*zsnowdz(ji,jj)
1661  zopticalpath3(ji) = zopticalpath3(ji) + zbeta3(ji,jj)*zsnowdz(ji,jj)
1662  !
1663  zcoef(ji,jj) = xvspec1*(1.0-pspectralalbedo(ji,1))*exp(-zopticalpath1(ji)*zprojlat(ji)) &
1664  + xvspec2*(1.0-pspectralalbedo(ji,2))*exp(-zopticalpath2(ji)*zprojlat(ji)) &
1665  + xvspec3*(1.0-pspectralalbedo(ji,3))*exp(-zopticalpath3(ji)*zprojlat(ji))
1666  !
1667  ENDDO
1668  ENDDO
1669 
1670  pswnetsnow(:) = psw_rad(:)*(1.-psnowalb(:))
1671  pswnetsnows(:) = pswnetsnow(:)*(1.0-zcoef(:,1))
1672 
1673 ENDIF
1674 !
1675 ! 3. Radiation at each level: (W/m2)
1676 ! ----------------------------------
1677 !
1678 ! NOTE: as this is a "sink" term for heat, it
1679 ! as assigned a negative value.
1680 DO jj=1,inlvls
1681  DO ji=1,ini
1682  pradsink(ji,jj) = -psw_rad(ji)*zcoef(ji,jj)
1683  ENDDO
1684 ENDDO
1685 !
1686 ! For thin snow packs, radiation might reach base of
1687 ! snowpack...so we influence this amount with sfc albedo
1688 ! and (outside of this routine) add any excess heat
1689 ! to underlying soil:
1690 !
1691 pradsink(:,inlvls) = pradsink(:,inlvls)*(1.0-palb(:))
1692 !
1693 ! 4. Excess radiation not absorbed by snowpack (W/m2):
1694 ! ----------------------------------------------------
1695 !
1696 zradtot(:) = pradsink(:,1) + (1.-psnowalb(:))*psw_rad(:)
1697 DO jj=2,inlvls
1698  DO ji=1,ini
1699  zradtot(ji) = zradtot(ji) + pradsink(ji,jj)-pradsink(ji,jj-1)
1700  ENDDO
1701 ENDDO
1702 !
1703 pradxs(:) = (1.-psnowalb(:))*psw_rad(:) - zradtot(:)
1704 !
1705 IF (lhook) CALL dr_hook('SNOW3LRAD',1,zhook_handle)
1706 !
1707 !-------------------------------------------------------------------------------
1708 !
1709 END SUBROUTINE snow3lrad
1710 !####################################################################
1711 !####################################################################
1712 !####################################################################
1713  SUBROUTINE snow3lebud(HSNOWRES, HIMPLICIT_WIND, &
1714  ppew_a_coef, ppew_b_coef, &
1715  ppet_a_coef, ppeq_a_coef, ppet_b_coef, ppeq_b_coef, &
1716  psnowdzmin, &
1717  pzref,pts,psnowrho,psnowliq,pscap,pscond1,pscond2, &
1718  puref,pexns,pexna,pdircoszw,pvmod, &
1719  plw_rad,psw_rad,pta,pqa,pps,ptstep, &
1720  psnowdz1,psnowdz2,palbt,pz0,pz0eff,pz0h, &
1721  psfcfrz,pradsink,phpsnow, &
1722  pct,pemist,prhoa,ptsterm1,ptsterm2,pra,pcdsnow,pchsnow, &
1723  pqsat,pdqsat,prsra,pustar2_ic,pri, &
1724  ppet_a_coef_t,ppeq_a_coef_t,ppet_b_coef_t,ppeq_b_coef_t )
1725 !
1726 !! PURPOSE
1727 !! -------
1728 ! Calculate surface energy budget linearization (terms) and turbulent
1729 ! exchange coefficients/resistance between surface and atmosphere.
1730 ! (Noilhan and Planton 1989; Giordani 1993; Noilhan and Mahfouf 1996)
1731 !
1732 USE modd_csts, ONLY : xcpd, xrholw, xstefan, xlvtt, xlstt
1733 USE modd_snow_par, ONLY : x_ri_max, xemissn
1734 !
1735 USE mode_thermos
1736 !
1737 USE modi_surface_ri
1738 USE modi_surface_aero_cond
1739 USE modi_surface_cd
1740 !
1741 !
1742 IMPLICIT NONE
1743 !
1744 !* 0.1 declarations of arguments
1745 !
1746 REAL, INTENT(IN) :: ptstep, psnowdzmin
1747 !
1748  CHARACTER(LEN=*), INTENT(IN) :: hsnowres ! type of sfc resistance
1749 ! DEFAULT=Louis (1979), standard ISBA
1750 ! method. Option to limit Ri number
1751 ! for very srtable conditions
1752 !
1753  CHARACTER(LEN=*), INTENT(IN) :: himplicit_wind ! wind implicitation option
1754 ! ! 'OLD' = direct
1755 ! ! 'NEW' = Taylor serie, order 1
1756 !
1757 REAL, DIMENSION(:), INTENT(IN) :: ppew_a_coef, ppew_b_coef, &
1758  ppet_a_coef, ppeq_a_coef, ppet_b_coef, &
1759  ppeq_b_coef
1760 ! PPEW_A_COEF = wind coefficient (m2s/kg)
1761 ! PPEW_B_COEF = wind coefficient (m/s)
1762 ! PPET_A_COEF = A-air temperature coefficient
1763 ! PPET_B_COEF = B-air temperature coefficient
1764 ! PPEQ_A_COEF = A-air specific humidity coefficient
1765 ! PPEQ_B_COEF = B-air specific humidity coefficient
1766 !
1767 REAL, DIMENSION(:), INTENT(IN) :: pzref, pts, psnowdz1, psnowdz2, &
1768  pradsink, psnowrho, psnowliq, pscap, &
1769  pscond1, pscond2, &
1770  pz0, phpsnow, &
1771  palbt, pz0eff, pz0h
1772 !
1773 REAL, DIMENSION(:), INTENT(IN) :: psw_rad, plw_rad, pta, pqa, pps, prhoa
1774 !
1775 REAL, DIMENSION(:), INTENT(IN) :: puref, pexns, pexna, pdircoszw, pvmod
1776 !
1777 REAL, DIMENSION(:), INTENT(OUT) :: ptsterm1, ptsterm2, pemist, pra, &
1778  pct, psfcfrz, pcdsnow, pchsnow, &
1779  pqsat, pdqsat, prsra
1780 !
1781 REAL, DIMENSION(:), INTENT(OUT) :: pustar2_ic, &
1782  ppet_a_coef_t, ppeq_a_coef_t, &
1783  ppet_b_coef_t, ppeq_b_coef_t
1784 !
1785 REAL, DIMENSION(:), INTENT(OUT) :: pri
1786 !
1787 !* 0.2 declarations of local variables
1788 !
1789 REAL, DIMENSION(SIZE(PTS)) :: zac, zri, zcond1, zcond2, &
1790  zsconda, za, zb, zc, &
1791  zcdn, zsnowdzm1, zsnowdzm2, &
1792  zvmod, zustar2, zts3, zlvt, &
1793  z_ccoef
1794 REAL(KIND=JPRB) :: zhook_handle
1795 !
1796 !-------------------------------------------------------------------------------
1797 !
1798 ! 1. New saturated specific humidity and derrivative:
1799 ! ---------------------------------------------------
1800 !
1801 IF (lhook) CALL dr_hook('SNOW3LEBUD',0,zhook_handle)
1802 !
1803 zri(:) = xundef
1804 !
1805 pqsat(:) = qsati(pts(:),pps(:))
1806 pdqsat(:) = dqsati(pts(:),pps(:),pqsat(:))
1807 !
1808 !
1809 ! 2. Surface properties:
1810 ! ----------------------
1811 ! It might be of interest to use snow-specific roughness
1812 ! or a temperature dependence on emissivity:
1813 ! but for now, use ISBA defaults.
1814 !
1815 pemist(:) = xemissn
1816 !
1817 ! 2. Computation of resistance and drag coefficient
1818 ! -------------------------------------------------
1819 !
1820  CALL surface_ri(pts, pqsat, pexns, pexna, pta, pqa, &
1821  pzref, puref, pdircoszw, pvmod, zri )
1822 !
1823 ! Simple adaptation of method by Martin and Lejeune (1998)
1824 ! to apply a lower limit to turbulent transfer coef
1825 ! by defining a maximum Richarson number for stable
1826 ! conditions:
1827 !
1828 IF(hsnowres=='RIL') THEN
1829  DO jj = 1,SIZE(zri)
1830  zri(jj) = min(x_ri_max,zri(jj))
1831  ENDDO
1832 ENDIF
1833 !
1834 pri(:)=zri(:)
1835 !
1836 ! Surface aerodynamic resistance for heat transfers
1837 !
1838  CALL surface_aero_cond(zri, pzref, puref, pvmod, pz0, pz0h, zac, pra, pchsnow)
1839 !
1840 ! For atmospheric model coupling:
1841 !
1842  CALL surface_cd(zri, pzref, puref, pz0eff, pz0h, pcdsnow, zcdn)
1843 !
1844 prsra(:) = prhoa(:) / pra(:)
1845 !
1846 !
1847 ! Modify flux-form implicit coupling coefficients:
1848 ! - wind components:
1849 !
1850 IF(himplicit_wind=='OLD')THEN
1851 ! old implicitation
1852  zustar2(:) = ( pcdsnow(:)*pvmod(:)*ppew_b_coef(:)) / &
1853  (1.0-prhoa(:)*pcdsnow(:)*pvmod(:)*ppew_a_coef(:))
1854 ELSE
1855 ! new implicitation
1856  zustar2(:) = (pcdsnow(:)*pvmod(:)*(2.*ppew_b_coef(:)-pvmod(:))) &
1857  / (1.0-2.0*prhoa(:)*pcdsnow(:)*pvmod(:)*ppew_a_coef(:))
1858 ENDIF
1859 !
1860 zvmod(:) = prhoa(:)*ppew_a_coef(:)*zustar2(:) + ppew_b_coef(:)
1861 zvmod(:) = max(zvmod(:),0.)
1862 !
1863 WHERE(ppew_a_coef(:)/= 0.)
1864  zustar2(:) = max( ( zvmod(:) - ppew_b_coef(:) ) / (prhoa(:)*ppew_a_coef(:)), 0.)
1865 ENDWHERE
1866 !
1867 ! implicit wind friction
1868 zustar2(:) = max(zustar2(:),0.)
1869 !
1870 pustar2_ic(:) = zustar2(:)
1871 !
1872 ! 3. Calculate linearized surface energy budget components:
1873 ! ---------------------------------------------------------
1874 ! To prevent numerical difficulties for very thin snow
1875 ! layers, limit the grid "thinness": this is important as
1876 ! layers become vanishing thin:
1877 !
1878 zsnowdzm1(:) = max(psnowdz1(:), psnowdzmin)
1879 zsnowdzm2(:) = max(psnowdz2(:), psnowdzmin)
1880 !
1881 ! Surface thermal inertia:
1882 !
1883 pct(:) = 1.0/(pscap(:)*zsnowdzm1(:))
1884 !
1885 ! Fraction of surface frozen (sublimation) with the remaining
1886 ! fraction being liquid (evaporation):
1887 ! NOTE: currently, for simplicity, assume no liquid water flux
1888 ! OFF: PSFCFRZ(:) = 1.0 - PSNOWLIQ(:)*XRHOLW/(ZSNOWDZM1(:)*PSNOWRHO(:))
1889 !
1890 psfcfrz(:) = 1.0
1891 !
1892 ! Thermal conductivity between uppermost and lower snow layers:
1893 !
1894 zcond1(:) = zsnowdzm1(:)/((zsnowdzm1(:)+zsnowdzm2(:))*pscond1(:))
1895 zcond2(:) = zsnowdzm2(:)/((zsnowdzm1(:)+zsnowdzm2(:))*pscond2(:))
1896 !
1897 zsconda(:) = 1.0/(zcond1(:)+zcond2(:))
1898 !
1899 ! Transform implicit coupling coefficients:
1900 ! Note, surface humidity is 100% over snow.
1901 !
1902 ! - specific humidity:
1903 !
1904 z_ccoef(:) = 1.0 - ppeq_a_coef(:)*prsra(:)
1905 !
1906 ppeq_a_coef_t(:) = - ppeq_a_coef(:)*prsra(:)*pdqsat(:)/z_ccoef(:)
1907 !
1908 ppeq_b_coef_t(:) = ( ppeq_b_coef(:) - ppeq_a_coef(:)*prsra(:)*(pqsat(:) - &
1909  pdqsat(:)*pts(:)) )/z_ccoef(:)
1910 !
1911 ! - air temperature:
1912 ! (assumes A and B correspond to potential T):
1913 !
1914 z_ccoef(:) = (1.0 - ppet_a_coef(:)*prsra(:))/pexna(:)
1915 !
1916 ppet_a_coef_t(:) = - ppet_a_coef(:)*prsra(:)/(pexns(:)*z_ccoef(:))
1917 !
1918 ppet_b_coef_t(:) = ppet_b_coef(:)/z_ccoef(:)
1919 !
1920 !
1921 ! Energy budget solution terms:
1922 !
1923 zts3(:) = pemist(:) * xstefan * pts(:)**3
1924 zlvt(:) = (1.-psfcfrz(:))*xlvtt + psfcfrz(:)*xlstt
1925 !
1926 za(:) = 1. / ptstep + pct(:) * (4. * zts3(:) + &
1927  prsra(:) * zlvt(:) * (pdqsat(:) - ppeq_a_coef_t(:)) &
1928  + prsra(:) * xcpd * ( (1./pexns(:))-(ppet_a_coef_t(:)/pexna(:)) ) &
1929  + (2.*zsconda(:)/(zsnowdzm2(:)+zsnowdzm1(:))) )
1930 !
1931 zb(:) = 1. / ptstep + pct(:) * (3. * zts3(:) + &
1932  prsra(:) * pdqsat(:) * zlvt(:) )
1933 !
1934 zc(:) = pct(:) * (prsra(:) * xcpd * ppet_b_coef_t(:)/pexna(:) + psw_rad(:) * &
1935  (1. - palbt(:)) + pemist(:)*plw_rad(:) - prsra(:) * &
1936  zlvt(:) * (pqsat(:)-ppeq_b_coef_t(:)) &
1937  + phpsnow(:) + pradsink(:) )
1938 !
1939 !
1940 ! Coefficients needed for implicit solution
1941 ! of linearized surface energy budget:
1942 !
1943 ptsterm2(:) = 2.*zsconda(:)*pct(:)/(za(:)*(zsnowdzm2(:)+zsnowdzm1(:)))
1944 !
1945 ptsterm1(:) = (pts(:)*zb(:) + zc(:))/za(:)
1946 IF (lhook) CALL dr_hook('SNOW3LEBUD',1,zhook_handle)
1947 !
1948 !-------------------------------------------------------------------------------
1949 !
1950 END SUBROUTINE snow3lebud
1951 !####################################################################
1952 !####################################################################
1953 !####################################################################
1954  SUBROUTINE snow3lsolvt(OMEB,PTSTEP,PSNOWDZMIN, &
1955  psnowdz,pscond,pscap,ptg, &
1956  psoilcond,pd_g, &
1957  pradsink,pct,pterm1,pterm2, &
1958  ppet_a_coef_t,ppeq_a_coef_t, &
1959  ppet_b_coef_t,ppeq_b_coef_t, &
1960  pta_ic, pqa_ic, &
1961  pgrndflux,pgrndfluxo,psnowtemp, &
1962  psnowflux )
1963 !
1964 !! PURPOSE
1965 !! -------
1966 ! This subroutine solves the 1-d diffusion of 'ZSNOWTEMP' using a
1967 ! layer averaged set of equations which are time differenced
1968 ! using the backward-difference scheme (implicit).
1969 ! The eqs are solved rapidly by taking advantage of the
1970 ! fact that the matrix is tridiagonal. This is a very
1971 ! general routine and can be used for the 1-d diffusion of any
1972 ! quantity as long as the diffusity is not a function of the
1973 ! quantity being diffused. Aaron Boone 8-98. Soln to the eq:
1974 !
1975 ! c dQ d k dQ dS
1976 ! -- = -- -- - --
1977 ! dt dx dx dx
1978 !
1979 ! where k = k(x) (thermal conductivity), c = c(x) (heat capacity)
1980 ! as an eg. for temperature/heat eq. S is a sink (-source) term.
1981 ! Diffusivity is k/c
1982 !
1983 !
1984 USE modd_csts,ONLY : xtt
1985 !
1986 USE modi_tridiag_ground
1987 !
1988 IMPLICIT NONE
1989 !
1990 !* 0.1 declarations of arguments
1991 !
1992 LOGICAL, INTENT(IN) :: omeb
1993 !
1994 REAL, INTENT(IN) :: ptstep, psnowdzmin
1995 !
1996 REAL, DIMENSION(:), INTENT(IN) :: ptg, psoilcond, pd_g, &
1997  pct, pterm1, pterm2
1998 
1999 !
2000 REAL, DIMENSION(:,:), INTENT(IN) :: psnowdz, pscond, pscap, &
2001  pradsink
2002 !
2003 REAL, DIMENSION(:), INTENT(IN) :: ppet_a_coef_t, ppeq_a_coef_t, &
2004  ppet_b_coef_t, ppeq_b_coef_t
2005 !
2006 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowtemp
2007 !
2008 REAL, DIMENSION(:), INTENT(OUT) :: pgrndflux, pgrndfluxo, psnowflux, &
2009  pta_ic, pqa_ic
2010 !
2011 !
2012 !* 0.2 declarations of local variables
2013 !
2014 !
2015 INTEGER :: jj, ji
2016 !
2017 INTEGER :: ini
2018 INTEGER :: inlvls
2019 !
2020 REAL, DIMENSION(SIZE(PTG)) :: zsnowtemp_delta
2021 !
2022 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)) :: zsnowtemp, zdterm, zcterm, &
2023  zfrcv, zamtrx, zbmtrx, &
2024  zcmtrx
2025 !
2026 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)) :: zwork1, zwork2, zdzdif, &
2027  zsnowdzm
2028 !
2029 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)-1) :: zsnowtemp_m, &
2030  zfrcv_m, zamtrx_m, &
2031  zbmtrx_m, zcmtrx_m
2032 REAL(KIND=JPRB) :: zhook_handle
2033 !-------------------------------------------------------------------------------
2034 !
2035 ! 0. Initialize:
2036 ! ------------------
2037 !
2038 IF (lhook) CALL dr_hook('SNOW3LSOLVT',0,zhook_handle)
2039 zsnowtemp(:,:) = psnowtemp(:,:)
2040 ini = SIZE(psnowdz(:,:),1)
2041 inlvls = SIZE(psnowdz(:,:),2)
2042 !
2043 !
2044 ! 1. Calculate tri-diagnoal matrix coefficients:
2045 ! ----------------------------------------------
2046 ! For heat transfer, assume a minimum grid
2047 ! thickness (to prevent numerical
2048 ! problems for very thin snow cover):
2049 !
2050 zsnowdzm(:,:) = max(psnowdz(:,:), psnowdzmin)
2051 !
2052 DO jj=1,inlvls-1
2053  DO ji=1,ini
2054  zdzdif(ji,jj) = 0.5*(zsnowdzm(ji,jj)+zsnowdzm(ji,jj+1))
2055  zwork1(ji,jj) = zsnowdzm(ji,jj )/(2.0*zdzdif(ji,jj)*pscond(ji,jj ))
2056  zwork2(ji,jj) = zsnowdzm(ji,jj+1)/(2.0*zdzdif(ji,jj)*pscond(ji,jj+1))
2057  ENDDO
2058 ENDDO
2059 !
2060 zdzdif(:,inlvls) = 0.5*(zsnowdzm(:,inlvls)+pd_g(:))
2061 zwork1(:,inlvls) = zsnowdzm(:,inlvls)/(2.0*zdzdif(:,inlvls)*pscond(:,inlvls))
2062 zwork2(:,inlvls) = pd_g(: )/(2.0*zdzdif(:,inlvls)*psoilcond(: ))
2063 !
2064 zdterm(:,:) = 1.0/(zdzdif(:,:)*(zwork1(:,:)+zwork2(:,:)))
2065 !
2066 zcterm(:,:) = pscap(:,:)*zsnowdzm(:,:)/ptstep
2067 !
2068 ! 2. Set up tri-diagonal matrix
2069 ! -----------------------------
2070 !
2071 ! Upper BC
2072 !
2073 zamtrx(:,1) = 0.0
2074 zbmtrx(:,1) = 1./(pct(:)*ptstep)
2075 zcmtrx(:,1) = -pterm2(:)*zbmtrx(:,1)
2076 zfrcv(:,1) = pterm1(:)*zbmtrx(:,1)
2077 !
2078 !
2079 ! Interior Grid
2080 !
2081 DO jj=2,inlvls-1
2082  DO ji=1,ini
2083  zamtrx(ji,jj) = -zdterm(ji,jj-1)
2084  zbmtrx(ji,jj) = zcterm(ji,jj) + zdterm(ji,jj-1) + zdterm(ji,jj)
2085  zcmtrx(ji,jj) = -zdterm(ji,jj)
2086  zfrcv(ji,jj) = zcterm(ji,jj)*psnowtemp(ji,jj) - (pradsink(ji,jj-1)-pradsink(ji,jj))
2087  ENDDO
2088 ENDDO
2089 !
2090 !Lower BC
2091 !
2092 zamtrx(:,inlvls) = -zdterm(:,inlvls-1)
2093 zbmtrx(:,inlvls) = zcterm(:,inlvls) + zdterm(:,inlvls-1) + &
2094  zdterm(:,inlvls)
2095 zcmtrx(:,inlvls) = 0.0
2096 zfrcv(:,inlvls) = zcterm(:,inlvls)*psnowtemp(:,inlvls) + &
2097  zdterm(:,inlvls)*ptg(:) &
2098  - (pradsink(:,inlvls-1)-pradsink(:,inlvls))
2099 !
2100 ! - - -------------------------------------------------
2101 !
2102 ! 4. Compute solution vector
2103 ! --------------------------
2104 !
2105  CALL tridiag_ground(zamtrx,zbmtrx,zcmtrx,zfrcv,zsnowtemp)
2106 !
2107 ! Heat flux between surface and 2nd snow layers: (W/m2)
2108 !
2109 psnowflux(:) = zdterm(:,1)*(zsnowtemp(:,1) - zsnowtemp(:,2))
2110 !
2111 !
2112 ! 5. Snow melt case
2113 ! -----------------
2114 ! If melting in uppermost layer, assume surface layer
2115 ! temperature at freezing point and re-evaluate lower
2116 ! snowpack temperatures. This is done as it is most likely
2117 ! most signigant melting will occur within a time step in surface layer.
2118 ! Surface energy budget (and fluxes) will
2119 ! be re-calculated (outside of this routine).
2120 !
2121 ! NOTE: if MEB is active, then surface fluxes have been defined outside
2122 ! of the snow routine and have been adjusted such that they are evaluated
2123 ! at a snow surface temperature no greater than Tf. Thus, the implicit surface temperature
2124 ! will likely never greatly exceed Tf (before melt computed and they are adjusted to Tf)
2125 ! so we can skip the next block of code when MEB is active.
2126 !
2127 IF(.NOT.omeb)THEN
2128 !
2129  zamtrx_m(:,1) = 0.0
2130  zbmtrx_m(:,1) = zcterm(:,2) + zdterm(:,1) + zdterm(:,2)
2131  zcmtrx_m(:,1) = -zdterm(:,2)
2132  zfrcv_m(:,1) = zcterm(:,2)*psnowtemp(:,2) + xtt*zdterm(:,1) - &
2133  (pradsink(:,1)-pradsink(:,2))
2134 !
2135  DO jj=2,inlvls-1
2136  DO ji=1,ini
2137  zamtrx_m(ji,jj) = zamtrx(ji,jj+1)
2138  zbmtrx_m(ji,jj) = zbmtrx(ji,jj+1)
2139  zcmtrx_m(ji,jj) = zcmtrx(ji,jj+1)
2140  zfrcv_m(ji,jj) = zfrcv(ji,jj+1)
2141  zsnowtemp_m(ji,jj) = psnowtemp(ji,jj+1)
2142  ENDDO
2143  ENDDO
2144 !
2145  CALL tridiag_ground(zamtrx_m,zbmtrx_m,zcmtrx_m,zfrcv_m,zsnowtemp_m)
2146 !
2147 ! If melting for 2 consecuative time steps, then replace current T-profile
2148 ! with one assuming T=Tf in surface layer:
2149 !
2150  zsnowtemp_delta(:) = 0.0
2151 !
2152  WHERE(zsnowtemp(:,1) > xtt .AND. psnowtemp(:,1) == xtt)
2153  psnowflux(:) = zdterm(:,1)*(xtt - zsnowtemp_m(:,1))
2154  zsnowtemp_delta(:) = 1.0
2155  END WHERE
2156 !
2157  DO jj=2,inlvls
2158  DO ji=1,ini
2159  zsnowtemp(ji,jj) = zsnowtemp_delta(ji)*zsnowtemp_m(ji,jj-1) &
2160  + (1.0-zsnowtemp_delta(ji))*zsnowtemp(ji,jj)
2161  ENDDO
2162  ENDDO
2163 !
2164 ENDIF
2165 !
2166 !
2167 ! 6. Lower boundary flux:
2168 ! -----------------------
2169 ! NOTE: evaluate this term assuming the snow layer
2170 ! can't exceed the freezing point as this adjustment
2171 ! is made in melting routine. Then must adjust temperature
2172 ! to conserve energy:
2173 !
2174 pgrndfluxo(:) = zdterm(:,inlvls)*(zsnowtemp(:,inlvls) -ptg(:))
2175 pgrndflux(:) = zdterm(:,inlvls)*(min(xtt,zsnowtemp(:,inlvls))-ptg(:))
2176 !
2177 zsnowtemp(:,inlvls) = zsnowtemp(:,inlvls) + (pgrndfluxo(:)-pgrndflux(:))/zcterm(:,inlvls)
2178 !
2179 ! 7. Update temperatute profile in time:
2180 ! --------------------------------------
2181 !
2182 psnowtemp(:,:) = zsnowtemp(:,:)
2183 !
2184 !
2185 ! 8. Compute new (implicit) air T and specific humidity
2186 ! -----------------------------------------------------
2187 !
2188 IF(.NOT.omeb)THEN
2189 !
2190  pta_ic(:) = ppet_b_coef_t(:) + ppet_a_coef_t(:)* psnowtemp(:,1)
2191 !
2192  pqa_ic(:) = ppeq_b_coef_t(:) + ppeq_a_coef_t(:)* psnowtemp(:,1)
2193 !
2194 ENDIF
2195 !
2196 IF (lhook) CALL dr_hook('SNOW3LSOLVT',1,zhook_handle)
2197 !
2198 !
2199 END SUBROUTINE snow3lsolvt
2200 !####################################################################
2201 !####################################################################
2202 !####################################################################
2203  SUBROUTINE snow3lmelt(PTSTEP,PSCAP,PSNOWTEMP,PSNOWDZ, &
2204  psnowrho,psnowliq,pmeltxs )
2205 !
2206 !
2207 !! PURPOSE
2208 !! -------
2209 ! Calculate snow melt (resulting from surface fluxes, ground fluxes,
2210 ! or internal shortwave radiation absorbtion). It is used to
2211 ! augment liquid water content, maintain temperatures
2212 ! at or below freezing, and possibly reduce the mass
2213 ! or compact the layer(s).
2214 !
2215 !
2216 USE modd_csts,ONLY : xtt, xlmtt, xrholw, xrholi
2217 !
2218 USE mode_snow3l
2219 !
2220 IMPLICIT NONE
2221 !
2222 !* 0.1 declarations of arguments
2223 !
2224 REAL, INTENT(IN) :: ptstep
2225 !
2226 REAL, DIMENSION(:,:), INTENT(IN) :: pscap
2227 !
2228 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowdz, psnowtemp, psnowrho, &
2229  psnowliq
2230 !
2231 REAL, DIMENSION(:), INTENT(OUT) :: pmeltxs
2232 !
2233 !
2234 !* 0.2 declarations of local variables
2235 !
2236 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zphase, zcmprsfact, &
2237  zsnowlwe, zwholdmax, &
2238  zsnowmelt, zsnowtemp, &
2239  zmeltxs
2240 !
2241 INTEGER :: jwrk, ji ! loop counter
2242 REAL(KIND=JPRB) :: zhook_handle
2243 !-------------------------------------------------------------------------------
2244 !
2245 ! 0. Initialize:
2246 ! ---------------------------
2247 !
2248 IF (lhook) CALL dr_hook('SNOW3LMELT',0,zhook_handle)
2249 zphase(:,:) = 0.0
2250 zcmprsfact(:,:) = 0.0
2251 zsnowlwe(:,:) = 0.0
2252 zwholdmax(:,:) = 0.0
2253 zsnowmelt(:,:) = 0.0
2254 zsnowtemp(:,:) = 0.0
2255 zmeltxs(:,:) = 0.0
2256 !
2257 ! 1. Determine amount of melt in each layer:
2258 ! ------------------------------------------
2259 !
2260 WHERE(psnowdz > 0.0)
2261 !
2262 ! Total Liquid equivalent water content of snow (m):
2263 !
2264  zsnowlwe(:,:) = psnowrho(:,:)*psnowdz(:,:)/xrholw
2265 !
2266 ! Melt snow if excess energy and snow available:
2267 ! Phase change (J/m2)
2268 !
2269  zphase(:,:) = min(pscap(:,:)*max(0.0, psnowtemp(:,:) - xtt)* &
2270  psnowdz(:,:), &
2271  max(0.0,zsnowlwe(:,:)-psnowliq(:,:))*xlmtt*xrholw)
2272 !
2273 !
2274 ! Update snow liq water content and temperature if melting:
2275 ! liquid water available for next layer from melting of snow
2276 ! which is assumed to be leaving the current layer (m):
2277 !
2278  zsnowmelt(:,:) = zphase(:,:)/(xlmtt*xrholw)
2279 !
2280 ! Cool off snow layer temperature due to melt:
2281 !
2282  zsnowtemp(:,:) = psnowtemp(:,:) - zphase(:,:)/(pscap(:,:)*psnowdz(:,:))
2283 !
2284  psnowtemp(:,:) = min(xtt, zsnowtemp(:,:))
2285 !
2286  zmeltxs(:,:) = (zsnowtemp(:,:)-psnowtemp(:,:))*pscap(:,:)*psnowdz(:,:)
2287 !
2288 END WHERE
2289 !
2290 ! Loss of snowpack depth: (m) and liquid equiv (m):
2291 ! Compression factor for melt loss: this decreases
2292 ! layer thickness and increases density thereby leaving
2293 ! total SWE constant. NOTE: All melt water
2294 ! in excess of the holding capacity is NOT used
2295 ! for compression, rather it decreases the layer
2296 ! thickness ONLY, causing a loss of SWE (outside
2297 ! of this routine).
2298 !
2299 zwholdmax(:,:) = snow3lhold(psnowrho,psnowdz)
2300 !
2301 WHERE(psnowdz > 0.0)
2302 !
2303  zcmprsfact(:,:) = (zsnowlwe(:,:)-min(psnowliq(:,:)+zsnowmelt(:,:), &
2304  zwholdmax(:,:)))/ &
2305  (zsnowlwe(:,:)-min(psnowliq(:,:),zwholdmax(:,:)))
2306 !
2307  psnowdz(:,:) = psnowdz(:,:)*zcmprsfact(:,:)
2308  psnowrho(:,:) = zsnowlwe(:,:)*xrholw/psnowdz(:,:)
2309 !
2310 ! Make sure maximum density is not surpassed! If it is, lower the density
2311 ! and increase the snow thickness accordingly:
2312 
2313  zcmprsfact(:,:) = max(xrholi, psnowrho(:,:))/xrholi
2314  psnowdz(:,:) = psnowdz(:,:)*zcmprsfact(:,:)
2315  psnowrho(:,:) = zsnowlwe(:,:)*xrholw/psnowdz(:,:)
2316 !
2317 !
2318 ! 2. Add snow melt to current snow liquid water content:
2319 ! ------------------------------------------------------
2320 !
2321  psnowliq(:,:) = psnowliq(:,:) + zsnowmelt(:,:)
2322 !
2323 END WHERE
2324 !
2325 ! 3. Excess heat from melting
2326 ! ---------------------------
2327 ! use it to warm underlying ground/vegetation layer to conserve energy
2328 !
2329 pmeltxs(:) = 0.
2330 DO jwrk = 1, SIZE(zmeltxs,2)
2331  DO ji = 1, SIZE(zmeltxs,1)
2332  pmeltxs(ji) = pmeltxs(ji) + zmeltxs(ji,jwrk)
2333  ENDDO
2334 ENDDO
2335 pmeltxs(:) = pmeltxs(:) / ptstep ! (W/m2)
2336 !
2337 IF (lhook) CALL dr_hook('SNOW3LMELT',1,zhook_handle)
2338 !
2339 !
2340 !
2341 END SUBROUTINE snow3lmelt
2342 !####################################################################
2343 !####################################################################
2344 !####################################################################
2345  SUBROUTINE snow3lrefrz(PTSTEP,PRR, &
2346  psnowrho,psnowtemp,psnowdz,psnowliq, &
2347  pthrufal )
2348 !
2349 !
2350 !! PURPOSE
2351 !! -------
2352 ! Calculate any freezing/refreezing of liquid water in the snowpack.
2353 ! Also, calculate liquid water transmission and snow runoff.
2354 ! Water flow causes layer thickness reduction and can cause
2355 ! rapid densification of a layer.
2356 !
2357 !
2358 USE modd_csts, ONLY : xtt, xlmtt, xrholw
2359 USE modd_snow_par, ONLY : xsnowdmin
2360 !
2361 USE mode_snow3l
2362 !
2363 IMPLICIT NONE
2364 !
2365 !* 0.1 declarations of arguments
2366 !
2367 REAL, INTENT(IN) :: ptstep
2368 !
2369 REAL, DIMENSION(:), INTENT(IN) :: prr
2370 !
2371 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowdz, psnowtemp, psnowliq, psnowrho
2372 !
2373 REAL, DIMENSION(:), INTENT(INOUT) :: pthrufal
2374 !
2375 !
2376 !
2377 !* 0.2 declarations of local variables
2378 !
2379 INTEGER :: jj, ji
2380 !
2381 INTEGER :: ini
2382 INTEGER :: inlvls
2383 !
2384 REAL, DIMENSION(SIZE(PRR)) :: zpcpxs, ztotwcap, zrainfall
2385 !
2386 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zflowliq, zwork, &
2387  zsnowliq, zsnowrho, &
2388  zwholdmax, zsnowdz, &
2389  zsnowtemp, zscap, &
2390  zsnowheat
2391 !
2392 REAL, DIMENSION(SIZE(PSNOWRHO,1),0:SIZE(PSNOWRHO,2)):: zflowliqt
2393 !
2394 REAL(KIND=JPRB) :: zhook_handle
2395 !
2396 !-------------------------------------------------------------------------------
2397 !
2398 ! 0. Initialize:
2399 ! --------------
2400 !
2401 IF (lhook) CALL dr_hook('SNOW3LREFRZ',0,zhook_handle)
2402 !
2403 zsnowrho(:,:) = psnowrho(:,:)
2404 zsnowliq(:,:) = psnowliq(:,:)
2405 zsnowtemp(:,:) = psnowtemp(:,:)
2406 ini = SIZE(psnowdz(:,:),1)
2407 inlvls = SIZE(psnowdz(:,:),2)
2408 !
2409 !
2410 ! 1. Refreeze due to heat transfer
2411 ! -----------------------------
2412 ! Freeze liquid water in any layers which cooled due
2413 ! to heat transfer. First, update H and re-diagnose (update) T and Wl:
2414 !
2415 zscap(:,:) = snow3lscap(zsnowrho)
2416 !
2417 zsnowheat(:,:) = psnowdz(:,:)*( zscap(:,:)*(zsnowtemp(:,:)-xtt) &
2418  - xlmtt*zsnowrho(:,:) ) + xlmtt*xrholw*zsnowliq(:,:)
2419 !
2420 zsnowtemp(:,:) = xtt + ( ((zsnowheat(:,:)/max(psnowdz(:,:),xsnowdmin/inlvls)) &
2421  + xlmtt*zsnowrho(:,:))/zscap(:,:) )
2422 !
2423 zsnowliq(:,:) = max(0.0,zsnowtemp(:,:)-xtt)*zscap(:,:)*psnowdz(:,:)/(xlmtt*xrholw)
2424 !
2425 zsnowtemp(:,:) = min(xtt,zsnowtemp(:,:))
2426 !
2427 !
2428 ! 2. Reduce thickness due to snowmelt in excess of holding capacity
2429 ! --------------------------------------------------------------
2430 ! Any water in excess of the
2431 ! Maximum holding space for liquid water
2432 ! amount is drained into next layer down.
2433 ! Loss of water due to snowmelt causes a reduction
2434 ! in snow layer mass by a reduction in thickness. Owing to a consistent
2435 ! decrease in both liq and thickness (and the fact that T=Tf when liquid present),
2436 ! enthalpy is conserved.
2437 !
2438 zwholdmax(:,:) = snow3lhold(psnowrho,psnowdz)
2439 !
2440 zflowliq(:,:) = max(0.,zsnowliq(:,:)-zwholdmax(:,:))
2441 !
2442 zsnowliq(:,:) = zsnowliq(:,:) - zflowliq(:,:)
2443 !
2444 zsnowdz(:,:) = psnowdz(:,:) - zflowliq(:,:)*xrholw/zsnowrho(:,:)
2445 !
2446 zsnowdz(:,:) = max(0.0, zsnowdz(:,:)) ! to prevent possible very small
2447 ! negative values (machine prescision
2448 ! as snow vanishes
2449 !
2450 !
2451 ! 3. Liquid water flow: liquid precipitation and meltwater
2452 ! -----------------------------------------------------
2453 !
2454 ! Rainfall flowing into uppermost snow layer:
2455 ! If rainfall is excessive enough (or layers thin enough)
2456 ! it is simply routed directly to runoff: First calculate
2457 ! the total snow pack available liquid water holding capacity:
2458 !
2459 ztotwcap(:) = 0.
2460 DO jj=1,inlvls
2461  DO ji=1,ini
2462  ztotwcap(ji) = ztotwcap(ji) + zwholdmax(ji,jj)
2463  ENDDO
2464 ENDDO
2465 !
2466 ! Rain entering snow (m):
2467 !
2468 zrainfall(:) = prr(:)*ptstep/xrholw ! rainfall (m)
2469 !
2470 zflowliqt(:,0)= min(zrainfall(:),ztotwcap(:))
2471 !
2472 ! Rain assumed to directly pass through the pack to runoff (m):
2473 !
2474 zpcpxs(:) = zrainfall(:) - zflowliqt(:,0)
2475 !
2476 DO jj=1,inlvls
2477  DO ji=1,ini
2478  zflowliqt(ji,jj) = zflowliq(ji,jj)
2479  ENDDO
2480 ENDDO
2481 !
2482 !
2483 ! Thickness is maintained during water through-flow,
2484 ! so that mass transfer is represented by
2485 ! density changes: NOTE a maximum density
2486 ! is assumed (XRHOSMAX_ES) so that all flow
2487 ! which would result in densities exceeding
2488 ! this limit are routed to next layer down.
2489 ! First test for saturation, then
2490 ! rout excess water down to next layer down
2491 ! and repeat calculation. Net gain in liquid (mass) is
2492 ! translated into a density increase:
2493 !
2494 zflowliq(:,:) = 0.0 ! clear this array for work
2495 psnowliq(:,:) = zsnowliq(:,:) ! reset liquid water content
2496 !
2497 DO jj=1,inlvls
2498  DO ji=1,ini
2499  zsnowliq(ji,jj) = zsnowliq(ji,jj) + zflowliqt(ji,jj-1)
2500  zflowliq(ji,jj) = max(0.0, zsnowliq(ji,jj)-zwholdmax(ji,jj))
2501  zsnowliq(ji,jj) = zsnowliq(ji,jj) - zflowliq(ji,jj)
2502  zflowliqt(ji,jj) = zflowliqt(ji,jj) + zflowliq(ji,jj)
2503  ENDDO
2504 ENDDO
2505 !
2506 zwork(:,:) = max(xsnowdmin/inlvls,zsnowdz(:,:))
2507 zsnowrho(:,:) = zsnowrho(:,:)+(zsnowliq(:,:)-psnowliq(:,:))*xrholw/zwork(:,:)
2508 zscap(:,:) = snow3lscap(zsnowrho(:,:))
2509 zsnowtemp(:,:) = xtt +(((zsnowheat(:,:)/zwork(:,:))+xlmtt*zsnowrho(:,:))/zscap(:,:))
2510 zsnowliq(:,:) = max(0.0,zsnowtemp(:,:)-xtt)*zscap(:,:)*zsnowdz(:,:)/(xlmtt*xrholw)
2511 zsnowtemp(:,:) = min(xtt,zsnowtemp(:,:))
2512 !
2513 ! Any remaining throughflow after freezing is available to
2514 ! the soil for infiltration or surface runoff (m).
2515 ! I.E. This is the amount of water leaving the snowpack:
2516 ! Rate water leaves the snowpack [kg/(m2 s)]:
2517 !
2518 pthrufal(:) = pthrufal(:) + zflowliqt(:,inlvls)
2519 !
2520 ! Add excess rain (rain which flowed directly through the snow
2521 ! due to saturation):
2522 !
2523 pthrufal(:) = (pthrufal(:) + zpcpxs(:))*xrholw/ptstep
2524 !
2525 ! 4. Update thickness and density and any freezing:
2526 ! ----------------------------------------------
2527 !
2528 psnowtemp(:,:)= zsnowtemp(:,:)
2529 psnowdz(:,:) = zsnowdz(:,:)
2530 psnowrho(:,:) = zsnowrho(:,:)
2531 psnowliq(:,:) = zsnowliq(:,:)
2532 !
2533 IF (lhook) CALL dr_hook('SNOW3LREFRZ',1,zhook_handle)
2534 !
2535 !-------------------------------------------------------------------------------
2536 !
2537 END SUBROUTINE snow3lrefrz
2538 !####################################################################
2539 !####################################################################
2540 !####################################################################
2541  SUBROUTINE snow3lflux(PSNOWTEMP,PSNOWDZ,PEXNS,PEXNA, &
2542  pustar2_ic, &
2543  ptstep,palbt,psw_rad,pemist,plwupsnow, &
2544  plw_rad,plwnetsnow, &
2545  pta,psfcfrz,pqa,phpsnow, &
2546  psnowtempo1,psnowflux,pct,pradsink, &
2547  pqsat,pdqsat,prsra, &
2548  prn,ph,pgflux,ples3l,plel3l,pevap, &
2549  pustar,osfcmelt )
2550 !
2551 !
2552 !! PURPOSE
2553 !! -------
2554 ! Calculate the surface fluxes (atmospheric/surface).
2555 ! (Noilhan and Planton 1989; Noilhan and Mahfouf 1996)
2556 !
2557 !
2558 USE modd_csts,ONLY : xstefan, xcpd, xlstt, xlvtt, xtt
2559 !
2560 USE mode_thermos
2561 !
2562 IMPLICIT NONE
2563 !
2564 !* 0.1 declarations of arguments
2565 !
2566 REAL, INTENT(IN) :: ptstep
2567 !
2568 REAL, DIMENSION(:), INTENT(IN) :: psnowdz, psnowtempo1, psnowflux, pct, &
2569  pradsink, pexns, pexna
2570 !
2571 REAL, DIMENSION(:), INTENT(IN) :: palbt, psw_rad, pemist, plw_rad, &
2572  pta, psfcfrz, pqa, &
2573  phpsnow, pqsat, pdqsat, prsra, &
2574  pustar2_ic
2575 !
2576 REAL, DIMENSION(:), INTENT(INOUT) :: psnowtemp
2577 !
2578 REAL, DIMENSION(:), INTENT(OUT) :: prn, ph, pgflux, ples3l, plel3l, &
2579  pevap, plwupsnow, pustar, &
2580  plwnetsnow
2581 !
2582 LOGICAL, DIMENSION(:), INTENT(OUT) :: osfcmelt
2583 !
2584 !
2585 !* 0.2 declarations of local variables
2586 !
2587 REAL, DIMENSION(SIZE(PSNOWDZ)) :: zevapc, zle, zsnowtemp, zsmsnow, zgflux, &
2588  zdeltat, zsnowto3
2589 REAL(KIND=JPRB) :: zhook_handle
2590 !
2591 !-------------------------------------------------------------------------------
2592 !
2593 ! 0. Initialize:
2594 ! --------------
2595 !
2596 IF (lhook) CALL dr_hook('SNOW3LFLUX',0,zhook_handle)
2597 zsnowtemp(:) = psnowtemp(:)
2598 zle(:) = 0.0
2599 zsmsnow(:) = 0.0
2600 zgflux(:) = 0.0
2601 !
2602 osfcmelt(:) = .false.
2603 !
2604 zsnowto3(:) = psnowtempo1(:) ** 3 ! to save some CPU time, store this
2605 !
2606 ! 1. Flux calculations when melt not occuring at surface (W/m2):
2607 ! --------------------------------------------------------------
2608 !
2609 !
2610 zdeltat(:) = psnowtemp(:) - psnowtempo1(:) ! surface T time change:
2611 !
2612 plwupsnow(:) = pemist(:) * xstefan * zsnowto3(:)*( psnowtempo1(:) + 4.* zdeltat(:) )
2613 !
2614 plwnetsnow(:)= pemist(:) * plw_rad(:) - plwupsnow(:)
2615 !
2616 prn(:) = (1. - palbt(:)) * psw_rad(:) + plwnetsnow(:)
2617 !
2618 ph(:) = prsra(:) * xcpd * (psnowtemp(:)/pexns(:) - pta(:)/pexna(:))
2619 !
2620 zevapc(:) = prsra(:) * ( (pqsat(:) - pqa(:)) + pdqsat(:)*zdeltat(:) )
2621 !
2622 ples3l(:) = psfcfrz(:) * xlstt * zevapc(:)
2623 !
2624 plel3l(:) = (1.-psfcfrz(:))* xlvtt * zevapc(:)
2625 !
2626 zle(:) = ples3l(:) + plel3l(:)
2627 !
2628 pgflux(:) = prn(:) - ph(:) - zle(:) + phpsnow(:)
2629 !
2630 !
2631 ! 2. Initial melt adjustment
2632 ! --------------------------
2633 ! If energy avalabile to melt snow, then recalculate fluxes
2634 ! at the freezing point and add residual heat to layer
2635 ! average heat.
2636 !
2637 ! A) If temperature change is > 0 and passes freezing point this timestep,
2638 ! then recalculate fluxes at freezing point and excess energy
2639 ! will be used outside of this routine to change snow heat content:
2640 !
2641 WHERE (psnowtemp > xtt .AND. psnowtempo1 < xtt)
2642 !
2643  osfcmelt(:)= .true.
2644 !
2645  zdeltat(:) = xtt - psnowtempo1(:)
2646 !
2647  plwupsnow(:) = pemist(:) * xstefan * zsnowto3(:)*( psnowtempo1(:) + 4.* zdeltat(:) )
2648 !
2649  plwnetsnow(:)= pemist(:) * plw_rad(:) - plwupsnow(:)
2650 !
2651  prn(:) = (1. - palbt(:)) * psw_rad(:) + plwnetsnow(:)
2652 !
2653  ph(:) = prsra(:) * xcpd * (xtt/pexns(:) - pta(:)/pexna(:))
2654 !
2655  zevapc(:) = prsra(:) * ( (pqsat(:) - pqa(:)) + pdqsat(:)*zdeltat(:) )
2656 !
2657  ples3l(:) = psfcfrz(:) * xlstt * zevapc(:)
2658 !
2659  plel3l(:) = (1.-psfcfrz(:))* xlvtt * zevapc(:)
2660 !
2661  zle(:) = ples3l(:) + plel3l(:)
2662 !
2663  zgflux(:) = prn(:) - ph(:) - zle(:) + phpsnow(:)
2664 !
2665  zsmsnow(:) = pgflux(:) - zgflux(:)
2666 !
2667  pgflux(:) = zgflux(:)
2668 !
2669 ! This will be used to change heat content of snow:
2670 !
2671  zsnowtemp(:) = psnowtemp(:) - zsmsnow(:)*ptstep*pct(:)
2672 !
2673 END WHERE
2674 !
2675 ! 3. Ongoing melt adjustment: explicit solution
2676 ! ---------------------------------------------
2677 ! If temperature change is 0 and at freezing point this timestep,
2678 ! then recalculate fluxes and surface temperature *explicitly*
2679 ! as this is *exact* for snow at freezing point (Brun, Martin)
2680 !
2681 WHERE(psnowtemp(:) > xtt .AND. psnowtempo1(:) >= xtt)
2682 !
2683  osfcmelt(:) = .true.
2684 !
2685  plwupsnow(:) = pemist(:) * xstefan * (xtt ** 4)
2686 !
2687  plwnetsnow(:)= pemist(:) * plw_rad(:) - plwupsnow(:)
2688 !
2689  prn(:) = (1. - palbt(:)) * psw_rad(:) + plwnetsnow(:)
2690 !
2691  ph(:) = prsra(:) * xcpd * (xtt/pexns(:) - pta(:)/pexna(:))
2692 !
2693  zevapc(:) = prsra(:) * (pqsat(:) - pqa(:))
2694 !
2695  ples3l(:) = psfcfrz(:) * xlstt * zevapc(:)
2696 !
2697  plel3l(:) = (1.-psfcfrz(:))* xlvtt * zevapc(:)
2698 !
2699  zle(:) = ples3l(:) + plel3l(:)
2700 !
2701  pgflux(:) = prn(:) - ph(:) - zle(:) + phpsnow(:)
2702 !
2703  zsnowtemp(:) = xtt + ptstep*pct(:)*(pgflux(:) + pradsink(:) - psnowflux(:))
2704 !
2705 END WHERE
2706 !
2707 ! 4. Update surface temperature:
2708 ! ------------------------------
2709 !
2710 psnowtemp(:) = zsnowtemp(:)
2711 !
2712 ! 5. Final evaporative flux (kg/m2/s)
2713 !
2714 pevap(:) = zevapc(:)
2715 !
2716 ! 6. Friction velocity
2717 ! --------------------
2718 !
2719 pustar(:) = sqrt(pustar2_ic(:))
2720 !
2721 IF (lhook) CALL dr_hook('SNOW3LFLUX',1,zhook_handle)
2722 !
2723 !-------------------------------------------------------------------------------
2724 !
2725 END SUBROUTINE snow3lflux
2726 !####################################################################
2727 !####################################################################
2728 !####################################################################
2729  SUBROUTINE snow3levapn(PPSN3L,PLES3L,PLEL3L,PTSTEP,PSNOWTEMP, &
2730  psnowrho,psnowdz,psnowliq,pta, &
2731  plvtt,plstt,psnowheat,psoilcor )
2732 !
2733 !
2734 !! PURPOSE
2735 !! -------
2736 ! Remove mass from uppermost snow layer in response to
2737 ! evaporation (liquid) and sublimation.
2738 !
2739 !
2740 USE modd_csts, ONLY : xrholw, xlmtt, xci, xtt
2741 USE modd_snow_par, ONLY : xrhosmin_es, xsnowdmin
2742 !
2743 IMPLICIT NONE
2744 !
2745 !* 0.1 declarations of arguments
2746 !
2747 REAL, INTENT(IN) :: ptstep
2748 !
2749 REAL, DIMENSION(:), INTENT(IN) :: ppsn3l
2750 !
2751 REAL, DIMENSION(:), INTENT(IN) :: ples3l, plel3l ! (W/m2)
2752 !
2753 REAL, DIMENSION(:), INTENT(IN) :: pta, plvtt, plstt
2754 !
2755 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowheat, psnowdz
2756 !
2757 REAL, DIMENSION(:), INTENT(INOUT) :: psnowrho, psnowliq, &
2758  psnowtemp
2759 !
2760 REAL, DIMENSION(:), INTENT(OUT) :: psoilcor
2761 !
2762 !* 0.2 declarations of local variables
2763 !
2764 INTEGER :: ini, inlvls, jj, ji
2765 !
2766 REAL, DIMENSION(SIZE(PLES3L)) :: zsnowevaps, zsnowevap, zsnowevapx, &
2767  zsnowdz, zsnowheat, zscap, zsnowtemp
2768 !
2769 REAL, DIMENSION(SIZE(PLES3L)) :: zxse, zisnowd
2770 
2771 !* 0.3 declarations of local parameters
2772 
2773 REAL, PARAMETER :: zsnowdemin = 1.e-4 ! m
2774 REAL, PARAMETER :: ztdif = 15. ! K To prevent a possible
2775  ! decoupling of sfc snow T
2776  ! when vanishingly thin, impose
2777  ! this max T-diff based on obs...
2778  ! between Ta-Ts
2779 !
2780 REAL(KIND=JPRB) :: zhook_handle
2781 !
2782 !
2783 !-------------------------------------------------------------------------------
2784 !
2785 ! 0. Initialize:
2786 ! --------------
2787 !
2788 IF (lhook) CALL dr_hook('SNOW3LEVAPN',0,zhook_handle)
2789 !
2790 psoilcor(:) = 0.0
2791 !
2792 zsnowevaps(:) = 0.0
2793 zsnowevap(:) = 0.0
2794 zsnowevapx(:) = 0.0
2795 zsnowdz(:) = 0.0
2796 zscap(:) = 0.0
2797 zsnowheat(:) = 0.0
2798 zsnowtemp(:) = 0.0
2799 zxse(:) = 0.0
2800 !
2801 ini = SIZE(psnowdz(:,:),1)
2802 inlvls = SIZE(psnowdz(:,:),2)
2803 !
2804 !
2805 ! 1. Evaporation of snow liquid water
2806 ! -----------------------------------
2807 ! Evaporation reduces liq water equivalent content
2808 ! of snow pack either by reducing the snow density
2809 ! (evaporation) or the layer thickness (sublimation).
2810 ! Condensation does the reverse.
2811 !
2812 ! Multiply evaporation components by snow fraction
2813 ! to be consistent with fluxes from the snow covered
2814 ! fraction of grid box
2815 !
2816 WHERE(psnowdz(:,1) > 0.0)
2817 !
2818 ! Evaporation:
2819 ! Reduce density and liquid water content:
2820 !
2821  zsnowevap(:) = ppsn3l(:)*plel3l(:)*ptstep/(plvtt(:)*xrholw)
2822  zsnowevapx(:) = min(zsnowevap(:),psnowliq(:))
2823 !
2824 ! This should not change enthalpy (since already accounted for
2825 ! in sfc e budget), so make density change insuring constant enthalpy:
2826 !
2827  psnowliq(:) = psnowliq(:) - zsnowevapx(:)
2828  psnowrho(:) = (psnowheat(:,1)-xlmtt*xrholw*psnowliq(:))/ &
2829  (psnowdz(:,1)*(xci*(psnowtemp(:)-xtt)-xlmtt))
2830 !
2831 ! Budget check: If density drops below minimum, then extract the
2832 ! corresponding water mass from soil (for vanishingly thin snow covers):
2833 !
2834  psoilcor(:) = max(0.0,xrhosmin_es-psnowrho(:))*psnowdz(:,1)/ptstep
2835  psnowrho(:) = max(xrhosmin_es,psnowrho(:))
2836 !
2837 END WHERE
2838 !
2839 ! 2. Update heat capacity:
2840 ! ---------------------------
2841 !
2842 zscap(:) = snow3lscap(psnowrho(:))
2843 !
2844 !
2845 ! 3. Sublimation of snow ice
2846 ! ---------------------------
2847 !
2848 WHERE(psnowdz(:,1) > 0.0)
2849 
2850 ! Budget check: as last traces of liquid in snow evaporates, it is possible
2851 ! evaporation could exceed liquid present in snow. If this is the case,
2852 ! remove residual mass from solid snow in order to maintain high
2853 ! accuracy water balance:
2854 
2855  zsnowevapx(:) = max(0.0, zsnowevap(:) - zsnowevapx(:))
2856  zsnowdz(:) = psnowdz(:,1) - zsnowevapx(:)*xrholw/psnowrho(:)
2857  psnowdz(:,1) = max(0.0, zsnowdz(:))
2858  psoilcor(:) = psoilcor(:) + max(0.0,-zsnowdz(:))*psnowrho(:)/ptstep
2859 
2860 ! Sublimation: Reduce layer thickness and total snow depth
2861 ! if sublimation: add to correction term if potential
2862 ! sublimation exceeds available snow cover.
2863 !
2864  zsnowevaps(:) = ppsn3l(:)*ples3l(:)*ptstep/(plstt(:)*psnowrho(:))
2865  zsnowdz(:) = psnowdz(:,1) - zsnowevaps(:)
2866  psnowdz(:,1) = max(0.0, zsnowdz(:))
2867  psoilcor(:) = psoilcor(:) + max(0.0,-zsnowdz(:))*psnowrho(:)/ptstep
2868 
2869 ! Note that the effect on enthalpy is already accounted for via the surface energy
2870 ! budget, so the change in enthalpy
2871 ! here owing to a mass loss must be offset by possible adjustments to T and Wl
2872 ! to conserve Enthalpy.
2873 
2874  psnowtemp(:) = xtt + ( ((psnowheat(:,1)/max(zsnowdemin,psnowdz(:,1))) &
2875  + xlmtt*psnowrho(:))/zscap(:) )
2876 
2877  psnowliq(:) = max(0.0,psnowtemp(:)-xtt)*zscap(:)* &
2878  psnowdz(:,1)/(xlmtt*xrholw)
2879 
2880  psnowtemp(:) = min(xtt,psnowtemp(:))
2881 
2882 ! For vanishingly thin snow, a decoupling can occur between forcing level T
2883 ! and snow sfc T (as snowpack vanishes owing to sublimation), so a simple limit
2884 ! imposed on this gradient:
2885 
2886  psnowtemp(:) = max(min(xtt,pta(:)-ztdif), psnowtemp(:))
2887 
2888 ! Update surface enthalpy:
2889 
2890  zsnowheat(:) = psnowheat(:,1)
2891  psnowheat(:,1) = psnowdz(:,1)*( zscap(:)*(psnowtemp(:)-xtt) &
2892  - xlmtt*psnowrho(:) ) + xlmtt*xrholw*psnowliq(:)
2893 
2894  zxse(:) = psnowheat(:,1) - zsnowheat(:) ! excess cooling (removed from sfc)
2895 
2896 END WHERE
2897 !
2898 ! If the sfc T-Gradient was limited (and thus the enthalpy in the uppermost layer changed),
2899 ! distribute this energy correction (cooling) in the layers below (thickness-based weighting)
2900 ! to conserve total snowpack enthalpy. Normally this only occurs for this snowpacks,
2901 ! but it can rarely occur otherwise also.
2902 !
2903 zisnowd(:) = 0.
2904 DO jj=2,inlvls
2905  DO ji=1,ini
2906  zisnowd(ji) = zisnowd(ji) + psnowdz(ji,jj) ! m
2907  ENDDO
2908 END DO
2909 zisnowd(:) = zxse(:)/max(zisnowd(:),zsnowdemin) ! J kg-1 m-1
2910 !
2911 DO jj=2,inlvls
2912  DO ji=1,ini
2913  psnowheat(ji,jj) = psnowheat(ji,jj) - psnowdz(ji,jj)*zisnowd(ji)
2914  ENDDO
2915 ENDDO
2916 !
2917 IF (lhook) CALL dr_hook('SNOW3LEVAPN',1,zhook_handle)
2918 !
2919 !-------------------------------------------------------------------------------
2920 !
2921 END SUBROUTINE snow3levapn
2922 !####################################################################
2923 !####################################################################
2924 !####################################################################
2925 SUBROUTINE snow3lgone(PTSTEP,PLEL3L,PLES3L,PSNOWRHO, &
2926  psnowheat,pradsink,pevapcor,pthrufal,pgrndflux, &
2927  pgfluxsnow,pgrndfluxo,psnowdz,psnowliq,psnowtemp, &
2928  plvtt,plstt,pradxs)
2929 !
2930 !
2931 !
2932 !! PURPOSE
2933 !! -------
2934 ! Account for the case when the last trace of snow melts
2935 ! during a time step: ensure mass and heat balance of
2936 ! snow AND underlying surface.
2937 !
2938 !
2939 USE modd_csts, ONLY : xtt, xlstt, xlvtt
2940 USE modd_snow_metamo, ONLY : xsnowdzmin
2941 !
2942 IMPLICIT NONE
2943 !
2944 !* 0.1 declarations of arguments
2945 !
2946 REAL, INTENT(IN) :: ptstep
2947 !
2948 REAL, DIMENSION(:), INTENT(IN) :: plel3l, ples3l, pgfluxsnow, &
2949  pradsink, pgrndfluxo, &
2950  plvtt, plstt
2951 !
2952 REAL, DIMENSION(:,:), INTENT(IN) :: psnowrho, psnowheat
2953 !
2954 REAL, DIMENSION(:), INTENT(INOUT) :: pgrndflux, pradxs
2955 !
2956 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowdz, psnowliq, psnowtemp
2957 !
2958 REAL, DIMENSION(:), INTENT(OUT) :: pthrufal ! melt water [kg/(m2 s)]
2959 !
2960 REAL, DIMENSION(:), INTENT(OUT) :: pevapcor ! [kg/(m2 s)]
2961 ! PEVAPCOR = for vanishingy thin snow cover,
2962 ! allow any excess evaporation
2963 ! to be extracted from the soil
2964 ! to maintain an accurate water
2965 ! balance.
2966 !
2967 !* 0.2 declarations of local variables
2968 !
2969 INTEGER :: jj, ji
2970 !
2971 INTEGER :: ini
2972 INTEGER :: inlvls
2973 !
2974 REAL, DIMENSION(SIZE(PLES3L)) :: zsnowheatc, zsnowgone_delta, zsnow
2975 REAL(KIND=JPRB) :: zhook_handle
2976 !
2977 !-------------------------------------------------------------------------------
2978 !
2979 ! 0. Initialize:
2980 ! --------------
2981 !
2982 IF (lhook) CALL dr_hook('SNOW3LGONE',0,zhook_handle)
2983 !
2984 ini = SIZE(psnowdz(:,:),1)
2985 inlvls = SIZE(psnowdz(:,:),2)
2986 !
2987 pevapcor(:) = 0.0
2988 pthrufal(:) = 0.0
2989 !
2990 zsnowheatc(:) = 0.
2991 zsnow(:) = 0.
2992 DO jj=1,inlvls
2993  DO ji=1,ini
2994  zsnowheatc(ji) = zsnowheatc(ji) + psnowheat(ji,jj) ! total heat content (J m-2)
2995  zsnow(ji) = zsnow(ji) + psnowdz(ji,jj) ! total snow depth (m)
2996  ENDDO
2997 ENDDO
2998 zsnowgone_delta(:) = 1.0
2999 !
3000 ! 1. Simple test to see if snow vanishes:
3001 ! ---------------------------------------
3002 ! If so, set thicknesses (and therefore mass and heat) and liquid content
3003 ! to zero, and adjust fluxes of water, evaporation and heat into underlying
3004 ! surface. Note, test with flux computed *before* correction since this represents
3005 ! actual inflow of heat from below (as heat content correction owing to a corrected
3006 ! flux has not yet been done: here we compare to pre-corrected heat content).
3007 !
3008 WHERE(pgfluxsnow(:) + pradsink(:) >= (-zsnowheatc(:)/ptstep) )
3009  pgrndflux(:) = pgfluxsnow(:) + (zsnowheatc(:)/ptstep)
3010  pevapcor(:) = (plel3l(:)/plvtt(:)) + (ples3l(:)/plstt(:))
3011  pradxs(:) = 0.0
3012  zsnowgone_delta(:) = 0.0 ! FLAG...if=0 then snow vanishes, else=1
3013 END WHERE
3014 !
3015 DO jj=1,inlvls
3016  DO ji=1,ini
3017  pthrufal(ji) = pthrufal(ji) + (1.0-zsnowgone_delta(ji))*psnowrho(ji,jj)*psnowdz(ji,jj)/ptstep
3018  ENDDO
3019 END DO
3020 !
3021 ! 2. Final update of snow state
3022 ! -----------------------------
3023 ! (either still present or not):
3024 !
3025 DO jj=1,inlvls
3026  DO ji=1,ini
3027  psnowdz(ji,jj) = psnowdz(ji,jj)*zsnowgone_delta(ji)
3028  psnowliq(ji,jj) = psnowliq(ji,jj)*zsnowgone_delta(ji)
3029  psnowtemp(ji,jj) = (1.0-zsnowgone_delta(ji))*xtt + psnowtemp(ji,jj)*zsnowgone_delta(ji)
3030  ENDDO
3031 ENDDO
3032 IF (lhook) CALL dr_hook('SNOW3LGONE',1,zhook_handle)
3033 !
3034 !
3035 END SUBROUTINE snow3lgone
3036 !####################################################################
3037 !####################################################################
3038 !####################################################################
3039 SUBROUTINE snow3levapgone(PSNOWHEAT,PSNOWDZ,PSNOWRHO,PSNOWTEMP,PSNOWLIQ)
3040 !
3041 !! PURPOSE
3042 !! -------
3043 !
3044 ! If all snow in uppermost layer evaporates/sublimates, re-distribute
3045 ! grid (below assumes very thin snowpacks so layer-thicknesses are
3046 ! constant).
3047 !
3048 !
3049 USE modd_csts, ONLY : xtt, xrholw, xlmtt
3050 USE modd_snow_par, ONLY : xrhosmin_es, xsnowdmin, xrhosmax_es
3051 !
3052 IMPLICIT NONE
3053 !
3054 !* 0.1 declarations of arguments
3055 !
3056 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowrho ! snow density profile (kg/m3)
3057 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowdz ! snow layer thickness profile (m)
3058 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowheat ! snow heat content/enthalpy (J/m2)
3059 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowtemp ! snow temperature profile (K)
3060 REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowliq ! snow liquid water profile (m)
3061 !
3062 !* 0.2 declarations of local variables
3063 !
3064 INTEGER :: jj, ji
3065 !
3066 INTEGER :: ini
3067 INTEGER :: inlvls
3068 !
3069 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: zsnowheat_1d ! total heat content (J/m2)
3070 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: zsnow ! total snow depth (m)
3071 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: zmass ! total mass
3072 !
3073 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)) :: zscap ! Snow layer heat capacity (J/K/m3)
3074 !
3075 REAL(KIND=JPRB) :: zhook_handle
3076 !
3077 !-------------------------------------------------------------------------------
3078 !
3079 ! Initialize:
3080 !
3081 IF (lhook) CALL dr_hook('SNOW3LEVAPGONE',0,zhook_handle)
3082 ini = SIZE(psnowdz,1)
3083 inlvls = SIZE(psnowdz,2)
3084 !
3085 ! First, determine where uppermost snow layer has completely
3086 ! evaporated/sublimated (as it becomes thin):
3087 !
3088 zsnowheat_1d(:) = 0.
3089 zsnow(:) = 0.
3090 zmass(:) = 0.
3091 !
3092 zscap(:,:) = snow3lscap(psnowrho(:,:))
3093 !
3094 DO jj=2,inlvls
3095  DO ji=1,ini
3096  IF(psnowdz(ji,1) == 0.0)THEN
3097  zsnowheat_1d(ji) = zsnowheat_1d(ji) + xlmtt*xrholw*psnowliq(ji,jj) &
3098  + psnowdz(ji,jj)*(zscap(ji,jj)*(psnowtemp(ji,jj)-xtt) &
3099  - xlmtt*psnowrho(ji,jj) )
3100  zsnow(ji) = zsnow(ji) + psnowdz(ji,jj)
3101  zmass(ji) = zmass(ji) + psnowdz(ji,jj)*psnowrho(ji,jj)
3102  ENDIF
3103  ENDDO
3104 ENDDO
3105 !
3106 ! Where uppermost snow layer has vanished, redistribute vertical
3107 ! snow mass and heat profiles (and associated quantities):
3108 !
3109 DO jj=1,inlvls
3110  DO ji=1,ini
3111  IF(zsnow(ji)/= 0.0)THEN
3112  zsnow(ji) = max(0.5*xsnowdmin,zsnow(ji))
3113  psnowdz(ji,jj) = zsnow(ji)/REAL(inlvls)
3114  psnowheat(ji,jj) = zsnowheat_1d(ji)/REAL(inlvls)
3115  psnowrho(ji,jj) = zmass(ji)/zsnow(ji)
3116  ENDIF
3117  ENDDO
3118 ENDDO
3119 !
3120 zscap(:,:) = snow3lscap(psnowrho(:,:))
3121 !
3122 DO jj=1,inlvls
3123  DO ji=1,ini
3124  IF(zsnow(ji)/= 0.0)THEN
3125  psnowtemp(ji,jj) = xtt + ( ((psnowheat(ji,jj)/psnowdz(ji,jj)) &
3126  + xlmtt*psnowrho(ji,jj))/zscap(ji,jj) )
3127  psnowliq(ji,jj) = max(0.0,psnowtemp(ji,jj)-xtt)*zscap(ji,jj)* &
3128  psnowdz(ji,jj)/(xlmtt*xrholw)
3129  psnowtemp(ji,jj) = min(xtt,psnowtemp(ji,jj))
3130  ENDIF
3131  ENDDO
3132 ENDDO
3133 IF (lhook) CALL dr_hook('SNOW3LEVAPGONE',1,zhook_handle)
3134 !
3135 END SUBROUTINE snow3levapgone
3136 !####################################################################
3137 !####################################################################
3138 !####################################################################
3139 SUBROUTINE snow3lebudmeb(PTSTEP, PSNOWDZMIN, &
3140  pts, psnowdz1, psnowdz2, pscond1, pscond2, pscap, &
3141  pswnetsnows, plwnetsnow, &
3142  phsnow, ples3l, plel3l, phpsnow, &
3143  pct, ptsterm1, ptsterm2, pgfluxsnow )
3144 !
3145 !! PURPOSE
3146 !! -------
3147 ! Calculate surface energy budget with surface fluxes imposed.
3148 !
3149 IMPLICIT NONE
3150 !
3151 !* 0.1 declarations of arguments
3152 !
3153 REAL, INTENT(IN) :: ptstep, psnowdzmin
3154 !
3155 REAL, DIMENSION(:), INTENT(IN) :: pts, psnowdz1, psnowdz2, pscond1, pscond2, pscap, &
3156  phsnow, ples3l, plel3l, phpsnow, &
3157  pswnetsnows, plwnetsnow
3158 !
3159 REAL, DIMENSION(:), INTENT(OUT) :: pct, ptsterm1, ptsterm2, pgfluxsnow
3160 !
3161 !* 0.2 declarations of local variables
3162 !
3163 REAL, DIMENSION(SIZE(PTS)) :: zsconda, za, zb, zc, &
3164  zsnowdzm1, zsnowdzm2
3165 !
3166 REAL(KIND=JPRB) :: zhook_handle
3167 !-------------------------------------------------------------------------------
3168 !
3169 IF (lhook) CALL dr_hook('SNOW3LEBUDMEB',0,zhook_handle)
3170 !
3171 ! Calculate surface energy budget components:
3172 ! ---------------------------------------------------------
3173 ! To prevent numerical difficulties for very thin snow
3174 ! layers, limit the grid "thinness": this is important as
3175 ! layers become vanishing thin:
3176 !
3177 zsnowdzm1(:) = max(psnowdz1(:), psnowdzmin)
3178 zsnowdzm2(:) = max(psnowdz2(:), psnowdzmin)
3179 !
3180 ! Surface thermal inertia:
3181 !
3182 pct(:) = 1.0/(pscap(:)*zsnowdzm1(:))
3183 !
3184 ! Surface fluxes entering the snowpack (radiative and turbulent):
3185 !
3186 pgfluxsnow(:) = pswnetsnows(:) + plwnetsnow(:) - phsnow(:) - ples3l(:) - plel3l(:)
3187 !
3188 ! Thermal conductivity between uppermost and lower snow layers:
3189 !
3190 zsconda(:) = (zsnowdzm1(:)+zsnowdzm2(:))/ &
3191  ((zsnowdzm1(:)/pscond1(:)) + (zsnowdzm2(:)/pscond2(:)))
3192 !
3193 !
3194 ! Energy budget solution terms (with surface flux imposed):
3195 !
3196 zb(:) = 1./ptstep
3197 !
3198 za(:) = zb(:) + pct(:)*(2*zsconda(:)/(zsnowdzm2(:)+zsnowdzm1(:)))
3199 !
3200 zc(:) = pct(:)*( pgfluxsnow(:) + phpsnow(:) )
3201 !
3202 ! Coefficients needed for implicit solution
3203 ! of linearized surface energy budget:
3204 !
3205 ptsterm2(:) = 2*zsconda(:)*pct(:)/(za(:)*(zsnowdzm2(:)+zsnowdzm1(:)))
3206 !
3207 ptsterm1(:) = (pts(:)*zb(:) + zc(:))/za(:)
3208 !
3209 !-------------------------------------------------------------------------------
3210 IF (lhook) CALL dr_hook('SNOW3LEBUDMEB',1,zhook_handle)
3211 !
3212 END SUBROUTINE snow3lebudmeb
3213 !####################################################################
3214 !####################################################################
3215 !####################################################################
3216 !
3217 !
3218 !
3219 END SUBROUTINE snow3l
snow3lrad
subroutine snow3lrad(OMEB, PSNOWDZMIN, PSW_RAD, PSNOWALB, PSPECTRALALBEDO, PSNOWDZ, PSNOWRHO, PALB, PPERMSNOWFRAC, PZENITH, PSWNETSNOW, PSWNETSNOWS, PRADSINK, PRADXS, PSNOWAGE)
Definition: snow3l.F90:1510
snow3ldrift
subroutine snow3ldrift(PTSTEP, PFORESTFRAC, PVMOD, PTA, PQA, PPS, PRHOA, PSNOWRHO, PSNOWDZ, PSNOW, OSNOWDRIFT_SUBLIM, PSNDRIFT)
Definition: snow3l.F90:1075
tridiag_ground
subroutine tridiag_ground(PA, PB, PC, PY, PX)
Definition: tridiag_ground.F90:6
mode_snow3l
Definition: mode_snow3l.F90:6
surface_ri
subroutine surface_ri(PTG, PQS, PEXNS, PEXNA, PTA, PQA, PZREF, PUREF, PDIRCOSZW, PVMOD, PRI)
Definition: surface_ri.F90:6
mode_snow3l::snow3lthrm
Definition: mode_snow3l.F90:122
snow3lrefrz
subroutine snow3lrefrz(PTSTEP, PRR, PSNOWRHO, PSNOWTEMP, PSNOWDZ, PSNOWLIQ, PTHRUFAL)
Definition: snow3l.F90:2345
mode_snow3l::snow3lhold
Definition: mode_snow3l.F90:67
snow3levapgone
subroutine snow3levapgone(PSNOWHEAT, PSNOWDZ, PSNOWRHO, PSNOWTEMP, PSNOWLIQ)
Definition: snow3l.F90:3039
snow3lebud
subroutine snow3lebud(HSNOWRES, HIMPLICIT_WIND, PPEW_A_COEF, PPEW_B_COEF, PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF, PSNOWDZMIN, PZREF, PTS, PSNOWRHO, PSNOWLIQ, PSCAP, PSCOND1, PSCOND2, PUREF, PEXNS, PEXNA, PDIRCOSZW, PVMOD, PLW_RAD, PSW_RAD, PTA, PQA, PPS, PTSTEP, PSNOWDZ1, PSNOWDZ2, PALBT, PZ0, PZ0EFF, PZ0H, PSFCFRZ, PRADSINK, PHPSNOW, PCT, PEMIST, PRHOA, PTSTERM1, PTSTERM2, PRA, PCDSNOW, PCHSNOW, PQSAT, PDQSAT, PRSRA, PUSTAR2_IC, PRI, PPET_A_COEF_T, PPEQ_A_COEF_T, PPET_B_COEF_T, PPEQ_B_COEF_T)
Definition: snow3l.F90:1713
snow3lgone
subroutine snow3lgone(PTSTEP, PLEL3L, PLES3L, PSNOWRHO, PSNOWHEAT, PRADSINK, PEVAPCOR, PTHRUFAL, PGRNDFLUX, PGFLUXSNOW, PGRNDFLUXO, PSNOWDZ, PSNOWLIQ, PSNOWTEMP, PLVTT, PLSTT, PRADXS)
Definition: snow3l.F90:2925
snow3levapn
subroutine snow3levapn(PPSN3L, PLES3L, PLEL3L, PTSTEP, PSNOWTEMP, PSNOWRHO, PSNOWDZ, PSNOWLIQ, PTA, PLVTT, PLSTT, PSNOWHEAT, PSOILCOR)
Definition: snow3l.F90:2729
mode_thermos
Definition: mode_thermos.F90:6
snow3lebudmeb
subroutine snow3lebudmeb(PTSTEP, PSNOWDZMIN, PTS, PSNOWDZ1, PSNOWDZ2, PSCOND1, PSCOND2, PSCAP, PSWNETSNOWS, PLWNETSNOW, PHSNOW, PLES3L, PLEL3L, PHPSNOW, PCT, PTSTERM1, PTSTERM2, PGFLUXSNOW)
Definition: snow3l.F90:3139
mode_snow3l::snow3lscap
Definition: mode_snow3l.F90:80
snow3lsolvt
subroutine snow3lsolvt(OMEB, PTSTEP, PSNOWDZMIN, PSNOWDZ, PSCOND, PSCAP, PTG, PSOILCOND, PD_G, PRADSINK, PCT, PTERM1, PTERM2, PPET_A_COEF_T, PPEQ_A_COEF_T, PPET_B_COEF_T, PPEQ_B_COEF_T, PTA_IC, PQA_IC, PGRNDFLUX, PGRNDFLUXO, PSNOWTEMP, PSNOWFLUX)
Definition: snow3l.F90:1954
snow3lmelt
subroutine snow3lmelt(PTSTEP, PSCAP, PSNOWTEMP, PSNOWDZ, PSNOWRHO, PSNOWLIQ, PMELTXS)
Definition: snow3l.F90:2203
snow3lcompactn
subroutine snow3lcompactn(PTSTEP, PSNOWDZMIN, PSNOWRHO, PSNOWDZ, PSNOWTEMP, PSNOW, PSNOWLIQ)
Definition: preps_for_meb_ebud_rad.F90:862
modd_meb_par
Definition: modd_meb_par.F90:6
surface_aero_cond
subroutine surface_aero_cond(PRI, PZREF, PUREF, PVMOD, PZ0, PZ0H, PAC, PRA, PCH)
Definition: surface_aero_cond.F90:6
mode_thermos::dqsati
Definition: mode_thermos.F90:75
snow3lflux
subroutine snow3lflux(PSNOWTEMP, PSNOWDZ, PEXNS, PEXNA, PUSTAR2_IC, PTSTEP, PALBT, PSW_RAD, PEMIST, PLWUPSNOW, PLW_RAD, PLWNETSNOW, PTA, PSFCFRZ, PQA, PHPSNOW, PSNOWTEMPO1, PSNOWFLUX, PCT, PRADSINK, PQSAT, PDQSAT, PRSRA, PRN, PH, PGFLUX, PLES3L, PLEL3L, PEVAP, PUSTAR, OSFCMELT)
Definition: snow3l.F90:2541
mode_snow3l::snow3lgrid
Definition: mode_snow3l.F90:91
surface_cd
subroutine surface_cd(PRI, PZREF, PUREF, PZ0EFF, PZ0H, PCD, PCDN)
Definition: surface_cd.F90:6
modd_type_date_surf::date_time
Definition: modd_type_date_surf.F90:46
mode_thermos::qsati
Definition: mode_thermos.F90:71
snow3l
subroutine snow3l(HSNOWRES, TPTIME, OMEB, HIMPLICIT_WIND, PPEW_A_COEF, PPEW_B_COEF, PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF, PSNOWSWE, PSNOWRHO, PSNOWHEAT, PSNOWALB, PSNOWGRAN1, PSNOWGRAN2, PSNOWHIST, PSNOWAGE, PTSTEP, PPS, PSR, PRR, PPSN3L, PTA, PTG, PSW_RAD, PQA, PVMOD, PLW_RAD, PRHOA, PUREF, PEXNS, PEXNA, PDIRCOSZW, PZREF, PZ0, PZ0EFF, PZ0H, PALB, PSOILCOND, PD_G, PLVTT, PLSTT, PSNOWLIQ, PSNOWTEMP, PSNOWDZ, PTHRUFAL, PGRNDFLUX, PEVAPCOR, PSOILCOR, PGFLXCOR, PSNOWSFCH, PDELHEATN, PDELHEATN_SFC, PSWNETSNOW, PSWNETSNOWS, PLWNETSNOW, PSNOWFLUX, PRNSNOW, PHSNOW, PGFLUXSNOW, PHPSNOW, PLES3L, PLEL3L, PEVAP, PSNDRIFT, PRI, PEMISNOW, PCDSNOW, PUSTAR, PCHSNOW, PSNOWHMASS, PQS, PPERMSNOWFRAC, PFORESTFRAC, PZENITH, PXLAT, PXLON, OSNOWDRIFT, OSNOWDRIFT_SUBLIM)
Definition: snow3l.F90:6
modd_surf_par
Definition: modd_surf_par.F90:6
modd_snow_metamo
Definition: modd_snow_metamo.F90:8
modd_snow_par
Definition: modd_snow_par.F90:6
snow3lfall
subroutine snow3lfall(PTSTEP, PSR, PTA, PVMOD, PSNOW, PSNOWRHO, PSNOWDZ, PSNOWHEAT, PSNOWHMASS, PSNOWAGE, PPERMSNOWFRAC)
Definition: preps_for_meb_ebud_rad.F90:172
modd_type_date_surf
Definition: modd_type_date_surf.F90:6
mode_snow3l::snow3ldopt
Definition: mode_snow3l.F90:126
mode_snow3l::snow3lalb
Definition: mode_snow3l.F90:132
modd_csts
Definition: modd_csts.F90:6
snow3ltransf
subroutine snow3ltransf(PSNOW, PSNOWDZ, PSNOWDZN, PSNOWRHO, PSNOWHEAT, PSNOWAGE)
Definition: preps_for_meb_ebud_rad.F90:334
Generated on Fri May 13 2016 09:31:05 for SURFEX by   doxygen 1.8.5