SURFEX v8.1
General documentation of Surfex
snow3L_isba.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_isba(IO, G, PK, PEK, DK, DEK, DMK, OMEB, HIMPLICIT_WIND, &
7  TPTIME, PTSTEP, PVEGTYPE, PTG, PCT, PSOILHCAPZ, &
8  PSOILCONDZ, PPS, PTA, PSW_RAD, PQA, PVMOD, PLW_RAD, PRR, &
9  PSR, PRHOA, PUREF, PEXNS, PEXNA, PDIRCOSZW, PZREF, &
10  PALB, PD_G, PDZG, PPEW_A_COEF, PPEW_B_COEF, PPET_A_COEF, &
11  PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF, PTHRUFAL, &
12  PGRNDFLUX, PFLSN_COR, PGSFCSNOW, PEVAPCOR, PLES3L, PLEL3L,&
13  PEVAP, PSNOWSFCH, PDELHEATN, PDELHEATN_SFC, PRI, PZENITH, &
14  PDELHEATG, PDELHEATG_SFC, PQS )
15 ! ######################################################################################
16 !
17 !!**** *SNOW3L_ISBA*
18 !!
19 !! PURPOSE
20 !! -------
21 !
22 ! 3-Layer snow scheme option (Boone and Etchevers 1999)
23 ! This routine is NOT called as snow depth goes below
24 ! a critical threshold which is vanishingly small.
25 ! This routine acts as an interface between SNOW3L and ISBA.
26 !
27 !!** METHOD
28 !! ------
29 !
30 ! Direct calculation
31 !
32 !! EXTERNAL
33 !! --------
34 !
35 ! None
36 !!
37 !! IMPLICIT ARGUMENTS
38 !! ------------------
39 !!
40 !!
41 !! REFERENCE
42 !! ---------
43 !!
44 !! Boone and Etchevers (1999)
45 !! Belair (1995)
46 !! Noilhan and Planton (1989)
47 !! Noilhan and Mahfouf (1996)
48 !!
49 !! AUTHOR
50 !! ------
51 !! A. Boone * Meteo-France *
52 !!
53 !! MODIFICATIONS
54 !! -------------
55 !! Original 7/99 Boone
56 !! Packing added 4/00 Masson & Boone
57 !! z0h and snow 2/06 LeMoigne
58 !!
59 !! Modified by B. Decharme (03/2009): Consistency with Arpege permanent
60 !! snow/ice treatment
61 !! Modified by A. Boone (04/2010): Implicit coupling with atmosphere permitted.
62 !!
63 !! Modified by B. Decharme (04/2010): check suspicious low temperature for ES and CROCUS
64 !! Modified by B. Decharme (08/2013): Qsat as argument (needed for coupling with atm)
65 !! Modified by A. Boone (10/2014): MEB: pass in fluxes when using MEB
66 !! Modified by B. Decharme (03/2016): No snowdrift under forest
67 !! Modified by M. Lafaysse (08/2015): MEB-Crocus coupling
68 !-------------------------------------------------------------------------------
69 !
70 USE modd_isba_options_n, ONLY : isba_options_t
71 USE modd_sfx_grid_n, ONLY : grid_t
72 USE modd_isba_n, ONLY : isba_pe_t, isba_p_t
73 USE modd_diag_n, ONLY : diag_t
74 USE modd_diag_evap_isba_n, ONLY : diag_evap_isba_t
75 USE modd_diag_misc_isba_n, ONLY : diag_misc_isba_t
76 !
77 USE modd_csts, ONLY : xtt, xpi, xday, xlmtt, xlstt
78 USE modd_snow_par, ONLY : xrhosmax_es, xsnowdmin, xrhosmin_es, xemissn
79 USE modd_surf_par, ONLY : xundef
80 USE modd_type_date_surf, ONLY: date_time
81 !
82 USE modd_data_cover_par, ONLY : nvt_snow, &
83  nvt_tebd, nvt_trbe, nvt_bone, &
84  nvt_trbd, nvt_tebe, nvt_tene, &
85  nvt_bobd, nvt_bond, nvt_shrb
86 !
87 USE modi_snow3l
88 USE modi_snowcro
89 USE modi_snowcro_diag
90 !
91 #ifdef SFX_OL
92 USE modn_io_offline, ONLY : xtstep_output
93 #endif
94 !
95 USE modi_abor1_sfx
96 !
97 USE yomhook ,ONLY : lhook, dr_hook
98 USE parkind1 ,ONLY : jprb
99 !
100 IMPLICIT NONE
101 !
102 !
103 !
104 !* 0.1 declarations of arguments
105 !
106 TYPE(isba_options_t), INTENT(INOUT) :: IO
107 TYPE(grid_t), INTENT(INOUT) :: G
108 TYPE(isba_p_t), INTENT(INOUT) :: PK
109 TYPE(isba_pe_t), INTENT(INOUT) :: PEK
110 TYPE(diag_t), INTENT(INOUT) :: DK
111 TYPE(diag_evap_isba_t), INTENT(INOUT) :: DEK
112 TYPE(diag_misc_isba_t), INTENT(INOUT) :: DMK
113 !
114 LOGICAL, INTENT(IN) :: OMEB ! True = coupled to MEB. This means surface fluxes ae IMPOSED
115 ! ! as an upper boundary condition to the explicit snow schemes.
116 ! ! If = False, then energy
117 ! ! budget and fluxes are computed herein.
118 !
119  CHARACTER(LEN=*), INTENT(IN) :: HIMPLICIT_WIND ! wind implicitation option
120 ! ! 'OLD' = direct
121 ! ! 'NEW' = Taylor serie, order 1
122 !
123 TYPE(date_time), INTENT(IN) :: TPTIME ! current date and time
124 REAL, INTENT(IN) :: PTSTEP
125 ! PTSTEP = time step of the integration
126 !
127 REAL, DIMENSION(:,:), INTENT(IN) :: PVEGTYPE ! fraction of each vegetation
128 !
129 !
130 REAL, DIMENSION(:,:), INTENT(INOUT) :: PTG
131 ! PTG = Soil temperature profile (K)
132 !
133 REAL, DIMENSION(:,:), INTENT(IN) :: PSOILHCAPZ, PD_G, PDZG
134 REAL, DIMENSION(:), INTENT(IN) :: PCT, PSOILCONDZ
135 ! PD_G = Depth to bottom of each soil layer (m)
136 ! PDZG = Soil layer thicknesses (m)
137 ! DMK%XCG = area-averaged soil heat capacity [(K m2)/J]
138 ! PCT = area-averaged surface heat capacity [(K m2)/J]
139 ! PSOILCONDZ= soil thermal conductivity (W m-1 K-1)
140 ! PSOILHCAPZ= soil heat capacity (J m-3 K-1)
141 !
142 REAL, DIMENSION(:), INTENT(IN) :: PPS, PTA, PSW_RAD, PQA, &
143  PVMOD, PLW_RAD, PSR, PRR
144 ! PSW_RAD = incoming solar radiation (W/m2)
145 ! PLW_RAD = atmospheric infrared radiation (W/m2)
146 ! PRR = rain rate [kg/(m2 s)]
147 ! PSR = snow rate (SWE) [kg/(m2 s)]
148 ! PTA = atmospheric temperature at level za (K)
149 ! PVMOD = modulus of the wind parallel to the orography (m/s)
150 ! PPS = surface pressure
151 ! PQA = atmospheric specific humidity
152 ! at level za
153 !
154 REAL, DIMENSION(:), INTENT(IN) :: PZREF, PUREF, PEXNS, PEXNA, PDIRCOSZW, PRHOA, PALB
155 ! PZREF = reference height of the first
156 ! atmospheric level
157 ! PUREF = reference height of the wind
158 ! PRHOA = air density
159 ! PEXNS = Exner function at surface
160 ! PEXNA = Exner function at lowest atmos level
161 ! PDIRCOSZW = Cosinus of the angle between the
162 ! normal to the surface and the vertical
163 ! PALB = soil/vegetation albedo
164 !
165 REAL, DIMENSION(:), INTENT(IN) :: PPEW_A_COEF, PPEW_B_COEF, &
166  PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, &
167  PPEQ_B_COEF
168 ! PPEW_A_COEF = wind coefficient
169 ! PPEW_B_COEF = wind coefficient
170 ! PPET_A_COEF = A-air temperature coefficient
171 ! PPET_B_COEF = B-air temperature coefficient
172 ! PPEQ_A_COEF = A-air specific humidity coefficient
173 ! PPEQ_B_COEF = B-air specific humidity coefficient !
174 !
175 REAL, DIMENSION(:), INTENT(INOUT) :: PLES3L, PLEL3L, PEVAP, PGRNDFLUX, PDELHEATG, PDELHEATG_SFC
176 ! PLEL3L = evaporation heat flux from snow (W/m2)
177 ! PLES3L = sublimation (W/m2)
178 ! PEVAP = total evaporative flux from snow (kg/m2/s)
179 ! PGRNDFLUX = soil/snow interface heat flux (W/m2)
180 ! PDELHEATG = ground heat content change (diagnostic) (W/m2)
181 ! note, modified if ground-snow flux adjusted
182 ! PDELHEATG_SFC = ground heat content change in sfc only (diagnostic) (W/m2)
183 ! note, modified if ground-snow flux adjusted
184 !
185 REAL, DIMENSION(:), INTENT(INOUT) :: PRI
186 ! PRI = Richardson number (-)
187 !
188 REAL, DIMENSION(:), INTENT(OUT) :: PTHRUFAL, PFLSN_COR, PEVAPCOR, PGSFCSNOW
189 ! PTHRUFAL = rate that liquid water leaves snow pack:
190 ! paritioned into soil infiltration/runoff
191 ! by ISBA [kg/(m2 s)]
192 ! PFLSN_COR = soil/snow correction heat flux (W/m2) (not MEB)
193 ! PEVAPCOR = evaporation/sublimation correction term:
194 ! extract any evaporation exceeding the
195 ! actual snow cover (as snow vanishes)
196 ! and apply it as a surface soil water
197 ! sink. [kg/(m2 s)]
198 ! PGSFCSNOW = heat flux between the surface and sub-surface
199 ! snow layers (for energy budget diagnostics) (W/m2)
200 !
201 REAL, DIMENSION(:), INTENT(OUT) :: PSNOWSFCH, PDELHEATN, PDELHEATN_SFC
202 !
203 REAL, DIMENSION(:), INTENT(OUT) :: PQS
204 ! PQS = surface humidity (kg/kg)
205 !
206 ! ajout_EB pour prendre en compte angle zenithal du soleil dans LRAD
207 ! puis plus tard dans LALB
208 REAL, DIMENSION(:), INTENT(IN) :: PZENITH ! solar zenith angle
209 !
210 !* 0.2 declarations of local variables
211 !
212 REAL, PARAMETER :: ZCHECK_TEMP = 50.0
213 ! Limit to check suspicious low temperature (K)
214 !
215 INTEGER :: JWRK, JJ ! Loop control
216 !
217 INTEGER :: INLVLS ! maximum number of snow layers
218 INTEGER :: INLVLG ! number of ground layers
219 !
220 REAL, DIMENSION(SIZE(PTA)) :: ZRRSNOW, ZSOILCOND, ZSNOW, ZSNOWFALL, &
221  ZSNOWABLAT_DELTA, ZSNOWSWE_1D, ZSNOWD, &
222  ZSNOWH, ZSNOWH1, ZGRNDFLUXN, ZPSN, &
223  ZSOILCOR, ZSNOWSWE_OUT, ZTHRUFAL, &
224  ZSNOW_MASS_BUDGET, ZWGHT, ZWORK, ZC2
225 ! ZSOILCOND = soil thermal conductivity [W/(m K)]
226 ! ZRRSNOW = rain rate over snow [kg/(m2 s)]
227 ! ZSNOW = snow depth (m)
228 ! ZSNOWFALL = minimum equivalent snow depth
229 ! for snow falling during the
230 ! current time step (m)
231 ! ZSNOWABLAT_DELTA = FLAG =1 if snow ablates completely
232 ! during current time step, else=0
233 ! ZSNOWSWE_1D = TOTAL snowpack SWE (kg m-2)
234 ! ZSNOWD = snow depth
235 ! ZSNOWH = snow total heat content (J m-2)
236 ! ZSNOWH1 = snow surface layer heat content (J m-2)
237 ! ZGRNDFLUXN = corrected snow-ground flux (if snow fully ablated during timestep)
238 ! ZPSN = snow fraction working array
239 ! ZSOILCOR = for vanishingy thin snow cover,
240 ! allow any excess evaporation
241 ! to be extracted from the soil
242 ! to maintain an accurate water
243 ! balance [kg/(m2 s)]
244 ! ZSNOW_MASS_BUDGET = snow water equivalent budget (kg/m2/s)
245 ! ZWGHT = MEB surface layer weight for distributing energy
246 ! between litter and ground layers for the case
247 ! of total ablation during a timestep (-).
248 ! ZWORK = local working variable (*)
249 ! ZC2 = sub-surface heat capacity [(K m2)/J]
250 !
251 !* 0.3 declarations of packed variables
252 !
253 INTEGER :: ISIZE_SNOW ! number of points where computations are done
254 INTEGER, DIMENSION(SIZE(PTA)) :: NMASK ! indices correspondance between arrays
255 !
256 LOGICAL, DIMENSION(SIZE(PTA)) :: LREMOVE_SNOW
257 !
258 REAL, DIMENSION(SIZE(PTA)) :: ZSWNET_N, ZSWNET_NS, ZLWNET_N
259 !
260 LOGICAL :: GCOMPUTECRODIAG ! flag to compute Crocus-MEPRA diagnostics
261 !
262 REAL(KIND=JPRB) :: ZHOOK_HANDLE
263 !
264 ! - - ---------------------------------------------------
265 !
266 IF (lhook) CALL dr_hook('SNOW3L_ISBA',0,zhook_handle)
267 !
268 !* 0. Initialize variables:
269 ! ---------------------
270 !
271 IF (SIZE(dmk%XSNOWDEND)>0) THEN
272  dmk%XSNOWDEND (:,:) = xundef
273  dmk%XSNOWSPHER(:,:) = xundef
274  dmk%XSNOWSIZE (:,:) = xundef
275  dmk%XSNOWSSA (:,:) = xundef
276  dmk%XSNOWRAM (:,:) = xundef
277  dmk%XSNOWSHEAR(:,:) = xundef
278  dmk%XSNOWTYPEMEPRA(:,:) = xundef
279 ENDIF
280 !
281 dek%XSNDRIFT(:) = 0.0
282 dmk%XSNOWHMASS(:) = 0.0
283 dmk%XSRSFC(:) = psr(:) ! these are snow and rain rates passed to ISBA,
284 dmk%XRRSFC(:) = prr(:) ! so initialize here if SNOW3L not used:
285 !
286 pflsn_cor(:) = 0.0
287 pthrufal(:) = 0.0
288 pevapcor(:) = 0.0
289 pqs(:) = xundef
290 !
291 zsnow(:) = 0.0
292 zsnowd(:) = 0.0
293 zgrndfluxn(:) = 0.0
294 zsnowh(:) = 0.0
295 zsnowh1(:) = 0.0
296 zsnowswe_1d(:) = 0.0
297 zsnowswe_out(:)= 0.0
298 zsoilcond(:) = 0.0
299 zrrsnow(:) = 0.0
300 zsnowfall(:) = 0.0
301 zsnowablat_delta(:) = 0.0
302 !
303 zwght(:) = 0.0
304 zwork(:) = 0.0
305 zc2(:) = pct(:)
306 !
307 dmk%XSNOWLIQ(:,:) = 0.0
308 dmk%XSNOWDZ(:,:) = 0.0
309 !
310 inlvls = SIZE(pek%TSNOW%WSNOW(:,:),2)
311 inlvlg = min(SIZE(pd_g(:,:),2),SIZE(ptg(:,:),2))
312 !
313 !
314 IF(.NOT.omeb)THEN
315  !
316  ! If MEB activated, these values are input, else initialize here:
317  !
318  pgrndflux(:) = 0.0
319  ples3l(:) = 0.0
320  plel3l(:) = 0.0
321  pevap(:) = 0.0
322  pri(:) = xundef
323  pek%TSNOW%EMIS(:) = xemissn
324  dmk%XRNSNOW(:) = 0.0
325  dmk%XHSNOW(:) = 0.0
326  dmk%XGFLUXSNOW(:) = 0.0
327  dmk%XHPSNOW(:) = 0.0
328  dmk%XUSTARSNOW(:) = 0.0
329  dmk%XCDSNOW(:) = 0.0
330  dmk%XCHSNOW(:) = 0.0
331  !
332 ENDIF
333 !
334 ! Use ISBA-SNOW3L or NOT: NOTE that if explicit soil diffusion method in use,
335 ! then *must* use explicit snow model:
336 !
337 IF (pek%TSNOW%SCHEME=='3-L' .OR. io%CISBA == 'DIF' .OR. pek%TSNOW%SCHEME == 'CRO') THEN
338 !
339  IF(.NOT.omeb)THEN
340  !
341  ! If MEB activated, these values are input, else initialize here:
342  !
343  zswnet_n(:) = 0.0
344  zswnet_ns(:) = 0.0
345  zlwnet_n(:) = 0.0
346  ELSE
347  zswnet_n(:) = dek%XSWNET_N(:)
348  zswnet_ns(:) = dek%XSWNET_NS(:)
349  zlwnet_n(:) = dek%XLWNET_N(:)
350  END IF
351 
352 ! - Snow and rain falling onto the 3-L grid space:
353 !
354  dmk%XSRSFC(:) = 0.0
355 !
356  DO jj=1,SIZE(psr)
357  zrrsnow(jj) = pek%XPSN(jj)*prr(jj)
358  dmk%XRRSFC(jj) = prr(jj) - zrrsnow(jj)
359  zsnowfall(jj) = psr(jj)*ptstep/xrhosmax_es ! maximum possible snowfall depth (m)
360  ENDDO
361 !
362 ! Calculate preliminary snow depth (m)
363 
364  zsnow(:) =0.
365  zsnowh(:) =0.
366  zsnowswe_1d(:)=0.
367  zsnowh1(:) = pek%TSNOW%HEAT(:,1)*pek%TSNOW%WSNOW(:,1)/pek%TSNOW%RHO(:,1) ! sfc layer only
368 !
369  DO jwrk=1,SIZE(pek%TSNOW%WSNOW(:,:),2)
370  DO jj=1,SIZE(pek%TSNOW%WSNOW(:,:),1)
371  zsnowswe_1d(jj) = zsnowswe_1d(jj) + pek%TSNOW%WSNOW(jj,jwrk)
372  zsnow(jj) = zsnow(jj) + pek%TSNOW%WSNOW(jj,jwrk)/pek%TSNOW%RHO(jj,jwrk)
373  zsnowh(jj) = zsnowh(jj) + pek%TSNOW%HEAT (jj,jwrk)*pek%TSNOW%WSNOW(jj,jwrk)/pek%TSNOW%RHO(jj,jwrk)
374  END DO
375  ENDDO
376 !
377  IF(io%CISBA == 'DIF')THEN
378  zsoilcond(:) = psoilcondz(:)
379  ELSE
380 !
381 ! - Soil thermal conductivity
382 ! is implicit in Force-Restore soil method, so it
383 ! must be backed-out of surface thermal coefficients
384 ! (Etchevers and Martin 1997):
385 !
386  zsoilcond(:) = 4.*xpi/( dmk%XCG(:)*dmk%XCG(:)*xday/(pd_g(:,1)*pct(:)) )
387 !
388  ENDIF
389 !
390 ! ===============================================================
391 ! === Packing: Only call snow model when there is snow on the surface
392 ! exceeding a minimum threshold OR if the equivalent
393 ! snow depth falling during the current time step exceeds
394 ! this limit.
395 !
396 ! counts the number of points where the computations will be made
397 !
398 !
399  isize_snow = 0
400  nmask(:) = 0
401 !
402  DO jj=1,SIZE(zsnow)
403  IF (zsnow(jj) >= xsnowdmin .OR. zsnowfall(jj) >= xsnowdmin) THEN
404  isize_snow = isize_snow + 1
405  nmask(isize_snow) = jj
406  ENDIF
407  ENDDO
408 !
409  IF (isize_snow>0) CALL call_model(isize_snow,inlvls,inlvlg,nmask)
410 !
411 ! ===============================================================
412 !
413 ! Remove trace amounts of snow and reinitialize snow prognostic variables
414 ! if snow cover is ablated.
415 ! If MEB used, soil T already computed, therefore correct heating/cooling
416 ! effect of updated snow-soil flux
417 !
418  zsnowd(:) = 0.
419  zsnowswe_out(:) = 0.
420  DO jwrk=1,SIZE(pek%TSNOW%WSNOW(:,:),2)
421  DO jj=1,SIZE(pek%TSNOW%WSNOW(:,:),1)
422  zsnowd(jj) = zsnowd(jj) + pek%TSNOW%WSNOW(jj,jwrk)/pek%TSNOW%RHO(jj,jwrk)
423  zsnowswe_out(jj) = zsnowswe_out(jj) + pek%TSNOW%WSNOW(jj,jwrk)
424  ENDDO
425  END DO
426 !
427  lremove_snow(:)=(zsnowd(:)<xsnowdmin*1.1)
428 !
429 !
430  IF(omeb)THEN
431  zpsn(:) = 1.0
432  IF(io%CISBA == 'DIF')THEN
433  zwght(:) = psoilhcapz(:,2)*pdzg(:,2)/(psoilhcapz(:,1)*pdzg(:,1) + psoilhcapz(:,2)*pdzg(:,2))
434  zc2(:) = 1/(psoilhcapz(:,2)*pdzg(:,2))
435  ELSE
436  zwght(:) = (pd_g(:,2)-pd_g(:,1))/pd_g(:,2)
437  ENDIF
438  ELSE
439 ! To Conserve mass in ISBA without MEB,
440 ! EVAP must be weignted by the snow fraction
441 ! in the calulation of THRUFAL
442  zpsn(:) = pek%XPSN(:)
443  ENDIF
444 !
445  zsnowablat_delta(:) = 0.0
446  zthrufal(:) = pthrufal(:)
447 !
448  WHERE(lremove_snow(:))
449  !
450  zsnowswe_out(:) = 0.0
451  ples3l(:) = min(ples3l(:), xlstt*(zsnowswe_1d(:)/ptstep + psr(:)))
452  plel3l(:) = 0.0
453  pevap(:) = ples3l(:)/pk%XLSTT(:)
454  pthrufal(:) = max(0.0, zsnowswe_1d(:)/ptstep + psr(:) - pevap(:)*zpsn(:) + zrrsnow(:)) ! kg m-2 s-1
455  zthrufal(:) = max(0.0, zsnowswe_1d(:)/ptstep + psr(:) - pevap(:) + zrrsnow(:)) ! kg m-2 s-1
456  !
457  dmk%XSRSFC(:) = 0.0
458  dmk%XRRSFC(:) = dmk%XRRSFC(:)
459  !
460  zsnowablat_delta(:) = 1.0
461  !
462  pek%TSNOW%ALB(:) = xundef
463  !
464  pevapcor(:) = 0.0
465  zsoilcor(:) = 0.0
466  !
467  dmk%XGFLUXSNOW(:) = dmk%XRNSNOW(:) - dmk%XHSNOW(:) - ples3l(:) - plel3l(:)
468  dmk%XSNOWHMASS(:) = -psr(:)*(xlmtt*ptstep)
469  !
470  pgsfcsnow(:) = 0.0
471  pdelheatn(:) = -zsnowh(:) /ptstep
472  pdelheatn_sfc(:) = -zsnowh1(:)/ptstep
473  psnowsfch(:) = pdelheatn_sfc(:) - (zswnet_ns(:) + zlwnet_n(:) &
474  - dmk%XHSNOW(:) - ples3l(:) - plel3l(:)) + pgsfcsnow(:) &
475  - dmk%XSNOWHMASS(:)/ptstep
476  zgrndfluxn(:) = (zsnowh(:)+dmk%XSNOWHMASS(:))/ptstep + dmk%XGFLUXSNOW(:)
477  zwork(:) = ptstep * zpsn(:) * (zgrndfluxn(:) - pgrndflux(:) - pflsn_cor(:))
478  ptg(:,1) = ptg(:,1) + zwork(:)*(1.-zwght(:))*pct(:)
479  ptg(:,2) = ptg(:,2) + zwork(:)* zwght(:) *zc2(:)
480  zwork(:) = zwork(:) / ptstep
481  pdelheatg(:) = pdelheatg(:) + zwork(:)
482  pdelheatg_sfc(:) = pdelheatg_sfc(:) + zwork(:)
483  pgrndflux(:) = zgrndfluxn(:)
484  pflsn_cor(:) = 0.0
485  !
486  END WHERE
487 !
488 !
489  DO jwrk=1,inlvls
490  DO jj=1,SIZE(pek%TSNOW%WSNOW(:,:),1)
491  pek%TSNOW%WSNOW(jj,jwrk) = (1.0-zsnowablat_delta(jj))*pek%TSNOW%WSNOW(jj,jwrk)
492  pek%TSNOW%HEAT (jj,jwrk) = (1.0-zsnowablat_delta(jj))*pek%TSNOW%HEAT (jj,jwrk)
493  pek%TSNOW%RHO (jj,jwrk) = (1.0-zsnowablat_delta(jj))*pek%TSNOW%RHO (jj,jwrk) + &
494  zsnowablat_delta(jj)*xrhosmin_es
495  pek%TSNOW%AGE(jj,jwrk) = (1.0-zsnowablat_delta(jj))*pek%TSNOW%AGE (jj,jwrk)
496  dmk%XSNOWTEMP(jj,jwrk) = (1.0-zsnowablat_delta(jj))*dmk%XSNOWTEMP(jj,jwrk) + &
497  zsnowablat_delta(jj)*xtt
498  dmk%XSNOWLIQ (jj,jwrk) = (1.0-zsnowablat_delta(jj))*dmk%XSNOWLIQ(jj,jwrk)
499  dmk%XSNOWDZ (jj,jwrk) = (1.0-zsnowablat_delta(jj))*dmk%XSNOWDZ (jj,jwrk)
500  ENDDO
501  ENDDO
502 !
503  IF (pek%TSNOW%SCHEME=='CRO') THEN
504  DO jwrk=1,inlvls
505  DO jj=1,SIZE(pek%TSNOW%GRAN1(:,:),1)
506  pek%TSNOW%GRAN1(jj,jwrk) = (1.0-zsnowablat_delta(jj))*pek%TSNOW%GRAN1(jj,jwrk)
507  pek%TSNOW%GRAN2(jj,jwrk) = (1.0-zsnowablat_delta(jj))*pek%TSNOW%GRAN2(jj,jwrk)
508  pek%TSNOW%HIST (jj,jwrk) = (1.0-zsnowablat_delta(jj))*pek%TSNOW%HIST (jj,jwrk)
509  ENDDO
510  ENDDO
511  ENDIF
512 !
513 ! ===============================================================
514 !
515 ! Compute snow mass budget
516 !
517  zsnow_mass_budget(:) = (zsnowswe_1d(:)-zsnowswe_out(:))/ptstep + psr(:)+zrrsnow(:) &
518  - pevap(:)-zthrufal(:) &
519  + pevapcor(:)+zsoilcor(:)
520 !
521 !
522 ! ===============================================================
523 !
524 ! To Conserve mass in ISBA, the latent heat flux part of
525 ! the EVAPCOR term must be weignted by the snow fraction
526 !
527  pevapcor(:) = pevapcor(:)*zpsn(:) + zsoilcor(:)
528 !
529 ! ===============================================================
530 !
531 ! check suspicious low temperature
532 !
533  DO jwrk=1,inlvls
534  !
535  DO jj=1,SIZE(pek%TSNOW%WSNOW,1)
536  !
537  IF (pek%TSNOW%WSNOW(jj,jwrk)>0.0) THEN
538  !
539  IF (dmk%XSNOWTEMP(jj,jwrk)<zcheck_temp) THEN
540  WRITE(*,*) 'Suspicious low temperature :',dmk%XSNOWTEMP(jj,jwrk)
541  WRITE(*,*) 'At point and location :',jj,'LAT=',g%XLAT(jj),'LON=',g%XLON(jj)
542  WRITE(*,*) 'At snow level / total layer:',jwrk,'/',inlvls
543  WRITE(*,*) 'SNOW MASS BUDGET (kg/m2/s) :',zsnow_mass_budget(jj)
544  WRITE(*,*) 'SWE BY LAYER (kg/m2) :',pek%TSNOW%WSNOW (jj,1:inlvls)
545  WRITE(*,*) 'DEKTH BY LAYER (m) :',dmk%XSNOWDZ (jj,1:inlvls)
546  WRITE(*,*) 'DENSITY BY LAYER (kg/m3) :',pek%TSNOW%RHO(jj,1:inlvls)
547  WRITE(*,*) 'TEMPERATURE BY LAYER (K) :',dmk%XSNOWTEMP(jj,1:inlvls)
548  CALL abor1_sfx('SNOW3L_ISBA: Suspicious low temperature')
549  ENDIF
550  !
551  ELSE
552  !
553  !Prognostic variables forced to XUNDEF for correct outputs
554  dmk%XSNOWDZ(jj,jwrk)=xundef
555  ! Careful : to compute average surface temperature in ISBA_SNOW_AGR
556  ! PSNOWTEMP(JJ,1) is required when PPSN(JJ)>0 even if PSNOWSWE(JJ,1)==0
557  ! (vanishing snowpack)
558  IF (.NOT.(pek%XPSN(jj)>0.0.AND.jwrk==1)) dmk%XSNOWTEMP(jj,jwrk) = xundef
559  dmk%XSNOWLIQ (jj,jwrk) = xundef
560  pek%TSNOW%HEAT(jj,jwrk) = xundef
561  pek%TSNOW%RHO (jj,jwrk) = xundef
562  pek%TSNOW%AGE (jj,jwrk) = xundef
563  IF (pek%TSNOW%SCHEME=='CRO') THEN
564  pek%TSNOW%GRAN1(jj,jwrk) = xundef
565  pek%TSNOW%GRAN2(jj,jwrk) = xundef
566  pek%TSNOW%HIST (jj,jwrk) = xundef
567  END IF
568  ENDIF
569  ENDDO
570  ENDDO
571 !
572  IF(omeb)THEN
573  dek%XSWNET_N(:) = zswnet_n(:)
574  dek%XSWNET_NS(:) = zswnet_ns(:)
575  dek%XLWNET_N(:) = zlwnet_n(:)
576  END IF
577 
578 ! ===============================================================
579 !
580 ENDIF
581 !
582 IF (lhook) CALL dr_hook('SNOW3L_ISBA',1,zhook_handle)
583 !
584  CONTAINS
585 !
586 !================================================================
587 SUBROUTINE call_model(KSIZE1,KSIZE2,KSIZE3,KMASK)
588 !
589 IMPLICIT NONE
590 !
591 INTEGER, INTENT(IN) :: KSIZE1
592 INTEGER, INTENT(IN) :: KSIZE2
593 INTEGER, INTENT(IN) :: KSIZE3
594 INTEGER, DIMENSION(:), INTENT(IN) :: KMASK
595 !
596 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWSWE
597 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWDZ
598 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWRHO
599 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWHEAT
600 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWTEMP
601 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWLIQ
602 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWGRAN1
603 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWGRAN2
604 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWHIST
605 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWAGE
606 REAL, DIMENSION(KSIZE1) :: ZP_SNOWALB
607 REAL, DIMENSION(KSIZE1) :: ZP_SWNETSNOW
608 REAL, DIMENSION(KSIZE1) :: ZP_SWNETSNOWS
609 REAL, DIMENSION(KSIZE1) :: ZP_LWNETSNOW
610 REAL, DIMENSION(KSIZE1) :: ZP_PS
611 REAL, DIMENSION(KSIZE1) :: ZP_SRSNOW
612 REAL, DIMENSION(KSIZE1) :: ZP_RRSNOW
613 REAL, DIMENSION(KSIZE1) :: ZP_PSN3L
614 REAL, DIMENSION(KSIZE1) :: ZP_TA
615 REAL, DIMENSION(KSIZE1) :: ZP_CT
616 REAL, DIMENSION(KSIZE1,KSIZE3) :: ZP_TG
617 REAL, DIMENSION(KSIZE1,KSIZE3) :: ZP_D_G
618 REAL, DIMENSION(KSIZE1,KSIZE3) :: ZP_DZG
619 REAL, DIMENSION(KSIZE1,KSIZE3) :: ZP_SOILHCAPZ
620 REAL, DIMENSION(KSIZE1) :: ZP_SOILD
621 REAL, DIMENSION(KSIZE1) :: ZP_DELHEATG
622 REAL, DIMENSION(KSIZE1) :: ZP_DELHEATG_SFC
623 REAL, DIMENSION(KSIZE1) :: ZP_SW_RAD
624 REAL, DIMENSION(KSIZE1) :: ZP_QA
625 REAL, DIMENSION(KSIZE1) :: ZP_LVTT
626 REAL, DIMENSION(KSIZE1) :: ZP_LSTT
627 REAL, DIMENSION(KSIZE1) :: ZP_VMOD
628 REAL, DIMENSION(KSIZE1) :: ZP_LW_RAD
629 REAL, DIMENSION(KSIZE1) :: ZP_RHOA
630 REAL, DIMENSION(KSIZE1) :: ZP_UREF
631 REAL, DIMENSION(KSIZE1) :: ZP_EXNS
632 REAL, DIMENSION(KSIZE1) :: ZP_EXNA
633 REAL, DIMENSION(KSIZE1) :: ZP_DIRCOSZW
634 REAL, DIMENSION(KSIZE1) :: ZP_ZREF
635 REAL, DIMENSION(KSIZE1) :: ZP_Z0NAT
636 REAL, DIMENSION(KSIZE1) :: ZP_Z0HNAT
637 REAL, DIMENSION(KSIZE1) :: ZP_Z0EFF
638 REAL, DIMENSION(KSIZE1) :: ZP_ALB
639 REAL, DIMENSION(KSIZE1) :: ZP_SOILCOND
640 REAL, DIMENSION(KSIZE1) :: ZP_THRUFAL
641 REAL, DIMENSION(KSIZE1) :: ZP_GRNDFLUX
642 REAL, DIMENSION(KSIZE1) :: ZP_FLSN_COR
643 REAL, DIMENSION(KSIZE1) :: ZP_GSFCSNOW
644 REAL, DIMENSION(KSIZE1) :: ZP_EVAPCOR
645 REAL, DIMENSION(KSIZE1) :: ZP_SOILCOR
646 REAL, DIMENSION(KSIZE1) :: ZP_GFLXCOR
647 REAL, DIMENSION(KSIZE1) :: ZP_RNSNOW
648 REAL, DIMENSION(KSIZE1) :: ZP_HSNOW
649 REAL, DIMENSION(KSIZE1) :: ZP_GFLUXSNOW
650 REAL, DIMENSION(KSIZE1) :: ZP_DELHEATN
651 REAL, DIMENSION(KSIZE1) :: ZP_DELHEATN_SFC
652 REAL, DIMENSION(KSIZE1) :: ZP_SNOWSFCH
653 REAL, DIMENSION(KSIZE1) :: ZP_HPSNOW
654 REAL, DIMENSION(KSIZE1) :: ZP_LES3L
655 REAL, DIMENSION(KSIZE1) :: ZP_LEL3L
656 REAL, DIMENSION(KSIZE1) :: ZP_EVAP
657 REAL, DIMENSION(KSIZE1) :: ZP_SNDRIFT
658 REAL, DIMENSION(KSIZE1) :: ZP_RI
659 REAL, DIMENSION(KSIZE1) :: ZP_QS
660 REAL, DIMENSION(KSIZE1) :: ZP_EMISNOW
661 REAL, DIMENSION(KSIZE1) :: ZP_CDSNOW
662 REAL, DIMENSION(KSIZE1) :: ZP_USTARSNOW
663 REAL, DIMENSION(KSIZE1) :: ZP_CHSNOW
664 REAL, DIMENSION(KSIZE1) :: ZP_SNOWHMASS
665 REAL, DIMENSION(KSIZE1) :: ZP_VEGTYPE
666 REAL, DIMENSION(KSIZE1) :: ZP_FOREST
667 REAL, DIMENSION(KSIZE1) :: ZP_PEW_A_COEF
668 REAL, DIMENSION(KSIZE1) :: ZP_PEW_B_COEF
669 REAL, DIMENSION(KSIZE1) :: ZP_PET_A_COEF
670 REAL, DIMENSION(KSIZE1) :: ZP_PET_B_COEF
671 REAL, DIMENSION(KSIZE1) :: ZP_PEQ_A_COEF
672 REAL, DIMENSION(KSIZE1) :: ZP_PEQ_B_COEF
673 REAL, DIMENSION(KSIZE1) :: ZP_ZENITH
674 REAL, DIMENSION(KSIZE1) :: ZP_LAT,ZP_LON
675 REAL, DIMENSION(KSIZE1) :: ZP_PSN_INV
676 REAL, DIMENSION(KSIZE1) :: ZP_PSN
677 REAL, DIMENSION(KSIZE1) :: ZP_PSN_GFLXCOR
678 REAL, DIMENSION(KSIZE1) :: ZP_WORK
679 !
680 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWDEND
681 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWSPHER
682 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWSIZE
683 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWSSA
684 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWTYPEMEPRA
685 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWRAM
686 REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWSHEAR
687 REAL, DIMENSION(KSIZE1) :: ZP_SNDPT_1DY
688 REAL, DIMENSION(KSIZE1) :: ZP_SNDPT_3DY
689 REAL, DIMENSION(KSIZE1) :: ZP_SNDPT_5DY
690 REAL, DIMENSION(KSIZE1) :: ZP_SNDPT_7DY
691 REAL, DIMENSION(KSIZE1) :: ZP_SNSWE_1DY
692 REAL, DIMENSION(KSIZE1) :: ZP_SNSWE_3DY
693 REAL, DIMENSION(KSIZE1) :: ZP_SNSWE_5DY
694 REAL, DIMENSION(KSIZE1) :: ZP_SNSWE_7DY
695 REAL, DIMENSION(KSIZE1) :: ZP_SNRAM_SONDE
696 REAL, DIMENSION(KSIZE1) :: ZP_SN_WETTHCKN
697 REAL, DIMENSION(KSIZE1) :: ZP_SN_REFRZNTHCKN
698 !
699 REAL, PARAMETER :: ZDEPTHABS = 0.60 ! m
700 !
701 INTEGER :: JWRK, JJ, JI
702 REAL(KIND=JPRB) :: ZHOOK_HANDLE
703 !
704 IF (lhook) CALL dr_hook('SNOW3L_ISBA:CALL_MODEL',0,zhook_handle)
705 !
706 ! Initialize:
707 !
708 zp_psn_gflxcor(:) = 0.
709 zp_work(:) = 0.
710 zp_soild(:) = 0.
711 !
712 ! pack the variables
713 !
714 DO jwrk=1,ksize2
715  DO jj=1,ksize1
716  ji = kmask(jj)
717  zp_snowswe(jj,jwrk) = pek%TSNOW%WSNOW (ji,jwrk)
718  zp_snowrho(jj,jwrk) = pek%TSNOW%RHO (ji,jwrk)
719  zp_snowheat(jj,jwrk) = pek%TSNOW%HEAT(ji,jwrk)
720  zp_snowage(jj,jwrk) = pek%TSNOW%AGE (ji,jwrk)
721  zp_snowtemp(jj,jwrk) = dmk%XSNOWTEMP(ji,jwrk)
722  zp_snowliq(jj,jwrk) = dmk%XSNOWLIQ (ji,jwrk)
723  zp_snowdz(jj,jwrk) = dmk%XSNOWDZ (ji,jwrk)
724  ENDDO
725 ENDDO
726 !
727 IF (pek%TSNOW%SCHEME=='CRO') THEN
728  DO jwrk=1,ksize2
729  DO jj=1,ksize1
730  ji = kmask(jj)
731  zp_snowgran1(jj,jwrk) = pek%TSNOW%GRAN1 (ji,jwrk)
732  zp_snowgran2(jj,jwrk) = pek%TSNOW%GRAN2 (ji,jwrk)
733  zp_snowhist(jj,jwrk) = pek%TSNOW%HIST (ji,jwrk)
734  ENDDO
735  ENDDO
736 ELSE
737  DO jwrk=1,ksize2
738  DO jj=1,ksize1
739  zp_snowgran1(jj,jwrk) = xundef
740  zp_snowgran2(jj,jwrk) = xundef
741  zp_snowhist(jj,jwrk) = xundef
742  ENDDO
743  ENDDO
744 ENDIF
745 !
746 DO jwrk=1,ksize3
747  DO jj=1,ksize1
748  ji = kmask(jj)
749  zp_tg(jj,jwrk) = ptg(ji,jwrk)
750  zp_d_g(jj,jwrk) = pd_g(ji,jwrk)
751  zp_soilhcapz(jj,jwrk) = psoilhcapz(ji,jwrk)
752  ENDDO
753 ENDDO
754 !
755 IF (omeb) THEN
756  DO jwrk=1,ksize3
757  DO jj=1,ksize1
758  ji = kmask(jj)
759  zp_dzg(jj,jwrk) = pdzg(ji,jwrk)
760  ENDDO
761  ENDDO
762 ENDIF
763 !
764 DO jj=1,ksize1
765  ji = kmask(jj)
766  zp_lvtt(jj) = pk%XLVTT (ji)
767  zp_lstt(jj) = pk%XLSTT (ji)
768  zp_emisnow(jj) = pek%TSNOW%EMIS(ji)
769  zp_snowalb(jj) = pek%TSNOW%ALB (ji)
770  zp_psn3l(jj) = pek%XPSN (ji)
771  zp_z0nat(jj) = dk%XZ0 (ji)
772  zp_z0hnat(jj) = dk%XZ0H (ji)
773  zp_z0eff(jj) = dk%XZ0EFF(ji)
774  zp_rnsnow(jj) = dmk%XRNSNOW (ji)
775  zp_hsnow(jj) = dmk%XHSNOW (ji)
776  zp_hpsnow(jj) = dmk%XHPSNOW (ji)
777 
778  zp_ps(jj) = pps(ji)
779  zp_srsnow(jj) = psr(ji)
780  zp_ct(jj) = pct(ji)
781  zp_ta(jj) = pta(ji)
782  zp_delheatg(jj) = pdelheatg(ji)
783  zp_delheatg_sfc(jj) = pdelheatg_sfc(ji)
784  zp_sw_rad(jj) = psw_rad(ji)
785  zp_qa(jj) = pqa(ji)
786  zp_vmod(jj) = pvmod(ji)
787  zp_lw_rad(jj) = plw_rad(ji)
788  zp_rhoa(jj) = prhoa(ji)
789  zp_uref(jj) = puref(ji)
790  zp_exns(jj) = pexns(ji)
791  zp_exna(jj) = pexna(ji)
792  zp_dircoszw(jj) = pdircoszw(ji)
793  zp_zref(jj) = pzref(ji)
794  zp_alb(jj) = palb(ji)
795 
796  zp_rrsnow(jj) = zrrsnow(ji)
797  zp_soilcond(jj) = zsoilcond(ji)
798 
799  !
800  zp_pew_a_coef(jj) = ppew_a_coef(ji)
801  zp_pew_b_coef(jj) = ppew_b_coef(ji)
802  zp_pet_a_coef(jj) = ppet_a_coef(ji)
803  zp_peq_a_coef(jj) = ppeq_a_coef(ji)
804  zp_pet_b_coef(jj) = ppet_b_coef(ji)
805  zp_peq_b_coef(jj) = ppeq_b_coef(ji)
806  !
807  zp_lat(jj) = g%XLAT(ji)
808  zp_lon(jj) = g%XLON(ji)
809 
810  zp_zenith(jj) = pzenith(ji)
811 !
812  zp_grndflux(jj) = pgrndflux(ji)
813  zp_delheatn(jj) = pdelheatn(ji)
814  zp_delheatn_sfc(jj) = pdelheatn_sfc(ji)
815  zp_snowsfch(jj) = psnowsfch(ji)
816  zp_les3l(jj) = ples3l(ji)
817  zp_lel3l(jj) = plel3l(ji)
818  zp_evap(jj) = pevap(ji)
819  !
820  zp_swnetsnow(jj) = zswnet_n(ji)
821  zp_swnetsnows(jj) = zswnet_ns(ji)
822  zp_lwnetsnow(jj) = zlwnet_n(ji)
823 ENDDO
824 !
825 DO jj=1,ksize1
826  ji = kmask(jj)
827  zp_vegtype(jj) = pvegtype(ji,nvt_snow)
828  zp_forest(jj) = pvegtype(ji,nvt_tebd) + pvegtype(ji,nvt_trbe) + pvegtype(ji,nvt_bone) &
829  + pvegtype(ji,nvt_trbd) + pvegtype(ji,nvt_tebe) + pvegtype(ji,nvt_tene) &
830  + pvegtype(ji,nvt_bobd) + pvegtype(ji,nvt_bond) + pvegtype(ji,nvt_shrb)
831 ENDDO
832 !
833 !
834 ! ===============================================================
835 ! conversion of snow heat from J/m3 into J/m2
836 WHERE(zp_snowswe(:,:)>0.) &
837  zp_snowheat(:,:) = zp_snowheat(:,:) / zp_snowrho(:,:) * zp_snowswe(:,:)
838 ! ===============================================================
839 !
840 zp_psn_inv(:) = 0.0
841 zp_psn(:) = zp_psn3l(:)
842 !
843 IF(omeb)THEN
844 !
845 ! MEB (case of imposed surface fluxes)
846 ! - Prepare inputs for explicit snow scheme(s):
847 ! If using MEB, these are INPUTs ONLY:
848 ! divide fluxes by snow fraction to make "snow-relative"
849 !
850  zp_psn(:) = max(1.e-4, zp_psn3l(:))
851  zp_psn_inv(:) = 1.0/zp_psn(:)
852 !
853  zp_rnsnow(:) = zp_rnsnow(:) *zp_psn_inv(:)
854  zp_swnetsnow(:) = zp_swnetsnow(:) *zp_psn_inv(:)
855  zp_swnetsnows(:) = zp_swnetsnows(:) *zp_psn_inv(:)
856  zp_lwnetsnow(:) = zp_lwnetsnow(:) *zp_psn_inv(:)
857  zp_hsnow(:) = zp_hsnow(:) *zp_psn_inv(:)
858  zp_gfluxsnow(:) = zp_gfluxsnow(:) *zp_psn_inv(:)
859  zp_gsfcsnow(:) = zp_gsfcsnow(:) *zp_psn_inv(:)
860  zp_snowhmass(:) = zp_snowhmass(:) *zp_psn_inv(:)
861  zp_les3l(:) = zp_les3l(:) *zp_psn_inv(:)
862  zp_lel3l(:) = zp_lel3l(:) *zp_psn_inv(:)
863  zp_grndflux(:) = zp_grndflux(:) *zp_psn_inv(:)
864  zp_evap(:) = zp_evap(:) *zp_psn_inv(:)
865  zp_hpsnow(:) = zp_hpsnow(:) *zp_psn_inv(:)
866  zp_delheatn(:) = zp_delheatn(:) *zp_psn_inv(:)
867  zp_delheatn_sfc(:)= zp_delheatn_sfc(:)*zp_psn_inv(:)
868  zp_snowsfch(:) = zp_snowsfch(:) *zp_psn_inv(:)
869 
870  zp_srsnow(:) = zp_srsnow(:) *zp_psn_inv(:)
871  zp_rrsnow(:) = zp_rrsnow(:) *zp_psn_inv(:)
872 
873  DO jj=1,ksize2
874  DO ji=1,ksize1
875  zp_snowswe(ji,jj) = zp_snowswe(ji,jj) *zp_psn_inv(ji)
876  zp_snowheat(ji,jj) = zp_snowheat(ji,jj)*zp_psn_inv(ji)
877  zp_snowdz(ji,jj) = zp_snowdz(ji,jj) *zp_psn_inv(ji)
878  ENDDO
879  ENDDO
880  !
881 ENDIF
882 !
883 ! Call ISBA-SNOW3L model:
884 !
885 IF (pek%TSNOW%SCHEME=='CRO') THEN
886  CALL snowcro(io%CSNOWRES, tptime, io%LGLACIER, himplicit_wind, &
887  zp_pew_a_coef, zp_pew_b_coef, zp_pet_a_coef, zp_peq_a_coef,&
888  zp_pet_b_coef, zp_peq_b_coef, zp_snowswe, zp_snowrho, &
889  zp_snowheat, zp_snowalb, zp_snowgran1, zp_snowgran2, &
890  zp_snowhist, zp_snowage, ptstep, zp_ps, zp_srsnow, &
891  zp_rrsnow, zp_psn3l, zp_ta, zp_tg(:,1), zp_sw_rad, zp_qa, &
892  zp_vmod, zp_lw_rad, zp_rhoa, zp_uref, zp_exns, zp_exna, &
893  zp_dircoszw, zp_zref, zp_z0nat, zp_z0eff, zp_z0hnat, &
894  zp_alb, zp_soilcond, zp_d_g(:,1), zp_snowliq, zp_snowtemp, &
895  zp_snowdz, zp_thrufal, zp_grndflux, zp_evapcor, zp_rnsnow, &
896  zp_hsnow, zp_gfluxsnow, zp_hpsnow, zp_les3l, zp_lel3l, &
897  zp_evap, zp_sndrift, zp_ri,zp_emisnow, zp_cdsnow, &
898  zp_ustarsnow, zp_chsnow, zp_snowhmass, zp_qs, zp_vegtype, &
899  zp_zenith, zp_lat, zp_lon, io%LSNOWDRIFT, &
900  io%LSNOWDRIFT_SUBLIM, io%LSNOW_ABS_ZENITH, io%CSNOWMETAMO, &
901  io%CSNOWRAD )
902 !
903  zp_gflxcor(:) = 0.0
904  zp_flsn_cor(:) = 0.0
905  zp_soilcor(:) = 0.0
906 !
907 #ifndef SFX_OL
908  ! En couplé il faudra voir si on veut virer les diagnostics, les calculer tout le temps, ou trouver une autre solution
909  gcomputecrodiag = (SIZE(dmk%XSNOWDEND)>0)
910 #else
911  gcomputecrodiag = (SIZE(dmk%XSNOWDEND)>0).AND.(mod(tptime%TIME,xtstep_output)==0.)
912 #endif
913  !
914  !Ajout test sur pas de temps de sortie
915  IF (gcomputecrodiag) THEN
916  CALL snowcro_diag(io%CSNOWMETAMO,&
917  zp_snowdz, zp_snowswe, zp_snowrho, zp_snowgran1, zp_snowgran2, zp_snowage, &
918  zp_snowhist, zp_snowtemp, zp_snowliq, zp_dircoszw, zp_snowdend, zp_snowspher, &
919  zp_snowsize, zp_snowssa, zp_snowtypemepra, zp_snowram, zp_snowshear, &
920  zp_sndpt_1dy, zp_sndpt_3dy, zp_sndpt_5dy, zp_sndpt_7dy, zp_snswe_1dy, &
921  zp_snswe_3dy, zp_snswe_5dy, zp_snswe_7dy, zp_snram_sonde, zp_sn_wetthckn, &
922  zp_sn_refrznthckn )
923  ENDIF
924  !
925 ELSE
926 !
927  CALL snow3l(io%CSNOWRES, tptime, omeb, himplicit_wind, &
928  zp_pew_a_coef, zp_pew_b_coef, &
929  zp_pet_a_coef, zp_peq_a_coef,zp_pet_b_coef, zp_peq_b_coef, &
930  zp_snowswe, zp_snowrho, zp_snowheat, zp_snowalb, &
931  zp_snowgran1, zp_snowgran2, zp_snowhist, zp_snowage, ptstep, &
932  zp_ps, zp_srsnow, zp_rrsnow, zp_psn3l, zp_ta, zp_tg(:,1), &
933  zp_sw_rad, zp_qa, zp_vmod, zp_lw_rad, zp_rhoa, zp_uref, &
934  zp_exns, zp_exna, zp_dircoszw, zp_zref, zp_z0nat, zp_z0eff, &
935  zp_z0hnat, zp_alb, zp_soilcond, zp_d_g(:,1), &
936  zp_lvtt, zp_lstt, zp_snowliq, &
937  zp_snowtemp, zp_snowdz, zp_thrufal, zp_grndflux , &
938  zp_evapcor, zp_soilcor, zp_gflxcor, zp_snowsfch, &
939  zp_delheatn, zp_delheatn_sfc, &
940  zp_swnetsnow, zp_swnetsnows, zp_lwnetsnow, zp_gsfcsnow, &
941  zp_rnsnow, zp_hsnow, zp_gfluxsnow, zp_hpsnow, zp_les3l, &
942  zp_lel3l, zp_evap, zp_sndrift, zp_ri, &
943  zp_emisnow, zp_cdsnow, zp_ustarsnow, &
944  zp_chsnow, zp_snowhmass, zp_qs, zp_vegtype, zp_forest, &
945  zp_zenith, zp_lat, zp_lon, io%LSNOWDRIFT, io%LSNOWDRIFT_SUBLIM )
946 !
947  IF(omeb)THEN
948 !
949 ! - reverse transform: back to surface-relative
950 !
951  zp_rnsnow(:) = zp_rnsnow(:) /zp_psn_inv(:)
952  zp_swnetsnow(:) = zp_swnetsnow(:) /zp_psn_inv(:)
953  zp_swnetsnows(:) = zp_swnetsnows(:) /zp_psn_inv(:)
954  zp_lwnetsnow(:) = zp_lwnetsnow(:) /zp_psn_inv(:)
955  zp_hsnow(:) = zp_hsnow(:) /zp_psn_inv(:)
956  zp_les3l(:) = zp_les3l(:) /zp_psn_inv(:)
957  zp_lel3l(:) = zp_lel3l(:) /zp_psn_inv(:)
958  zp_grndflux(:) = zp_grndflux(:) /zp_psn_inv(:)
959  zp_evap(:) = zp_evap(:) /zp_psn_inv(:)
960  zp_hpsnow(:) = zp_hpsnow(:) /zp_psn_inv(:)
961  zp_gfluxsnow(:) = zp_gfluxsnow(:) /zp_psn_inv(:)
962  zp_delheatn(:) = zp_delheatn(:) /zp_psn_inv(:)
963  zp_delheatn_sfc(:)= zp_delheatn_sfc(:)/zp_psn_inv(:)
964  zp_snowsfch(:) = zp_snowsfch(:) /zp_psn_inv(:)
965  zp_gsfcsnow(:) = zp_gsfcsnow(:) /zp_psn_inv(:)
966 
967  zp_srsnow(:) = zp_srsnow(:) /zp_psn_inv(:)
968  zp_rrsnow(:) = zp_rrsnow(:) /zp_psn_inv(:)
969  DO jj=1,ksize2
970  DO ji=1,ksize1
971  zp_snowswe(ji,jj) = zp_snowswe(ji,jj) /zp_psn_inv(ji)
972  zp_snowheat(ji,jj) = zp_snowheat(ji,jj)/zp_psn_inv(ji)
973  zp_snowdz(ji,jj) = zp_snowdz(ji,jj) /zp_psn_inv(ji)
974  ENDDO
975  ENDDO
976 
977  zp_snowhmass(:) = zp_snowhmass(:)/zp_psn_inv(:)
978  zp_thrufal(:) = zp_thrufal(:) /zp_psn_inv(:)
979 !
980 ! Final Adjustments:
981 ! ------------------
982 ! Add cooling/heating flux correction to underlying soil.
983 ! This term is usually active for vanishingly thin snowpacks..
984 ! it is put outside of the snow scheme owing to it's dependence on
985 ! snow fraction. It is related to a possible correction to the ground-snow
986 ! heat flux when it is imposed (using MEB).
987 ! Also, it is added as a heat sink/source here since
988 ! fluxes have already be computed and should not be adjusted at this point:
989 ! applying it to the soil has the same impact as soil freeze-thaw, in the
990 ! sense it is computed after the fluxes have been updated.
991 ! (and update heat storage diagnostic in a consistent manner)
992 !
993 ! Energy is thickness weighted, thus thicker layers receive more energy and energy
994 ! is evenly distributed to depth ZDEPTHABS. An
995 ! alternate method is to weight near surface layers more and diminish weights
996 ! (thus eenrgy received by each layer) with depth. Both methods conserve energy as
997 ! long as vertical weights are normalized.
998 
999 ! i) Determine soil depth for energy absorption:
1000 
1001  zp_soild(:) = zp_dzg(:,1)
1002  DO jj=2,ksize3
1003  DO ji=1,ksize1
1004  IF(zp_dzg(ji,jj) <= zdepthabs)THEN
1005  zp_soild(ji) = zp_dzg(ji,jj)
1006  ENDIF
1007  ENDDO
1008  ENDDO
1009 
1010 ! ii) Distribute (possible) energy to absorb vertically over some layer (defined above):
1011 
1012  zp_psn_gflxcor(:) = zp_psn(:)*zp_gflxcor(:) ! (W/m2)
1013  zp_work(:) = zp_psn_gflxcor(:)*ptstep/zp_soild(:)
1014 
1015  zp_tg(:,1) = zp_tg(:,1) + zp_work(:)*zp_ct(:)*zp_d_g(:,1) ! (K)
1016  DO jj=2,ksize3
1017  DO ji=1,ksize1
1018  IF (zp_soild(ji) <= zdepthabs) THEN
1019  zp_tg(ji,jj) = zp_tg(ji,jj) + zp_work(ji)/zp_soilhcapz(ji,jj) ! K
1020  ENDIF
1021  ENDDO
1022  ENDDO
1023 
1024  zp_delheatg(:) = zp_delheatg(:) + zp_psn_gflxcor(:) ! (W/m2)
1025  zp_delheatg_sfc(:) = zp_delheatg_sfc(:) + zp_psn_gflxcor(:) ! (W/m2)
1026 !
1027  zp_flsn_cor(:) = 0.0
1028 !
1029  ELSE
1030 !
1031 ! To conserve energy in ISBA, the correction flux must be distributed at least
1032 ! over the first 60cm depth. This method prevent numerical oscillations
1033 ! especially when explicit snow vanishes. Final Adjustments are done in ISBA_CEB
1034 !
1035  zp_flsn_cor(:) = zp_gflxcor(:) ! (W/m2)
1036 !
1037  ENDIF
1038 !
1039 ENDIF
1040 !
1041 !===============================================================
1042 !conversion of snow heat from J/m2 into J/m3
1043 WHERE(zp_snowswe(:,:)>0.)
1044  zp_snowheat(:,:)=zp_snowheat(:,:)*zp_snowrho(:,:)/zp_snowswe(:,:)
1045 ENDWHERE
1046 !===============================================================
1047 !
1048 ! === Packing:
1049 !
1050 ! unpack variables
1051 !
1052 DO jwrk=1,ksize2
1053  DO jj=1,ksize1
1054  ji = kmask(jj)
1055  pek%TSNOW%WSNOW(ji,jwrk) = zp_snowswe(jj,jwrk)
1056  pek%TSNOW%RHO (ji,jwrk) = zp_snowrho(jj,jwrk)
1057  pek%TSNOW%HEAT (ji,jwrk) = zp_snowheat(jj,jwrk)
1058  pek%TSNOW%AGE (ji,jwrk) = zp_snowage(jj,jwrk)
1059  dmk%XSNOWTEMP(ji,jwrk) = zp_snowtemp(jj,jwrk)
1060  dmk%XSNOWLIQ (ji,jwrk) = zp_snowliq(jj,jwrk)
1061  dmk%XSNOWDZ (ji,jwrk) = zp_snowdz(jj,jwrk)
1062  ENDDO
1063 ENDDO
1064 !
1065 IF (pek%TSNOW%SCHEME=='CRO') THEN
1066  DO jwrk=1,ksize2
1067  DO jj=1,ksize1
1068  ji = kmask(jj)
1069  pek%TSNOW%GRAN1(ji,jwrk) = zp_snowgran1(jj,jwrk)
1070  pek%TSNOW%GRAN2(ji,jwrk) = zp_snowgran2(jj,jwrk)
1071  pek%TSNOW%HIST (ji,jwrk) = zp_snowhist(jj,jwrk)
1072  ENDDO
1073  ENDDO
1074 
1075  IF (SIZE(dmk%XSNOWDEND)>0) THEN
1076  ! This is equivalent to test the value of DGMI%LPROSNOW which does not enter in ISBA
1077  DO jwrk = 1,ksize2
1078  DO jj=1,ksize1
1079  ji = kmask(jj)
1080  dmk%XSNOWDEND (ji,jwrk) = zp_snowdend(jj,jwrk)
1081  dmk%XSNOWSPHER (ji,jwrk) = zp_snowspher(jj,jwrk)
1082  dmk%XSNOWSIZE (ji,jwrk) = zp_snowsize(jj,jwrk)
1083  dmk%XSNOWSSA (ji,jwrk) = zp_snowssa(jj,jwrk)
1084  dmk%XSNOWTYPEMEPRA(ji,jwrk) = zp_snowtypemepra(jj,jwrk)
1085  dmk%XSNOWRAM (ji,jwrk) = zp_snowram(jj,jwrk)
1086  dmk%XSNOWSHEAR (ji,jwrk) = zp_snowshear(jj,jwrk)
1087  ENDDO
1088  ENDDO
1089  ENDIF
1090 
1091 ENDIF
1092 !
1093 DO jwrk=1,ksize3
1094  DO jj=1,ksize1
1095  ji = kmask(jj)
1096  ptg(ji,jwrk)= zp_tg(jj,jwrk)
1097  ENDDO
1098 ENDDO
1099 !
1100 DO jj=1,ksize1
1101  ji = kmask(jj)
1102  pek%TSNOW%ALB(ji) = zp_snowalb(jj)
1103  pek%TSNOW%EMIS(ji) = zp_emisnow(jj)
1104  dmk%XCDSNOW (ji) = zp_cdsnow(jj)
1105  dmk%XUSTARSNOW(ji) = zp_ustarsnow(jj)
1106  dmk%XCHSNOW (ji) = zp_chsnow(jj)
1107  dmk%XSNOWHMASS(ji) = zp_snowhmass(jj)
1108  dmk%XRNSNOW (ji) = zp_rnsnow(jj)
1109  dmk%XHSNOW (ji) = zp_hsnow(jj)
1110  dmk%XHPSNOW (ji) = zp_hpsnow(jj)
1111  dmk%XGFLUXSNOW(ji) = zp_gfluxsnow(jj)
1112  !
1113  pdelheatg(ji) = zp_delheatg(jj)
1114  pdelheatg_sfc(ji) = zp_delheatg_sfc(jj)
1115  pthrufal(ji) = zp_thrufal(jj)
1116  pevapcor(ji) = zp_evapcor(jj)
1117  pri(ji) = zp_ri(jj)
1118  pqs(ji) = zp_qs(jj)
1119  pgrndflux(ji) = zp_grndflux(jj)
1120  pflsn_cor(ji) = zp_flsn_cor(jj)
1121  pdelheatn(ji) = zp_delheatn(jj)
1122  pdelheatn_sfc(ji) = zp_delheatn_sfc(jj)
1123  psnowsfch(ji) = zp_snowsfch(jj)
1124  pgsfcsnow(ji) = zp_gsfcsnow(jj)
1125  ples3l(ji) = zp_les3l(jj)
1126  plel3l(ji) = zp_lel3l(jj)
1127  pevap(ji) = zp_evap(jj)
1128  zsoilcor(ji) = zp_soilcor(jj)
1129  !
1130  zswnet_n(ji) = zp_swnetsnow(jj)
1131  zswnet_ns(ji) = zp_swnetsnows(jj)
1132  zlwnet_n(ji) = zp_lwnetsnow(jj)
1133 ENDDO
1134 !
1135 IF ( SIZE(dmk%XSNOWDEND)>0 ) THEN
1136  ! This is equivalent to test the value of DGMI%LPROSNOW which does not enter in ISBATHEN
1137  dmk%XSNDPT_1DY(:) = xundef
1138  dmk%XSNDPT_3DY(:) = xundef
1139  dmk%XSNDPT_5DY(:) = xundef
1140  dmk%XSNDPT_7DY(:) = xundef
1141  dmk%XSNSWE_1DY(:) = xundef
1142  dmk%XSNSWE_3DY(:) = xundef
1143  dmk%XSNSWE_5DY(:) = xundef
1144  dmk%XSNSWE_7DY(:) = xundef
1145  dmk%XSNRAM_SONDE (:) = xundef
1146  dmk%XSN_WETTHCKN (:) = xundef
1147  dmk%XSN_REFRZNTHCKN(:) = xundef
1148  DO jj=1,ksize1
1149  ji = kmask(jj)
1150  dmk%XSNDPT_1DY(ji) = zp_sndpt_1dy(jj)
1151  dmk%XSNDPT_3DY(ji) = zp_sndpt_3dy(jj)
1152  dmk%XSNDPT_5DY(ji) = zp_sndpt_5dy(jj)
1153  dmk%XSNDPT_7DY(ji) = zp_sndpt_7dy(jj)
1154  dmk%XSNSWE_1DY(ji) = zp_snswe_1dy(jj)
1155  dmk%XSNSWE_3DY(ji) = zp_snswe_3dy(jj)
1156  dmk%XSNSWE_5DY(ji) = zp_snswe_5dy(jj)
1157  dmk%XSNSWE_7DY(ji) = zp_snswe_7dy(jj)
1158  dmk%XSNRAM_SONDE (ji) = zp_snram_sonde(jj)
1159  dmk%XSN_WETTHCKN (ji) = zp_sn_wetthckn(jj)
1160  dmk%XSN_REFRZNTHCKN(ji) = zp_sn_refrznthckn(jj)
1161  ENDDO
1162 ENDIF
1163 !
1164 IF (lhook) CALL dr_hook('SNOW3L_ISBA:CALL_MODEL',1,zhook_handle)
1165 !
1166 END SUBROUTINE call_model
1167 !
1168 END SUBROUTINE snow3l_isba
modd_surf_par
Definition: modd_surf_par.F90:6
modd_isba_n
Definition: modd_isban.F90:6
yomhook::dr_hook
Definition: yomhook.F90:20
call_model
subroutine call_model(KSIZE1, KSIZE2, KSIZE3, KMASK)
Definition: snow3L_isba.F90:588
modd_csts::xpi
real, save xpi
Definition: modd_csts.F90:43
snowcro_diag
subroutine snowcro_diag(HSNOWMETAMO, PSNOWDZ, PSNOWSWE, PSNOWRHO, PSNOWGRAN1, PSNOWGRAN2, PSNOWAGE, PSNOWHIST, PSNOWTEMP, PSNOWLIQ, PDIRCOSZW, PSNOWDEND, PSNOWSPHER, PSNOWSIZE, PSNOWSSA, PSNOWTYPEMEPRA, PSNOWRAM, PSNOWSHEAR, PSNOWDEPTH_1DAYS, PSNOWDEPTH_3DAYS, PSNOWDEPTH_5DAYS, PSNOWDEPTH_7DAYS, PSNOWSWE_1DAYS, PSNOWSWE_3DAYS, PSNOWSWE_5DAYS, PSNOWSWE_7DAYS, PSNOWRAM_SONDE, PSNOW_WETTHICKNESS, PSNOW_REFROZENTHICKNESS)
Definition: snowcro_diag.F90:10
modd_csts::xlstt
real, save xlstt
Definition: modd_csts.F90:71
modd_diag_evap_isba_n
Definition: modd_diag_evap_isban.F90:6
modd_diag_misc_isba_n
Definition: modd_diag_misc_isban.F90:6
abor1_sfx
subroutine abor1_sfx(YTEXT)
Definition: abor1_sfx.F90:7
modd_surf_par::xundef
real, parameter xundef
Definition: modd_surf_par.F90:43
parkind1::jprb
integer, parameter jprb
Definition: parkind1.F90:32
modd_sfx_grid_n::grid_t
Definition: modd_sfx_gridn.F90:39
modn_io_offline
Definition: modn_io_offline.F90:6
modd_isba_options_n::isba_options_t
Definition: modd_isba_optionsn.F90:49
snow3l_isba
subroutine snow3l_isba(IO, G, PK, PEK, DK, DEK, DMK, OMEB, HIMPLICIT_WIND, TPTIME, PTSTEP, PVEGTYPE, PTG, PCT, PSOILHCAPZ, PSOILCONDZ, PPS, PTA, PSW_RAD, PQA, PVMOD, PLW_RAD, PRR, PSR, PRHOA, PUREF, PEXNS, PEXNA, PDIRCOSZW, PZREF, PALB, PD_G, PDZG, PPEW_A_COEF, PPEW_B_COEF, PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF, PTHRUFAL, PGRNDFLUX, PFLSN_COR, PGSFCSNOW, PEVAPCOR, PLES3L, PLEL3L, PEVAP, PSNOWSFCH, PDELHEATN, PDELHEATN_SFC, PRI, PZENITH, PDELHEATG, PDELHEATG_SFC, PQS)
Definition: snow3L_isba.F90:15
modd_csts::xday
real, save xday
Definition: modd_csts.F90:45
modd_isba_options_n
Definition: modd_isba_optionsn.F90:6
modd_diag_n::diag_t
Definition: modd_diagn.F90:81
modd_diag_evap_isba_n::diag_evap_isba_t
Definition: modd_diag_evap_isban.F90:40
modd_sfx_grid_n
Definition: modd_sfx_gridn.F90:6
modd_type_date_surf::date_time
Definition: modd_type_date_surf.F90:46
modd_diag_misc_isba_n::diag_misc_isba_t
Definition: modd_diag_misc_isban.F90:53
yomhook::lhook
logical lhook
Definition: yomhook.F90:15
snow3l
subroutine snow3l(HSNOWRES, TPTIME, OMEB, HIMPLICIT_WIND,
Definition: snow3l.F90:7
snowcro
subroutine snowcro(HSNOWRES, TPTIME, OGLACIER, HIMPLICIT_WIND,
Definition: snowcro.F90:7
modn_io_offline::xtstep_output
real xtstep_output
Definition: modn_io_offline.F90:96
modd_isba_n::isba_pe_t
Definition: modd_isban.F90:341
modd_csts
Definition: modd_csts.F90:6
modd_csts::xtt
real, save xtt
Definition: modd_csts.F90:66
modd_csts::xlmtt
real, save xlmtt
Definition: modd_csts.F90:72
modd_type_date_surf
Definition: modd_type_date_surf.F90:6
modd_diag_n
Definition: modd_diagn.F90:6
modd_isba_n::isba_p_t
Definition: modd_isban.F90:219
Generated on Tue Jan 16 2018 16:23:30 for SURFEX v8.1 by   doxygen 1.8.13