SURFEX v8.1
General documentation of Surfex
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 REAL, DIMENSION(:), INTENT(INOUT) :: PSNOWALB
188 ! PSNOWALB = Prognostic surface snow albedo
189 ! (does not include anything but
190 ! the actual snow cover)
191 !
192 REAL, DIMENSION(:,:), INTENT(INOUT):: PSNOWHEAT, PSNOWRHO, PSNOWSWE
193 ! PSNOWHEAT = Snow layer(s) heat content (J/m2)
194 ! PSNOWRHO = Snow layer(s) averaged density (kg/m3)
195 ! PSNOWSWE = Snow layer(s) liquid Water Equivalent (SWE:kg m-2)
196 !
197 REAL, DIMENSION(:,:), INTENT(INOUT):: PSNOWGRAN1, PSNOWGRAN2, PSNOWHIST
198 ! PSNOWGRAN1 = Snow layers grain feature 1
199 ! PSNOWGRAN2 = Snow layer grain feature 2
200 ! PSNOWHIST = Snow layer grain historical
201 ! parameter (only for non
202 ! dendritic snow)
203 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWAGE ! Snow grain age
204 !
205 REAL, DIMENSION(:), INTENT(INOUT) :: PRNSNOW, PHSNOW, PLES3L, PLEL3L, &
206  PHPSNOW, PEVAP, PGRNDFLUX, PEMISNOW
207 ! PLES3L = evaporation heat flux from snow (W/m2)
208 ! PLEL3L = sublimation (W/m2)
209 ! PHPSNOW = heat release from rainfall (W/m2)
210 ! PRNSNOW = net radiative flux from snow (W/m2)
211 ! PHSNOW = sensible heat flux from snow (W/m2)
212 ! PEVAP = total evaporative flux (kg/m2/s)
213 ! PGRNDFLUX = soil/snow interface heat flux (W/m2)
214 ! PEMISNOW = snow surface emissivity
215 !
216 REAL, DIMENSION(:), INTENT(OUT) :: PGFLUXSNOW
217 ! PGFLUXSNOW = net heat flux from snow (W/m2)
218 !
219 REAL, DIMENSION(:), INTENT(INOUT) :: PSWNETSNOW, PLWNETSNOW, PSWNETSNOWS
220 ! PSWNETSNOW = net shortwave radiation entering top of snowpack
221 ! (W m-2) Imposed if MEB=T, diagnosed herein if MEB=F
222 ! PSWNETSNOWS= net shortwave radiation in uppermost layer of snowpack
223 ! (W m-2) Imposed if MEB=T, diagnosed herein if MEB=F
224 ! Used for surface energy budget diagnostics
225 ! PLWNETSNOW = net longwave radiation entering top of snowpack
226 ! (W m-2) Imposed if MEB=T, diagnosed herein if MEB=F
227 !
228 REAL, DIMENSION(:), INTENT(INOUT) :: PUSTAR, PCDSNOW, PCHSNOW, PRI
229 ! PCDSNOW = drag coefficient for momentum over snow (-)
230 ! PUSTAR = friction velocity over snow (m/s)
231 ! PCHSNOW = drag coefficient for heat over snow (-)
232 ! PRI = Richardson number (-)
233 !
234 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWTEMP
235 REAL, DIMENSION(:,:), INTENT(OUT) :: PSNOWLIQ, PSNOWDZ
236 ! PSNOWLIQ = Snow layer(s) liquid water content (m)
237 ! PSNOWTEMP = Snow layer(s) temperature (m)
238 ! PSNOWDZ = Snow layer(s) thickness (m)
239 !
240 REAL, DIMENSION(:), INTENT(OUT) :: PTHRUFAL, PEVAPCOR, PSOILCOR, PGFLXCOR, &
241  PSNOWFLUX, PSNOWSFCH, PDELHEATN, PDELHEATN_SFC
242 ! PTHRUFAL = rate that liquid water leaves snow pack:
243 ! paritioned into soil infiltration/runoff
244 ! by ISBA [kg/(m2 s)]
245 ! PEVAPCOR = evaporation/sublimation correction term:
246 ! extract any evaporation exceeding the
247 ! actual snow cover (as snow vanishes)
248 ! and apply it as a surface soil water
249 ! sink. [kg/(m2 s)]
250 ! PSOILCOR = for vanishingy thin snow cover,
251 ! allow any excess evaporation
252 ! to be extracted from the soil
253 ! to maintain an accurate water
254 ! balance [kg/(m2 s)]
255 ! PGFLXCOR = flux correction to underlying soil for vanishing snowpack
256 ! (to put any energy excess from snow to soil) (W/m2)
257 ! PSNOWFLUX = heat flux between the surface and sub-surface
258 ! snow layers (W/m2)
259 ! PSNOWSFCH = snow surface layer pseudo-heating term owing to
260 ! changes in grid thickness (W m-2)
261 ! PDELHEATN = total snow heat content change in the surface layer (W m-2)
262 ! PDELHEATN_SFC = total snow heat content change during the timestep (W m-2)
263 !
264 REAL, DIMENSION(:), INTENT(OUT) :: PSNDRIFT
265 ! PSNDRIFT = blowing snow sublimation (kg/m2/s)
266 !
267 REAL, DIMENSION(:), INTENT(OUT) :: PSNOWHMASS
268 ! PSNOWHMASS = heat content change due to mass
269 ! changes in snowpack (J/m2): for budget
270 ! calculations only.
271 !
272 REAL, DIMENSION(:), INTENT(OUT) :: PQS
273 ! PQS = surface humidity
274 !
275 REAL, DIMENSION(:), INTENT(IN) :: PZENITH ! solar zenith angle
276 REAL, DIMENSION(:), INTENT(IN) :: PXLAT,PXLON ! LAT/LON after packing
277 !
278 LOGICAL, INTENT(IN) :: OSNOWDRIFT, OSNOWDRIFT_SUBLIM ! activate snowdrift, sublimation during drift
279 !
280 !* 0.2 declarations of local variables
281 !
282 INTEGER :: JJ, JI ! Loop control
283 !
284 INTEGER :: INI ! number of point
285 INTEGER :: INLVLS ! number of snow layers
286 !
287 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWTEMP, ZSCAP, ZSNOWDZN, ZSCOND, &
288  ZRADSINK, ZWORK2D, ZSNOWTEMPO
289 ! ZSNOWTEMP = Snow layer(s) averaged temperature (K)
290 ! ZSCAP = Snow layer(s) heat capacity [J/(K m3)]
291 ! ZSNOWDZN = Updated snow layer thicknesses (m)
292 ! ZSCOND = Snow layer(s) thermal conducivity [W/(m K)]
293 ! ZRADSINK = Snow solar Radiation source terms (W/m2)
294 ! ZWORK2D = working variable (*)
295 !
296 REAL, DIMENSION(SIZE(PTA)) :: ZSNOW, ZSFCFRZ, ZTSTERM1, ZTSTERM2, &
297  ZCT, ZRA, ZSNOWTEMPO1
298 ! ZSNOW = Total snow depth (m)
299 ! ZCT = inverse of the product of snow heat capacity
300 ! and layer thickness [(m2 K)/J]
301 ! ZRA = Surface aerodynamic resistance
302 ! ZTSTERM1,ZTSTERM2 = Surface energy budget coefficients
303 ! ZSNOWTEMPO1= value of uppermost snow temperature
304 ! before time integration (K)
305 !
306 LOGICAL, DIMENSION(SIZE(PTA)) :: GSFCMELT
307 ! GSFCMELT = FLAG if surface melt is occurring, used
308 ! for surface albedo calculation.
309 !
310 REAL, DIMENSION(SIZE(PTA)) :: ZRSRA, ZDQSAT, ZQSAT, ZRADXS, ZMELTXS, ZLIQHEATXS, &
311  ZLWUPSNOW, ZGRNDFLUX, ZGRNDFLUXO, ZGRNDFLUXI, ZPSN3L
312 ! ZRSRA = air density over aerodynamic resistance
313 ! ZDQSAT = derrivative of saturation specific humidity
314 ! ZQSAT = saturation specific humidity
315 ! ZRADXS = shortwave radiation absorbed by soil surface
316 ! (for thin snow sover) (W m-2)
317 ! ZMELTXS = excess energy for snowmelt for vanishingly
318 ! thin snowpacks: used to heat underlying surface
319 ! in order to maintain energy conservation (W m-2)
320 ! ZLIQHEATXS = excess snowpack heating for vanishingly thin
321 ! snow cover: add energy to snow/ground heat
322 ! flux (W m-2)
323 ! ZLWUPSNOW = upwelling longwave radiative flux (W m-2)
324 ! ZGRNDFLUXO= snow-ground flux before correction (W m-2)
325 ! ZGRNDFLUX = snow-ground flux after correction (W m-2)
326 ! ZGRNDFLUXI= for the case where the ground flux is imposed,
327 ! this is the actual imposed value.
328 ! ZPSN3L = snow fraction: different use if MEB "on".
329 ! In this case, it is only used for Tg update
330 ! since only this variable has a sub-grid relevance.
331 !
332 REAL, DIMENSION(SIZE(PTA)) :: ZUSTAR2_IC, ZTA_IC, ZQA_IC, ZWORK, ZWORK2, ZWORK3, &
333  ZPET_A_COEF_T, ZPEQ_A_COEF_T, ZPET_B_COEF_T, ZPEQ_B_COEF_T
334 ! ZUSTAR2_IC = implicit lowest atmospheric level friction (m2/s2)
335 ! ZTA_IC = implicit lowest atmospheric level air temperature
336 ! ZQA_IC = implicit lowest atmospheric level specific humidity
337 ! ZPET_A_COEF_T = transformed A-air temperature coefficient
338 ! ZPET_B_COEF_T = transformed B-air temperature coefficient
339 ! ZPEQ_A_COEF_T = transformed A-air specific humidity coefficient
340 ! ZPEQ_B_COEF_T = transformed B-air specific humidity coefficient
341 !
342 REAL, DIMENSION(SIZE(PSNOWRHO,1),NSPEC_BAND_SNOW) :: ZSPECTRALALBEDO, ZSPECTRALWORK
343 ! ZSPECTRALALBEDO=spectral albedo
344 !
345 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWHEAT0
346 !
347 REAL(KIND=JPRB) :: ZHOOK_HANDLE
348 !
349 ! - - ---------------------------------------------------
350 !
351 ! 0. Initialization
352 ! --------------
353 ! NOTE that snow layer thickness is used throughout this code: SWE
354 ! is only used to diagnose the thickness at the beginning of this routine
355 ! and it is updated at the end of this routine.
356 !
357 IF (lhook) CALL dr_hook('SNOW3L',0,zhook_handle)
358 !
359 psnowdz(:,:) = psnowswe(:,:)/psnowrho(:,:)
360 !
361 ini = SIZE(psnowswe(:,:),1)
362 inlvls = SIZE(psnowswe(:,:),2) ! total snow layers
363 !
364 zustar2_ic = 0.0
365 zta_ic = 0.0
366 zqa_ic = 0.0
367 !
368 zgrndfluxo = 0.0
369 zgrndflux = pgrndflux
370 !
371 zsnowheat0(:,:) = psnowheat(:,:) ! save initial heat content
372 !
373 zsnowtempo(:,:) = psnowtemp(:,:) ! for MEB, this is the updated T profile
374 !
375 !* 1. Snow total depth
376 ! ----------------
377 !
378 zsnow(:) = 0.
379 DO jj=1,inlvls
380  DO ji=1,ini
381  zsnow(ji) = zsnow(ji) + psnowdz(ji,jj) ! m
382  ENDDO
383 END DO
384 !
385 zwork(:)=zsnow(:)
386 !
387 ! Caluclate new snow albedo at time t
388 !
389 zwork2(:)=psnowalb(:)
390 !
391  CALL snow3lalb(zwork2,zspectralalbedo,psnowrho(:,1),psnowage(:,1),ppermsnowfrac,pps)
392 zwork3(:) = psnowalb(:)/zwork2(:)
393 DO jj=1,SIZE(zspectralalbedo,2)
394  DO ji=1,ini
395  zspectralalbedo(ji,jj)=zspectralalbedo(ji,jj)*zwork3(ji)
396  ENDDO
397 ENDDO
398 !
399 !* 2. Snowfall
400 ! --------
401 !
402 ! Caluclate uppermost density and thickness changes due to snowfall,
403 ! and add heat content of falling snow
404 !
405  CALL snow3lfall(ptstep,psr,pta,pvmod,zsnow,psnowrho,psnowdz, &
406  psnowheat,psnowhmass,psnowage,ppermsnowfrac )
407 !
408 ! Caluclate new snow albedo at time t if snowfall
409 !
410  CALL snow3lalb(zwork2,zspectralwork,psnowrho(:,1),psnowage(:,1),ppermsnowfrac,pps)
411 !
412 DO jj=1,SIZE(zspectralalbedo,2)
413  DO ji=1,ini
414  IF(zwork(ji)==0.0.AND.psr(ji)>0.0)THEN
415  zspectralalbedo(ji,jj) = zspectralwork(ji,jj)
416  ENDIF
417  ENDDO
418 ENDDO
419 WHERE(zwork(:)==0.0.AND.psr(:)>0.0)
420  psnowalb(:) = zwork2(:)
421 ENDWHERE
422 !
423 !* 3. Update grid/discretization
424 ! --------------------------
425 ! Reset grid to conform to model specifications:
426 !
427  CALL snow3lgrid(zsnowdzn,zsnow,psnowdz_old=psnowdz)
428 !
429 ! Mass/Heat redistribution:
430 !
431  CALL snow3ltransf(zsnow,psnowdz,zsnowdzn,psnowrho,psnowheat,psnowage)
432 !
433 !
434 !* 4. Liquid water content and snow temperature
435 ! -----------------------------------------
436 !
437 ! First diagnose snow temperatures and liquid
438 ! water portion of the snow from snow heat content:
439 !
440 zscap(:,:) = snow3lscap(psnowrho)
441 !
442 zsnowtemp(:,:) = xtt + ( ((psnowheat(:,:)/psnowdz(:,:)) &
443  + xlmtt*psnowrho(:,:))/zscap(:,:) )
444 !
445 psnowliq(:,:) = max(0.0,zsnowtemp(:,:)-xtt)*zscap(:,:)* &
446  psnowdz(:,:)/(xlmtt*xrholw)
447 !
448 zsnowtemp(:,:) = min(xtt,zsnowtemp(:,:))
449 !
450 !
451 !* 5. Snow Compaction
452 ! ---------------
453 !
454 ! Calculate snow density: compaction/aging: density increases
455 !
456  CALL snow3lcompactn(ptstep,xsnowdzmin,psnowrho,psnowdz,zsnowtemp,zsnow,psnowliq)
457 !
458 ! Snow compaction and metamorphism due to drift
459 !
460 psndrift(:) = 0.0
461 IF (osnowdrift) THEN
462  CALL snow3ldrift(ptstep,pforestfrac,pvmod,pta,pqa,pps,prhoa,&
463  psnowrho,psnowdz,zsnow,osnowdrift_sublim,psndrift)
464 ENDIF
465 !
466 ! Update snow heat content (J/m2):
467 !
468 zscap(:,:) = snow3lscap(psnowrho)
469 psnowheat(:,:) = psnowdz(:,:)*( zscap(:,:)*(zsnowtemp(:,:)-xtt) &
470  - xlmtt*psnowrho(:,:) ) + xlmtt*xrholw*psnowliq(:,:)
471 !
472 !
473 !* 6. Solar radiation transmission
474 ! -----------------------------
475 !
476 ! Heat source (-sink) term due to shortwave
477 ! radiation transmission within the snowpack:
478 !
479  CALL snow3lrad(omeb,xsnowdzmin,psw_rad,psnowalb, &
480  zspectralalbedo,psnowdz,psnowrho,palb, &
481  ppermsnowfrac,pzenith,pswnetsnow, &
482  pswnetsnows,zradsink,zradxs,psnowage)
483 !
484 !
485 !* 7. Heat transfer and surface energy budget
486 ! ---------------------------------------
487 ! Snow thermal conductivity:
488 !
489  CALL snow3lthrm(psnowrho,zscond,zsnowtemp,pps)
490 !
491 ! Precipitation heating term:
492 ! Rainfall renders it's heat to the snow when it enters
493 ! the snowpack:
494 ! NOTE: for now set to zero because we'd need to remove this heat from
495 ! the atmosphere to conserve energy. This can be done, but should
496 ! be tested within an atmos model first probably...NOTE, this
497 ! could be actiavted for use in OFFLINE mode however.
498 !
499 !PHPSNOW(:) = PRR(:)*XCL*(MAX(XTT,PTA(:))-XTT) ! (W/m2)
500 phpsnow(:) = 0.0
501 !
502 ! Surface Energy Budget calculations using ISBA linearized form
503 ! and standard ISBA turbulent transfer formulation
504 !
505 IF(omeb)THEN
506 
507 ! - surface fluxes prescribed:
508 
509  CALL snow3lebudmeb(ptstep,xsnowdzmin, &
510  zsnowtemp(:,1),psnowdz(:,1),psnowdz(:,2), &
511  zscond(:,1),zscond(:,2),zscap(:,1), &
512  pswnetsnows,plwnetsnow, &
513  phsnow,ples3l,plel3l,phpsnow, &
514  zct,ztsterm1,ztsterm2,pgfluxsnow)
515 
516 ! - for mass fluxes, snow fraction notion removed
517 
518  zpsn3l(:) = 1.0
519 
520 ELSE
521 
522  zpsn3l(:) = ppsn3l(:)
523 
524 ! - compute surface energy budget and fluxes:
525 
526  CALL snow3lebud(hsnowres, himplicit_wind, &
527  ppew_a_coef, ppew_b_coef, &
528  ppet_a_coef, ppeq_a_coef, ppet_b_coef, ppeq_b_coef, &
529  xsnowdzmin, &
530  pzref,zsnowtemp(:,1),psnowrho(:,1),psnowliq(:,1),zscap(:,1), &
531  zscond(:,1),zscond(:,2), &
532  puref,pexns,pexna,pdircoszw,pvmod, &
533  plw_rad,psw_rad,pta,pqa,pps,ptstep, &
534  psnowdz(:,1),psnowdz(:,2),psnowalb,pz0,pz0eff,pz0h, &
535  zsfcfrz,zradsink(:,1),phpsnow, &
536  zct,pemisnow,prhoa,ztsterm1,ztsterm2,zra,pcdsnow,pchsnow, &
537  zqsat, zdqsat, zrsra, zustar2_ic, pri, &
538  zpet_a_coef_t,zpeq_a_coef_t,zpet_b_coef_t,zpeq_b_coef_t )
539 !
540 ENDIF
541 !
542 ! Heat transfer: simple diffusion along the thermal gradient
543 !
544 zsnowtempo1(:) = zsnowtemp(:,1) ! save surface snow temperature before update
545 !
546 zgrndfluxi(:) = zgrndflux(:)
547 !
548  CALL snow3lsolvt(omeb,ptstep,xsnowdzmin,psnowdz,zscond,zscap,ptg, &
549  psoilcond,pd_g,zradsink,zct,ztsterm1,ztsterm2, &
550  zpet_a_coef_t,zpeq_a_coef_t,zpet_b_coef_t,zpeq_b_coef_t, &
551  zta_ic,zqa_ic,zgrndflux,zgrndfluxo,zsnowtemp,psnowflux )
552 !
553 !
554 !* 8. Surface fluxes
555 ! --------------
556 !
557 ! Since surface fluxes already computed under MEB option (and already
558 ! recomputed for the case when T>Tf), Only call if MEB not in use:
559 !
560 IF(.NOT.omeb)THEN
561 
562  CALL snow3lflux(zsnowtemp(:,1),psnowdz(:,1),pexns,pexna, &
563  zustar2_ic, &
564  ptstep,psnowalb,psw_rad, &
565  pemisnow,zlwupsnow,plw_rad,plwnetsnow, &
566  zta_ic,zsfcfrz,zqa_ic,phpsnow, &
567  zsnowtempo1,psnowflux,zct,zradsink(:,1), &
568  zqsat,zdqsat,zrsra, &
569  prnsnow,phsnow,pgfluxsnow,ples3l,plel3l,pevap, &
570  pustar,gsfcmelt )
571 !
572 ENDIF
573 !
574 !
575 !* 9. Snow melt
576 ! ---------
577 !
578 ! First Test to see if snow pack vanishes during this time step:
579 !
580  CALL snow3lgone(ptstep,plel3l,ples3l,psnowrho, &
581  psnowheat,zradsink(:,inlvls),pevapcor,pthrufal,zgrndflux, &
582  pgfluxsnow,zgrndfluxo,psnowdz,psnowliq,zsnowtemp, &
583  plvtt,plstt,zradxs )
584 !
585 ! For "normal" melt: transform excess heat content into snow liquid:
586 !
587  CALL snow3lmelt(ptstep,zscap,zsnowtemp,psnowdz,psnowrho,psnowliq,zmeltxs)
588 !
589 !
590 !* 10. Snow water flow and refreezing
591 ! ------------------------------
592 ! Liquid water vertical transfer and possible snowpack runoff
593 ! And refreezing/freezing of meltwater/rainfall (ripening of the snow)
594 !
595  CALL snow3lrefrz(ptstep,prr,psnowrho,zsnowtemp,psnowdz,psnowliq,pthrufal)
596 !
597 zscap(:,:) = snow3lscap(psnowrho)
598 psnowheat(:,:) = psnowdz(:,:)*( zscap(:,:)*(zsnowtemp(:,:)-xtt) &
599  - xlmtt*psnowrho(:,:) ) + xlmtt*xrholw*psnowliq(:,:)
600 !
601 !
602 !* 11. Snow Evaporation/Sublimation mass updates:
603 ! ------------------------------------------
604 !
605  CALL snow3levapn(zpsn3l,ples3l,plel3l,ptstep,zsnowtemp(:,1),psnowrho(:,1), &
606  psnowdz,psnowliq(:,1),pta,plvtt,plstt,psnowheat,psoilcor )
607 !
608 ! Update snow temperatures and liquid
609 ! water portion of the snow from snow heat content
610 ! since vertical distribution of enthalpy might have been
611 ! modified in the previous routine:
612 !
613 zscap(:,:) = snow3lscap(psnowrho)
614 zwork2d(:,:) = min(1.0, psnowdz(:,:)/xsnowdmin)
615 zsnowtemp(:,:) = xtt + zwork2d(:,:)*( ((psnowheat(:,:)/max(xsnowdmin,psnowdz(:,:))) &
616  + xlmtt*psnowrho(:,:))/zscap(:,:) )
617 psnowliq(:,:) = max(0.0,zsnowtemp(:,:)-xtt)*zscap(:,:)*psnowdz(:,:)/(xlmtt*xrholw)
618 zsnowtemp(:,:) = min(xtt,zsnowtemp(:,:))
619 !
620 !
621 ! If all snow in uppermost layer evaporates/sublimates, re-distribute
622 ! grid (below could be evoked for vanishingly thin snowpacks):
623 !
624  CALL snow3levapgone(psnowheat,psnowdz,psnowrho,zsnowtemp,psnowliq)
625 !
626 !
627 !* 12. Update surface albedo:
628 ! ----------------------
629 ! Snow clear sky albedo:
630 !
631  CALL snow3lalb(psnowalb,zspectralalbedo,psnowrho(:,1),psnowage(:,1),ppermsnowfrac,pps)
632 !
633 !
634 !* 13. Update snow heat content:
635 ! -------------------------
636 ! Update the heat content (variable stored each time step)
637 ! using current snow temperature and liquid water content:
638 !
639 ! First, make check to make sure heat content not too large
640 ! (this can result due to signifigant heating of thin snowpacks):
641 !
642 DO jj=1,inlvls
643  DO ji=1,ini
644  zliqheatxs(ji) = max(0.0,psnowliq(ji,jj)*xrholw-psnowdz(ji,jj)*psnowrho(ji,jj))*xlmtt/ptstep
645  psnowliq(ji,jj)= psnowliq(ji,jj) - zliqheatxs(ji)*ptstep/(xrholw*xlmtt)
646  psnowliq(ji,jj)= max(0.0, psnowliq(ji,jj))
647  ENDDO
648 ENDDO
649 !
650 psnowtemp(:,:) = zsnowtemp(:,:)
651 !
652 zscap(:,:) = snow3lscap(psnowrho)
653 !
654 psnowheat(:,:) = psnowdz(:,:)*( zscap(:,:)*(psnowtemp(:,:)-xtt) &
655  - xlmtt*psnowrho(:,:) ) + xlmtt*xrholw*psnowliq(:,:)
656 !
657 !
658 !* 14. Snow/Ground heat flux:
659 ! ----------------------
660 !
661 !Add radiation not absorbed by snow to soil/vegetation interface flux
662 !
663 pgrndflux(:) = zgrndfluxo(:)+zradxs(:)
664 !
665 !Comput soil/snow correction heat flux to prevent TG1 numerical oscillations
666 !Add any excess heat for melting and liquid water :
667 !
668 pgflxcor(:) = (zgrndflux(:)-zgrndfluxo(:))+zmeltxs(:)+zliqheatxs(:)
669 !
670 !
671 !* 15. Update snow mass and age:
672 ! -------------------------
673 !
674 psnowswe(:,:) = psnowdz(:,:)*psnowrho(:,:)
675 !
676 WHERE (psnowswe(:,:)>0)
677  psnowage(:,:)=psnowage(:,:)+(ptstep/xday)
678 ELSEWHERE
679  psnowage(:,:)= xundef
680 ENDWHERE
681 !
682 !
683 !* 16. Update surface specific humidity:
684 ! ---------------------------------
685 !
686 IF(omeb)THEN
687  pqs(:) = qsati(psnowtemp(:,1),pps) ! purely diagnostic
688 ELSE
689  pqs(:) = zqsat(:)
690 ENDIF
691 !
692 !
693 !* 17. MEB related checks
694 ! ------------------
695 !
696 IF(omeb)THEN
697 
698 ! Final adjustment: If using MEB, then all excess heat (above initial guess ground flux)
699 ! is used herein to adjust the ground temperature since for the MEB case, boundary
700 ! fluxes are imposed (the energy budget of the ground has already been computed).
701 ! This correction is needed since snow-soil flux here is implicit, while guess
702 ! used in surface energy soil budget was semi-implicit. This forces the flux
703 ! seen by the ground to be *the same* as that leaving the snow (energy conservation).
704 ! Also, add any excessing cooling from sublimation as snowpack becomes vanishingly thin.
705 ! This is added back to total heat content of the snowpack (and distributed among
706 ! all snow layers while conserving total heat content plus correction)
707 
708  zwork(:) = 0.
709  DO jj=1,inlvls
710  DO ji=1,ini
711  zwork(ji) = zwork(ji) + psnowheat(ji,jj)
712  ENDDO
713  ENDDO
714  zwork2(:) = min(0.0, zwork(:) + pgrndflux(:) - zgrndfluxi(:))
715  pgflxcor(:) = pgflxcor(:) + max(0., zwork2(:)) ! add any possible (rare!) excess to soil
716 
717  WHERE(zwork(:) > -1.e-10)
718  zwork(:) = 1. ! i.e. no modifs to H profile
719  ELSEWHERE
720  zwork(:) = zwork2(:)/zwork(:)
721  END WHERE
722 
723  DO jj=1,inlvls
724  DO ji=1,ini
725  psnowheat(ji,jj) = psnowheat(ji,jj)*zwork(ji)
726  ENDDO
727  ENDDO
728 
729 ! these are just updated diagnostics at this point:
730 
731  zwork2d(:,:) = min(1.0, psnowdz(:,:)/xsnowdmin)/max(xsnowdmin,psnowdz(:,:))
732  psnowtemp(:,:) = xtt + zwork2d(:,:)*( (psnowheat(:,:) + xlmtt*psnowswe(:,:))/zscap(:,:) )
733  psnowliq(:,:) = max(0.0,psnowtemp(:,:)-xtt)*zscap(:,:)*psnowdz(:,:)/(xlmtt*xrholw)
734  psnowtemp(:,:) = min(xtt,psnowtemp(:,:))
735 
736 ENDIF
737 !
738 !
739 !* 18. Energy Budget Diagnostics:
740 ! --------------------------
741 !
742 ! NOTE: To check enthalpy conservation, the error in W m-2 is defined HERE as
743 !
744 ! error = (SUM(PSNOWHEAT(:,:),2)-SUM(ZSNOWHEAT0(:,:),2))/PTSTEP &
745 ! - PSNOWHMASS(:)/PTSTEP - (PGFLUXSNOW(:) - PGRNDFLUX(:))
746 !
747 ! where ZSNOWHEAT0 is the value of PSNOWHEAT that enters this routine before any adjustments/processes.
748 ! Errors should be VERY small ( ~0) at each time step.
749 ! The above is "snow relative"...the actual tile or grid box error is multiplied by PPSN3L.
750 !
751 ! Snow heat storage change over the time step (W m-2):
752 
753 pdelheatn(:) = 0.0
754 DO jj=1,inlvls
755  DO ji=1,ini
756  pdelheatn(ji) = pdelheatn(ji) + (psnowheat(ji,jj)-zsnowheat0(ji,jj))
757  ENDDO
758 END DO
759 pdelheatn(:) = pdelheatn(:) /ptstep
760 pdelheatn_sfc(:) = (psnowheat(:,1)-zsnowheat0(:,1))/ptstep
761 
762 !
763 ! The pseudo-heating (W m-2) is derrived as a diagnostic from the surface energy budget
764 ! here...NOTE that the total snow energy budget is explicitly computed, and since
765 ! it balances, then simply compute this term as a residue (which ensures sfc energy balance):
766 !
767 psnowsfch(:) = pdelheatn_sfc(:) - (pswnetsnows(:) +plwnetsnow(:) - phsnow(:) -ples3l(:)-plel3l(:)) &
768  + psnowflux(:) - psnowhmass(:)/ptstep
769 !
770 IF (lhook) CALL dr_hook('SNOW3L',1,zhook_handle)
771 !
772 !
773  CONTAINS
774 !
775 !
776 !
777 !####################################################################
778 !####################################################################
779 !####################################################################
780 SUBROUTINE snow3ldrift(PTSTEP,PFORESTFRAC,PVMOD,PTA,PQA,PPS,PRHOA,&
781  PSNOWRHO,PSNOWDZ,PSNOW,OSNOWDRIFT_SUBLIM,PSNDRIFT)
782 !
783 USE modd_surf_par, ONLY : xundef
784 !
785 USE modd_snow_par, ONLY : xvtime, xvromax, xvromin, xvmob1, &
786  xvdrift1, xvdrift2, xvdrift3, &
787  xcoef_ff, xcoef_effect, xqs_ref
788 !
789 USE mode_thermos
790 !
791 !
792 !! PURPOSE
793 !! -------
794 ! Snow compaction and metamorphism due to drift
795 ! Mass is unchanged: layer thickness is reduced
796 ! in proportion to density increases. Method inspired from
797 ! Brun et al. (1997) and Guyomarch
798 !
799 ! - computes a mobility index of each snow layer from its density
800 ! - computes a drift index of each layer from its mobility and wind speed
801 ! - computes a transport index
802 ! - increases density in case of transport
803 !
804 ! HISTORY:
805 ! Basic parameterization from Crocus/ARPEGE Coupling (1997)
806 ! Adaptation from Crocus SNOWDRIFT subroutine
807 !
808 !! PURPOSE
809 !! -------
810 ! Snow compaction due to overburden and settling.
811 ! Mass is unchanged: layer thickness is reduced
812 ! in proportion to density increases.
813 !
814 IMPLICIT NONE
815 !
816 !* 0.1 declarations of arguments
817 !
818 REAL, INTENT(IN) :: PTSTEP
819 !
820 REAL, DIMENSION(:), INTENT(IN) :: PFORESTFRAC
821 REAL, DIMENSION(:), INTENT(IN) :: PVMOD
822 REAL, DIMENSION(:), INTENT(IN) :: PTA
823 REAL, DIMENSION(:), INTENT(IN) :: PQA
824 REAL, DIMENSION(:), INTENT(IN) :: PPS
825 REAL, DIMENSION(:), INTENT(IN) :: PRHOA
826 !
827 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWRHO, PSNOWDZ
828 !
829 REAL, DIMENSION(:), INTENT(OUT) :: PSNOW
830 !
831 LOGICAL, INTENT(IN) :: OSNOWDRIFT_SUBLIM
832 REAL, DIMENSION(:), INTENT(OUT) :: PSNDRIFT !blowing snow sublimation (kg/m2/s)
833 !
834 !* 0.2 declarations of local variables
835 !
836 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWRHO2
837 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZRMOB
838 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZRDRIFT
839 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZRT
840 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZDRO
841 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZQS_EFFECT
842 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZDRIFT_EFFECT
843 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: ZPROFEQU
844 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: ZWIND
845 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: ZQSATI
846 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: ZVT
847 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: ZQS
848 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: ZW
849 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: ZT
850 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: ZSNOWDZ1
851 REAL, DIMENSION(SIZE(PSNOWRHO,1) ) :: ZFOREST_EFFECT
852 !
853 LOGICAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: GDRIFT
854 !
855 INTEGER :: JJ, JI
856 !
857 INTEGER :: INI
858 INTEGER :: INLVLS
859 !
860 REAL(KIND=JPRB) :: ZHOOK_HANDLE
861 !
862 !-------------------------------------------------------------------------------
863 !
864 IF (lhook) CALL dr_hook('SNOW3LDRIFT',0,zhook_handle)
865 !
866 ! 0. Initialization:
867 ! ------------------
868 !
869 ini = SIZE(psnowdz(:,:),1)
870 inlvls = SIZE(psnowdz(:,:),2)
871 !
872 zsnowrho2(:,:) = psnowrho(:,:)
873 zrdrift(:,:) = xundef
874 zrt(:,:) = xundef
875 zdro(:,:) = xundef
876 zqs_effect(:,:) = 0.0
877 zdrift_effect(:,:) = 0.0
878 gdrift(:,:) = .false.
879 !
880 zprofequ(:) = 0.0
881 !
882 zsnowdz1(:) = psnowdz(:,1)
883 !
884 ! 1. Initialazation of drift and induced settling
885 ! -----------------------------------------------
886 !
887 !Limitation of wind speed in Forest : ~ 15% of wind in open area
888 zforest_effect(:) = 1.0 - 0.85 * pforestfrac(:)
889 !
890 zwind(:) = xcoef_ff * zforest_effect(:) * pvmod(:)
891 !
892 DO jj=1,inlvls
893  DO ji=1,ini
894 !
895 ! mobility index computation of a layer as a function of its density
896  zrmob(ji,jj)= 1.25-1.25e-3*(max(psnowrho(ji,jj),xvromin)-xvromin) / xvmob1
897 !
898 ! computation of the drift index inclunding the decay by overburden snow
899  zrdrift(ji,jj) = zrmob(ji,jj)-(xvdrift1*exp(-xvdrift2*zwind(ji))-1.0)
900 !
901  ENDDO
902 ENDDO
903 !
904 !Compaction from the surface to layers beneath until
905 !a snow layer having negative drift index
906 gdrift(:,:) = (zrdrift(:,:)>0.0)
907 DO ji=1,ini
908  DO jj=1,inlvls
909  IF(.NOT.gdrift(ji,jj))THEN
910  gdrift(ji,jj:inlvls)=.false.
911  EXIT
912  ENDIF
913  ENDDO
914 ENDDO
915 
916 ! 2. Specific case for blowing snow sublimation
917 ! ---------------------------------------------
918 !
919 IF(osnowdrift_sublim)THEN
920 !
921  zqsati(:)= qsati(pta(:),pps(:))
922 !
923 ! computation of wind speed threshold QSATI and RH withe respect to ice
924 !
925  zw(:)=max(-0.99,zrmob(:,1))
926  zvt(:)=log(xvdrift1/(1.0+zw(:)))/xvdrift2
927 !
928  WHERE(zwind(:)>0.0)
929  zw(:)=log(zwind(:)/zvt(:))
930  zw(:)=exp(3.6*zw(:))
931  ELSEWHERE
932  zw(:)=0.0
933  ENDWHERE
934 !
935 ! computation of sublimation rate according to Gordon's PhD
936 !
937  zt(:)=xtt/pta(:)
938  zt(:)=0.0018*(zt(:)**4)
939 !
940  zqs(:)=zt(:)*zvt(:)*prhoa(:)*zqsati(:)*(1.-pqa(:)/zqsati(:))*zw(:)
941 !
942 ! surface depth decrease in case of blowing snow sublimation
943 !
944  psnowdz(:,1)=max(0.5*psnowdz(:,1),psnowdz(:,1)-max(0.,zqs(:))*ptstep/(xcoef_ff*psnowrho(:,1)))
945 !
946  psndrift(:) = (zsnowdz1(:)-psnowdz(:,1))*psnowrho(:,1)/ptstep
947 !
948  zqs_effect(:,1)=min(3.,max(0.,zqs(:))/xqs_ref)
949 !
950 ENDIF
951 !
952 ! 3. Computation of drift and induced settling
953 ! --------------------------------------------
954 !
955 DO jj=1,inlvls
956  DO ji=1,ini
957 
958 ! update the decay coeff by half the current layer
959  zprofequ(ji) = zprofequ(ji) + 0.5 * psnowdz(ji,jj) * 0.1 * (xvdrift3-zrdrift(ji,jj))
960 !
961  IF(gdrift(ji,jj).AND.psnowrho(ji,jj)<xvromax)THEN
962 !
963 ! computation of the drift index inclunding the decay by overburden snow
964  zrt(ji,jj) = max(0.0,zrdrift(ji,jj)*exp(-zprofequ(ji)*100.0))
965 !
966  zdrift_effect(ji,jj) = (zqs_effect(ji,jj)+xcoef_effect)*zrt(ji,jj)/(xvtime*xcoef_ff)
967 !
968 ! settling by wind transport only in case of not too dense snow
969  zdro(ji,jj) = (xvromax - psnowrho(ji,jj)) * zdrift_effect(ji,jj) * ptstep
970 !
971 ! Calculate new snow density:
972  zsnowrho2(ji,jj) = min(xvromax,psnowrho(ji,jj)+zdro(ji,jj))
973 
974 ! Conserve mass by decreasing grid thicknesses in response to density increases
975  psnowdz(ji,jj) = psnowdz(ji,jj)*(psnowrho(ji,jj)/zsnowrho2(ji,jj))
976 !
977  ENDIF
978 !
979 ! update the decay coeff by half the current layer
980  zprofequ(ji) = zprofequ(ji) + 0.5 * psnowdz(ji,jj) * 0.1 * (xvdrift3-zrdrift(ji,jj))
981 !
982  ENDDO
983 ENDDO
984 !
985 ! 4. Update total snow depth and density profile:
986 ! -----------------------------------------------
987 !
988 ! Compaction/augmentation of total snowpack depth
989 !
990 psnow(:) = 0.
991 DO jj=1,inlvls
992  DO ji=1,ini
993  psnow(ji) = psnow(ji) + psnowdz(ji,jj)
994  ENDDO
995 ENDDO
996 !
997 ! Update density (kg m-3):
998 !
999 psnowrho(:,:) = zsnowrho2(:,:)
1000 !
1001 IF (lhook) CALL dr_hook('SNOW3LDRIFT',1,zhook_handle)
1002 !
1003 !-------------------------------------------------------------------------------
1004 !
1005 END SUBROUTINE snow3ldrift
1006 !####################################################################
1007 !####################################################################
1008 !####################################################################
1009  SUBROUTINE snow3lrad(OMEB, PSNOWDZMIN, PSW_RAD, PSNOWALB, &
1010  PSPECTRALALBEDO, PSNOWDZ, PSNOWRHO, PALB, &
1011  PPERMSNOWFRAC, PZENITH, PSWNETSNOW, &
1012  PSWNETSNOWS, PRADSINK, PRADXS, PSNOWAGE )
1013 !
1014 !! PURPOSE
1015 !! -------
1016 ! Calculate the transmission of shortwave (solar) radiation
1017 ! through the snowpack (using a form of Beer's Law: exponential
1018 ! decay of radiation with increasing snow depth).
1019 !
1020 USE modd_surf_par, ONLY : xundef
1021 !
1022 USE modd_meb_par, ONLY : xsw_wght_vis, xsw_wght_nir
1023 !
1024 USE mode_snow3l, ONLY : snow3ldopt, snow3lradabs_sfc
1025 !
1026 IMPLICIT NONE
1027 !
1028 !* 0.1 declarations of arguments
1029 !
1030 LOGICAL, INTENT(IN) :: OMEB ! if=T, then uppermost abs is diagnosed
1031 ! ! since fluxes known
1032 !
1033 REAL, INTENT(IN) :: PSNOWDZMIN
1034 !
1035 REAL, DIMENSION(:), INTENT(IN) :: PSW_RAD
1036 REAL, DIMENSION(:), INTENT(IN) :: PSNOWALB
1037 REAL, DIMENSION(:), INTENT(IN) :: PALB
1038 REAL, DIMENSION(:), INTENT(IN) :: PPERMSNOWFRAC
1039 REAL, DIMENSION(:), INTENT(IN) :: PZENITH
1040 !
1041 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWRHO, PSNOWDZ, PSNOWAGE
1042 REAL, DIMENSION(:,:), INTENT(IN) :: PSPECTRALALBEDO
1043 !
1044 REAL, DIMENSION(:), INTENT(INOUT) :: PSWNETSNOW, PSWNETSNOWS
1045 !
1046 REAL, DIMENSION(:), INTENT(OUT) :: PRADXS
1047 !
1048 REAL, DIMENSION(:,:), INTENT(OUT) :: PRADSINK
1049 !
1050 !
1051 !* 0.2 declarations of local variables
1052 !
1053 INTEGER :: JJ, JI
1054 !
1055 INTEGER :: INI
1056 INTEGER :: INLVLS
1057 !
1058 REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: ZRADTOT
1059 !
1060 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZDSGRAIN, ZCOEF, ZSNOWDZ, ZAGE
1061 REAL, DIMENSION(SIZE(PSPECTRALALBEDO,1),SIZE(PSPECTRALALBEDO,2)) :: ZSPECTRALALBEDO
1062 !
1063 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1064 !-------------------------------------------------------------------------------
1065 !
1066 ! 0. Initialization:
1067 ! ------------------
1068 !
1069 IF (lhook) CALL dr_hook('SNOW3LRAD',0,zhook_handle)
1070 !
1071 ini = SIZE(psnowdz(:,:),1)
1072 inlvls = SIZE(psnowdz(:,:),2)
1073 !
1074 zspectralalbedo(:,:) = xundef ! Init
1075 !
1076 ! 1. Vanishingly thin snowpack check:
1077 ! -----------------------------------
1078 ! For vanishingly thin snowpacks, much of the radiation
1079 ! can pass through snowpack into underlying soil, making
1080 ! a large (albeit temporary) thermal gradient: by imposing
1081 ! a minimum thickness, this increases the radiation absorbtion
1082 ! for vanishingly thin snowpacks.
1083 !
1084 zsnowdz(:,:) = max(psnowdzmin, psnowdz(:,:))
1085 !
1086 !
1087 ! 2. Extinction of net shortwave radiation
1088 ! ----------------------------------------
1089 ! Fn of snow depth and density (Loth and Graf 1993:
1090 ! SNOWCVEXT => from Bohren and Barkstrom 1974
1091 ! SNOWAGRAIN and SNOWBGRAIN=> from Jordan 1976)
1092 !
1093 !
1094 ! Snow optical grain diameter (no age dependency over polar regions):
1095 !
1096 DO jj=1,inlvls
1097  DO ji=1,ini
1098  zage(ji,jj) = (1.0-ppermsnowfrac(ji))*psnowage(ji,jj)
1099  ENDDO
1100 ENDDO
1101 !
1102 zdsgrain(:,:) = snow3ldopt(psnowrho,zage)
1103 !
1104 !
1105 ! Calculate the transmission of shortwave radiation within the snowpack:
1106 !
1107 IF(omeb)THEN
1108 
1109 ! Only 2 bands currently considered (for vegetation and soil)
1110 ! thus eliminate a band and renormalize (as was done for the surface snow albedo
1111 ! for the MEB snow surface energy budget computations)
1112 
1113  zspectralalbedo(:,1) = pspectralalbedo(:,1)
1114  zspectralalbedo(:,2) = (psnowalb(:) - xsw_wght_vis*zspectralalbedo(:,1))/xsw_wght_nir
1115 !
1116  zcoef(:,:) = snow3lradabs_sfc(psnowrho,zsnowdz,zspectralalbedo,pzenith,ppermsnowfrac,zdsgrain)
1117 
1118 ! Diagnose surface layer coef (should be very close/identical to ZCOEF(:,1) computed above)
1119 
1120  zcoef(:,1) = 1.0 - pswnetsnows(:)/max(1.e-4,pswnetsnow(:))
1121 
1122 ELSE
1123 
1124 ! Consider 3 bands:
1125 
1126  zcoef(:,:) = snow3lradabs_sfc(psnowrho,zsnowdz,pspectralalbedo,pzenith,ppermsnowfrac,zdsgrain)
1127 
1128  pswnetsnow(:) = psw_rad(:)*(1.-psnowalb(:))
1129  pswnetsnows(:) = pswnetsnow(:)*(1.0-zcoef(:,1))
1130 
1131 ENDIF
1132 !
1133 ! 3. Radiation at each level: (W/m2)
1134 ! ----------------------------------
1135 !
1136 ! NOTE: as this is a "sink" term for heat, it
1137 ! as assigned a negative value.
1138 DO jj=1,inlvls
1139  DO ji=1,ini
1140  pradsink(ji,jj) = -psw_rad(ji)*zcoef(ji,jj)
1141  ENDDO
1142 ENDDO
1143 !
1144 ! For thin snow packs, radiation might reach base of
1145 ! snowpack...so we influence this amount with sfc albedo
1146 ! and (outside of this routine) add any excess heat
1147 ! to underlying soil:
1148 !
1149 pradsink(:,inlvls) = pradsink(:,inlvls)*(1.0-palb(:))
1150 !
1151 ! 4. Excess radiation not absorbed by snowpack (W/m2):
1152 ! ----------------------------------------------------
1153 !
1154 zradtot(:) = pradsink(:,1) + (1.-psnowalb(:))*psw_rad(:)
1155 DO jj=2,inlvls
1156  DO ji=1,ini
1157  zradtot(ji) = zradtot(ji) + pradsink(ji,jj)-pradsink(ji,jj-1)
1158  ENDDO
1159 ENDDO
1160 !
1161 pradxs(:) = (1.-psnowalb(:))*psw_rad(:) - zradtot(:)
1162 !
1163 IF (lhook) CALL dr_hook('SNOW3LRAD',1,zhook_handle)
1164 !
1165 !-------------------------------------------------------------------------------
1166 !
1167 END SUBROUTINE snow3lrad
1168 !####################################################################
1169 !####################################################################
1170 !####################################################################
1171  SUBROUTINE snow3lebud(HSNOWRES, HIMPLICIT_WIND, &
1172  PPEW_A_COEF, PPEW_B_COEF, &
1173  PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF, &
1174  PSNOWDZMIN, &
1175  PZREF,PTS,PSNOWRHO,PSNOWLIQ,PSCAP,PSCOND1,PSCOND2, &
1176  PUREF,PEXNS,PEXNA,PDIRCOSZW,PVMOD, &
1177  PLW_RAD,PSW_RAD,PTA,PQA,PPS,PTSTEP, &
1178  PSNOWDZ1,PSNOWDZ2,PALBT,PZ0,PZ0EFF,PZ0H, &
1179  PSFCFRZ,PRADSINK,PHPSNOW, &
1180  PCT,PEMIST,PRHOA,PTSTERM1,PTSTERM2,PRA,PCDSNOW,PCHSNOW, &
1181  PQSAT,PDQSAT,PRSRA,PUSTAR2_IC,PRI, &
1182  PPET_A_COEF_T,PPEQ_A_COEF_T,PPET_B_COEF_T,PPEQ_B_COEF_T )
1183 !
1184 !! PURPOSE
1185 !! -------
1186 ! Calculate surface energy budget linearization (terms) and turbulent
1187 ! exchange coefficients/resistance between surface and atmosphere.
1188 ! (Noilhan and Planton 1989; Giordani 1993; Noilhan and Mahfouf 1996)
1189 !
1190 USE modd_csts, ONLY : xcpd, xrholw, xstefan, xlvtt, xlstt
1191 USE modd_snow_par, ONLY : x_ri_max, xemissn
1192 !
1193 USE mode_thermos
1194 !
1195 USE modi_surface_ri
1196 USE modi_surface_aero_cond
1197 USE modi_surface_cd
1198 !
1199 !
1200 IMPLICIT NONE
1201 !
1202 !* 0.1 declarations of arguments
1203 !
1204 REAL, INTENT(IN) :: PTSTEP, PSNOWDZMIN
1205 !
1206  CHARACTER(LEN=*), INTENT(IN) :: HSNOWRES ! type of sfc resistance
1207 ! DEFAULT=Louis (1979), standard ISBA
1208 ! method. Option to limit Ri number
1209 ! for very srtable conditions
1210 !
1211  CHARACTER(LEN=*), INTENT(IN) :: HIMPLICIT_WIND ! wind implicitation option
1212 ! ! 'OLD' = direct
1213 ! ! 'NEW' = Taylor serie, order 1
1214 !
1215 REAL, DIMENSION(:), INTENT(IN) :: PPEW_A_COEF, PPEW_B_COEF, &
1216  PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, &
1217  PPEQ_B_COEF
1218 ! PPEW_A_COEF = wind coefficient (m2s/kg)
1219 ! PPEW_B_COEF = wind coefficient (m/s)
1220 ! PPET_A_COEF = A-air temperature coefficient
1221 ! PPET_B_COEF = B-air temperature coefficient
1222 ! PPEQ_A_COEF = A-air specific humidity coefficient
1223 ! PPEQ_B_COEF = B-air specific humidity coefficient
1224 !
1225 REAL, DIMENSION(:), INTENT(IN) :: PZREF, PTS, PSNOWDZ1, PSNOWDZ2, &
1226  PRADSINK, PSNOWRHO, PSNOWLIQ, PSCAP, &
1227  PSCOND1, PSCOND2, &
1228  PZ0, PHPSNOW, &
1229  PALBT, PZ0EFF, PZ0H
1230 !
1231 REAL, DIMENSION(:), INTENT(IN) :: PSW_RAD, PLW_RAD, PTA, PQA, PPS, PRHOA
1232 !
1233 REAL, DIMENSION(:), INTENT(IN) :: PUREF, PEXNS, PEXNA, PDIRCOSZW, PVMOD
1234 !
1235 REAL, DIMENSION(:), INTENT(OUT) :: PTSTERM1, PTSTERM2, PEMIST, PRA, &
1236  PCT, PSFCFRZ, PCDSNOW, PCHSNOW, &
1237  PQSAT, PDQSAT, PRSRA
1238 !
1239 REAL, DIMENSION(:), INTENT(OUT) :: PUSTAR2_IC, &
1240  PPET_A_COEF_T, PPEQ_A_COEF_T, &
1241  PPET_B_COEF_T, PPEQ_B_COEF_T
1242 !
1243 REAL, DIMENSION(:), INTENT(OUT) :: PRI
1244 !
1245 !* 0.2 declarations of local variables
1246 !
1247 REAL, DIMENSION(SIZE(PTS)) :: ZAC, ZRI, ZCOND1, ZCOND2, &
1248  ZSCONDA, ZA, ZB, ZC, &
1249  ZCDN, ZSNOWDZM1, ZSNOWDZM2, &
1250  ZVMOD, ZUSTAR2, ZTS3, ZLVT, &
1251  Z_CCOEF
1252 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1253 !
1254 !-------------------------------------------------------------------------------
1255 !
1256 ! 1. New saturated specific humidity and derrivative:
1257 ! ---------------------------------------------------
1258 !
1259 IF (lhook) CALL dr_hook('SNOW3LEBUD',0,zhook_handle)
1260 !
1261 zri(:) = xundef
1262 !
1263 pqsat(:) = qsati(pts(:),pps(:))
1264 pdqsat(:) = dqsati(pts(:),pps(:),pqsat(:))
1265 !
1266 !
1267 ! 2. Surface properties:
1268 ! ----------------------
1269 ! It might be of interest to use snow-specific roughness
1270 ! or a temperature dependence on emissivity:
1271 ! but for now, use ISBA defaults.
1272 !
1273 pemist(:) = xemissn
1274 !
1275 ! 2. Computation of resistance and drag coefficient
1276 ! -------------------------------------------------
1277 !
1278  CALL surface_ri(pts, pqsat, pexns, pexna, pta, pqa, &
1279  pzref, puref, pdircoszw, pvmod, zri )
1280 !
1281 ! Simple adaptation of method by Martin and Lejeune (1998)
1282 ! to apply a lower limit to turbulent transfer coef
1283 ! by defining a maximum Richarson number for stable
1284 ! conditions:
1285 !
1286 IF(hsnowres=='RIL') THEN
1287  DO jj = 1,SIZE(zri)
1288  zri(jj) = min(x_ri_max,zri(jj))
1289  ENDDO
1290 ENDIF
1291 !
1292 pri(:)=zri(:)
1293 !
1294 ! Surface aerodynamic resistance for heat transfers
1295 !
1296  CALL surface_aero_cond(zri, pzref, puref, pvmod, pz0, pz0h, zac, pra, pchsnow)
1297 !
1298 ! For atmospheric model coupling:
1299 !
1300  CALL surface_cd(zri, pzref, puref, pz0eff, pz0h, pcdsnow, zcdn)
1301 !
1302 prsra(:) = prhoa(:) / pra(:)
1303 !
1304 !
1305 ! Modify flux-form implicit coupling coefficients:
1306 ! - wind components:
1307 !
1308 IF(himplicit_wind=='OLD')THEN
1309 ! old implicitation
1310  zustar2(:) = ( pcdsnow(:)*pvmod(:)*ppew_b_coef(:)) / &
1311  (1.0-prhoa(:)*pcdsnow(:)*pvmod(:)*ppew_a_coef(:))
1312 ELSE
1313 ! new implicitation
1314  zustar2(:) = (pcdsnow(:)*pvmod(:)*(2.*ppew_b_coef(:)-pvmod(:))) &
1315  / (1.0-2.0*prhoa(:)*pcdsnow(:)*pvmod(:)*ppew_a_coef(:))
1316 ENDIF
1317 !
1318 zvmod(:) = prhoa(:)*ppew_a_coef(:)*zustar2(:) + ppew_b_coef(:)
1319 zvmod(:) = max(zvmod(:),0.)
1320 !
1321 WHERE(ppew_a_coef(:)/= 0.)
1322  zustar2(:) = max( ( zvmod(:) - ppew_b_coef(:) ) / (prhoa(:)*ppew_a_coef(:)), 0.)
1323 ENDWHERE
1324 !
1325 ! implicit wind friction
1326 zustar2(:) = max(zustar2(:),0.)
1327 !
1328 pustar2_ic(:) = zustar2(:)
1329 !
1330 ! 3. Calculate linearized surface energy budget components:
1331 ! ---------------------------------------------------------
1332 ! To prevent numerical difficulties for very thin snow
1333 ! layers, limit the grid "thinness": this is important as
1334 ! layers become vanishing thin:
1335 !
1336 zsnowdzm1(:) = max(psnowdz1(:), psnowdzmin)
1337 zsnowdzm2(:) = max(psnowdz2(:), psnowdzmin)
1338 !
1339 ! Surface thermal inertia:
1340 !
1341 pct(:) = 1.0/(pscap(:)*zsnowdzm1(:))
1342 !
1343 ! Fraction of surface frozen (sublimation) with the remaining
1344 ! fraction being liquid (evaporation):
1345 ! NOTE: currently, for simplicity, assume no liquid water flux
1346 ! OFF: PSFCFRZ(:) = 1.0 - PSNOWLIQ(:)*XRHOLW/(ZSNOWDZM1(:)*PSNOWRHO(:))
1347 !
1348 psfcfrz(:) = 1.0
1349 !
1350 ! Thermal conductivity between uppermost and lower snow layers:
1351 !
1352 zcond1(:) = zsnowdzm1(:)/((zsnowdzm1(:)+zsnowdzm2(:))*pscond1(:))
1353 zcond2(:) = zsnowdzm2(:)/((zsnowdzm1(:)+zsnowdzm2(:))*pscond2(:))
1354 !
1355 zsconda(:) = 1.0/(zcond1(:)+zcond2(:))
1356 !
1357 ! Transform implicit coupling coefficients:
1358 ! Note, surface humidity is 100% over snow.
1359 !
1360 ! - specific humidity:
1361 !
1362 z_ccoef(:) = 1.0 - ppeq_a_coef(:)*prsra(:)
1363 !
1364 ppeq_a_coef_t(:) = - ppeq_a_coef(:)*prsra(:)*pdqsat(:)/z_ccoef(:)
1365 !
1366 ppeq_b_coef_t(:) = ( ppeq_b_coef(:) - ppeq_a_coef(:)*prsra(:)*(pqsat(:) - &
1367  pdqsat(:)*pts(:)) )/z_ccoef(:)
1368 !
1369 ! - air temperature:
1370 ! (assumes A and B correspond to potential T):
1371 !
1372 z_ccoef(:) = (1.0 - ppet_a_coef(:)*prsra(:))/pexna(:)
1373 !
1374 ppet_a_coef_t(:) = - ppet_a_coef(:)*prsra(:)/(pexns(:)*z_ccoef(:))
1375 !
1376 ppet_b_coef_t(:) = ppet_b_coef(:)/z_ccoef(:)
1377 !
1378 !
1379 ! Energy budget solution terms:
1380 !
1381 zts3(:) = pemist(:) * xstefan * pts(:)**3
1382 zlvt(:) = (1.-psfcfrz(:))*xlvtt + psfcfrz(:)*xlstt
1383 !
1384 za(:) = 1. / ptstep + pct(:) * (4. * zts3(:) + &
1385  prsra(:) * zlvt(:) * (pdqsat(:) - ppeq_a_coef_t(:)) &
1386  + prsra(:) * xcpd * ( (1./pexns(:))-(ppet_a_coef_t(:)/pexna(:)) ) &
1387  + (2.*zsconda(:)/(zsnowdzm2(:)+zsnowdzm1(:))) )
1388 !
1389 zb(:) = 1. / ptstep + pct(:) * (3. * zts3(:) + &
1390  prsra(:) * pdqsat(:) * zlvt(:) )
1391 !
1392 zc(:) = pct(:) * (prsra(:) * xcpd * ppet_b_coef_t(:)/pexna(:) + psw_rad(:) * &
1393  (1. - palbt(:)) + pemist(:)*plw_rad(:) - prsra(:) * &
1394  zlvt(:) * (pqsat(:)-ppeq_b_coef_t(:)) &
1395  + phpsnow(:) + pradsink(:) )
1396 !
1397 !
1398 ! Coefficients needed for implicit solution
1399 ! of linearized surface energy budget:
1400 !
1401 ptsterm2(:) = 2.*zsconda(:)*pct(:)/(za(:)*(zsnowdzm2(:)+zsnowdzm1(:)))
1402 !
1403 ptsterm1(:) = (pts(:)*zb(:) + zc(:))/za(:)
1404 IF (lhook) CALL dr_hook('SNOW3LEBUD',1,zhook_handle)
1405 !
1406 !-------------------------------------------------------------------------------
1407 !
1408 END SUBROUTINE snow3lebud
1409 !####################################################################
1410 !####################################################################
1411 !####################################################################
1412  SUBROUTINE snow3lsolvt(OMEB,PTSTEP,PSNOWDZMIN, &
1413  PSNOWDZ,PSCOND,PSCAP,PTG, &
1414  PSOILCOND,PD_G, &
1415  PRADSINK,PCT,PTERM1,PTERM2, &
1416  PPET_A_COEF_T,PPEQ_A_COEF_T, &
1417  PPET_B_COEF_T,PPEQ_B_COEF_T, &
1418  PTA_IC, PQA_IC, &
1419  PGRNDFLUX,PGRNDFLUXO,PSNOWTEMP, &
1420  PSNOWFLUX )
1421 !
1422 !! PURPOSE
1423 !! -------
1424 ! This subroutine solves the 1-d diffusion of 'ZSNOWTEMP' using a
1425 ! layer averaged set of equations which are time differenced
1426 ! using the backward-difference scheme (implicit).
1427 ! The eqs are solved rapidly by taking advantage of the
1428 ! fact that the matrix is tridiagonal. This is a very
1429 ! general routine and can be used for the 1-d diffusion of any
1430 ! quantity as long as the diffusity is not a function of the
1431 ! quantity being diffused. Aaron Boone 8-98. Soln to the eq:
1432 !
1433 ! c dQ d k dQ dS
1434 ! -- = -- -- - --
1435 ! dt dx dx dx
1436 !
1437 ! where k = k(x) (thermal conductivity), c = c(x) (heat capacity)
1438 ! as an eg. for temperature/heat eq. S is a sink (-source) term.
1439 ! Diffusivity is k/c
1440 !
1441 !
1442 USE modd_csts,ONLY : xtt
1443 !
1444 USE modi_tridiag_ground
1445 !
1446 IMPLICIT NONE
1447 !
1448 !* 0.1 declarations of arguments
1449 !
1450 LOGICAL, INTENT(IN) :: OMEB
1451 !
1452 REAL, INTENT(IN) :: PTSTEP, PSNOWDZMIN
1453 !
1454 REAL, DIMENSION(:), INTENT(IN) :: PTG, PSOILCOND, PD_G, &
1455  PCT, PTERM1, PTERM2
1456 
1457 !
1458 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWDZ, PSCOND, PSCAP, &
1459  PRADSINK
1460 !
1461 REAL, DIMENSION(:), INTENT(IN) :: PPET_A_COEF_T, PPEQ_A_COEF_T, &
1462  PPET_B_COEF_T, PPEQ_B_COEF_T
1463 !
1464 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWTEMP
1465 !
1466 REAL, DIMENSION(:), INTENT(OUT) :: PGRNDFLUX, PGRNDFLUXO, PSNOWFLUX, &
1467  PTA_IC, PQA_IC
1468 !
1469 !
1470 !* 0.2 declarations of local variables
1471 !
1472 !
1473 INTEGER :: JJ, JI
1474 !
1475 INTEGER :: INI
1476 INTEGER :: INLVLS
1477 !
1478 REAL, DIMENSION(SIZE(PTG)) :: ZSNOWTEMP_DELTA
1479 !
1480 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)) :: ZSNOWTEMP, ZDTERM, ZCTERM, &
1481  ZFRCV, ZAMTRX, ZBMTRX, &
1482  ZCMTRX
1483 !
1484 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)) :: ZWORK1, ZWORK2, ZDZDIF, &
1485  ZSNOWDZM
1486 !
1487 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)-1) :: ZSNOWTEMP_M, &
1488  ZFRCV_M, ZAMTRX_M, &
1489  ZBMTRX_M, ZCMTRX_M
1490 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1491 !-------------------------------------------------------------------------------
1492 !
1493 ! 0. Initialize:
1494 ! ------------------
1495 !
1496 IF (lhook) CALL dr_hook('SNOW3LSOLVT',0,zhook_handle)
1497 zsnowtemp(:,:) = psnowtemp(:,:)
1498 ini = SIZE(psnowdz(:,:),1)
1499 inlvls = SIZE(psnowdz(:,:),2)
1500 !
1501 !
1502 ! 1. Calculate tri-diagnoal matrix coefficients:
1503 ! ----------------------------------------------
1504 ! For heat transfer, assume a minimum grid
1505 ! thickness (to prevent numerical
1506 ! problems for very thin snow cover):
1507 !
1508 zsnowdzm(:,:) = max(psnowdz(:,:), psnowdzmin)
1509 !
1510 DO jj=1,inlvls-1
1511  DO ji=1,ini
1512  zdzdif(ji,jj) = 0.5*(zsnowdzm(ji,jj)+zsnowdzm(ji,jj+1))
1513  zwork1(ji,jj) = zsnowdzm(ji,jj )/(2.0*zdzdif(ji,jj)*pscond(ji,jj ))
1514  zwork2(ji,jj) = zsnowdzm(ji,jj+1)/(2.0*zdzdif(ji,jj)*pscond(ji,jj+1))
1515  ENDDO
1516 ENDDO
1517 !
1518 zdzdif(:,inlvls) = 0.5*(zsnowdzm(:,inlvls)+pd_g(:))
1519 zwork1(:,inlvls) = zsnowdzm(:,inlvls)/(2.0*zdzdif(:,inlvls)*pscond(:,inlvls))
1520 zwork2(:,inlvls) = pd_g(: )/(2.0*zdzdif(:,inlvls)*psoilcond(: ))
1521 !
1522 zdterm(:,:) = 1.0/(zdzdif(:,:)*(zwork1(:,:)+zwork2(:,:)))
1523 !
1524 zcterm(:,:) = pscap(:,:)*zsnowdzm(:,:)/ptstep
1525 !
1526 ! 2. Set up tri-diagonal matrix
1527 ! -----------------------------
1528 !
1529 ! Upper BC
1530 !
1531 zamtrx(:,1) = 0.0
1532 zbmtrx(:,1) = 1./(pct(:)*ptstep)
1533 zcmtrx(:,1) = -pterm2(:)*zbmtrx(:,1)
1534 zfrcv(:,1) = pterm1(:)*zbmtrx(:,1)
1535 !
1536 !
1537 ! Interior Grid
1538 !
1539 DO jj=2,inlvls-1
1540  DO ji=1,ini
1541  zamtrx(ji,jj) = -zdterm(ji,jj-1)
1542  zbmtrx(ji,jj) = zcterm(ji,jj) + zdterm(ji,jj-1) + zdterm(ji,jj)
1543  zcmtrx(ji,jj) = -zdterm(ji,jj)
1544  zfrcv(ji,jj) = zcterm(ji,jj)*psnowtemp(ji,jj) - (pradsink(ji,jj-1)-pradsink(ji,jj))
1545  ENDDO
1546 ENDDO
1547 !
1548 !Lower BC
1549 !
1550 zamtrx(:,inlvls) = -zdterm(:,inlvls-1)
1551 zbmtrx(:,inlvls) = zcterm(:,inlvls) + zdterm(:,inlvls-1) + &
1552  zdterm(:,inlvls)
1553 zcmtrx(:,inlvls) = 0.0
1554 zfrcv(:,inlvls) = zcterm(:,inlvls)*psnowtemp(:,inlvls) + &
1555  zdterm(:,inlvls)*ptg(:) &
1556  - (pradsink(:,inlvls-1)-pradsink(:,inlvls))
1557 !
1558 ! - - -------------------------------------------------
1559 !
1560 ! 4. Compute solution vector
1561 ! --------------------------
1562 !
1563  CALL tridiag_ground(zamtrx,zbmtrx,zcmtrx,zfrcv,zsnowtemp)
1564 !
1565 ! Heat flux between surface and 2nd snow layers: (W/m2)
1566 !
1567 psnowflux(:) = zdterm(:,1)*(zsnowtemp(:,1) - zsnowtemp(:,2))
1568 !
1569 !
1570 ! 5. Snow melt case
1571 ! -----------------
1572 ! If melting in uppermost layer, assume surface layer
1573 ! temperature at freezing point and re-evaluate lower
1574 ! snowpack temperatures. This is done as it is most likely
1575 ! most signigant melting will occur within a time step in surface layer.
1576 ! Surface energy budget (and fluxes) will
1577 ! be re-calculated (outside of this routine).
1578 !
1579 ! NOTE: if MEB is active, then surface fluxes have been defined outside
1580 ! of the snow routine and have been adjusted such that they are evaluated
1581 ! at a snow surface temperature no greater than Tf. Thus, the implicit surface temperature
1582 ! will likely never greatly exceed Tf (before melt computed and they are adjusted to Tf)
1583 ! so we can skip the next block of code when MEB is active.
1584 !
1585 IF(.NOT.omeb)THEN
1586 !
1587  zamtrx_m(:,1) = 0.0
1588  zbmtrx_m(:,1) = zcterm(:,2) + zdterm(:,1) + zdterm(:,2)
1589  zcmtrx_m(:,1) = -zdterm(:,2)
1590  zfrcv_m(:,1) = zcterm(:,2)*psnowtemp(:,2) + xtt*zdterm(:,1) - &
1591  (pradsink(:,1)-pradsink(:,2))
1592 !
1593  DO jj=2,inlvls-1
1594  DO ji=1,ini
1595  zamtrx_m(ji,jj) = zamtrx(ji,jj+1)
1596  zbmtrx_m(ji,jj) = zbmtrx(ji,jj+1)
1597  zcmtrx_m(ji,jj) = zcmtrx(ji,jj+1)
1598  zfrcv_m(ji,jj) = zfrcv(ji,jj+1)
1599  zsnowtemp_m(ji,jj) = psnowtemp(ji,jj+1)
1600  ENDDO
1601  ENDDO
1602 !
1603  CALL tridiag_ground(zamtrx_m,zbmtrx_m,zcmtrx_m,zfrcv_m,zsnowtemp_m)
1604 !
1605 ! If melting for 2 consecuative time steps, then replace current T-profile
1606 ! with one assuming T=Tf in surface layer:
1607 !
1608  zsnowtemp_delta(:) = 0.0
1609 !
1610  WHERE(zsnowtemp(:,1) > xtt .AND. psnowtemp(:,1) == xtt)
1611  psnowflux(:) = zdterm(:,1)*(xtt - zsnowtemp_m(:,1))
1612  zsnowtemp_delta(:) = 1.0
1613  END WHERE
1614 !
1615  DO jj=2,inlvls
1616  DO ji=1,ini
1617  zsnowtemp(ji,jj) = zsnowtemp_delta(ji)*zsnowtemp_m(ji,jj-1) &
1618  + (1.0-zsnowtemp_delta(ji))*zsnowtemp(ji,jj)
1619  ENDDO
1620  ENDDO
1621 !
1622 ENDIF
1623 !
1624 !
1625 ! 6. Lower boundary flux:
1626 ! -----------------------
1627 ! NOTE: evaluate this term assuming the snow layer
1628 ! can't exceed the freezing point as this adjustment
1629 ! is made in melting routine. Then must adjust temperature
1630 ! to conserve energy.
1631 ! NOTE: if MEB used, surface fluxes are imposed and a test semi-implicit
1632 ! ground-snow flux has already been estimated and snow-free ground fluxes
1633 ! are generally weaker (and snow fractions are generally set to go to unity
1634 ! faster than for the composite soil-veg case), thus this correction
1635 ! is not as essential and is off.
1636 !
1637 pgrndfluxo(:) = zdterm(:,inlvls)*(zsnowtemp(:,inlvls) -ptg(:))
1638 IF(omeb)THEN
1639  pgrndflux(:) = pgrndfluxo(:)
1640 ELSE
1641  pgrndflux(:) = zdterm(:,inlvls)*(min(xtt,zsnowtemp(:,inlvls))-ptg(:))
1642 ENDIF
1643 !
1644 zsnowtemp(:,inlvls) = zsnowtemp(:,inlvls) + (pgrndfluxo(:)-pgrndflux(:))/zcterm(:,inlvls)
1645 !
1646 ! 7. Update temperatute profile in time:
1647 ! --------------------------------------
1648 !
1649 psnowtemp(:,:) = zsnowtemp(:,:)
1650 !
1651 !
1652 ! 8. Compute new (implicit) air T and specific humidity
1653 ! -----------------------------------------------------
1654 !
1655 IF(.NOT.omeb)THEN
1656 !
1657  pta_ic(:) = ppet_b_coef_t(:) + ppet_a_coef_t(:)* psnowtemp(:,1)
1658 !
1659  pqa_ic(:) = ppeq_b_coef_t(:) + ppeq_a_coef_t(:)* psnowtemp(:,1)
1660 !
1661 ENDIF
1662 !
1663 IF (lhook) CALL dr_hook('SNOW3LSOLVT',1,zhook_handle)
1664 !
1665 !
1666 END SUBROUTINE snow3lsolvt
1667 !####################################################################
1668 !####################################################################
1669 !####################################################################
1670  SUBROUTINE snow3lmelt(PTSTEP,PSCAP,PSNOWTEMP,PSNOWDZ, &
1671  PSNOWRHO,PSNOWLIQ,PMELTXS )
1672 !
1673 !
1674 !! PURPOSE
1675 !! -------
1676 ! Calculate snow melt (resulting from surface fluxes, ground fluxes,
1677 ! or internal shortwave radiation absorbtion). It is used to
1678 ! augment liquid water content, maintain temperatures
1679 ! at or below freezing, and possibly reduce the mass
1680 ! or compact the layer(s).
1681 !
1682 !
1683 USE modd_csts,ONLY : xtt, xlmtt, xrholw, xrholi
1684 !
1685 USE mode_snow3l
1686 !
1687 IMPLICIT NONE
1688 !
1689 !* 0.1 declarations of arguments
1690 !
1691 REAL, INTENT(IN) :: PTSTEP
1692 !
1693 REAL, DIMENSION(:,:), INTENT(IN) :: PSCAP
1694 !
1695 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ, PSNOWTEMP, PSNOWRHO, &
1696  PSNOWLIQ
1697 !
1698 REAL, DIMENSION(:), INTENT(OUT) :: PMELTXS
1699 !
1700 !
1701 !* 0.2 declarations of local variables
1702 !
1703 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZPHASE, ZCMPRSFACT, &
1704  ZSNOWLWE, ZWHOLDMAX, &
1705  ZSNOWMELT, ZSNOWTEMP, &
1706  ZMELTXS
1707 !
1708 INTEGER :: JWRK, JI ! loop counter
1709 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1710 !-------------------------------------------------------------------------------
1711 !
1712 ! 0. Initialize:
1713 ! ---------------------------
1714 !
1715 IF (lhook) CALL dr_hook('SNOW3LMELT',0,zhook_handle)
1716 zphase(:,:) = 0.0
1717 zcmprsfact(:,:) = 0.0
1718 zsnowlwe(:,:) = 0.0
1719 zwholdmax(:,:) = 0.0
1720 zsnowmelt(:,:) = 0.0
1721 zsnowtemp(:,:) = 0.0
1722 zmeltxs(:,:) = 0.0
1723 !
1724 ! 1. Determine amount of melt in each layer:
1725 ! ------------------------------------------
1726 !
1727 WHERE(psnowdz > 0.0)
1728 !
1729 ! Total Liquid equivalent water content of snow (m):
1730 !
1731  zsnowlwe(:,:) = psnowrho(:,:)*psnowdz(:,:)/xrholw
1732 !
1733 ! Melt snow if excess energy and snow available:
1734 ! Phase change (J/m2)
1735 !
1736  zphase(:,:) = min(pscap(:,:)*max(0.0, psnowtemp(:,:) - xtt)* &
1737  psnowdz(:,:), &
1738  max(0.0,zsnowlwe(:,:)-psnowliq(:,:))*xlmtt*xrholw)
1739 !
1740 !
1741 ! Update snow liq water content and temperature if melting:
1742 ! liquid water available for next layer from melting of snow
1743 ! which is assumed to be leaving the current layer (m):
1744 !
1745  zsnowmelt(:,:) = zphase(:,:)/(xlmtt*xrholw)
1746 !
1747 ! Cool off snow layer temperature due to melt:
1748 !
1749  zsnowtemp(:,:) = psnowtemp(:,:) - zphase(:,:)/(pscap(:,:)*psnowdz(:,:))
1750 !
1751  psnowtemp(:,:) = min(xtt, zsnowtemp(:,:))
1752 !
1753  zmeltxs(:,:) = (zsnowtemp(:,:)-psnowtemp(:,:))*pscap(:,:)*psnowdz(:,:)
1754 !
1755 END WHERE
1756 !
1757 ! Loss of snowpack depth: (m) and liquid equiv (m):
1758 ! Compression factor for melt loss: this decreases
1759 ! layer thickness and increases density thereby leaving
1760 ! total SWE constant. NOTE: All melt water
1761 ! in excess of the holding capacity is NOT used
1762 ! for compression, rather it decreases the layer
1763 ! thickness ONLY, causing a loss of SWE (outside
1764 ! of this routine).
1765 !
1766 zwholdmax(:,:) = snow3lhold(psnowrho,psnowdz)
1767 !
1768 WHERE(psnowdz > 0.0)
1769 !
1770  zcmprsfact(:,:) = (zsnowlwe(:,:)-min(psnowliq(:,:)+zsnowmelt(:,:), &
1771  zwholdmax(:,:)))/ &
1772  (zsnowlwe(:,:)-min(psnowliq(:,:),zwholdmax(:,:)))
1773 !
1774  psnowdz(:,:) = psnowdz(:,:)*zcmprsfact(:,:)
1775  psnowrho(:,:) = zsnowlwe(:,:)*xrholw/psnowdz(:,:)
1776 !
1777 ! Make sure maximum density is not surpassed! If it is, lower the density
1778 ! and increase the snow thickness accordingly:
1779 
1780  zcmprsfact(:,:) = max(xrholi, psnowrho(:,:))/xrholi
1781  psnowdz(:,:) = psnowdz(:,:)*zcmprsfact(:,:)
1782  psnowrho(:,:) = zsnowlwe(:,:)*xrholw/psnowdz(:,:)
1783 !
1784 !
1785 ! 2. Add snow melt to current snow liquid water content:
1786 ! ------------------------------------------------------
1787 !
1788  psnowliq(:,:) = psnowliq(:,:) + zsnowmelt(:,:)
1789 !
1790 END WHERE
1791 !
1792 ! 3. Excess heat from melting
1793 ! ---------------------------
1794 ! use it to warm underlying ground/vegetation layer to conserve energy
1795 !
1796 pmeltxs(:) = 0.
1797 DO jwrk = 1, SIZE(zmeltxs,2)
1798  DO ji = 1, SIZE(zmeltxs,1)
1799  pmeltxs(ji) = pmeltxs(ji) + zmeltxs(ji,jwrk)
1800  ENDDO
1801 ENDDO
1802 pmeltxs(:) = pmeltxs(:) / ptstep ! (W/m2)
1803 !
1804 IF (lhook) CALL dr_hook('SNOW3LMELT',1,zhook_handle)
1805 !
1806 !
1807 !
1808 END SUBROUTINE snow3lmelt
1809 !####################################################################
1810 !####################################################################
1811 !####################################################################
1812  SUBROUTINE snow3lrefrz(PTSTEP,PRR, &
1813  PSNOWRHO,PSNOWTEMP,PSNOWDZ,PSNOWLIQ, &
1814  PTHRUFAL )
1815 !
1816 !
1817 !! PURPOSE
1818 !! -------
1819 ! Calculate any freezing/refreezing of liquid water in the snowpack.
1820 ! Also, calculate liquid water transmission and snow runoff.
1821 ! Water flow causes layer thickness reduction and can cause
1822 ! rapid densification of a layer.
1823 !
1824 !
1825 USE modd_csts, ONLY : xtt, xlmtt, xrholw
1826 USE modd_snow_par, ONLY : xsnowdmin
1827 !
1828 USE mode_snow3l
1829 !
1830 IMPLICIT NONE
1831 !
1832 !* 0.1 declarations of arguments
1833 !
1834 REAL, INTENT(IN) :: PTSTEP
1835 !
1836 REAL, DIMENSION(:), INTENT(IN) :: PRR
1837 !
1838 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ, PSNOWTEMP, PSNOWLIQ, PSNOWRHO
1839 !
1840 REAL, DIMENSION(:), INTENT(INOUT) :: PTHRUFAL
1841 !
1842 !
1843 !
1844 !* 0.2 declarations of local variables
1845 !
1846 INTEGER :: JJ, JI
1847 !
1848 INTEGER :: INI
1849 INTEGER :: INLVLS
1850 !
1851 REAL, DIMENSION(SIZE(PRR)) :: ZPCPXS, ZTOTWCAP, ZRAINFALL
1852 !
1853 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZFLOWLIQ, ZWORK, &
1854  ZSNOWLIQ, ZSNOWRHO, &
1855  ZWHOLDMAX, ZSNOWDZ, &
1856  ZSNOWTEMP, ZSCAP, &
1857  ZSNOWHEAT
1858 !
1859 REAL, DIMENSION(SIZE(PSNOWRHO,1),0:SIZE(PSNOWRHO,2)):: ZFLOWLIQT
1860 !
1861 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1862 !
1863 !-------------------------------------------------------------------------------
1864 !
1865 ! 0. Initialize:
1866 ! --------------
1867 !
1868 IF (lhook) CALL dr_hook('SNOW3LREFRZ',0,zhook_handle)
1869 !
1870 zsnowrho(:,:) = psnowrho(:,:)
1871 zsnowliq(:,:) = psnowliq(:,:)
1872 zsnowtemp(:,:) = psnowtemp(:,:)
1873 ini = SIZE(psnowdz(:,:),1)
1874 inlvls = SIZE(psnowdz(:,:),2)
1875 !
1876 !
1877 ! 1. Refreeze due to heat transfer
1878 ! -----------------------------
1879 ! Freeze liquid water in any layers which cooled due
1880 ! to heat transfer. First, update H and re-diagnose (update) T and Wl:
1881 !
1882 zscap(:,:) = snow3lscap(zsnowrho)
1883 !
1884 zsnowheat(:,:) = psnowdz(:,:)*( zscap(:,:)*(zsnowtemp(:,:)-xtt) &
1885  - xlmtt*zsnowrho(:,:) ) + xlmtt*xrholw*zsnowliq(:,:)
1886 !
1887 zsnowtemp(:,:) = xtt + ( ((zsnowheat(:,:)/max(psnowdz(:,:),xsnowdmin/inlvls)) &
1888  + xlmtt*zsnowrho(:,:))/zscap(:,:) )
1889 !
1890 zsnowliq(:,:) = max(0.0,zsnowtemp(:,:)-xtt)*zscap(:,:)*psnowdz(:,:)/(xlmtt*xrholw)
1891 !
1892 zsnowtemp(:,:) = min(xtt,zsnowtemp(:,:))
1893 !
1894 !
1895 ! 2. Reduce thickness due to snowmelt in excess of holding capacity
1896 ! --------------------------------------------------------------
1897 ! Any water in excess of the
1898 ! Maximum holding space for liquid water
1899 ! amount is drained into next layer down.
1900 ! Loss of water due to snowmelt causes a reduction
1901 ! in snow layer mass by a reduction in thickness. Owing to a consistent
1902 ! decrease in both liq and thickness (and the fact that T=Tf when liquid present),
1903 ! enthalpy is conserved.
1904 !
1905 zwholdmax(:,:) = snow3lhold(psnowrho,psnowdz)
1906 !
1907 zflowliq(:,:) = max(0.,zsnowliq(:,:)-zwholdmax(:,:))
1908 !
1909 zsnowliq(:,:) = zsnowliq(:,:) - zflowliq(:,:)
1910 !
1911 zsnowdz(:,:) = psnowdz(:,:) - zflowliq(:,:)*xrholw/zsnowrho(:,:)
1912 !
1913 zsnowdz(:,:) = max(0.0, zsnowdz(:,:)) ! to prevent possible very small
1914 ! negative values (machine prescision
1915 ! as snow vanishes
1916 !
1917 !
1918 ! 3. Liquid water flow: liquid precipitation and meltwater
1919 ! -----------------------------------------------------
1920 !
1921 ! Rainfall flowing into uppermost snow layer:
1922 ! If rainfall is excessive enough (or layers thin enough)
1923 ! it is simply routed directly to runoff: First calculate
1924 ! the total snow pack available liquid water holding capacity:
1925 !
1926 ztotwcap(:) = 0.
1927 DO jj=1,inlvls
1928  DO ji=1,ini
1929  ztotwcap(ji) = ztotwcap(ji) + zwholdmax(ji,jj)
1930  ENDDO
1931 ENDDO
1932 !
1933 ! Rain entering snow (m):
1934 !
1935 zrainfall(:) = prr(:)*ptstep/xrholw ! rainfall (m)
1936 !
1937 zflowliqt(:,0)= min(zrainfall(:),ztotwcap(:))
1938 !
1939 ! Rain assumed to directly pass through the pack to runoff (m):
1940 !
1941 zpcpxs(:) = zrainfall(:) - zflowliqt(:,0)
1942 !
1943 DO jj=1,inlvls
1944  DO ji=1,ini
1945  zflowliqt(ji,jj) = zflowliq(ji,jj)
1946  ENDDO
1947 ENDDO
1948 !
1949 !
1950 ! Thickness is maintained during water through-flow,
1951 ! so that mass transfer is represented by
1952 ! density changes: NOTE a maximum density
1953 ! is assumed (XRHOSMAX_ES) so that all flow
1954 ! which would result in densities exceeding
1955 ! this limit are routed to next layer down.
1956 ! First test for saturation, then
1957 ! rout excess water down to next layer down
1958 ! and repeat calculation. Net gain in liquid (mass) is
1959 ! translated into a density increase:
1960 !
1961 zflowliq(:,:) = 0.0 ! clear this array for work
1962 psnowliq(:,:) = zsnowliq(:,:) ! reset liquid water content
1963 !
1964 DO jj=1,inlvls
1965  DO ji=1,ini
1966  zsnowliq(ji,jj) = zsnowliq(ji,jj) + zflowliqt(ji,jj-1)
1967  zflowliq(ji,jj) = max(0.0, zsnowliq(ji,jj)-zwholdmax(ji,jj))
1968  zsnowliq(ji,jj) = zsnowliq(ji,jj) - zflowliq(ji,jj)
1969  zflowliqt(ji,jj) = zflowliqt(ji,jj) + zflowliq(ji,jj)
1970  ENDDO
1971 ENDDO
1972 !
1973 zwork(:,:) = max(xsnowdmin/inlvls,zsnowdz(:,:))
1974 zsnowrho(:,:) = zsnowrho(:,:)+(zsnowliq(:,:)-psnowliq(:,:))*xrholw/zwork(:,:)
1975 zscap(:,:) = snow3lscap(zsnowrho(:,:))
1976 zsnowtemp(:,:) = xtt +(((zsnowheat(:,:)/zwork(:,:))+xlmtt*zsnowrho(:,:))/zscap(:,:))
1977 zsnowliq(:,:) = max(0.0,zsnowtemp(:,:)-xtt)*zscap(:,:)*zsnowdz(:,:)/(xlmtt*xrholw)
1978 zsnowtemp(:,:) = min(xtt,zsnowtemp(:,:))
1979 !
1980 ! Any remaining throughflow after freezing is available to
1981 ! the soil for infiltration or surface runoff (m).
1982 ! I.E. This is the amount of water leaving the snowpack:
1983 ! Rate water leaves the snowpack [kg/(m2 s)]:
1984 !
1985 pthrufal(:) = pthrufal(:) + zflowliqt(:,inlvls)
1986 !
1987 ! Add excess rain (rain which flowed directly through the snow
1988 ! due to saturation):
1989 !
1990 pthrufal(:) = (pthrufal(:) + zpcpxs(:))*xrholw/ptstep
1991 !
1992 ! 4. Update thickness and density and any freezing:
1993 ! ----------------------------------------------
1994 !
1995 psnowtemp(:,:)= zsnowtemp(:,:)
1996 psnowdz(:,:) = zsnowdz(:,:)
1997 psnowrho(:,:) = zsnowrho(:,:)
1998 psnowliq(:,:) = zsnowliq(:,:)
1999 !
2000 IF (lhook) CALL dr_hook('SNOW3LREFRZ',1,zhook_handle)
2001 !
2002 !-------------------------------------------------------------------------------
2003 !
2004 END SUBROUTINE snow3lrefrz
2005 !####################################################################
2006 !####################################################################
2007 !####################################################################
2008  SUBROUTINE snow3lflux(PSNOWTEMP,PSNOWDZ,PEXNS,PEXNA, &
2009  PUSTAR2_IC, &
2010  PTSTEP,PALBT,PSW_RAD,PEMIST,PLWUPSNOW, &
2011  PLW_RAD,PLWNETSNOW, &
2012  PTA,PSFCFRZ,PQA,PHPSNOW, &
2013  PSNOWTEMPO1,PSNOWFLUX,PCT,PRADSINK, &
2014  PQSAT,PDQSAT,PRSRA, &
2015  PRN,PH,PGFLUX,PLES3L,PLEL3L,PEVAP, &
2016  PUSTAR,OSFCMELT )
2017 !
2018 !
2019 !! PURPOSE
2020 !! -------
2021 ! Calculate the surface fluxes (atmospheric/surface).
2022 ! (Noilhan and Planton 1989; Noilhan and Mahfouf 1996)
2023 !
2024 !
2025 USE modd_csts,ONLY : xstefan, xcpd, xlstt, xlvtt, xtt
2026 !
2027 USE mode_thermos
2028 !
2029 IMPLICIT NONE
2030 !
2031 !* 0.1 declarations of arguments
2032 !
2033 REAL, INTENT(IN) :: PTSTEP
2034 !
2035 REAL, DIMENSION(:), INTENT(IN) :: PSNOWDZ, PSNOWTEMPO1, PSNOWFLUX, PCT, &
2036  PRADSINK, PEXNS, PEXNA
2037 !
2038 REAL, DIMENSION(:), INTENT(IN) :: PALBT, PSW_RAD, PEMIST, PLW_RAD, &
2039  PTA, PSFCFRZ, PQA, &
2040  PHPSNOW, PQSAT, PDQSAT, PRSRA, &
2041  PUSTAR2_IC
2042 !
2043 REAL, DIMENSION(:), INTENT(INOUT) :: PSNOWTEMP
2044 !
2045 REAL, DIMENSION(:), INTENT(OUT) :: PRN, PH, PGFLUX, PLES3L, PLEL3L, &
2046  PEVAP, PLWUPSNOW, PUSTAR, &
2047  PLWNETSNOW
2048 !
2049 LOGICAL, DIMENSION(:), INTENT(OUT) :: OSFCMELT
2050 !
2051 !
2052 !* 0.2 declarations of local variables
2053 !
2054 REAL, DIMENSION(SIZE(PSNOWDZ)) :: ZEVAPC, ZLE, ZSNOWTEMP, ZSMSNOW, ZGFLUX, &
2055  ZDELTAT, ZSNOWTO3
2056 REAL(KIND=JPRB) :: ZHOOK_HANDLE
2057 !
2058 !-------------------------------------------------------------------------------
2059 !
2060 ! 0. Initialize:
2061 ! --------------
2062 !
2063 IF (lhook) CALL dr_hook('SNOW3LFLUX',0,zhook_handle)
2064 zsnowtemp(:) = psnowtemp(:)
2065 zle(:) = 0.0
2066 zsmsnow(:) = 0.0
2067 zgflux(:) = 0.0
2068 !
2069 osfcmelt(:) = .false.
2070 !
2071 zsnowto3(:) = psnowtempo1(:) ** 3 ! to save some CPU time, store this
2072 !
2073 ! 1. Flux calculations when melt not occuring at surface (W/m2):
2074 ! --------------------------------------------------------------
2075 !
2076 !
2077 zdeltat(:) = psnowtemp(:) - psnowtempo1(:) ! surface T time change:
2078 !
2079 plwupsnow(:) = pemist(:) * xstefan * zsnowto3(:)*( psnowtempo1(:) + 4.* zdeltat(:) )
2080 !
2081 plwnetsnow(:)= pemist(:) * plw_rad(:) - plwupsnow(:)
2082 !
2083 prn(:) = (1. - palbt(:)) * psw_rad(:) + plwnetsnow(:)
2084 !
2085 ph(:) = prsra(:) * xcpd * (psnowtemp(:)/pexns(:) - pta(:)/pexna(:))
2086 !
2087 zevapc(:) = prsra(:) * ( (pqsat(:) - pqa(:)) + pdqsat(:)*zdeltat(:) )
2088 !
2089 ples3l(:) = psfcfrz(:) * xlstt * zevapc(:)
2090 !
2091 plel3l(:) = (1.-psfcfrz(:))* xlvtt * zevapc(:)
2092 !
2093 zle(:) = ples3l(:) + plel3l(:)
2094 !
2095 pgflux(:) = prn(:) - ph(:) - zle(:) + phpsnow(:)
2096 !
2097 !
2098 ! 2. Initial melt adjustment
2099 ! --------------------------
2100 ! If energy avalabile to melt snow, then recalculate fluxes
2101 ! at the freezing point and add residual heat to layer
2102 ! average heat.
2103 !
2104 ! A) If temperature change is > 0 and passes freezing point this timestep,
2105 ! then recalculate fluxes at freezing point and excess energy
2106 ! will be used outside of this routine to change snow heat content:
2107 !
2108 WHERE (psnowtemp > xtt .AND. psnowtempo1 < xtt)
2109 !
2110  osfcmelt(:)= .true.
2111 !
2112  zdeltat(:) = xtt - psnowtempo1(:)
2113 !
2114  plwupsnow(:) = pemist(:) * xstefan * zsnowto3(:)*( psnowtempo1(:) + 4.* zdeltat(:) )
2115 !
2116  plwnetsnow(:)= pemist(:) * plw_rad(:) - plwupsnow(:)
2117 !
2118  prn(:) = (1. - palbt(:)) * psw_rad(:) + plwnetsnow(:)
2119 !
2120  ph(:) = prsra(:) * xcpd * (xtt/pexns(:) - pta(:)/pexna(:))
2121 !
2122  zevapc(:) = prsra(:) * ( (pqsat(:) - pqa(:)) + pdqsat(:)*zdeltat(:) )
2123 !
2124  ples3l(:) = psfcfrz(:) * xlstt * zevapc(:)
2125 !
2126  plel3l(:) = (1.-psfcfrz(:))* xlvtt * zevapc(:)
2127 !
2128  zle(:) = ples3l(:) + plel3l(:)
2129 !
2130  zgflux(:) = prn(:) - ph(:) - zle(:) + phpsnow(:)
2131 !
2132  zsmsnow(:) = pgflux(:) - zgflux(:)
2133 !
2134  pgflux(:) = zgflux(:)
2135 !
2136 ! This will be used to change heat content of snow:
2137 !
2138  zsnowtemp(:) = psnowtemp(:) - zsmsnow(:)*ptstep*pct(:)
2139 !
2140 END WHERE
2141 !
2142 ! 3. Ongoing melt adjustment: explicit solution
2143 ! ---------------------------------------------
2144 ! If temperature change is 0 and at freezing point this timestep,
2145 ! then recalculate fluxes and surface temperature *explicitly*
2146 ! as this is *exact* for snow at freezing point (Brun, Martin)
2147 !
2148 WHERE(psnowtemp(:) > xtt .AND. psnowtempo1(:) >= xtt)
2149 !
2150  osfcmelt(:) = .true.
2151 !
2152  plwupsnow(:) = pemist(:) * xstefan * (xtt ** 4)
2153 !
2154  plwnetsnow(:)= pemist(:) * plw_rad(:) - plwupsnow(:)
2155 !
2156  prn(:) = (1. - palbt(:)) * psw_rad(:) + plwnetsnow(:)
2157 !
2158  ph(:) = prsra(:) * xcpd * (xtt/pexns(:) - pta(:)/pexna(:))
2159 !
2160  zevapc(:) = prsra(:) * (pqsat(:) - pqa(:))
2161 !
2162  ples3l(:) = psfcfrz(:) * xlstt * zevapc(:)
2163 !
2164  plel3l(:) = (1.-psfcfrz(:))* xlvtt * zevapc(:)
2165 !
2166  zle(:) = ples3l(:) + plel3l(:)
2167 !
2168  pgflux(:) = prn(:) - ph(:) - zle(:) + phpsnow(:)
2169 !
2170  zsnowtemp(:) = xtt + ptstep*pct(:)*(pgflux(:) + pradsink(:) - psnowflux(:))
2171 !
2172 END WHERE
2173 !
2174 ! 4. Update surface temperature:
2175 ! ------------------------------
2176 !
2177 psnowtemp(:) = zsnowtemp(:)
2178 !
2179 ! 5. Final evaporative flux (kg/m2/s)
2180 !
2181 pevap(:) = zevapc(:)
2182 !
2183 ! 6. Friction velocity
2184 ! --------------------
2185 !
2186 pustar(:) = sqrt(pustar2_ic(:))
2187 !
2188 IF (lhook) CALL dr_hook('SNOW3LFLUX',1,zhook_handle)
2189 !
2190 !-------------------------------------------------------------------------------
2191 !
2192 END SUBROUTINE snow3lflux
2193 !####################################################################
2194 !####################################################################
2195 !####################################################################
2196  SUBROUTINE snow3levapn(PPSN3L,PLES3L,PLEL3L,PTSTEP,PSNOWTEMP, &
2197  PSNOWRHO,PSNOWDZ,PSNOWLIQ,PTA, &
2198  PLVTT,PLSTT,PSNOWHEAT,PSOILCOR )
2199 !
2200 !
2201 !! PURPOSE
2202 !! -------
2203 ! Remove mass from uppermost snow layer in response to
2204 ! evaporation (liquid) and sublimation.
2205 !
2206 !
2207 USE modd_csts, ONLY : xrholw, xlmtt, xci, xtt
2208 USE modd_snow_par, ONLY : xrhosmin_es, xsnowdmin
2209 !
2210 IMPLICIT NONE
2211 !
2212 !* 0.1 declarations of arguments
2213 !
2214 REAL, INTENT(IN) :: PTSTEP
2215 !
2216 REAL, DIMENSION(:), INTENT(IN) :: PPSN3L
2217 !
2218 REAL, DIMENSION(:), INTENT(IN) :: PLES3L, PLEL3L ! (W/m2)
2219 !
2220 REAL, DIMENSION(:), INTENT(IN) :: PTA, PLVTT, PLSTT
2221 !
2222 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWHEAT, PSNOWDZ
2223 !
2224 REAL, DIMENSION(:), INTENT(INOUT) :: PSNOWRHO, PSNOWLIQ, &
2225  PSNOWTEMP
2226 !
2227 REAL, DIMENSION(:), INTENT(OUT) :: PSOILCOR
2228 !
2229 !* 0.2 declarations of local variables
2230 !
2231 INTEGER :: INI, INLVLS, JJ, JI
2232 !
2233 REAL, DIMENSION(SIZE(PLES3L)) :: ZSNOWEVAPS, ZSNOWEVAP, ZSNOWEVAPX, &
2234  ZSNOWDZ, ZSNOWHEAT, ZSCAP, ZSNOWTEMP
2235 !
2236 REAL, DIMENSION(SIZE(PLES3L)) :: ZXSE, ZISNOWD
2237 
2238 !* 0.3 declarations of local parameters
2239 
2240 REAL, PARAMETER :: ZSNOWDEMIN = 1.e-4 ! m
2241 REAL, PARAMETER :: ZTDIF = 15. ! K To prevent a possible
2242  ! decoupling of sfc snow T
2243  ! when vanishingly thin, impose
2244  ! this max T-diff based on obs...
2245  ! between Ta-Ts
2246 !
2247 REAL(KIND=JPRB) :: ZHOOK_HANDLE
2248 !
2249 !
2250 !-------------------------------------------------------------------------------
2251 !
2252 ! 0. Initialize:
2253 ! --------------
2254 !
2255 IF (lhook) CALL dr_hook('SNOW3LEVAPN',0,zhook_handle)
2256 !
2257 psoilcor(:) = 0.0
2258 !
2259 zsnowevaps(:) = 0.0
2260 zsnowevap(:) = 0.0
2261 zsnowevapx(:) = 0.0
2262 zsnowdz(:) = 0.0
2263 zscap(:) = 0.0
2264 zsnowheat(:) = 0.0
2265 zsnowtemp(:) = 0.0
2266 zxse(:) = 0.0
2267 !
2268 ini = SIZE(psnowdz(:,:),1)
2269 inlvls = SIZE(psnowdz(:,:),2)
2270 !
2271 !
2272 ! 1. Evaporation of snow liquid water
2273 ! -----------------------------------
2274 ! Evaporation reduces liq water equivalent content
2275 ! of snow pack either by reducing the snow density
2276 ! (evaporation) or the layer thickness (sublimation).
2277 ! Condensation does the reverse.
2278 !
2279 ! Multiply evaporation components by snow fraction
2280 ! to be consistent with fluxes from the snow covered
2281 ! fraction of grid box
2282 !
2283 WHERE(psnowdz(:,1) > 0.0)
2284 !
2285 ! Evaporation:
2286 ! Reduce density and liquid water content:
2287 !
2288  zsnowevap(:) = ppsn3l(:)*plel3l(:)*ptstep/(plvtt(:)*xrholw)
2289  zsnowevapx(:) = min(zsnowevap(:),psnowliq(:))
2290 !
2291 ! This should not change enthalpy (since already accounted for
2292 ! in sfc e budget), so make density change insuring constant enthalpy:
2293 !
2294  psnowliq(:) = psnowliq(:) - zsnowevapx(:)
2295  psnowrho(:) = (psnowheat(:,1)-xlmtt*xrholw*psnowliq(:))/ &
2296  (psnowdz(:,1)*(xci*(psnowtemp(:)-xtt)-xlmtt))
2297 !
2298 ! Budget check: If density drops below minimum, then extract the
2299 ! corresponding water mass from soil (for vanishingly thin snow covers):
2300 !
2301  psoilcor(:) = max(0.0,xrhosmin_es-psnowrho(:))*psnowdz(:,1)/ptstep
2302  psnowrho(:) = max(xrhosmin_es,psnowrho(:))
2303 !
2304 END WHERE
2305 !
2306 ! 2. Update heat capacity:
2307 ! ---------------------------
2308 !
2309 zscap(:) = snow3lscap(psnowrho(:))
2310 !
2311 !
2312 ! 3. Sublimation of snow ice
2313 ! ---------------------------
2314 !
2315 WHERE(psnowdz(:,1) > 0.0)
2316 
2317 ! Budget check: as last traces of liquid in snow evaporates, it is possible
2318 ! evaporation could exceed liquid present in snow. If this is the case,
2319 ! remove residual mass from solid snow in order to maintain high
2320 ! accuracy water balance:
2321 
2322  zsnowevapx(:) = max(0.0, zsnowevap(:) - zsnowevapx(:))
2323  zsnowdz(:) = psnowdz(:,1) - zsnowevapx(:)*xrholw/psnowrho(:)
2324  psnowdz(:,1) = max(0.0, zsnowdz(:))
2325  psoilcor(:) = psoilcor(:) + max(0.0,-zsnowdz(:))*psnowrho(:)/ptstep
2326 
2327 ! Sublimation: Reduce layer thickness and total snow depth
2328 ! if sublimation: add to correction term if potential
2329 ! sublimation exceeds available snow cover.
2330 !
2331  zsnowevaps(:) = ppsn3l(:)*ples3l(:)*ptstep/(plstt(:)*psnowrho(:))
2332  zsnowdz(:) = psnowdz(:,1) - zsnowevaps(:)
2333  psnowdz(:,1) = max(0.0, zsnowdz(:))
2334  psoilcor(:) = psoilcor(:) + max(0.0,-zsnowdz(:))*psnowrho(:)/ptstep
2335 
2336 ! Note that the effect on enthalpy is already accounted for via the surface energy
2337 ! budget, so the change in enthalpy
2338 ! here owing to a mass loss must be offset by possible adjustments to T and Wl
2339 ! to conserve Enthalpy.
2340 
2341  psnowtemp(:) = xtt + ( ((psnowheat(:,1)/max(zsnowdemin,psnowdz(:,1))) &
2342  + xlmtt*psnowrho(:))/zscap(:) )
2343 
2344  psnowliq(:) = max(0.0,psnowtemp(:)-xtt)*zscap(:)* &
2345  psnowdz(:,1)/(xlmtt*xrholw)
2346 
2347  psnowtemp(:) = min(xtt,psnowtemp(:))
2348 
2349 ! For vanishingly thin snow, a decoupling can occur between forcing level T
2350 ! and snow sfc T (as snowpack vanishes owing to sublimation), so a simple limit
2351 ! imposed on this gradient:
2352 
2353  psnowtemp(:) = max(min(xtt,pta(:)-ztdif), psnowtemp(:))
2354 
2355 ! Update surface enthalpy:
2356 
2357  zsnowheat(:) = psnowheat(:,1)
2358  psnowheat(:,1) = psnowdz(:,1)*( zscap(:)*(psnowtemp(:)-xtt) &
2359  - xlmtt*psnowrho(:) ) + xlmtt*xrholw*psnowliq(:)
2360 
2361  zxse(:) = psnowheat(:,1) - zsnowheat(:) ! excess cooling (removed from sfc)
2362 
2363 END WHERE
2364 !
2365 ! If the sfc T-Gradient was limited (and thus the enthalpy in the uppermost layer changed),
2366 ! distribute this energy correction (cooling) in the layers below (thickness-based weighting)
2367 ! to conserve total snowpack enthalpy. Normally this only occurs for this snowpacks,
2368 ! but it can rarely occur otherwise also.
2369 !
2370 zisnowd(:) = 0.
2371 DO jj=2,inlvls
2372  DO ji=1,ini
2373  zisnowd(ji) = zisnowd(ji) + psnowdz(ji,jj) ! m
2374  ENDDO
2375 END DO
2376 zisnowd(:) = zxse(:)/max(zisnowd(:),zsnowdemin) ! J kg-1 m-1
2377 !
2378 DO jj=2,inlvls
2379  DO ji=1,ini
2380  psnowheat(ji,jj) = psnowheat(ji,jj) - psnowdz(ji,jj)*zisnowd(ji)
2381  ENDDO
2382 ENDDO
2383 !
2384 IF (lhook) CALL dr_hook('SNOW3LEVAPN',1,zhook_handle)
2385 !
2386 !-------------------------------------------------------------------------------
2387 !
2388 END SUBROUTINE snow3levapn
2389 !####################################################################
2390 !####################################################################
2391 !####################################################################
2392 SUBROUTINE snow3lgone(PTSTEP,PLEL3L,PLES3L,PSNOWRHO, &
2393  PSNOWHEAT,PRADSINK,PEVAPCOR,PTHRUFAL,PGRNDFLUX, &
2394  PGFLUXSNOW,PGRNDFLUXO,PSNOWDZ,PSNOWLIQ,PSNOWTEMP, &
2395  PLVTT,PLSTT,PRADXS)
2396 !
2397 !
2398 !
2399 !! PURPOSE
2400 !! -------
2401 ! Account for the case when the last trace of snow melts
2402 ! during a time step: ensure mass and heat balance of
2403 ! snow AND underlying surface.
2404 !
2405 !
2406 USE modd_csts, ONLY : xtt, xlstt, xlvtt
2407 USE modd_snow_metamo, ONLY : xsnowdzmin
2408 !
2409 IMPLICIT NONE
2410 !
2411 !* 0.1 declarations of arguments
2412 !
2413 REAL, INTENT(IN) :: PTSTEP
2414 !
2415 REAL, DIMENSION(:), INTENT(IN) :: PLEL3L, PLES3L, PGFLUXSNOW, &
2416  PRADSINK, PGRNDFLUXO, &
2417  PLVTT, PLSTT
2418 !
2419 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWRHO, PSNOWHEAT
2420 !
2421 REAL, DIMENSION(:), INTENT(INOUT) :: PGRNDFLUX, PRADXS
2422 !
2423 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ, PSNOWLIQ, PSNOWTEMP
2424 !
2425 REAL, DIMENSION(:), INTENT(OUT) :: PTHRUFAL ! melt water [kg/(m2 s)]
2426 !
2427 REAL, DIMENSION(:), INTENT(OUT) :: PEVAPCOR ! [kg/(m2 s)]
2428 ! PEVAPCOR = for vanishingy thin snow cover,
2429 ! allow any excess evaporation
2430 ! to be extracted from the soil
2431 ! to maintain an accurate water
2432 ! balance.
2433 !
2434 !* 0.2 declarations of local variables
2435 !
2436 INTEGER :: JJ, JI
2437 !
2438 INTEGER :: INI
2439 INTEGER :: INLVLS
2440 !
2441 REAL, DIMENSION(SIZE(PLES3L)) :: ZSNOWHEATC, ZSNOWGONE_DELTA, ZSNOW
2442 REAL(KIND=JPRB) :: ZHOOK_HANDLE
2443 !
2444 !-------------------------------------------------------------------------------
2445 !
2446 ! 0. Initialize:
2447 ! --------------
2448 !
2449 IF (lhook) CALL dr_hook('SNOW3LGONE',0,zhook_handle)
2450 !
2451 ini = SIZE(psnowdz(:,:),1)
2452 inlvls = SIZE(psnowdz(:,:),2)
2453 !
2454 pevapcor(:) = 0.0
2455 pthrufal(:) = 0.0
2456 !
2457 zsnowheatc(:) = 0.
2458 zsnow(:) = 0.
2459 DO jj=1,inlvls
2460  DO ji=1,ini
2461  zsnowheatc(ji) = zsnowheatc(ji) + psnowheat(ji,jj) ! total heat content (J m-2)
2462  zsnow(ji) = zsnow(ji) + psnowdz(ji,jj) ! total snow depth (m)
2463  ENDDO
2464 ENDDO
2465 zsnowgone_delta(:) = 1.0
2466 !
2467 ! 1. Simple test to see if snow vanishes:
2468 ! ---------------------------------------
2469 ! If so, set thicknesses (and therefore mass and heat) and liquid content
2470 ! to zero, and adjust fluxes of water, evaporation and heat into underlying
2471 ! surface. Note, test with flux computed *before* correction since this represents
2472 ! actual inflow of heat from below (as heat content correction owing to a corrected
2473 ! flux has not yet been done: here we compare to pre-corrected heat content).
2474 !
2475 WHERE(pgfluxsnow(:) + pradsink(:) >= (-zsnowheatc(:)/ptstep) )
2476  pgrndflux(:) = pgfluxsnow(:) + (zsnowheatc(:)/ptstep)
2477  pevapcor(:) = (plel3l(:)/plvtt(:)) + (ples3l(:)/plstt(:))
2478  pradxs(:) = 0.0
2479  zsnowgone_delta(:) = 0.0 ! FLAG...if=0 then snow vanishes, else=1
2480 END WHERE
2481 !
2482 DO jj=1,inlvls
2483  DO ji=1,ini
2484  pthrufal(ji) = pthrufal(ji) + (1.0-zsnowgone_delta(ji))*psnowrho(ji,jj)*psnowdz(ji,jj)/ptstep
2485  ENDDO
2486 END DO
2487 !
2488 ! 2. Final update of snow state
2489 ! -----------------------------
2490 ! (either still present or not):
2491 !
2492 DO jj=1,inlvls
2493  DO ji=1,ini
2494  psnowdz(ji,jj) = psnowdz(ji,jj)*zsnowgone_delta(ji)
2495  psnowliq(ji,jj) = psnowliq(ji,jj)*zsnowgone_delta(ji)
2496  psnowtemp(ji,jj) = (1.0-zsnowgone_delta(ji))*xtt + psnowtemp(ji,jj)*zsnowgone_delta(ji)
2497  ENDDO
2498 ENDDO
2499 IF (lhook) CALL dr_hook('SNOW3LGONE',1,zhook_handle)
2500 !
2501 !
2502 END SUBROUTINE snow3lgone
2503 !####################################################################
2504 !####################################################################
2505 !####################################################################
2506 SUBROUTINE snow3levapgone(PSNOWHEAT,PSNOWDZ,PSNOWRHO,PSNOWTEMP,PSNOWLIQ)
2507 !
2508 !! PURPOSE
2509 !! -------
2510 !
2511 ! If all snow in uppermost layer evaporates/sublimates, re-distribute
2512 ! grid (below assumes very thin snowpacks so layer-thicknesses are
2513 ! constant).
2514 !
2515 !
2516 USE modd_csts, ONLY : xtt, xrholw, xlmtt
2517 USE modd_snow_par, ONLY : xrhosmin_es, xsnowdmin, xrhosmax_es
2518 !
2519 IMPLICIT NONE
2520 !
2521 !* 0.1 declarations of arguments
2522 !
2523 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWRHO ! snow density profile (kg/m3)
2524 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWDZ ! snow layer thickness profile (m)
2525 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWHEAT ! snow heat content/enthalpy (J/m2)
2526 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWTEMP ! snow temperature profile (K)
2527 REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWLIQ ! snow liquid water profile (m)
2528 !
2529 !* 0.2 declarations of local variables
2530 !
2531 INTEGER :: JJ, JI
2532 !
2533 INTEGER :: INI
2534 INTEGER :: INLVLS
2535 !
2536 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: ZSNOWHEAT_1D ! total heat content (J/m2)
2537 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: ZSNOW ! total snow depth (m)
2538 REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: ZMASS ! total mass
2539 !
2540 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)) :: ZSCAP ! Snow layer heat capacity (J/K/m3)
2541 !
2542 REAL(KIND=JPRB) :: ZHOOK_HANDLE
2543 !
2544 !-------------------------------------------------------------------------------
2545 !
2546 ! Initialize:
2547 !
2548 IF (lhook) CALL dr_hook('SNOW3LEVAPGONE',0,zhook_handle)
2549 ini = SIZE(psnowdz,1)
2550 inlvls = SIZE(psnowdz,2)
2551 !
2552 ! First, determine where uppermost snow layer has completely
2553 ! evaporated/sublimated (as it becomes thin):
2554 !
2555 zsnowheat_1d(:) = 0.
2556 zsnow(:) = 0.
2557 zmass(:) = 0.
2558 !
2559 zscap(:,:) = snow3lscap(psnowrho(:,:))
2560 !
2561 DO jj=2,inlvls
2562  DO ji=1,ini
2563  IF(psnowdz(ji,1) == 0.0)THEN
2564  zsnowheat_1d(ji) = zsnowheat_1d(ji) + xlmtt*xrholw*psnowliq(ji,jj) &
2565  + psnowdz(ji,jj)*(zscap(ji,jj)*(psnowtemp(ji,jj)-xtt) &
2566  - xlmtt*psnowrho(ji,jj) )
2567  zsnow(ji) = zsnow(ji) + psnowdz(ji,jj)
2568  zmass(ji) = zmass(ji) + psnowdz(ji,jj)*psnowrho(ji,jj)
2569  ENDIF
2570  ENDDO
2571 ENDDO
2572 !
2573 ! Where uppermost snow layer has vanished, redistribute vertical
2574 ! snow mass and heat profiles (and associated quantities):
2575 !
2576 DO jj=1,inlvls
2577  DO ji=1,ini
2578  IF(zsnow(ji)/= 0.0)THEN
2579  zsnow(ji) = max(0.5*xsnowdmin,zsnow(ji))
2580  psnowdz(ji,jj) = zsnow(ji)/REAL(inlvls)
2581  psnowheat(ji,jj) = zsnowheat_1d(ji)/REAL(inlvls)
2582  psnowrho(ji,jj) = zmass(ji)/zsnow(ji)
2583  ENDIF
2584  ENDDO
2585 ENDDO
2586 !
2587 zscap(:,:) = snow3lscap(psnowrho(:,:))
2588 !
2589 DO jj=1,inlvls
2590  DO ji=1,ini
2591  IF(zsnow(ji)/= 0.0)THEN
2592  psnowtemp(ji,jj) = xtt + ( ((psnowheat(ji,jj)/psnowdz(ji,jj)) &
2593  + xlmtt*psnowrho(ji,jj))/zscap(ji,jj) )
2594  psnowliq(ji,jj) = max(0.0,psnowtemp(ji,jj)-xtt)*zscap(ji,jj)* &
2595  psnowdz(ji,jj)/(xlmtt*xrholw)
2596  psnowtemp(ji,jj) = min(xtt,psnowtemp(ji,jj))
2597  ENDIF
2598  ENDDO
2599 ENDDO
2600 IF (lhook) CALL dr_hook('SNOW3LEVAPGONE',1,zhook_handle)
2601 !
2602 END SUBROUTINE snow3levapgone
2603 !####################################################################
2604 !####################################################################
2605 !####################################################################
2606 SUBROUTINE snow3lebudmeb(PTSTEP, PSNOWDZMIN, &
2607  PTS, PSNOWDZ1, PSNOWDZ2, PSCOND1, PSCOND2, PSCAP, &
2608  PSWNETSNOWS, PLWNETSNOW, &
2609  PHSNOW, PLES3L, PLEL3L, PHPSNOW, &
2610  PCT, PTSTERM1, PTSTERM2, PGFLUXSNOW )
2611 !
2612 !! PURPOSE
2613 !! -------
2614 ! Calculate surface energy budget with surface fluxes imposed.
2615 !
2616 IMPLICIT NONE
2617 !
2618 !* 0.1 declarations of arguments
2619 !
2620 REAL, INTENT(IN) :: PTSTEP, PSNOWDZMIN
2621 !
2622 REAL, DIMENSION(:), INTENT(IN) :: PTS, PSNOWDZ1, PSNOWDZ2, PSCOND1, PSCOND2, PSCAP, &
2623  PHSNOW, PLES3L, PLEL3L, PHPSNOW, &
2624  PSWNETSNOWS, PLWNETSNOW
2625 !
2626 REAL, DIMENSION(:), INTENT(OUT) :: PCT, PTSTERM1, PTSTERM2, PGFLUXSNOW
2627 !
2628 !* 0.2 declarations of local variables
2629 !
2630 REAL, DIMENSION(SIZE(PTS)) :: ZSCONDA, ZA, ZB, ZC, &
2631  ZSNOWDZM1, ZSNOWDZM2
2632 !
2633 REAL(KIND=JPRB) :: ZHOOK_HANDLE
2634 !-------------------------------------------------------------------------------
2635 !
2636 IF (lhook) CALL dr_hook('SNOW3LEBUDMEB',0,zhook_handle)
2637 !
2638 ! Calculate surface energy budget components:
2639 ! ---------------------------------------------------------
2640 ! To prevent numerical difficulties for very thin snow
2641 ! layers, limit the grid "thinness": this is important as
2642 ! layers become vanishing thin:
2643 !
2644 zsnowdzm1(:) = max(psnowdz1(:), psnowdzmin)
2645 zsnowdzm2(:) = max(psnowdz2(:), psnowdzmin)
2646 !
2647 ! Surface thermal inertia:
2648 !
2649 pct(:) = 1.0/(pscap(:)*zsnowdzm1(:))
2650 !
2651 ! Surface fluxes entering the snowpack (radiative and turbulent):
2652 !
2653 pgfluxsnow(:) = pswnetsnows(:) + plwnetsnow(:) - phsnow(:) - ples3l(:) - plel3l(:)
2654 !
2655 ! Thermal conductivity between uppermost and lower snow layers:
2656 !
2657 zsconda(:) = (zsnowdzm1(:)+zsnowdzm2(:))/ &
2658  ((zsnowdzm1(:)/pscond1(:)) + (zsnowdzm2(:)/pscond2(:)))
2659 !
2660 !
2661 ! Energy budget solution terms (with surface flux imposed):
2662 !
2663 zb(:) = 1./ptstep
2664 !
2665 za(:) = zb(:) + pct(:)*(2*zsconda(:)/(zsnowdzm2(:)+zsnowdzm1(:)))
2666 !
2667 zc(:) = pct(:)*( pgfluxsnow(:) + phpsnow(:) )
2668 !
2669 ! Coefficients needed for implicit solution
2670 ! of linearized surface energy budget:
2671 !
2672 ptsterm2(:) = 2*zsconda(:)*pct(:)/(za(:)*(zsnowdzm2(:)+zsnowdzm1(:)))
2673 !
2674 ptsterm1(:) = (pts(:)*zb(:) + zc(:))/za(:)
2675 !
2676 !-------------------------------------------------------------------------------
2677 IF (lhook) CALL dr_hook('SNOW3LEBUDMEB',1,zhook_handle)
2678 !
2679 END SUBROUTINE snow3lebudmeb
2680 !####################################################################
2681 !####################################################################
2682 !####################################################################
2683 !
2684 !
2685 !
2686 END SUBROUTINE snow3l
mode_snow3l::snow3ltransf
Definition: mode_snow3l.F90:144
snow3lrad
subroutine snow3lrad(OMEB, PSNOWDZMIN, PSW_RAD, PSNOWALB, PSPECTRALALBEDO, PSNOWDZ, PSNOWRHO, PALB, PPERMSNOWFRAC, PZENITH, PSWNETSNOW, PSWNETSNOWS, PRADSINK, PRADXS, PSNOWAGE)
Definition: snow3l.F90:1013
modd_surf_par
Definition: modd_surf_par.F90:6
tridiag_ground
subroutine tridiag_ground(PA, PB, PC, PY, PX)
Definition: tridiag_ground.F90:7
modd_snow_metamo
Definition: modd_snow_metamo.F90:8
snow3ldrift
subroutine snow3ldrift(PTSTEP, PFORESTFRAC, PVMOD, PTA, PQA, PPS, PRHOA, PSNOWRHO, PSNOWDZ, PSNOW, OSNOWDRIFT_SUBLIM, PSNDRIFT
Definition: snow3l.F90:782
modd_meb_par::xsw_wght_vis
real, parameter xsw_wght_vis
Definition: modd_meb_par.F90:55
modd_snow_metamo::xsnowdzmin
real, parameter xsnowdzmin
Definition: modd_snow_metamo.F90:45
modd_csts::xcpd
real, save xcpd
Definition: modd_csts.F90:63
yomhook::dr_hook
Definition: yomhook.F90:20
surface_ri
subroutine surface_ri(PTG, PQS, PEXNS, PEXNA, PTA, PQA, PZREF, PUREF, PDIRCOSZW, PVMOD, PRI)
Definition: surface_ri.F90:8
mode_snow3l::snow3lthrm
Definition: mode_snow3l.F90:126
modd_csts::xstefan
real, save xstefan
Definition: modd_csts.F90:59
snow3lrefrz
subroutine snow3lrefrz(PTSTEP, PRR, PSNOWRHO, PSNOWTEMP, PSNOWDZ, PSNOWLIQ, PTHRUFAL)
Definition: snow3l.F90:1815
mode_snow3l::snow3lhold
Definition: mode_snow3l.F90:67
snow3levapgone
subroutine snow3levapgone(PSNOWHEAT, PSNOWDZ, PSNOWRHO, PSNOWTEMP, PSNOWLIQ)
Definition: snow3l.F90:2507
snow3lgone
subroutine snow3lgone(PTSTEP, PLEL3L, PLES3L, PSNOWRHO, PSNOWHEAT, PRADSINK, PEVAPCOR, PTHRUFAL, PGRNDFLUX, PGFLUXSNOW, PGRNDFLUXO, PSNOWDZ, PSNOWLIQ, PSNOWTEMP, PLVTT, PLSTT, PRADXS)
Definition: snow3l.F90:2396
snow3levapn
subroutine snow3levapn(PPSN3L, PLES3L, PLEL3L, PTSTEP, PSNOWTEMP, PSNOWRHO, PSNOWDZ, PSNOWLIQ, PTA, PLVTT, PLSTT, PSNOWHEAT, PSOILCOR)
Definition: snow3l.F90:2199
modd_csts::xlvtt
real, save xlvtt
Definition: modd_csts.F90:70
snow3lebudmeb
subroutine snow3lebudmeb(PTSTEP, PSNOWDZMIN, PTS, PSNOWDZ1, PSNOWDZ2, PSCOND1, PSCOND2, PSCAP, PSWNETSNOWS, PLWNETSNOW, PHSNOW, PLES3L, PLEL3L, PHPSNOW, PCT, PTSTERM1, PTSTERM2, PGFLUXSNOW)
Definition: snow3l.F90:2611
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:1421
modd_csts::xlstt
real, save xlstt
Definition: modd_csts.F90:71
snow3lmelt
subroutine snow3lmelt(PTSTEP, PSCAP, PSNOWTEMP, PSNOWDZ, PSNOWRHO, PSNOWLIQ, PMELTXS)
Definition: snow3l.F90:1672
modd_surf_par::xundef
real, parameter xundef
Definition: modd_surf_par.F90:43
surface_aero_cond
subroutine surface_aero_cond(PRI, PZREF, PUREF, PVMOD, PZ0, PZ0H, PAC, PRA, PCH)
Definition: surface_aero_cond.F90:8
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:2017
parkind1::jprb
integer, parameter jprb
Definition: parkind1.F90:32
mode_snow3l::snow3lgrid
Definition: mode_snow3l.F90:91
modd_csts::xci
real, save xci
Definition: modd_csts.F90:65
mode_thermos
Definition: mode_thermos.F90:6
modd_csts::xday
real, save xday
Definition: modd_csts.F90:45
surface_cd
subroutine surface_cd(PRI, PZREF, PUREF, PZ0EFF, PZ0H, PCD, PCDN)
Definition: surface_cd.F90:8
modd_csts::xcl
real, save xcl
Definition: modd_csts.F90:65
mode_snow3l
Definition: mode_snow3l.F90:6
modd_meb_par
Definition: modd_meb_par.F90:6
modd_type_date_surf::date_time
Definition: modd_type_date_surf.F90:46
mode_thermos::qsati
Definition: mode_thermos.F90:71
yomhook::lhook
logical lhook
Definition: yomhook.F90:15
modd_csts::xrholi
real, save xrholi
Definition: modd_csts.F90:81
snow3l
subroutine snow3l(HSNOWRES, TPTIME, OMEB, HIMPLICIT_WIND,
Definition: snow3l.F90:7
snow3lebud
subroutine snow3lebud(HSNOWRES, HIMPLICIT_WIND,
Definition: snow3l.F90:1172
modd_csts::xrholw
real, save xrholw
Definition: modd_csts.F90:64
modd_csts
Definition: modd_csts.F90:6
mode_snow3l::snow3lfall
Definition: mode_snow3l.F90:140
modd_csts::xtt
real, save xtt
Definition: modd_csts.F90:66
mode_snow3l::snow3lradabs_sfc
Definition: mode_snow3l.F90:122
mode_snow3l::snow3ldopt
Definition: mode_snow3l.F90:130
modd_csts::xlmtt
real, save xlmtt
Definition: modd_csts.F90:72
modd_type_date_surf
Definition: modd_type_date_surf.F90:6
mode_snow3l::snow3lalb
Definition: mode_snow3l.F90:136
mode_snow3l::snow3lcompactn
Definition: mode_snow3l.F90:148
modd_meb_par::xsw_wght_nir
real, parameter xsw_wght_nir
Definition: modd_meb_par.F90:57
Generated on Tue Jan 16 2018 16:23:30 for SURFEX v8.1 by   doxygen 1.8.13