SURFEX v8.1
General documentation of Surfex
isba_meb.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 isba_meb(IO, KK, PK, PEK, DK, DEK, DMK, G, AG, &
7  TPTIME, OMEB, OSHADE, HIMPLICIT_WIND, PTSTEP, &
8  PSOILHCAPZ, PSOILCONDZ, PFROZEN1, PPS, PZENITH, &
9  PSCA_SW, PSW_RAD, PVMOD, PRR, PSR, PRHOA, PTA, PQA, &
10  PDIRCOSZW, PEXNS, PEXNA, PPET_A_COEF, PPET_B_COEF, &
11  PPEQ_A_COEF, PPEQ_B_COEF, PPEW_A_COEF, PPEW_B_COEF, &
12  PZREF, PUREF, PZ0G_WITHOUT_SNOW, PZ0_MEBV, PZ0H_MEBV, &
13  PZ0EFF_MEBV, PZ0_MEBN, PZ0H_MEBN, PZ0EFF_MEBN, &
14  PALBNIR_TVEG, PALBVIS_TVEG, PALBNIR_TSOIL, PALBVIS_TSOIL, &
15  PABC, PIACAN, PPOI, PCSP, PRESP_BIOMASS_INST, PPALPHAN, &
16  PF2, PLW_RAD, PGRNDFLUX, PFLSN_COR, PUSTAR, PEMIST, &
17  PHU_AGG, PAC_AGG, PDELHEATV_SFC, PDELHEATG_SFC, PDELHEATG,&
18  PDELHEATN, PDELHEATN_SFC, PRESTOREN, PTDEEP_A, PDEEP_FLUX,&
19  PRISNOW, PSNOW_THRUFAL, PSNOW_THRUFAL_SOIL, PEVAPCOR, &
20  PSUBVCOR, PLITCOR, PSNOWSFCH, PQSNOW)
21 ! ##########################################################################
22 !
23 !
24 !!**** *isba_meb*
25 !!
26 !! PURPOSE
27 !! -------
28 ! Monitor for the calculation of the surface fluxes and of the
29 ! prognostic variables of the surface over natural areas
30 ! with an explicit vegetation layer
31 !
32 ! NOTE...currently MEB can be coupled with
33 ! IO%CISBA='DIF' or '3-L' soil options
34 ! HSNOW='3-L' snow scheme
35 ! Soon, HSNOW=CRO and IO%CPHOTO/=NON (i.e. Ags will be added)
36 !
37 !!** METHOD
38 !! ------
39 !
40 !! EXTERNAL
41 !! --------
42 !!
43 !! IMPLICIT ARGUMENTS
44 !! ------------------
45 !!
46 !!
47 !! REFERENCE
48 !! ---------
49 !!
50 !! Noilhan and Planton (1989)
51 !!
52 !! AUTHOR
53 !! ------
54 !! A. Boone * Meteo-France *
55 !! P. Samuelsson * SMHI *
56 !!
57 !! MODIFICATIONS
58 !! -------------
59 !! Original 10/2014
60 !! (A. Napoly) 09/2015 Add Litter layer option code
61 !! (A. Boone) 02/2017 Owing to fix to FAPAIR.F90 routine (called by
62 !! RAIDATIVE_TRANSFERT.F90 herein), edited slightly
63 !! SWnet computations to be compatible.
64 !!
65 !-------------------------------------------------------------------------------
66 !
67 !* 0. DECLARATIONS
68 ! ------------
69 !
72 USE modd_sfx_grid_n, ONLY : grid_t
73 USE modd_agri_n, ONLY : agri_t
74 USE modd_diag_n, ONLY : diag_t
77 !
78 USE modd_surf_par, ONLY : xundef
79 USE modd_csts, ONLY : xcpd, xday, xrholw, xlvtt, xlstt
81 USE modd_isba_par, ONLY : xrs_max, xlimh
82 USE modd_data_cover_par, ONLY : nvt_snow
83 !
85 !
86 USE mode_thermos
87 USE mode_meb, ONLY : snow_intercept_eff
88 !
89 USE modi_wet_leaves_frac
90 USE modi_veg
91 USE modi_snow_leaves_frac_meb
92 USE modi_preps_for_meb_ebud_rad
93 USE modi_isba_lwnet_meb
94 USE modi_drag_meb
95 USE modi_e_budget_meb
96 USE modi_isba_fluxes_meb
97 USE modi_snow_load_meb
98 USE modi_hydro_veg
99 USE modi_snow3l_isba
100 USE modi_radiative_transfert
101 USE modi_cotwores
102 !
103 !
104 USE yomhook ,ONLY : lhook, dr_hook
105 USE parkind1 ,ONLY : jprb
106 !
107 IMPLICIT NONE
108 !
109 !* 0.1 declarations of arguments
110 ! -------------------------
111 !
112 !
113 !* general variables
114 ! -----------------
115 !
116 TYPE(isba_options_t), INTENT(INOUT) :: IO
117 TYPE(isba_k_t), INTENT(INOUT) :: KK
118 TYPE(isba_p_t), INTENT(INOUT) :: PK
119 TYPE(isba_pe_t), INTENT(INOUT) :: PEK
120 TYPE(grid_t), INTENT(INOUT) :: G
121 TYPE(agri_t), INTENT(INOUT) :: AG
122 TYPE(diag_t), INTENT(INOUT) :: DK
123 TYPE(diag_evap_isba_t), INTENT(INOUT) :: DEK
124 TYPE(diag_misc_isba_t), INTENT(INOUT) :: DMK
125 !
126 TYPE(date_time), INTENT(IN) :: TPTIME ! current date and time
127 !
128 LOGICAL, INTENT(IN) :: OMEB ! True = patch with multi-energy balance
129 ! ! False = patch with classical ISBA
130 LOGICAL, DIMENSION(:),INTENT(INOUT) :: OSHADE ! where vegetation evolution occurs
131  CHARACTER(LEN=*), INTENT(IN) :: HIMPLICIT_WIND! wind implicitation option
132 ! ! 'OLD' = direct
133 ! ! 'NEW' = Taylor serie, order 1
134 REAL, INTENT(IN) :: PTSTEP ! Model time step (s)
135 !
136 REAL, DIMENSION(:), INTENT(IN) :: PPS ! Pressure [Pa]
137 REAL, DIMENSION(:), INTENT(IN) :: PZENITH ! solar zenith angle
138 REAL, DIMENSION(:), INTENT(IN) :: PSW_RAD ! solar (shortwave) incoming radiation [W/m2]
139 REAL, DIMENSION(:), INTENT(IN) :: PLW_RAD ! thermal (longwave) incoming radiation [W/m2]
140 REAL, DIMENSION(:), INTENT(IN) :: PSCA_SW ! solar diffuse incoming radiation [W/m2]
141 REAL, DIMENSION(:), INTENT(IN) :: PEXNA ! Exner function: forcing level (-)
142 REAL, DIMENSION(:), INTENT(IN) :: PEXNS ! Exner function: surface (-)
143 REAL, DIMENSION(:), INTENT(IN) :: PRR ! Rain rate (kg/m2/s)
144 REAL, DIMENSION(:), INTENT(IN) :: PSR ! Snow rate (kg/m2/s)
145 REAL, DIMENSION(:), INTENT(IN) :: PRHOA ! air density (kg/m3)
146 REAL, DIMENSION(:), INTENT(IN) :: PVMOD ! modulus of the wind
147 ! ! parallel to the orography (m/s)
148 REAL, DIMENSION(:), INTENT(IN) :: PTA ! Temperature of atmosphere (K)
149 REAL, DIMENSION(:), INTENT(IN) :: PQA ! specific humidity of atmosphere (kg/kg)
150 REAL, DIMENSION(:), INTENT(IN) :: PZREF ! normal distance of the first
151 ! ! atmospheric level to the
152 ! ! orography (m)
153 REAL, DIMENSION(:), INTENT(IN) :: PUREF ! reference height of the wind (m)
154 ! ! NOTE this is different from ZZREF
155 ! ! ONLY in stand-alone/forced mode,
156 ! ! NOT when coupled to a model (MesoNH)
157 REAL, DIMENSION(:), INTENT(IN) :: PDIRCOSZW ! Director Cosinus along the z
158 ! ! direction at the surface w-point
159 REAL, DIMENSION(:,:), INTENT(IN) :: PSOILHCAPZ ! ISBA-DF Soil heat capacity
160 ! ! profile [J/(m3 K)]
161 REAL, DIMENSION(:,:), INTENT(IN) :: PSOILCONDZ ! ISBA-DF Soil conductivity
162 ! ! profile [W/(m K)]
163 REAL, DIMENSION(:), INTENT(IN) :: PFROZEN1 ! surface frozen fraction (-)
164 !
165 REAL, DIMENSION(:), INTENT(IN) :: PPALPHAN ! snow/canopy transition coefficient
166 REAL, DIMENSION(:), INTENT(IN) :: PALBNIR_TVEG ! albedo of vegetation in NIR
167 ! ! (needed for LM_TR or MEB)
168 REAL, DIMENSION(:), INTENT(IN) :: PALBVIS_TVEG ! albedo of vegetation in VIS
169 ! ! (needed for LM_TR or MEB)
170 REAL, DIMENSION(:), INTENT(IN) :: PALBNIR_TSOIL ! albedo of bare soil in NIR
171 ! ! (needed for LM_TR or MEB)
172 REAL, DIMENSION(:), INTENT(IN) :: PALBVIS_TSOIL ! albedo of bare soil in VIS
173 REAL, DIMENSION(:), INTENT(IN) :: PF2 ! Soil water stress factor for transpiration (-)
174 REAL, DIMENSION(:), INTENT(IN) :: PZ0G_WITHOUT_SNOW ! roughness length for momentum at snow-free canopy floor (m)
175 REAL, DIMENSION(:), INTENT(IN) :: PZ0_MEBV ! roughness length for momentum over MEB vegetation part of patch (m)
176 REAL, DIMENSION(:), INTENT(IN) :: PZ0H_MEBV ! roughness length for heat over MEB vegetation part of path (m)
177 REAL, DIMENSION(:), INTENT(IN) :: PZ0EFF_MEBV ! roughness length for momentum over MEB vegetation part of patch (m)
178 REAL, DIMENSION(:), INTENT(IN) :: PZ0_MEBN ! roughness length for momentum over MEB snow part of patch (m)
179 REAL, DIMENSION(:), INTENT(IN) :: PZ0H_MEBN ! roughness length for heat over MEB snow part of path (m)
180 REAL, DIMENSION(:), INTENT(IN) :: PZ0EFF_MEBN ! roughness length for momentum over MEB snow part of patch (m)
181 !
182 ! implicit atmospheric coupling coefficients:
183 !
184 REAL, DIMENSION(:), INTENT(IN) :: PPET_A_COEF, PPET_B_COEF, &
185  PPEQ_A_COEF, PPEQ_B_COEF, &
186  PPEW_A_COEF, PPEW_B_COEF
187 ! ! PPEW_A_COEF A-wind coefficient
188 ! ! PPEW_B_COEF B-wind coefficient
189 ! ! PPET_A_COEF A-air temperature coefficient
190 ! ! PPET_B_COEF B-air temperature coefficient
191 ! ! PPEQ_A_COEF A-air specific humidity coefficient
192 ! ! PPEQ_B_COEF B-air specific humidity coefficient
193 REAL, DIMENSION(:), INTENT(IN) :: PTDEEP_A ! Deep soil temperature boundary condition
194 ! ! (prescribed)
195 ! PTDEEP_A = Deep soil temperature
196 ! coefficient depending on flux
197 !
198 ! ISBA-Ags parameters
199 ! (see also parameters with 'Ags:' in comments)
200 !
201 REAL, DIMENSION(:), INTENT(IN) :: PCSP ! atmospheric CO2 concentration
202 ! [ppmm]=[kg CO2 / kg air]
203 REAL, DIMENSION(:), INTENT(IN) :: PPOI ! Gaussian weights (as above)
204 ! ! temperature
205 ! - - - - - - - - - - - - - - - - - - - -
206 !
207 REAL, DIMENSION(:), INTENT(INOUT) :: PABC ! Ags: abscissa needed for integration
208 ! ! of net assimilation and stomatal
209 ! ! conductance over canopy depth
210 REAL, DIMENSION(:,:), INTENT(OUT) :: PIACAN ! PAR in the canopy at different gauss levels
211 ! ! when using the DIF soil option (W/m2)
212 REAL, DIMENSION(:), INTENT(OUT) :: PUSTAR ! friction velocity
213 !
214 REAL, DIMENSION(:), INTENT(OUT) :: PGRNDFLUX ! snow/soil-biomass interface flux (W/m2)
215 REAL, DIMENSION(:), INTENT(OUT) :: PFLSN_COR ! soil/snow interface correction flux to conserve energy (W/m2)
216 REAL, DIMENSION(:), INTENT(OUT) :: PEMIST ! total effective surface emissivity...LWUP = EMIST*TS_RAD**4 (-)
217 REAL, DIMENSION(:), INTENT(OUT) :: PAC_AGG ! aggregated aerodynamic conductance
218  ! for evaporative flux calculations
219 REAL, DIMENSION(:), INTENT(OUT) :: PHU_AGG ! aggregated relative humidity
220  ! for evaporative flux calculations
221 REAL, DIMENSION(:), INTENT(OUT) :: PDELHEATV_SFC ! change in heat storage of the vegetation canopy layer over the current time step (W/m2)
222 REAL, DIMENSION(:), INTENT(OUT) :: PDELHEATG_SFC ! change in heat storage of the ground sfc layer over the current time step (W/m2)
223 REAL, DIMENSION(:), INTENT(OUT) :: PDELHEATG ! change in heat storage of the entire soil column over the current time step (W/m2)
224 REAL, DIMENSION(:), INTENT(OUT) :: PRESTOREN ! conductive heat flux between the surface and sub-surface soil layers
225 ! ! for the multi-layer snow schemes..for composite snow, it is
226 ! ! equal to DEK%XRESTORE (W/m2)
227 REAL, DIMENSION(:), INTENT(OUT) :: PDELHEATN ! change in heat storage of the entire snow column over the current time step (W/m2)
228 REAL, DIMENSION(:), INTENT(OUT) :: PDELHEATN_SFC ! change in heat storage of the surface snow layer over the current time step (W/m2)
229 REAL, DIMENSION(:), INTENT(OUT) :: PDEEP_FLUX ! Heat flux at bottom of ISBA (W/m2)
230 REAL, DIMENSION(:), INTENT(OUT) :: PRISNOW ! Richarson number over ground-based snowpack (-)
231 REAL, DIMENSION(:), INTENT(OUT) :: PSNOW_THRUFAL ! rate that liquid water leaves (explicit) snow pack:
232 ! ! ISBA-ES or CROCUS [kg/(m2 s)]
233 REAL, DIMENSION(:), INTENT(OUT) :: PSNOW_THRUFAL_SOIL !liquid water leaving the snowpack directly to the
234 ! !soil, ISBA-ES: [kg/(m2 s)] (equal to ZSNOW_THRUFAL
235 ! !if OMEB_LITTER=False and zero if OMEB_LITTER=True)
236 ! ! ISBA-ES or CROCUS [kg/(m2 s)]
237 REAL, DIMENSION(:), INTENT(OUT) :: PEVAPCOR ! evaporation correction as last traces of snow
238 ! ! cover ablate..if sublimation exceeds trace amounts
239  ! of snow during time step, required residual mass taken
240  ! from sfc soil layer [kg/(m2 s)]
241 REAL, DIMENSION(:), INTENT(OUT) :: PSUBVCOR ! A possible snow mass correction (to be potentially
242 ! ! removed from soil) (kg/m2/s)
243 REAL, DIMENSION(:), INTENT(OUT) :: PLITCOR ! A possible ice mass correction in litter layer (to be potentially
244 ! ! removed from soil) (kg/m2/s)
245 REAL, DIMENSION(:), INTENT(OUT) :: PSNOWSFCH ! snow surface layer pseudo-heating term owing to
246 ! ! changes in grid thickness (W/m2)
247 REAL, DIMENSION(:), INTENT(OUT) :: PQSNOW ! snow surface specific humidity (kg/kg)
248 !
249 ! diagnostic variables for Carbon assimilation:
250 !
251 REAL, DIMENSION(:,:), INTENT(OUT) :: PRESP_BIOMASS_INST ! instantaneous biomass respiration (kgCO2/kgair m/s)
252 !
253 !* 0.2 declarations of local variables
254 !
255 !
256 REAL, PARAMETER :: ZTSTEP_EB = 300. ! s Minimum time tstep required
257 ! ! to time-split MEB energy budget
258 REAL, PARAMETER :: ZSWRAD_MIN = 1.e-6! W/m2 Threshold SWdown for which Sun is up...approx
259  ! (No need to do nonsense SW computations during the night)
260 !
261 INTEGER :: JTSPLIT_EB ! number of time splits
262 INTEGER :: JDT ! time split loop index
263 !
264 REAL :: ZTSTEP ! Local time split timestep (s)
265 REAL, DIMENSION(SIZE(PPS)) :: ZWORK,ZWORK2,ZWORK3,ZWORK4 ! Working variables [*]
266 REAL, DIMENSION(SIZE(PEK%TSNOW%WSNOW,1),SIZE(PEK%TSNOW%WSNOW,2)) :: ZSNOWCOND ! snow thermal conductivity [W/(m K)]
267 REAL, DIMENSION(SIZE(PEK%TSNOW%WSNOW,1),SIZE(PEK%TSNOW%WSNOW,2)) :: ZSNOWHCAP ! snow heat capacity [J/(m3 K)]
268 REAL, DIMENSION(SIZE(PEK%TSNOW%WSNOW,1),SIZE(PEK%TSNOW%WSNOW,2)) :: ZSNOWRHO ! snow layer density (kg/m3)
269 REAL, DIMENSION(SIZE(PEK%TSNOW%WSNOW,1),SIZE(PEK%TSNOW%WSNOW,2)) :: ZSNOWAGE ! snow layer grain age
270 REAL, DIMENSION(SIZE(PEK%TSNOW%WSNOW,1),SIZE(PEK%TSNOW%WSNOW,2)) :: ZSNOWSWE ! snow layer liquid water equivalent (kg/m2)
271 REAL, DIMENSION(SIZE(PEK%TSNOW%WSNOW,1),SIZE(PEK%TSNOW%WSNOW,2)) :: ZSNOWLIQ ! snow layer liquid water (m)
272 REAL, DIMENSION(SIZE(PEK%TSNOW%WSNOW,1),SIZE(PEK%TSNOW%WSNOW,2)) :: ZTAU_N ! snow rad transmission coef at layer base (-)
273 REAL, DIMENSION(SIZE(PPS)) :: ZCHIP !
274 REAL, DIMENSION(SIZE(PPS)) :: ZALBS ! Effective surface (ground) albedo
275 REAL, DIMENSION(SIZE(PPS)) :: ZSIGMA_F ! LW transmission factor
276 REAL, DIMENSION(SIZE(PPS)) :: ZSIGMA_FN ! LW transmission factor - including buried (snow)
277 ! ! vegetation effect
278 REAL, DIMENSION(SIZE(PPS)) :: ZDLWNET_V_DTV ! LW Jacobian: flux derrivative d LWnet_v/dTv [W/(m K2)]
279 REAL, DIMENSION(SIZE(PPS)) :: ZDLWNET_V_DTG ! LW Jacobian: flux derrivative d LWnet_v/dTg [W/(m K2)]
280 REAL, DIMENSION(SIZE(PPS)) :: ZDLWNET_V_DTN ! LW Jacobian: flux derrivative d LWnet_v/dTn [W/(m K2)]
281 REAL, DIMENSION(SIZE(PPS)) :: ZDLWNET_G_DTV ! LW Jacobian: flux derrivative d LWnet_g/dTv [W/(m K2)]
282 REAL, DIMENSION(SIZE(PPS)) :: ZDLWNET_G_DTG ! LW Jacobian: flux derrivative d LWnet_g/dTg [W/(m K2)]
283 REAL, DIMENSION(SIZE(PPS)) :: ZDLWNET_G_DTN ! LW Jacobian: flux derrivative d LWnet_g/dTn [W/(m K2)]
284 REAL, DIMENSION(SIZE(PPS)) :: ZDLWNET_N_DTV ! LW Jacobian: flux derrivative d LWnet_n/dTv [W/(m K2)]
285 REAL, DIMENSION(SIZE(PPS)) :: ZDLWNET_N_DTG ! LW Jacobian: flux derrivative d LWnet_n/dTg [W/(m K2)]
286 REAL, DIMENSION(SIZE(PPS)) :: ZDLWNET_N_DTN ! LW Jacobian: flux derrivative d LWnet_n/dTn [W/(m K2)]
287 REAL, DIMENSION(SIZE(PPS)) :: ZWRMAX ! maximum canopy water equivalent interception capacity [kg/m2]
288 REAL, DIMENSION(SIZE(PPS)) :: ZWRLMAX ! maximum litter water equivalent interception capacity [kg/m2]
289 REAL, DIMENSION(SIZE(PPS)) :: ZRS ! stomatal resistance (s/m)
290 REAL, DIMENSION(SIZE(PPS)) :: ZRSN ! stomatal resistance of non-snow-buried canopy (s/m)
291 ! ! Etv=>0 as F2=>0 (-)
292 REAL, DIMENSION(SIZE(PPS)) :: ZWRVNMAX ! maximum snow water equivalent interception capacity (kg/m2)
293 REAL, DIMENSION(SIZE(PPS)) :: ZPSNCV ! intercepted canopy snow fraction (-) NOTE! Not the same as the
294 ! ! ground-based snowpack
295 REAL, DIMENSION(SIZE(PPS)) :: ZMELTVN ! intercepted canopy snow net freeze/melt rate (kg/m2/s)
296 ! ! (if it is < 0, this signifies freezing)
297 REAL, DIMENSION(SIZE(PPS)) :: ZTHRMA_TA ! linear transform energy budget coefficient for Ta
298 REAL, DIMENSION(SIZE(PPS)) :: ZTHRMB_TA ! linear transform energy budget coefficient for Ta
299 REAL, DIMENSION(SIZE(PPS)) :: ZTHRMA_TC ! linear transform energy budget coefficient for Tc
300 REAL, DIMENSION(SIZE(PPS)) :: ZTHRMB_TC ! linear transform energy budget coefficient for Tc
301 REAL, DIMENSION(SIZE(PPS)) :: ZTHRMA_TN ! linear transform energy budget coefficient for Tn
302 REAL, DIMENSION(SIZE(PPS)) :: ZTHRMB_TN ! linear transform energy budget coefficient for Tn
303 REAL, DIMENSION(SIZE(PPS)) :: ZTHRMA_TG ! linear transform energy budget coefficient for Tg
304 REAL, DIMENSION(SIZE(PPS)) :: ZTHRMB_TG ! linear transform energy budget coefficient for Tg
305 REAL, DIMENSION(SIZE(PPS)) :: ZTHRMA_TV ! linear transform energy budget coefficient for Tv
306 REAL, DIMENSION(SIZE(PPS)) :: ZTHRMB_TV ! linear transform energy budget coefficient for Tv
307 REAL, DIMENSION(SIZE(PPS)) :: ZPET_A_COEF ! atmospheric coupling coefficient: Ta
308 REAL, DIMENSION(SIZE(PPS)) :: ZPET_B_COEF ! atmospheric coupling coefficient: Ta
309 REAL, DIMENSION(SIZE(PPS)) :: ZKVN ! snow interception efficiency
310 REAL, DIMENSION(SIZE(PPS)) :: ZVELC ! wind speed at the top of the canopy (m/s)
311 REAL, DIMENSION(SIZE(PPS)) :: ZDELTA ! fraction of the foliage
312 ! ! covered with intercepted water (-)
313 REAL, DIMENSION(SIZE(PPS)) :: ZHUGI ! humidity over frozen bare ground (-)
314 REAL, DIMENSION(SIZE(PPS)) :: ZHVN ! Halstead coefficient vegetation canopy above snow (-)
315 REAL, DIMENSION(SIZE(PPS)) :: ZHVG ! Halstead coefficient vegetation canopy above snow-free ground (-)
316 REAL, DIMENSION(SIZE(PPS)) :: ZLEG_DELTA ! soil evaporation delta fn (-)
317 REAL, DIMENSION(SIZE(PPS)) :: ZLEGI_DELTA ! soil sublimation delta fn (-)
318 REAL, DIMENSION(SIZE(PPS)) :: ZHSGL ! surface halstead cofficient for bare soil (-)
319 REAL, DIMENSION(SIZE(PPS)) :: ZHSGF ! surface halstead cofficient for bare soil ice (-)
320 REAL, DIMENSION(SIZE(PPS)) :: ZFLXC_CA ! turb transfer coef between vegetation canopy air and atmosphere (kg/m2/s)
321 REAL, DIMENSION(SIZE(PPS)) :: ZFLXC_N_A ! ...between the snow on the ground and atmosphere (kg/m2/s)
322 REAL, DIMENSION(SIZE(PPS)) :: ZFLXC_GV ! ...between snow-free ground and canopy air (kg/m2/s)
323 REAL, DIMENSION(SIZE(PPS)) :: ZFLXC_GN ! ...between snow on the ground and canopy air (kg/m2/s)
324 REAL, DIMENSION(SIZE(PPS)) :: ZFLXC_VG_C ! ...between vegetation canopy over snow-free ground and canopy air (kg/m2/s)
325 REAL, DIMENSION(SIZE(PPS)) :: ZFLXC_VN_C ! ...between vegetation canopy over the snow on the ground and canopy air (kg/m2/s)
326 REAL, DIMENSION(SIZE(PPS)) :: ZFLXC_CV ! ...between vegetation canopy and canopy air (kg/m2/s)
327 REAL, DIMENSION(SIZE(PPS)) :: ZFLXC_MOM ! Effective drag coefficient for momentum [kg/(m2 s)]
328 REAL, DIMENSION(SIZE(PPS)) :: ZQSATG ! saturation specific humidity for PEK%XTG (ground surface: kg kg-1)
329 REAL, DIMENSION(SIZE(PPS)) :: ZQSATV ! saturation specific humidity for PEK%XTV (vegetation canopy: kg kg-1)
330 REAL, DIMENSION(SIZE(PPS)) :: ZQSATC ! saturation specific humidity for PEK%XTC (canopy air: kg kg-1)
331 REAL, DIMENSION(SIZE(PPS)) :: ZQSATN ! saturation specific humidity for DMK%XSNOWTEMP (snow surface: kg kg-1)
332 REAL, DIMENSION(SIZE(PPS)) :: ZDELTAVK ! canopy interception capacity fraction including K-factor (-)
333 REAL, DIMENSION(SIZE(PPS)) :: ZCHEATV ! Vegetation canopy *effective surface* heat capacity (J m-2 K-1)
334 REAL, DIMENSION(SIZE(PPS)) :: ZCHEATG ! Understory-ground *effective surface* heat capacity (J m-2 K-1)
335 REAL, DIMENSION(SIZE(PPS)) :: ZCHEATN ! Ground-based snow *effective surface* heat capacity (J m-2 K-1)
336 REAL, DIMENSION(SIZE(PPS)) :: ZHVGS ! Dimensionless pseudo humidity factor for computing
337 ! ! vapor fluxes from the non-buried part of the canopy
338 ! ! to the canopy air (-)
339 REAL, DIMENSION(SIZE(PPS)) :: ZHVNS ! Dimensionless pseudo humidity factor for computing
340 ! ! vapor fluxes from the partly-buried part of the canopy
341 ! ! to the canopy air (-)
342 REAL, DIMENSION(SIZE(PPS)) :: ZDQSAT_G ! saturation specific humidity derivative for understory (kg kg-1 K-1)
343 REAL, DIMENSION(SIZE(PPS)) :: ZDQSAT_V ! saturation specific humidity derivative for the
344 ! ! vegetation canopy (kg kg-1 K-1)
345 REAL, DIMENSION(SIZE(PPS)) :: ZDQSATI_N ! saturation specific humidity derivative over ice for
346 ! ! the ground-based snowpack (kg kg-1 K-1)
347 REAL, DIMENSION(SIZE(PPS)) :: ZDELTAT_G ! Time change in soil surface temperature (K)
348 REAL, DIMENSION(SIZE(PPS)) :: ZDELTAT_V ! Time change in vegetation canopy temperature (K)
349 REAL, DIMENSION(SIZE(PPS)) :: ZDELTAT_N ! Time change in snowpack surface temperature (K)
350 REAL, DIMENSION(SIZE(PPS)) :: ZRNET_V ! Net vegetation canopy radiation (W/m2)
351 REAL, DIMENSION(SIZE(PPS)) :: ZRNET_G ! Net understory-ground radiation (W/m2)
352 REAL, DIMENSION(SIZE(PPS)) :: ZFLXC_C_A_F ! Exchange coefficient between the snow on the ground and
353 ! ! atmosphere modified by a partially to fully buried
354 ! ! vegetation canopy [kg/(m2 s)]
355 REAL, DIMENSION(SIZE(PPS)) :: ZFLXC_N_A_F ! Exchange coefficient between vegetation canopy air and
356 ! ! atmosphere modified by a partially to fully buried
357 ! ! vegetation canopy [kg/(m2 s)]
358 REAL, DIMENSION(SIZE(PPS)) :: ZEVAP_C_A ! Total canopy evapotranspiration and sublimation
359 ! ! of intercepted snow (kg/m2/s)
360 REAL, DIMENSION(SIZE(PPS)) :: ZEVAP_N_A ! Vapor flux from the ground-based snowpack (part burying
361 ! ! the canopy vegetation) to the atmosphere (kg/m2/s)
362 REAL, DIMENSION(SIZE(PPS)) :: ZH_N_A ! Sensible heat flux from the ground-based snowpack (part
363 ! ! burying the canopy vegetation) to the atmosphere (W/m2)
364 REAL, DIMENSION(SIZE(PPS)) :: ZVEGFACT ! Fraction of canopy vegetation possibly receiving
365 ! ! rainfall (-)
366 REAL, DIMENSION(SIZE(PPS)) :: ZRRSFC ! The sum of all non-intercepted rain and canopy drip (kg/m2/s)
367 REAL, DIMENSION(SIZE(PPS)) :: ZRRSFCL ! The sum of all non-intercepted rain and drip from litter (kg/m2/s)
368 REAL, DIMENSION(SIZE(PPS)) :: ZLES3L ! latent heat flux - sublimation of ice from the ground
369 ! ! based snowpack (W/m2)
370 REAL, DIMENSION(SIZE(PPS)) :: ZLEL3L ! latent heat flux - evaporation of liquid water from the
371 ! ! ground based snowpack (W/m2))
372 REAL, DIMENSION(SIZE(PPS)) :: ZEVAP3L ! total mass loss via evap & sublm from the ground based snowpack (kg/m2/s)
373 REAL, DIMENSION(SIZE(PPS)) :: ZUSTAR2_IC ! friction velocity (possibly implicitly coupled) (m/s)
374 REAL, DIMENSION(SIZE(PPS)) :: ZTA_IC ! atmospheric temperature (possibly implicitly coupled) (m/s)
375 REAL, DIMENSION(SIZE(PPS)) :: ZQA_IC ! atmospheric specific humidity (possibly implicitly coupled) (m/s)
376 REAL, DIMENSION(SIZE(PPS)) :: ZSWUP ! net upwelling shortwave radiation [W/m2]
377 REAL, DIMENSION(SIZE(PPS)) :: ZLWUP ! net upwelling longwave radiation [W/m2]
378 REAL, DIMENSION(SIZE(PPS)) :: ZUSTAR2SNOW ! snow fraciton velocity squared (m2/s2)
379 REAL, DIMENSION(SIZE(PPS)) :: ZVMOD ! lowest level atmospheric wind speed update estimate (K)
380 REAL, DIMENSION(SIZE(PPS)) :: ZRR ! combined rain rate (above canopy) and irrigation need (kg/m2/s)
381 REAL, DIMENSION(SIZE(PPS)) :: ZFLSN_COR ! snow/soil-biomass correction flux (W/m2) (not MEB)
382 REAL, DIMENSION(SIZE(PPS)) :: ZWSFC ! surface liquid water content for resistances (m3/m3)
383 REAL, DIMENSION(SIZE(PPS)) :: ZWISFC ! surface frozen water content for resistances (m3/m3)
384 REAL, DIMENSION(SIZE(PPS)) :: ZLESFC ! evaporation from the surface (soil or litter) (W/m2)
385 REAL, DIMENSION(SIZE(PPS)) :: ZLESFCI ! sublimation from the surface (soil or litter) (W/m2)
386 REAL, DIMENSION(SIZE(PPS)) :: ZPERMSNOWFRAC ! fraction of permanent snow/ice
387 !
388 ! - TR_ML radiation option: NOTE...always used by MEB
389 !
390 REAL, DIMENSION(SIZE(PPS),SIZE(PABC)) :: ZIACAN_SUNLIT ! Absorbed PAR of each level within the
391 REAL, DIMENSION(SIZE(PPS),SIZE(PABC)) :: ZIACAN_SHADE ! canopy - Split into SHADEd and SUNLIT
392 REAL, DIMENSION(SIZE(PPS),SIZE(PABC)) :: ZFRAC_SUN ! fraction of sunlit leaves
393 !
394 REAL, DIMENSION(SIZE(PPS)) :: ZLAI ! Potentially covered/buried canopy LAI (m2/m2)
395 REAL, DIMENSION(SIZE(PPS)) :: ZALBVIS_TSOIL ! average snow-free ground VIS albedo (soil plus flooded fraction)
396 REAL, DIMENSION(SIZE(PPS)) :: ZALBNIR_TSOIL ! average snow-free ground NIR albedo (soil plus flooded fraction)
397 REAL, DIMENSION(SIZE(PPS)) :: ZSWNET_S ! Net SW radiation at the surface (below canopy snow/ground/flooded zone)
398 REAL, DIMENSION(SIZE(PPS)) :: ZLTT ! Average latent heat (normalization factor) (J/kg)
399 REAL, DIMENSION(SIZE(PPS)) :: ZLSTTC ! Working coefficient to compute ZLTT: frozen part (-)
400 REAL, DIMENSION(SIZE(PPS)) :: ZLVTTC ! Working coefficient to compute ZLTT: non-frozen part (-)
401 REAL, DIMENSION(SIZE(PPS)) :: ZZREF
402 REAL, DIMENSION(SIZE(PPS)) :: ZUREF
403 !
404 !
405 ! - CPHOTO/=NON (Ags Option(s)):
406 !
407 REAL, DIMENSION(SIZE(PPS)) :: ZQSAT ! CPHOTO/=NON (Ags Option(s))diagnosed (past time step) Qsat relative to canopy (for Ags)
408 REAL, DIMENSION(SIZE(PPS)) :: ZFFV ! submerged vegetation (by flooding) fraction (-)
409 REAL, DIMENSION(SIZE(PPS),SIZE(PABC)) :: ZIACAN ! PAR in the canopy at different gauss levels: local working needed if
410 ! ! Ags if off (i.e. CPHOTO==NON)
411 !
412 REAL, DIMENSION(:,:), ALLOCATABLE :: ZTGL ! Temporary temperature of litter + soil
413 REAL, DIMENSION(:,:), ALLOCATABLE :: ZSOILHCAPZ ! Temporary heat capacity of litter + soil
414 REAL, DIMENSION(:,:), ALLOCATABLE :: ZSOILCONDZ ! Temporary heat conductivity of litter + soil
415 REAL, DIMENSION(:,:), ALLOCATABLE :: ZD_G ! Temporary depth of bottom litter + soil layers
416 REAL, DIMENSION(:,:), ALLOCATABLE :: ZDZG ! Temporary thickness of litter + soil layers
417 REAL, DIMENSION(:,:), ALLOCATABLE :: ZWFC ! Temporary Wfc of bottom litter + soil layers
418 REAL, DIMENSION(:,:), ALLOCATABLE :: ZWSAT ! Temporary Wsat of bottom litter + soil layers
419 !
420 ! Working sums for flux averaging over MEB time split
421 !
422 REAL, DIMENSION(SIZE(PPS)) :: ZH_SUM, ZH_CA_SUM, ZH_N_A_SUM, ZH_CV_SUM, ZH_GV_SUM, &
423  ZH_GN_SUM, ZHSNOW_SUM, ZHPSNOW_SUM
424 REAL, DIMENSION(SIZE(PPS)) :: ZHU_AGG_SUM, ZAC_AGG_SUM
425 
426 REAL, DIMENSION(SIZE(PPS)) :: ZLE_SUM, ZLE_CA_SUM, ZLE_CV_SUM, ZLE_GV_SUM, &
427  ZLE_GN_SUM, ZLETR_CV_SUM, ZLEG_SUM,ZLEGI_SUM,ZLESFC_SUM, &
428  ZLESFCI_SUM, ZLER_CV_SUM, ZLE_FLOOD_SUM, ZLEI_FLOOD_SUM, &
429  ZLES_CV_SUM, ZLETR_SUM, ZLER_SUM, ZLEV_SUM, &
430  ZLEI_SUM, ZLES3L_SUM, ZLEL3L_SUM, ZEVAP3L_SUM, &
431  ZUSTAR2_SUM, ZUSTAR2SNOW_SUM, ZCDSNOW_SUM, &
432  ZCHSNOW_SUM, ZRISNOW_SUM, ZEVAP_SUM
433 
434 REAL, DIMENSION(SIZE(PPS)) :: ZGRNDFLUX_SUM, ZRESTORE_SUM
435 
436 REAL, DIMENSION(SIZE(PPS)) :: ZSWNET_V_SUM, ZSWNET_G_SUM, ZSWNET_N_SUM, ZLWNET_V_SUM, &
437  ZLWNET_G_SUM, ZLWNET_N_SUM, ZEMIST_SUM, ZSWUP_SUM, &
438  ZLWUP_SUM
439 REAL, DIMENSION(SIZE(PPS)) :: ZDELHEATG_SFC_SUM, ZDELHEATV_SFC_SUM, ZDELHEATG_SUM
440 !
441 REAL, DIMENSION(SIZE(PPS)) :: ZALBV, ZALBG, ZTR
442 !
443 REAL, DIMENSION(SIZE(PEK%XWR,1)) :: ZPHASEL ! Phase changement in litter (W/m2)
444 REAL, DIMENSION(SIZE(PEK%XWR,1)) :: ZCTSFC
445 REAL, DIMENSION(SIZE(PEK%XWR,1)) :: ZFROZEN1SFC
446 !
447 INTEGER :: INJ, INL, JJ, JL
448 !
449 REAL(KIND=JPRB) :: ZHOOK_HANDLE
450 !-------------------------------------------------------------------------------
451 !
452 !* 1.0 Preliminaries
453 ! -------------
454 !
455 IF (lhook) CALL dr_hook('ISBA_MEB',0,zhook_handle)
456 !
457 piacan(:,:) = 0.0
458 !
459 dmk%XFAPAR(:) = 0.0
460 dmk%XFAPIR(:) = 0.0
461 dmk%XFAPAR_BS(:) = 0.0
462 dmk%XFAPIR_BS(:) = 0.0
463 !
464 dek%XRRLIT(:) = 0.0
465 dek%XDRIPLIT(:) = 0.0
466 !
467 dek%XLEGI(:) = 0.0
468 dek%XLEG(:) = 0.0
469 !
470 zlesfci(:) = 0.0
471 zlesfc(:) = 0.0
472 !
473 ziacan_sunlit(:,:) = xundef
474 ziacan_shade(:,:) = xundef
475 zfrac_sun(:,:) = xundef
476 zlai(:) = xundef
477 zalbvis_tsoil(:) = xundef
478 zalbnir_tsoil(:) = xundef
479 zswnet_s(:) = xundef
480 zqsat(:) = xundef
481 !
482 zwork(:) = xundef
483 zwork2(:) = xundef
484 zwork3(:) = xundef
485 zwork4(:) = xundef
486 !
487 ztr(:) = 0.0
488 !
489 !* 1.1 Preliminaries for litter parameters
490 ! -----------------------------------
491 !
492 inj=SIZE(pek%XWG,1)
493 inl=SIZE(pek%XWG,2)
494 !
496 !
497 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
498 !
499 !* 1.2 Preliminaries for litter temperature
500 ! ------------------------------------
501 !
502 ! Concatenate PEK%XTL and PEK%XTG and the parameters linked to heat transfer into the soil
503 !
504 
505  CALL prep_meb_soil(io%LMEB_LITTER, psoilhcapz, psoilcondz, kk, pk, pek, &
506  zd_g, zdzg,ztgl, zsoilhcapz, zsoilcondz, zwsat, zwfc,&
507  zwsfc, zwisfc, zctsfc, dmk%XCT, pfrozen1, zfrozen1sfc )
508 !
509 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
510 !
511 !* 2.0 Preliminaries for energy and radiation budget
512 ! ---------------------------------------------
513 !
514 zpermsnowfrac(:) = pk%XVEGTYPE_PATCH(:,nvt_snow)
515 !
516 ! Local working:
517 ! - possibly adjust these prognostic variables locally, but do not save
518 !
519 zsnowrho(:,:) = pek%TSNOW%RHO (:,:)
520 zsnowage(:,:) = pek%TSNOW%AGE (:,:)
521 zsnowswe(:,:) = pek%TSNOW%WSNOW(:,:)
522 !
523  CALL preps_for_meb_ebud_rad(pps, pek%XLAI, zsnowrho, zsnowswe, pek%TSNOW%HEAT, zsnowliq, &
524  dmk%XSNOWTEMP, dmk%XSNOWDZ, zsnowcond, zsnowhcap, &
525  pek%TSNOW%EMIS, zsigma_f, zchip, ptstep, psr, pta, &
526  pvmod, zsnowage, zpermsnowfrac )
527 !
528 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
529 !
530 !* 3.0 Shortwave radiative transfer
531 ! ----------------------------
532 !
533 ! Calculate snow albedo: split into spectral bands:
534 !
535  CALL snowalb_spectral_bands_meb(pk%XVEGTYPE_PATCH, pek, zsnowrho, zsnowage, pps,&
536  dmk%XSNOWDZ,pzenith, ztau_n)
537 !
538 !
539 ! NOTE, currently MEB only uses 2 of 3 potential snow albedo spectral bands
540 !
541 !
542 WHERE(pek%TSNOW%ALB(:) /= xundef)
543  zlai(:) = pek%XLAI(:)*(1.0-ppalphan(:))
544  zalbvis_tsoil(:) = palbvis_tsoil(:)*(1.-pek%XPSN(:)) + pek%XPSN(:)*pek%TSNOW%ALBVIS(:)
545  zalbnir_tsoil(:) = palbnir_tsoil(:)*(1.-pek%XPSN(:)) + pek%XPSN(:)*pek%TSNOW%ALBNIR(:)
546 ELSEWHERE
547  zlai(:) = pek%XLAI(:)
548  zalbvis_tsoil(:) = palbvis_tsoil(:)
549  zalbnir_tsoil(:) = palbnir_tsoil(:)
550 END WHERE
551 !
552  CALL radiative_transfert(io%LAGRI_TO_GRASS, pk%XVEGTYPE_PATCH, &
553  palbvis_tveg, zalbvis_tsoil, palbnir_tveg, zalbnir_tsoil, &
554  psw_rad, zlai, pzenith, pabc, pek%XFAPARC, pek%XFAPIRC, &
555  pek%XMUS, pek%XLAI_EFFC, oshade, ziacan, ziacan_sunlit, &
556  ziacan_shade, zfrac_sun, dmk%XFAPAR, dmk%XFAPIR, &
557  dmk%XFAPAR_BS, dmk%XFAPIR_BS )
558 
559 ! Compute all-wavelength effective ground (soil+snow) surface,
560 ! soil and veg albedos, respectively:
561 
562 zalbs(:) = xsw_wght_nir*zalbnir_tsoil(:) + xsw_wght_vis*zalbvis_tsoil(:)
563 
564 zalbg(:) = xsw_wght_nir*palbnir_tsoil(:) + xsw_wght_vis*palbvis_tsoil(:)
565 
566 zalbv(:) = xsw_wght_nir*palbnir_tveg(:) + xsw_wght_vis*palbvis_tveg(:)
567 !
568 WHERE(psw_rad(:) > zswrad_min) ! Sun is up...approx
569 !
570 ! Total effective surface (canopy, ground/flooded zone, snow) all-wavelength
571 ! albedo: diagnosed from shortwave energy budget closure
572 
573  dk%XALBT(:) = 1. - (xsw_wght_vis*(dmk%XFAPAR(:)+dmk%XFAPAR_BS(:)) + &
574  xsw_wght_nir*(dmk%XFAPIR(:)+dmk%XFAPIR_BS(:)))
575  zswup(:) = psw_rad(:)*dk%XALBT(:)
576  dk%XALBT(:) = zswup(:)/max(1.e-5, psw_rad(:))
577 
578 ! Diagnose all-wavelength SW radiative budget
579 ! for the canopy and below canopy (surface) components;
580 
581  dek%XSWNET_V(:) = psw_rad(:)*(xsw_wght_vis*dmk%XFAPAR (:) + xsw_wght_nir*dmk%XFAPIR (:))
582  zswnet_s(:) = psw_rad(:)*(xsw_wght_vis*dmk%XFAPAR_BS(:) + xsw_wght_nir*dmk%XFAPIR_BS(:))
583 !
584 ! Compute total all wavelength SW transmission:
585 ! A solution of a quadradic equation based on ground energy budget Eq in FAPAIR.F90
586 ! Tr = [ -b - sqrt(b^2 - 4ac) ]/(2a) (this is good root for this computation)
587 ! Here we derrive Eq so that a=1
588 
589  zwork4(:) = zalbs(:)*zalbv(:)
590  zwork2(:) = -(1.-zalbs(:)*(1.-zalbv(:)))/zwork4(:) ! b
591  zwork3(:) = zswnet_s(:)/(psw_rad(:)*zwork4(:)) ! c
592  zwork(:) = sqrt(max(0.0, zwork2(:)**2 - 4*zwork3(:))) ! sqrt(b**2 - 4c)
593  ztr(:) = 0.5*(-zwork2(:) - zwork(:)) ! -b - sqrt(b^2 - 4c) ]/2
594  ztr(:) = min(1.,max(0., ztr(:) ))
595 !
596 ! Downwelling SW radiation arriving at ground/snow surface
597 
598  dek%XSWDOWN_GN(:) = psw_rad(:)*ztr(:)
599 
600 ! Get snow and ground components:
601 
602  dek%XSWNET_G(:) = (1.-pek%XPSN(:))*dek%XSWDOWN_GN(:)*(1.-zalbg(:)+zalbs(:)*(1.-ztr(:))*zalbv(:))
603  dek%XSWNET_N(:) = zswnet_s(:)-dek%XSWNET_G(:) ! conservative computation
604 
605 ! Quantity of net shortwave radiation absorbed in surface snow layer
606 
607  dek%XSWNET_NS(:) = dek%XSWNET_N(:)*(1.0 - ztau_n(:,1))
608 !
609 ! Any SW radiation reaching the base of the lowest snow layer can pass
610 ! into the soil:
611 
612  ztau_n(:,SIZE(pek%TSNOW%WSNOW,2)) = ztau_n(:,SIZE(pek%TSNOW%WSNOW,2))*(1.-zalbg(:))
613 !
614 ELSEWHERE
615 
616 ! Sun is down:
617 
618  dk%XALBT(:) = 1.0
619  zswup(:) = psw_rad(:)
620  dek%XSWDOWN_GN(:) = 0.
621  dek%XSWNET_G(:) = 0.
622  dek%XSWNET_V(:) = 0.
623  dek%XSWNET_N(:) = 0.
624  dek%XSWNET_NS(:) = 0.
625  ztau_n(:,SIZE(pek%TSNOW%WSNOW,2)) = 0.
626 
627 END WHERE
628 !
629 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
630 !
631 !* 4.0 Longwave radiative transfer
632 ! ---------------------------
633 !
634  CALL isba_lwnet_meb(pek%XLAI, pek%XPSN, ppalphan,pek%TSNOW%EMIS, kk%XEMISF, kk%XFF, &
635  pek%XTV, ztgl(:,1), dmk%XSNOWTEMP(:,1), plw_rad, dek%XLWNET_N, &
636  dek%XLWNET_V, dek%XLWNET_G, zdlwnet_v_dtv, zdlwnet_v_dtg, zdlwnet_v_dtn, &
637  zdlwnet_g_dtv, zdlwnet_g_dtg, zdlwnet_g_dtn, zdlwnet_n_dtv, &
638  zdlwnet_n_dtg, zdlwnet_n_dtn, zsigma_f, zsigma_fn, dek%XLWDOWN_GN )
639 !
640 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
641 !
642 !* 5.0 Fraction of leaves occupied by intercepted water
643 ! ------------------------------------------------
644 !
645 ! Vegetation canopy:
646 !
647 ! First, compute an effective veg fraction: it can only be < unity if vegetation is buried by snowpack...
648 !
649 zwork(:) = (1.0 - pek%XPSN(:) + pek%XPSN(:)*(1.0 - ppalphan(:)))
650 !
651  CALL wet_leaves_frac(pek%XWR(:), zwork, pek%XWRMAX_CF(:), pz0_mebv, pek%XLAI(:), zwrmax, zdelta)
652 !
653 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
654 !
655 !* 6.0 Plant stress, stomatal resistance and, possibly, CO2 assimilatio
656 ! --------------------------------------------------------------------
657 !
658 ! MEB-NOTE here assumed IO%CPHOTO=='DEF' or 'AST' for now
659 ! More Ags options to be added later
660 !
661 IF (io%CPHOTO=='NON') THEN
662 !
663 ! Canopy vegetation (no snow, or snow below the main part of the canopy):
664 !
665  CALL veg(psw_rad, pek%XTC(:), pek%XQC(:), pps, pek%XRGL(:), pek%XLAI(:), &
666  pek%XRSMIN(:), pek%XGAMMA(:), pf2, zrs)
667 !
668 !
669 ELSE IF (maxval(pek%XGMES) /= xundef .OR. minval(pek%XGMES) /= xundef) THEN
670 !
671 ! NOTE: For now we assume that forest canopy can be flooded.
672 ! However, we need to likely compute a fraction like PALPHAN (for snow vertical extent)
673 ! for floods for grasses/crops/shrubs...i.e. low vegetation
674 
675  zffv(:) = 0.0
676 
677  zqsat(:) = qsat(pek%XTV(:),pps)
678 
679  zwork(:) = pek%XLE(:) ! passed in as LE: Since Qc corresponds to the effective
680  pek%XLE(:) = 0. ! surface specific humidity in ISBA-MEB,
681  ! so no need for LE correction (only required for composite ISBA)
682  CALL cotwores(ptstep, io, oshade, pk, pek, pk%XDMAX, ppoi, pcsp, pek%XTV, &
683  pf2, psw_rad, pek%XQC, zqsat, ppalphan, zdelta, prhoa, &
684  pzenith, zffv, ziacan_sunlit, ziacan_shade, zfrac_sun, &
685  ziacan, pabc, zrs, dek%XGPP, presp_biomass_inst(:,1))
686  pek%XLE(:) = zwork(:)
687 !
688  piacan(:,:) = ziacan(:,:)
689 !
690 ELSE
691  presp_biomass_inst(:,1) = 0.0
692  dek%XGPP(:) = 0.0
693 ENDIF
694 !
695 ! Additional resistance for possibly snow-buried canopy vegetation:
696 !
697 zrsn(:) = zrs(:)/( 1.0 - min(ppalphan(:), 1.0 - (zrs(:)/xrs_max)) )
698 !
699 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
700 !
701 !* 6.0 Canopy snow (intercepted) needed diagnostics:
702 ! ---------------------------------------------
703 !
704  CALL snow_leaves_frac_meb(pek%XPSN, ppalphan, pek%XWRVN, pek%XTV, zchip, &
705  pek%XLAI, zwrvnmax, zpsncv, zmeltvn)
706 !
707 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
708 !
709 !* 7.0 Aerodynamic drag and heat/mass transfer/fluxes
710 ! and energy budget solution
711 ! ----------------------------------------------
712 !
713 ! NOTE, this assumes thermodynamic variable herein is potential T
714 
715 zpet_a_coef(:) = -ppet_a_coef(:)*xcpd
716 zpet_b_coef(:) = ppet_b_coef(:)*xcpd
717 zthrma_ta(:) = xcpd/pexna(:)
718 zthrmb_ta(:) = 0.0
719 zwork(:) = xcpd/pexns(:)
720 zthrma_tc(:) = zwork(:)
721 zthrmb_tc(:) = 0.0
722 zthrma_tn(:) = zwork(:)
723 zthrmb_tn(:) = 0.0
724 zthrma_tg(:) = zwork(:)
725 zthrmb_tg(:) = 0.0
726 zthrma_tv(:) = zwork(:)
727 zthrmb_tv(:) = 0.0
728 !
729 !
730 ! For turbulence computations:
731 ! Adjust (shift upward) local reference heights "seen by the turbulence scheme"
732 ! if they are below the canopy:
733 ! NOTE, this approach is an alternative to LFORC_MEASURE=F, which shifts
734 ! vegetation downward by the displacement height. Both approaches essentially
735 ! conceptually assume that the vegetation is part of the terrain.
736 ! Also, here, conserve any UREF and ZREF differences.
737 !
738 zzref(:) = pzref(:)
739 zuref(:) = puref(:)
740 IF(io%LFORC_MEASURE)THEN
741  WHERE(pzref(:) - pek%XH_VEG(:) < xlimh)
742  zzref(:) = pek%XH_VEG(:) + xlimh
743  zuref(:) = pek%XH_VEG(:) + xlimh + max(0.,puref(:)-pzref(:))
744  END WHERE
745 ENDIF
746 !
747 ! Compute the average latent heat (normalization factor) (J kg-1):
748 ! NOTE that we could use a function which depends on the different resistances,
749 ! but this can make the average latent heat relatively noisy(leading to a slightly less
750 ! tightly closed budgets owing to numerical noise): here we opt for a more
751 ! temporally smooth approximation based on fractional coverage:
752 !
753 zlvttc(:) = ( zsigma_f(:)*(1.-zpsncv(:)) + (1.0-pek%XPSN(:)-kk%XFF(:))*(1.0-zfrozen1sfc(:)) )* &
754  (1.0 - pek%XPSN(:)*ppalphan(:))
755 zlsttc(:) = ( zsigma_f(:)* zpsncv(:) + (1.0-pek%XPSN(:)-kk%XFF(:))* zfrozen1sfc(:) + pek%XPSN(:) )* &
756  (1.0 - pek%XPSN(:)*ppalphan(:)) + pek%XPSN(:)*ppalphan(:)
757 zltt(:) = (zlvttc(:)*xlvtt + zlsttc(:)*xlstt)/max(1.e-12, zlvttc(:) + zlsttc(:))
758 zltt(:) = min(xlstt, max(zltt(:), xlvtt)) ! numerical check
759 !
760 ! Possibly split time step if large:
761 ! Although the energy budget is fully implicit, a very small canopy heat capacity
762 ! (and neglect of canopy air space heat capacity) can possibly lead to
763 ! numerical shocks, especially during transition periods between stable and unstable
764 ! regimes. Thus, for relatively large steps, a simple time split scheme is activated.
765 ! Note that soil moisture is held constant, while turbulent exchange coefficients are updated during the split.
766 ! Also, experience shows that splitting at least once for moderately sized time steps
767 ! is quite effective in removing any lingering small but possible oscillations.
768 ! Finally, for *very* small time steps (such as those for high res runs), no split is performed.
769 ! Fluxes are averaged over the time split for conservation.
770 !
771 jtsplit_eb = 1 + int(ptstep/ztstep_eb) ! number of split-time steps
772 ztstep = ptstep/jtsplit_eb ! split time step...for relatively small time steps, no split
773 !
774 ! initialize time split sums for fluxes:
775 !
777 !
778 !
779 ! Note, when implicitly coupled to the atmosphere, these
780 ! 3 variables will evolve during the split...we used updated
781 ! values for turbulent exchange computations (drag_meb).
782 ! NOTE...when explicit coupling used, these 3 variables do NOT vary
783 ! during the split.
784 !
785 zvmod(:) = pvmod(:)
786 zta_ic(:) = pta(:)
787 zqa_ic(:) = pqa(:)
788 !
789 !
790 loop_time_split_eb: DO jdt=1,jtsplit_eb
791 !* 7.1 Aerodynamic drag and heat transfer coefficients
792 ! -----------------------------------------------
793 !
794  CALL drag_meb(io, pek, dmk, dk, &
795  ztgl(:,1), zta_ic, zqa_ic, zvmod, zwsfc, zwisfc, &
796  zwsat(:,1), zwfc(:,1), pexns, pexna, pps, prr, psr, &
797  prhoa, pz0g_without_snow, pz0_mebv, pz0h_mebv, &
798  pz0eff_mebv, pz0_mebn, pz0h_mebn, pz0eff_mebn, &
799  zsnowswe(:,1), zchip, ztstep, zrs, zrsn, ppalphan, &
800  zzref, zuref, pdircoszw, zpsncv, zdelta, zvelc, &
801  prisnow, zustar2snow, zhugi, zhvg, &
802  zhvn, zleg_delta, zlegi_delta, zhsgl, zhsgf, &
803  zflxc_ca, zflxc_n_a, zflxc_gv, zflxc_gn, &
804  zflxc_vg_c, zflxc_vn_c, zflxc_mom, zqsatg, zqsatv, &
805  zqsatc, zqsatn, zdeltavk )
806 !
807  zkvn(:) = snow_intercept_eff(zchip,zvelc,zwrvnmax)
808 
809 !* 7.2 Resolution of the surface energy budgets
810 ! ----------------------------------------
811 !
812  CALL e_budget_meb(io, kk, pk, pek, dk, dek, dmk, &
813  ztstep, zltt, pps, zctsfc, ptdeep_a, zd_g, zsoilcondz, zsoilhcapz,&
814  zsnowrho, zsnowcond, zsnowhcap, ztau_n, zdlwnet_v_dtv, zdlwnet_v_dtg, &
815  zdlwnet_v_dtn, zdlwnet_g_dtv, zdlwnet_g_dtg, zdlwnet_g_dtn, &
816  zdlwnet_n_dtv, zdlwnet_n_dtg, zdlwnet_n_dtn, ppew_a_coef, &
817  ppew_b_coef, zpet_a_coef, ppeq_a_coef, zpet_b_coef, &
818  ppeq_b_coef, zthrma_ta, zthrmb_ta, zthrma_tc, zthrmb_tc, &
819  zthrma_tg, zthrmb_tg, zthrma_tv, zthrmb_tv, zthrma_tn, &
820  zthrmb_tn, zqsatg, zqsatv, zqsatn, ppalphan, zpsncv, &
821  zcheatv, zcheatg, zcheatn, zleg_delta, zlegi_delta, zhugi, &
822  zhvg, zhvn, zfrozen1sfc, zflxc_ca, zflxc_gv, zflxc_vg_c, &
823  zflxc_vn_c, zflxc_gn, zflxc_n_a, zflxc_mom, ztgl, zsnowliq, &
824  zflxc_cv, zhvgs, zhvns, zdqsat_g,zdqsat_v,zdqsati_n, &
825  zta_ic, zqa_ic, zustar2_ic, zvmod, zdeltat_g, zdeltat_v, &
826  zdeltat_n, pgrndflux, pdeep_flux, pdelheatv_sfc, &
827  pdelheatg_sfc, pdelheatg )
828 !
829 !* 7.3 Energy and momentum fluxes and radiative temperature and emissivity
830 ! -------------------------------------------------------------------
831 !
832  CALL isba_fluxes_meb(kk, pk, pek, dk, dek, dmk, prhoa, zltt, zsigma_f, zsigma_fn, &
833  zrnet_v, zrnet_g, zdlwnet_v_dtv, zdlwnet_v_dtg, &
834  zdlwnet_v_dtn, zdlwnet_g_dtv, zdlwnet_g_dtg, &
835  zdlwnet_g_dtn, zdlwnet_n_dtv, zdlwnet_n_dtg, &
836  zdlwnet_n_dtn, zthrma_ta, zthrmb_ta, zthrma_tc, &
837  zthrmb_tc, zthrma_tg, zthrmb_tg, zthrma_tv, zthrmb_tv, &
838  zthrma_tn, zthrmb_tn, zqsatg, zqsatv, zqsatn, ppalphan, &
839  zpsncv, zfrozen1sfc, zleg_delta, zlegi_delta, zhugi, &
840  zhvg, zhvn, zflxc_ca, zflxc_gv, zflxc_vg_c, zflxc_vn_c, &
841  zflxc_gn, zflxc_n_a, zflxc_mom, zflxc_cv, zhvgs, &
842  zhvns, ztgl, zdqsat_g, zdqsat_v, zdqsati_n, zta_ic, &
843  zqa_ic, zdeltavk, zdeltat_g, zdeltat_v, zdeltat_n, &
844  zswup, psw_rad, plw_rad, zlwup, zh_n_a, zevap_c_a, &
845  zevap_n_a, zlesfc, zlesfci, zles3l, zlel3l, zevap3l, &
846  pemist )
847 !
848 ! Compute aggregated coefficients for evaporation
849 ! Sum(LEC+LES+LEL) = ACagg * Lv * RHOA * (HUagg.Qsat - Qa)
850 !
851  zflxc_c_a_f(:) = zflxc_ca(:)*(1.0-pek%XPSN(:)*ppalphan(:))
852  zflxc_n_a_f(:) = zflxc_n_a(:)* pek%XPSN(:)*ppalphan(:)
853 
854  phu_agg(:) = (zflxc_c_a_f(:)*pek%XQC(:)+ zflxc_n_a_f(:)*zqsatn(:))/ &
855  (zflxc_c_a_f(:)*zqsatc(:) + zflxc_n_a_f(:)*zqsatn(:))
856 
857  pac_agg(:) = zflxc_c_a_f(:) + zflxc_n_a_f(:) ! kg/m2/s
858 !
859 ! Sum fluxes over time split:
860 
862 
863 ENDDO loop_time_split_eb
864 !
865  CALL avg_fluxes_meb_tsplit ! average fluxes over time split
866 !
867 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
868 !
869 !* 8.0 Snow explicit canopy loading/interception
870 ! ------------------------------------------
871 !
872  CALL snow_load_meb(pk, pek, dek, ptstep, psr, zwrvnmax, zkvn, zcheatv, zmeltvn, &
873  zvelc, psubvcor)
874 !
875 
876 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
877 !
878 !* 9.0 Snow explicit canopy loading/interception
879 ! ------------------------------------------
880 !
881 zrr(:) = prr(:)
882 dek%XIRRIG_FLUX(:) = 0.0
883 !
884 !* add irrigation over vegetation to liquid precipitation (rr)
885 ! Water "need" treated as if sprayed from above (over vegetation and soil):
886 !
887 IF (SIZE(ag%LIRRIGATE,1)>0) THEN
888  WHERE (ag%LIRRIGATE(:) .AND. pek%XIRRIG(:)>0. .AND. pek%XIRRIG(:) /= xundef .AND. (pf2(:)<ag%XTHRESHOLDSPT(:)) )
889  dek%XIRRIG_FLUX(:) = pek%XWATSUP(:) / xday
890  zrr(:) = prr(:) + pek%XWATSUP(:)/xday
891  ag%LIRRIDAY (:) = .true.
892  END WHERE
893 ENDIF
894 !
895 ! Call canopy interception...here because meltwater should be allowed to fall
896 ! on understory snowpack
897 !
898 ! Fraction of canopy vegetation possibly receiving rainfall/irrigation
899 !
900 zvegfact(:) = zsigma_f(:)*(1.0-ppalphan(:)*pek%XPSN(:))
901 !
902 ! The sum of all non-intercepted rain and drip is "ZRRSFC" (kg/m2/s):
903 ! this is then partitioned by snow scheme into part falling on
904 ! snowpack and part falling onto snow-free understory.
905 !
906 !
907  CALL hydro_veg(io%CRAIN, ptstep, kk%XMUF, zrr, dek%XLEV_CV, dek%XLETR_CV, &
908  zvegfact, zpsncv, pek%XWR, zwrmax, zrrsfc, dek%XDRIP, dek%XRRVEG, &
909  pk%XLVTT )
910 !
911 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
912 !
913 !* 10.0 Explicit snow scheme (MEB: impose surface fluxes as upper BC)
914 ! ----------------------------------------------------------------
915 !
916  CALL snow3l_isba(io, g, pk, pek, dk, dek, dmk, omeb, himplicit_wind, &
917  tptime, ptstep, pk%XVEGTYPE_PATCH, ztgl, zctsfc, &
918  zsoilhcapz, zsoilcondz(:,1), pps, pta, psw_rad, pqa, &
919  pvmod, plw_rad, zrrsfc, dek%XSR_GN, prhoa, zuref, pexns, &
920  pexna, pdircoszw, zzref, zalbg, zd_g, zdzg, ppew_a_coef, &
921  ppew_b_coef, ppet_a_coef, ppeq_a_coef, ppet_b_coef, &
922  ppeq_b_coef, psnow_thrufal, pgrndflux, pflsn_cor, &
923  prestoren, pevapcor, dek%XLES, dek%XLESL, zevap3l, psnowsfch, &
924  pdelheatn, pdelheatn_sfc, prisnow, pzenith, pdelheatg, &
925  pdelheatg_sfc, pqsnow )
926 !
927 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
928 !* 11.0 Litter layer hydrology:
929 ! -----------------------
930 !
931 IF(io%LMEB_LITTER)THEN
932 !
933  zwork(:) = 0.
934  zwork2(:) = pek%XWRL(:)
935  zwork3(:) = 1.
936  zwork4(:) = psnow_thrufal(:) + zrrsfc(:)*(1-pek%XPSN)
937  zwrlmax(:) = pek%XGNDLITTER(:)*zwfc(:,1)*xrholw
938 
939  CALL hydro_veg(io%CRAIN, ptstep, kk%XMUF, zwork4(:), zlesfc, zwork, zwork3, zwork, &
940  pek%XWRL, zwrlmax, zrrsfcl, dek%XDRIPLIT, dek%XRRLIT, pk%XLVTT )
941 
942  dmk%XRRSFC(:) = zrrsfcl(:)
943  psnow_thrufal_soil(:) = 0.0
944 
945 ELSE
946 
947  psnow_thrufal_soil(:) = psnow_thrufal(:)
948 
949 ENDIF
950 !
951 !* 11.0 Separate litter and soil temperature
952 ! ------------------------------------
953 !
954  CALL reshift_meb_soil(io%LMEB_LITTER, ztgl, zlesfc, zlesfci, pek, dek)
955 !
957 !
958 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
959 !* 13.0 Ice litter effect
960 ! ------------------
961 !
962 IF(io%LMEB_LITTER)THEN
963 !
964  CALL ice_litter(ptstep, dek%XLELITTERI, psoilhcapz, pek, pk%NWG_LAYER, &
965  pk%XDZG, zphasel,zctsfc,pk%XLSTT,plitcor )
966 !
967 ENDIF
968 !
969 IF (lhook) CALL dr_hook('ISBA_MEB',1,zhook_handle)
970 !
971 !-------------------------------------------------------------------------------
972 !
973 CONTAINS
974 !
975 !===============================================================================
976 !
977 SUBROUTINE init_sum_fluxes_meb_tsplit
978 !
979 IMPLICIT NONE
980 !
981 !
982 !* 0.2 declarations of local variables
983 !
984 REAL(KIND=JPRB) :: ZHOOK_HANDLE
985 !
986 !-------------------------------------------------------------------------------
987 !
988 IF (lhook) CALL dr_hook('ISBA_MEB:INIT_SUM_FLUXES_MEB_TSPLIT ',0,zhook_handle)
989 !
990 ! sensible heat fluxes:
991 !
992 zh_sum(:) = 0.0
993 zh_n_a_sum(:) = 0.0
994 zh_ca_sum(:) = 0.0
995 zh_cv_sum(:) = 0.0
996 zh_gv_sum(:) = 0.0
997 zh_gn_sum(:) = 0.0
998 zhsnow_sum(:) = 0.0
999 !
1000 ! latent heat/water vapor fluxes:
1001 !
1002 zle_sum(:) = 0.0
1003 !
1004 zle_ca_sum(:) = 0.0
1005 zle_cv_sum(:) = 0.0
1006 zle_gv_sum(:) = 0.0
1007 zle_gn_sum(:) = 0.0
1008 !
1009 zletr_cv_sum(:) = 0.0
1010 zler_cv_sum(:) = 0.0
1011 zles_cv_sum(:) = 0.0
1012 !
1013 zleg_sum(:) = 0.0
1014 zlegi_sum(:) = 0.0
1015 zlesfc_sum(:) = 0.0
1016 zlesfci_sum(:) = 0.0
1017 zle_flood_sum(:) = 0.0
1018 zlei_flood_sum(:)= 0.0
1019 zletr_sum(:) = 0.0
1020 zler_sum(:) = 0.0
1021 zlev_sum(:) = 0.0
1022 zlei_sum(:) = 0.0
1023 zles3l_sum(:) = 0.0
1024 zlel3l_sum(:) = 0.0
1025 zevap3l_sum(:) = 0.0
1026 zevap_sum(:) = 0.0
1027 !
1028 zhu_agg_sum(:) = 0.0
1029 zac_agg_sum(:) = 0.0
1030 !
1031 ! momentum/turb:
1032 !
1033 zustar2_sum(:) = 0.0
1034 zustar2snow_sum(:) = 0.
1035 zcdsnow_sum(:) = 0.
1036 zchsnow_sum(:) = 0.
1037 zrisnow_sum(:) = 0.
1038 !
1039 ! surface interfacial/sub-surface heat fluxes:
1040 !
1041 zgrndflux_sum(:) = 0.0
1042 zrestore_sum(:) = 0.0
1043 zhpsnow_sum(:) = 0.0
1044 !
1045 ! radiative fluxes:
1046 !
1047 zswnet_v_sum(:) = 0.0
1048 zswnet_g_sum(:) = 0.0
1049 zswnet_n_sum(:) = 0.0
1050 zlwnet_v_sum(:) = 0.0
1051 zlwnet_g_sum(:) = 0.0
1052 zlwnet_n_sum(:) = 0.0
1053 zemist_sum(:) = 0.0
1054 zswup_sum(:) = 0.0
1055 zlwup_sum(:) = 0.0
1056 !
1057 zdelheatv_sfc_sum(:) = 0.0
1058 zdelheatg_sfc_sum(:) = 0.0
1059 zdelheatg_sum(:) = 0.0
1060 !
1061 IF (lhook) CALL dr_hook('ISBA_MEB:INIT_SUM_FLUXES_MEB_TSPLIT ',1,zhook_handle)
1062 !
1063 END SUBROUTINE init_sum_fluxes_meb_tsplit
1064 !
1065 !===============================================================================
1066 !
1067 SUBROUTINE sum_fluxes_meb_tsplit
1069 IMPLICIT NONE
1070 !
1071 !* 0.2 declarations of local variables
1072 !
1073 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1074 !
1075 !-------------------------------------------------------------------------------
1076 !
1077 IF (lhook) CALL dr_hook('ISBA_MEB:SUM_FLUXES_MEB_TSPLIT ',0,zhook_handle)
1078 !
1079 ! Sum fluxes over MEB TIME SPLIT:
1080 !
1081 ! sensible heat fluxes:
1082 !
1083 zh_n_a_sum(:) = zh_n_a_sum(:) + zh_n_a(:)
1084 !
1085 zh_sum(:) = zh_sum(:) + dk%XH(:)
1086 !
1087 zh_ca_sum(:) = zh_ca_sum(:) + dek%XH_CA(:)
1088 zh_cv_sum(:) = zh_cv_sum(:) + dek%XH_CV(:)
1089 zh_gv_sum(:) = zh_gv_sum(:) + dek%XH_GV(:)
1090 zh_gn_sum(:) = zh_gn_sum(:) + dek%XH_GN(:)
1091 !
1092 zhsnow_sum(:) = zhsnow_sum(:) + dmk%XHSNOW(:)
1093 !
1094 ! latent heat/water vapor fluxes:
1095 !
1096 zle_sum(:) = zle_sum(:) + pek%XLE(:)
1097 !
1098 zlei_sum(:) = zlei_sum(:) + dk%XLEI(:)
1099 zevap_sum(:) = zevap_sum(:) + dk%XEVAP(:)
1100 !
1101 zle_ca_sum(:) = zle_ca_sum(:) + dek%XLE_CA(:)
1102 zle_cv_sum(:) = zle_cv_sum(:) + dek%XLE_CV(:)
1103 zle_gv_sum(:) = zle_gv_sum(:) + dek%XLE_GV(:)
1104 zle_gn_sum(:) = zle_gn_sum(:) + dek%XLE_GN(:)
1105 !
1106 zletr_cv_sum(:) = zletr_cv_sum(:) + dek%XLETR_CV(:)
1107 zler_cv_sum(:) = zler_cv_sum(:) + dek%XLER_CV(:)
1108 zles_cv_sum(:) = zles_cv_sum(:) + dek%XLES_CV(:)
1109 !
1110 zletr_sum(:) = zletr_sum(:) + dek%XLETR(:)
1111 zler_sum(:) = zler_sum(:) + dek%XLER(:)
1112 zlev_sum(:) = zlev_sum(:) + dek%XLEV(:)
1113 !
1114 zleg_sum(:) = zleg_sum(:) + dek%XLEG(:)
1115 zlegi_sum(:) = zlegi_sum(:) + dek%XLEGI(:)
1116 
1117 zle_flood_sum(:) = zle_flood_sum(:) + dek%XLE_FLOOD(:)
1118 zlei_flood_sum(:) = zlei_flood_sum(:) + dek%XLEI_FLOOD(:)
1119 !
1120 zlesfc_sum(:) = zlesfc_sum(:) + zlesfc(:)
1121 zlesfci_sum(:) = zlesfci_sum(:) + zlesfci(:)
1122 !
1123 zles3l_sum(:) = zles3l_sum(:) + zles3l(:)
1124 zlel3l_sum(:) = zlel3l_sum(:) + zlel3l(:)
1125 zevap3l_sum(:) = zevap3l_sum(:) + zevap3l(:)
1126 !
1127 zhu_agg_sum(:) = zhu_agg_sum(:) + phu_agg(:)
1128 zac_agg_sum(:) = zac_agg_sum(:) + pac_agg(:)
1129 !
1130 ! momentum/turb:
1131 !
1132 zcdsnow_sum(:) = zcdsnow_sum(:) + dmk%XCDSNOW(:)
1133 zchsnow_sum(:) = zchsnow_sum(:) + dmk%XCHSNOW(:)
1134 !
1135 zustar2_sum(:) = zustar2_sum(:) + zustar2_ic(:)
1136 zustar2snow_sum(:) = zustar2snow_sum(:) + zustar2snow(:)
1137 zrisnow_sum(:) = zrisnow_sum(:) + prisnow(:)
1138 !
1139 ! surface interfacial/sub-surface heat fluxes:
1140 !
1141 zgrndflux_sum(:) = zgrndflux_sum(:) + pgrndflux(:)
1142 !
1143 zrestore_sum(:) = zrestore_sum(:) + dek%XRESTORE(:)
1144 !
1145 zhpsnow_sum(:) = zhpsnow_sum(:) + dmk%XHPSNOW(:)
1146 !
1147 ! radiative fluxes:
1148 !
1149 zswnet_v_sum(:) = zswnet_v_sum(:) + dek%XSWNET_V(:)
1150 zswnet_g_sum(:) = zswnet_g_sum(:) + dek%XSWNET_G(:)
1151 zswnet_n_sum(:) = zswnet_n_sum(:) + dek%XSWNET_N(:)
1152 zlwnet_v_sum(:) = zlwnet_v_sum(:) + dek%XLWNET_V(:)
1153 zlwnet_g_sum(:) = zlwnet_g_sum(:) + dek%XLWNET_G(:)
1154 zlwnet_n_sum(:) = zlwnet_n_sum(:) + dek%XLWNET_N(:)
1155 !
1156 zemist_sum(:) = zemist_sum(:) + pemist(:)
1157 zswup_sum(:) = zswup_sum(:) + zswup(:)
1158 zlwup_sum(:) = zlwup_sum(:) + zlwup(:)
1159 !
1160 zdelheatv_sfc_sum(:) = zdelheatv_sfc_sum(:) + pdelheatv_sfc(:)
1161 zdelheatg_sfc_sum(:) = zdelheatg_sfc_sum(:) + pdelheatg_sfc(:)
1162 zdelheatg_sum(:) = zdelheatg_sum(:) + pdelheatg(:)
1163 !
1164 IF (lhook) CALL dr_hook('ISBA_MEB:SUM_FLUXES_MEB_TSPLIT ',1,zhook_handle)
1165 !
1166 END SUBROUTINE sum_fluxes_meb_tsplit
1167 !
1168 !===============================================================================
1169 !
1170 SUBROUTINE avg_fluxes_meb_tsplit
1172 USE modd_csts, ONLY : xstefan
1173 !
1174 IMPLICIT NONE
1175 !
1176 !* 0.2 declarations of local variables
1177 !
1178 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1179 !
1180 !-------------------------------------------------------------------------------
1181 !
1182 IF (lhook) CALL dr_hook('ISBA_MEB:AVG_FLUXES_MEB_TSPLIT ',0,zhook_handle)
1183 !
1184 ! Average fluxes over MEB TIME SPLIT:
1185 !
1186 ! sensible heat fluxes:
1187 !
1188 zh_n_a(:) = zh_n_a_sum(:) /jtsplit_eb
1189 !
1190 dk%XH(:) = zh_sum(:) /jtsplit_eb
1191 !
1192 dek%XH_CA(:) = zh_ca_sum(:) /jtsplit_eb
1193 dek%XH_CV(:) = zh_cv_sum(:) /jtsplit_eb
1194 dek%XH_GV(:) = zh_gv_sum(:) /jtsplit_eb
1195 dek%XH_GN(:) = zh_gn_sum(:) /jtsplit_eb
1196 !
1197 dmk%XHSNOW(:) = zhsnow_sum(:) /jtsplit_eb
1198 !
1199 ! latent heat/water vapor fluxes:
1200 !
1201 zlesfc(:) = zlesfc_sum(:) /jtsplit_eb
1202 zlesfci(:) = zlesfci_sum(:) /jtsplit_eb
1203 !
1204 dk%XLEI(:) = zlei_sum(:) /jtsplit_eb
1205 dk%XEVAP(:) = zevap_sum(:) /jtsplit_eb
1206 !
1207 pek%XLE(:) = zle_sum(:) /jtsplit_eb
1208 !
1209 dek%XLE_CA(:) = zle_ca_sum(:) /jtsplit_eb
1210 dek%XLE_CV(:) = zle_cv_sum(:) /jtsplit_eb
1211 dek%XLE_GV(:) = zle_gv_sum(:) /jtsplit_eb
1212 dek%XLE_GN(:) = zle_gn_sum(:) /jtsplit_eb
1213 !
1214 dek%XLETR_CV(:) = zletr_cv_sum(:) /jtsplit_eb
1215 dek%XLER_CV(:) = zler_cv_sum(:) /jtsplit_eb
1216 dek%XLES_CV(:) = zles_cv_sum(:) /jtsplit_eb
1217 !
1218 dek%XLETR(:) = zletr_sum(:) /jtsplit_eb
1219 dek%XLER(:) = zler_sum(:) /jtsplit_eb
1220 dek%XLEV(:) = zlev_sum(:) /jtsplit_eb
1221 !
1222 dek%XLEG(:) = zleg_sum(:) /jtsplit_eb
1223 dek%XLEGI(:) = zlegi_sum(:) /jtsplit_eb
1224 dek%XLE_FLOOD(:) = zle_flood_sum(:) /jtsplit_eb
1225 dek%XLEI_FLOOD(:) = zlei_flood_sum(:)/jtsplit_eb
1226 dek%XLES(:) = zles3l_sum(:) /jtsplit_eb
1227 dek%XLESL(:) = zlel3l_sum(:) /jtsplit_eb
1228 !
1229 zevap3l(:) = zevap3l_sum(:) /jtsplit_eb
1230 !
1231 phu_agg(:) = zhu_agg_sum(:) /jtsplit_eb
1232 pac_agg(:) = zac_agg_sum(:) /jtsplit_eb
1233 !
1234 ! momentum/turb:
1235 !
1236 pustar(:) = sqrt( zustar2_sum(:) /jtsplit_eb )
1237 prisnow(:) = zrisnow_sum(:) /jtsplit_eb
1238 !
1239 dmk%XUSTARSNOW(:) = sqrt( zustar2snow_sum(:)/jtsplit_eb )
1240 dmk%XCDSNOW(:) = zcdsnow_sum(:) /jtsplit_eb
1241 dmk%XCHSNOW(:) = zchsnow_sum(:) /jtsplit_eb
1242 !
1243 ! surface interfacial/sub-surface heat fluxes:
1244 !
1245 pgrndflux(:) = zgrndflux_sum(:) /jtsplit_eb
1246 !
1247 dek%XRESTORE(:) = zrestore_sum(:) /jtsplit_eb
1248 dmk%XHPSNOW(:) = zhpsnow_sum(:) /jtsplit_eb
1249 !
1250 ! radiative fluxes:
1251 !
1252 dek%XSWNET_V(:) = zswnet_v_sum(:) /jtsplit_eb
1253 dek%XSWNET_G(:) = zswnet_g_sum(:) /jtsplit_eb
1254 dek%XSWNET_N(:) = zswnet_n_sum(:) /jtsplit_eb
1255 dek%XLWNET_V(:) = zlwnet_v_sum(:) /jtsplit_eb
1256 dek%XLWNET_G(:) = zlwnet_g_sum(:) /jtsplit_eb
1257 dek%XLWNET_N(:) = zlwnet_n_sum(:) /jtsplit_eb
1258 !
1259 pemist(:) = zemist_sum(:) /jtsplit_eb
1260 zswup(:) = zswup_sum(:) /jtsplit_eb
1261 zlwup(:) = zlwup_sum(:) /jtsplit_eb
1262 !
1263 pdelheatv_sfc(:) = zdelheatv_sfc_sum(:) /jtsplit_eb
1264 pdelheatg_sfc(:) = zdelheatg_sfc_sum(:) /jtsplit_eb
1265 pdelheatg(:) = zdelheatg_sum(:) /jtsplit_eb
1266 !
1267 ! Additional diagnostics depending on AVG quantities:
1268 !
1269 dk%XTSRAD(:) = ((zlwup(:) - plw_rad(:)*(1.0-pemist(:)))/(xstefan*pemist(:)))**0.25
1270 !
1271 zrnet_v(:) = dek%XSWNET_V(:) + dek%XLWNET_V(:)
1272 !
1273 zrnet_g(:) = dek%XSWNET_G(:) + dek%XLWNET_G(:)
1274 !
1275 dmk%XRNSNOW(:) = dek%XSWNET_N(:) + dek%XLWNET_N(:)
1276 !
1277 dk%XRN(:) = zrnet_v(:) + zrnet_g(:) + dmk%XRNSNOW(:)
1278 !
1279 dek%XLEV_CV(:) = dek%XLE_CV(:) - dek%XLES_CV(:)
1280 !
1281 IF (lhook) CALL dr_hook('ISBA_MEB:AVG_FLUXES_MEB_TSPLIT ',1,zhook_handle)
1282 !
1283 END SUBROUTINE avg_fluxes_meb_tsplit
1284 !
1285 !===============================================================================
1286 SUBROUTINE snowalb_spectral_bands_meb(PVEGTYPE,PEK,PSNOWRHO,PSNOWAGE,PPS, &
1287  PSNOWDZ,PZENITH,PTAU_N)
1288 !
1289 ! Split Total snow albedo into N-spectral bands
1290 ! NOTE currently MEB only uses 2 bands of the 3 possible.
1291 !
1292 USE modd_surf_par, ONLY : xundef
1293 USE modd_data_cover_par, ONLY : nvt_snow
1295 USE modd_snow_par, ONLY : nspec_band_snow
1296 USE modd_snow_metamo, ONLY : xsnowdzmin
1297 !
1298 USE mode_snow3l, ONLY : snow3lalb, snow3ldopt
1299 !
1300 IMPLICIT NONE
1301 !
1302 !* 0.1 declarations of arguments
1303 !
1304 TYPE(isba_pe_t), INTENT(INOUT) :: PEK
1305 REAL, DIMENSION(:,:), INTENT(IN) :: PVEGTYPE ! fraction of each vegetation (-)
1306 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWRHO ! Snow layer average density (kg/m3)
1307 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWDZ ! Snow layer thickness (m)
1308 REAL, DIMENSION(:), INTENT(IN) :: PZENITH ! Zenith angle (rad)
1309 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWAGE ! Snow grain age
1310 REAL, DIMENSION(:), INTENT(IN) :: PPS ! Pressure [Pa]
1311 REAL, DIMENSION(:,:), INTENT(OUT) :: PTAU_N ! SW absorption (coef) in uppermost snow layer (-)
1312 !
1313 !* 0.2 declarations of local variables
1314 !
1315 INTEGER :: JJ, JI, INJ, INLVLS
1316 REAL, DIMENSION(SIZE(PPS)) :: ZWORK, ZWORKA, ZAGE
1317 REAL, DIMENSION(SIZE(PPS)) :: ZPROJLAT, ZDSGRAIN, ZBETA1, ZBETA2, ZBETA3, &
1318  ZOPTICALPATH1, ZOPTICALPATH2, ZOPTICALPATH3
1319 REAL, DIMENSION(SIZE(PPS)) :: ZPERMSNOWFRAC
1320 REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)) :: ZSNOWDZ
1321 REAL, DIMENSION(SIZE(PPS),NSPEC_BAND_SNOW) :: ZSPECTRALALBEDO
1322 ! ZSPECTRALALBEDO = spectral albedo (3 bands in algo:
1323 ! MEB currently uses 2)
1324 ! 1=VIS, 2=NIR, 3=UV
1325 !
1326 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1327 !
1328 !-------------------------------------------------------------------------------
1329 !
1330 IF (lhook) CALL dr_hook('ISBA_MEB:SNOWALB_SPECTRAL_BANDS_MEB',0,zhook_handle)
1331 !
1332 inj = SIZE(psnowdz,1)
1333 inlvls = SIZE(psnowdz,2)
1334 !
1335 ! 1) Spectral albedo
1336 ! ------------------
1337 !
1338 zwork(:) = 0.0
1339 zworka(:) = pek%TSNOW%ALB(:)
1340 zpermsnowfrac(:) = pvegtype(:,nvt_snow)
1341 !
1342  CALL snow3lalb(zworka,zspectralalbedo,psnowrho(:,1),psnowage(:,1),zpermsnowfrac,pps)
1343 !
1344 ! Since we only consider VIS and NIR bands for soil and veg in MEB currently:
1345 ! (also note, here PSNOWALB doesn't evolve...we just diagnose spectral components).
1346 !
1347 WHERE(pek%TSNOW%ALB(:)/=xundef)
1348 !
1349  pek%TSNOW%ALBVIS(:) = zspectralalbedo(:,1)
1350 !
1351 ! We diagnose NIR albedo such that total albedo is conserved
1352 ! (using just 2 spectral bands in MEB)
1353 !
1354  pek%TSNOW%ALBNIR(:) = (pek%TSNOW%ALB(:) - xsw_wght_vis*pek%TSNOW%ALBVIS(:))/xsw_wght_nir
1355 !
1356 ! currently NOT used by MEB
1357 !
1358  pek%TSNOW%ALBFIR(:) = xundef
1359 !
1360 ! For the surface layer absorbtion computation:
1361 !
1362  zspectralalbedo(:,1) = pek%TSNOW%ALBVIS(:)
1363  zspectralalbedo(:,2) = pek%TSNOW%ALBNIR(:)
1364  zspectralalbedo(:,3) = pek%TSNOW%ALBFIR(:)
1365 !
1366 ELSEWHERE
1367 !
1368  pek%TSNOW%ALBVIS(:) = xundef
1369  pek%TSNOW%ALBNIR(:) = xundef
1370  pek%TSNOW%ALBFIR(:) = xundef
1371 !
1372 END WHERE
1373 !
1374 ! Snow optical grain diameter (no age dependency over polar regions):
1375 !
1376 zage(:) = (1.0-zpermsnowfrac(:))*psnowage(:,1)
1377 !
1378 zdsgrain(:) = snow3ldopt(psnowrho(:,1),zage)
1379 !
1380 ! 2) SW absorption in uppermost snow layer
1381 ! ----------------------------------------
1382 ! For now, consider just 2 bands with MEB, so renormalize:
1383 
1384 zspectralalbedo(:,2) = (pek%TSNOW%ALB(:) - xsw_wght_vis*zspectralalbedo(:,1))/xsw_wght_nir
1385 zspectralalbedo(:,3) = xundef
1386 !
1387 ! Adjust thickness to be as in snow computations:
1388 !
1389 DO jj=1,inlvls
1390  DO ji=1,inj
1391  zsnowdz(ji,jj) = psnowdz(ji,jj)/max(1.e-4,pek%XPSN(ji))
1392  ENDDO
1393 ENDDO
1394 !
1395  CALL snow3lradtrans(xsnowdzmin, zspectralalbedo, zsnowdz, psnowrho, &
1396  zpermsnowfrac, pzenith, psnowage, ptau_n)
1397 !
1398 ! Note that because we force a snow thickness to compute tramission,
1399 ! a bogus value ( < 0) can be computed despite the non-estance of snow.
1400 ! To check/prevent any problems, simply make a simple check:
1401 !
1402 ptau_n(:,:) = max(0., ptau_n(:,:))
1403 !
1404 IF (lhook) CALL dr_hook('ISBA_MEB:SNOWALB_SPECTRAL_BANDS_MEB',1,zhook_handle)
1405 !
1406 END SUBROUTINE snowalb_spectral_bands_meb
1407 !===============================================================================
1408  SUBROUTINE snow3lradtrans(PSNOWDZMIN, PSPECTRALALBEDO, PSNOWDZ, PSNOWRHO, &
1409  PPERMSNOWFRAC, PZENITH, PSNOWAGE, PRADTRANS)
1410 !
1411 !! PURPOSE
1412 !! -------
1413 ! Calculate the transmission of shortwave (solar) radiation
1414 ! through the snowpack (using a form of Beer's Law: exponential
1415 ! decay of radiation with increasing snow depth).
1416 !
1417 USE modd_surf_par, ONLY : xundef
1419 !
1421 !
1422 IMPLICIT NONE
1423 !
1424 !* 0.1 declarations of arguments
1425 !
1426 REAL, INTENT(IN) :: PSNOWDZMIN
1427 !
1428 REAL, DIMENSION(:), INTENT(IN) :: PPERMSNOWFRAC
1429 REAL, DIMENSION(:), INTENT(IN) :: PZENITH
1430 REAL, DIMENSION(:,:), INTENT(IN) :: PSNOWRHO, PSNOWDZ, PSNOWAGE
1431 REAL, DIMENSION(:,:), INTENT(IN) :: PSPECTRALALBEDO
1432 !
1433 REAL, DIMENSION(:,:), INTENT(OUT) :: PRADTRANS
1434 !
1435 !
1436 !* 0.2 declarations of local variables
1437 !
1438 INTEGER :: JJ, JI
1439 !
1440 INTEGER :: INJ
1441 INTEGER :: INLVLS
1442 !
1443 REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZDSGRAIN, ZCOEF, ZSNOWDZ, ZAGE
1444 !
1445 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1446 !-------------------------------------------------------------------------------
1447 !
1448 ! 0. Initialization:
1449 ! ------------------
1450 !
1451 IF (lhook) CALL dr_hook('SNOW3LRADTRANS',0,zhook_handle)
1452 !
1453 inj = SIZE(psnowdz(:,:),1)
1454 inlvls = SIZE(psnowdz(:,:),2)
1455 !
1456 !
1457 ! 1. Vanishingly thin snowpack check:
1458 ! -----------------------------------
1459 ! For vanishingly thin snowpacks, much of the radiation
1460 ! can pass through snowpack into underlying soil, making
1461 ! a large (albeit temporary) thermal gradient: by imposing
1462 ! a minimum thickness, this increases the radiation absorbtion
1463 ! for vanishingly thin snowpacks.
1464 !
1465 zsnowdz(:,:) = max(psnowdzmin, psnowdz(:,:))
1466 !
1467 !
1468 ! 2. Extinction of net shortwave radiation
1469 ! ----------------------------------------
1470 ! Fn of snow depth and density (Loth and Graf 1993:
1471 ! SNOWCVEXT => from Bohren and Barkstrom 1974
1472 ! SNOWAGRAIN and SNOWBGRAIN=> from Jordan 1976)
1473 !
1474 ! Snow optical grain diameter (no age dependency over polar regions):
1475 !
1476 zage(:,:) = 0.
1477 DO jj=1,inlvls
1478  DO ji=1,inj
1479  IF(psnowage(ji,jj)/=xundef)THEN
1480  zage(ji,jj) = (1.0-ppermsnowfrac(ji))*psnowage(ji,jj)
1481  ENDIF
1482  ENDDO
1483 ENDDO
1484 !
1485 zdsgrain(:,:) = snow3ldopt(psnowrho,zage)
1486 !
1487 zcoef(:,:) = snow3lradabs_sfc(psnowrho,zsnowdz,pspectralalbedo, &
1488  pzenith,ppermsnowfrac,zdsgrain)
1489 !
1490 ! 3. Radiation trans at base of each layer
1491 ! ----------------------------------
1492 ! NOTE, at level=0, rad = Swnet*(1-alb) so radcoef(0)=1 implicitly
1493 !
1494 pradtrans(:,:) = zcoef(:,:)
1495 !
1496 IF (lhook) CALL dr_hook('SNOW3LRADTRANS',1,zhook_handle)
1497 !
1498 !-------------------------------------------------------------------------------
1499 !
1500 END SUBROUTINE snow3lradtrans
1501 !===============================================================================
1502 
1505 IMPLICIT NONE
1506 !
1507 !* 0.2 declarations of local variables
1508 !
1509 INTEGER :: INLL
1510 !
1511 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1512 !
1513 !-------------------------------------------------------------------------------
1514 
1515 IF (lhook) CALL dr_hook('ISBA_MEB:ALLOCATE_LOCAL_VARS_PREP_GRID_SOIL ',0,zhook_handle)
1516 
1517 inll = inl
1518 IF(io%LMEB_LITTER)inll = inl + 1
1519 
1520 ALLOCATE ( ztgl(inj, inll ))
1521 ALLOCATE ( zsoilhcapz(inj, inll ))
1522 ALLOCATE ( zsoilcondz(inj, inll ))
1523 ALLOCATE ( zd_g(inj, inll ))
1524 ALLOCATE ( zdzg(inj, inll ))
1525 ALLOCATE ( zwfc(inj, inll ))
1526 ALLOCATE ( zwsat(inj, inll ))
1527 
1528 IF (lhook) CALL dr_hook('ISBA_MEB:ALLOCATE_LOCAL_VARS_PREP_GRID_SOIL ',1,zhook_handle)
1529 
1531 !===============================================================================
1534 IMPLICIT NONE
1535 !
1536 !* 0.2 declarations of local variables
1537 !
1538 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1539 !
1540 !-------------------------------------------------------------------------------
1541 
1542 IF (lhook) CALL dr_hook('ISBA_MEB:DEALLOCATE_LOCAL_VARS_PREP_GRID_SOIL ',0,zhook_handle)
1543 
1544 DEALLOCATE ( ztgl )
1545 DEALLOCATE ( zsoilhcapz )
1546 DEALLOCATE ( zsoilcondz )
1547 DEALLOCATE ( zd_g )
1548 DEALLOCATE ( zdzg )
1549 DEALLOCATE ( zwsat )
1550 DEALLOCATE ( zwfc )
1551 
1552 IF (lhook) CALL dr_hook('ISBA_MEB:DEALLOCATE_LOCAL_VARS_PREP_GRID_SOIL ',1,zhook_handle)
1553 
1555 !===============================================================================
1556 SUBROUTINE reshift_meb_soil(OMEB_LITTER,PTGL,PLESFC,PLESFCI,PEK,DEK)
1558 IMPLICIT NONE
1559 !
1560 !* 0.1 declarations of arguments
1561 !
1562 LOGICAL, INTENT(IN) :: OMEB_LITTER
1563 REAL, DIMENSION(:,:), INTENT(IN) :: PTGL
1564 REAL, DIMENSION(:), INTENT(IN) :: PLESFC
1565 REAL, DIMENSION(:), INTENT(IN) :: PLESFCI
1566 TYPE(isba_pe_t), INTENT(INOUT) :: PEK
1567 TYPE(diag_evap_isba_t), INTENT(INOUT) :: DEK
1568 !
1569 !* 0.2 declarations of local variables
1570 !
1571 INTEGER :: JJ, JL
1572 !
1573 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1574 !
1575 !-------------------------------------------------------------------------------
1576 
1577 inj = SIZE(pek%XTG(:,:),1)
1578 inl = SIZE(pek%XTG(:,:),2)
1579 
1580 IF (lhook) CALL dr_hook('ISBA_MEB:FINISH_MEB_SOIL ',0,zhook_handle)
1581 
1582 IF (omeb_litter)THEN
1583 
1584  pek%XTL(:) = ptgl(:,1)
1585 
1586  DO jl=1,inl
1587  DO jj=1,inj
1588  pek%XTG(jj,jl) = ptgl(jj,jl+1)
1589  ENDDO
1590  ENDDO
1591 
1592  dek%XLEG(:) = 0.0
1593  dek%XLEGI(:) = 0.0
1594  dek%XLELITTER(:) = plesfc(:)
1595  dek%XLELITTERI(:) = plesfci(:)
1596 ELSE
1597 
1598  pek%XTG(:,:) = ptgl(:,:)
1599 
1600  dek%XLEG(:) = plesfc(:)
1601  dek%XLEGI(:) = plesfci(:)
1602  dek%XLELITTER(:) = 0.
1603  dek%XLELITTERI(:) = 0.
1604 
1605 ENDIF
1606 
1607 
1608 IF (lhook) CALL dr_hook('ISBA_MEB:FINISH_MEB_SOIL ',1,zhook_handle)
1609 
1610 END SUBROUTINE reshift_meb_soil
1611 !===============================================================================
1612 SUBROUTINE prep_meb_soil(OMEB_LITTER,PSOILHCAPZ,PSOILCONDZ,KK,PK,PEK,PD_GL,&
1613  PDZGL,PTGL,PSOILHCAPL,PSOILCONDL,PWSATL,PWFCL,PWSFC,&
1614  PWISFC,PCTSFC,PCT,PFROZEN1,PFROZEN1SFC)
1615 !
1616 USE modd_csts, ONLY : xrholw,xrholi, xcl, xci
1617 USE modd_isba_par, ONLY : xwgmin, xomsph
1618 !
1619 IMPLICIT NONE
1620 !
1621 !* 0.1 declarations of arguments
1622 !
1623 TYPE(isba_k_t), INTENT(INOUT) :: KK
1624 TYPE(isba_p_t), INTENT(INOUT) :: PK
1625 TYPE(isba_pe_t), INTENT(INOUT) :: PEK
1626 LOGICAL, INTENT(IN) :: OMEB_LITTER
1627 REAL, DIMENSION(:,:), INTENT(IN) :: PSOILHCAPZ
1628 REAL, DIMENSION(:,:), INTENT(IN) :: PSOILCONDZ
1629 !
1630 REAL, DIMENSION(:), INTENT(IN) :: PCT
1631 REAL, DIMENSION(:), INTENT(IN) :: PFROZEN1
1632 REAL, DIMENSION(:,:), INTENT(OUT) :: PD_GL
1633 REAL, DIMENSION(:,:), INTENT(OUT) :: PDZGL
1634 REAL, DIMENSION(:,:), INTENT(OUT) :: PTGL
1635 REAL, DIMENSION(:,:), INTENT(OUT) :: PSOILHCAPL
1636 REAL, DIMENSION(:,:), INTENT(OUT) :: PSOILCONDL
1637 REAL, DIMENSION(:,:), INTENT(OUT) :: PWSATL
1638 REAL, DIMENSION(:,:), INTENT(OUT) :: PWFCL
1639 REAL, DIMENSION(:), INTENT(OUT) :: PWSFC
1640 REAL, DIMENSION(:), INTENT(OUT) :: PWISFC
1641 REAL, DIMENSION(:), INTENT(OUT) :: PCTSFC
1642 REAL, DIMENSION(:), INTENT(OUT) :: PFROZEN1SFC
1643 !
1644 !* 0.2 declarations of local variables
1645 !
1646 INTEGER :: INJ, INL, JJ, JL
1647 !
1648 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1649 !
1650 !* 0.3 declarations of local parameters
1651 !
1652 REAL, PARAMETER :: Z1 = 45.0 !litter bulk density (kg/m3)
1653 REAL, PARAMETER :: Z2 = 0.1 !coeff for litter conductivity (W/(mK))
1654 REAL, PARAMETER :: Z3 = 0.03 !coeff for litter conductivity
1655 REAL, PARAMETER :: Z4 = 0.95 !litter porosity (m3/m3)
1656 REAL, PARAMETER :: Z5 = 0.12 !litter field capacity (m3/m3)
1657 !
1658 REAL, DIMENSION(SIZE(PEK%XWG,1)) :: ZWORK
1659 !
1660 !-------------------------------------------------------------------------------
1661 !
1662 IF (lhook) CALL dr_hook('ISBA_MEB:PREP_MEB_SOIL',0,zhook_handle)
1663 !
1664 inj = SIZE(pk%XDG,1)
1665 inl = SIZE(pk%XDG,2)
1666 !
1667 zwork(:) = 0.0
1668 IF(omeb_litter)THEN
1669  ptgl(:,1) = pek%XTL(:)
1670  zwork(:) = pek%XWRL(:)/(xrholw*pek%XGNDLITTER(:))
1671  psoilhcapl(:,1) = xomsph*z1 + (xcl*xrholw)*zwork(:) + (xci*xrholi/xrholw)*pek%XWRLI(:)/pek%XGNDLITTER(:)
1672  psoilcondl(:,1) = z2 + z3 * zwork(:)
1673  pwsatl(:,1) = z4
1674  pwfcl(:,1) = z5
1675  pd_gl(:,1) = pek%XGNDLITTER(:)
1676  pdzgl(:,1) = pek%XGNDLITTER(:)
1677  pctsfc(:) = 1. / (psoilhcapl(:,1) * pek%XGNDLITTER(:))
1678  pfrozen1sfc(:) = pek%XWRLI(:) / ( pek%XWRLI(:) + max(pek%XWRL(:), (xwgmin*xrholw)*pek%XGNDLITTER(:) ))
1679 
1680  DO jl=1,inl
1681  DO jj=1,inj
1682  ptgl(jj,jl+1) = pek%XTG(jj,min(jl,SIZE(pek%XTG,2)))
1683  psoilhcapl(jj,jl+1) = psoilhcapz(jj,jl)
1684  psoilcondl(jj,jl+1) = psoilcondz(jj,jl)
1685  pwsatl(jj,jl+1) = kk%XWSAT(jj,jl)
1686  pwfcl(jj,jl+1) = kk%XWFC(jj,jl)
1687  pd_gl(jj,jl+1) = pek%XGNDLITTER(jj) + pk%XDG(jj,jl)
1688  pdzgl(jj,jl+1) = pk%XDZG(jj,jl)
1689  ENDDO
1690  ENDDO
1691  pwsfc(:) = pek%XWRL(:) /(xrholw*pek%XGNDLITTER(:)) ! (m3/m3)
1692  pwisfc(:) = pek%XWRLI(:)/(xrholw*pek%XGNDLITTER(:)) ! (m3/m3)
1693 
1694 ELSE
1695  ptgl(:,:) = pek%XTG(:,:)
1696  psoilhcapl(:,:) = psoilhcapz(:,:)
1697  psoilcondl(:,:) = psoilcondz(:,:)
1698  pwsatl(:,:) = kk%XWSAT(:,:)
1699  pwfcl(:,:) = kk%XWFC(:,:)
1700  pd_gl(:,:) = pk%XDG(:,:)
1701  pdzgl(:,:) = pk%XDZG(:,:)
1702  pctsfc(:) = pct(:)
1703  pwsfc(:) = pek%XWG(:,1)
1704  pwisfc(:) = pek%XWGI(:,1)
1705  pfrozen1sfc(:) = pfrozen1(:)
1706 ENDIF
1707 IF (lhook) CALL dr_hook('ISBA_MEB:PREP_MEB_SOIL',1,zhook_handle)
1708 
1709 END SUBROUTINE prep_meb_soil
1710 !===============================================================================
1711 SUBROUTINE ice_litter(PTSTEP, PLELITTERI, PSOILHCAPZ, PEK, &
1712  KWG_LAYER, PDZG, PPHASEL, PCTSFC, PLSTT, PLITCOR )
1714 USE modd_csts, ONLY : xlmtt, xtt, xci, xrholi, xrholw
1715 !
1716 IMPLICIT NONE
1717 !
1718 !* 0.1 declarations of arguments
1719 !
1720 TYPE(isba_pe_t), INTENT(INOUT) :: PEK
1721 REAL, INTENT(IN) :: PTSTEP
1722 ! PTSTEP = Model time step (s)
1723 !
1724 REAL, DIMENSION(:), INTENT(IN) :: PLELITTERI
1725 ! PLELITTERI = ice sublimation (m s-1)
1726 REAL, DIMENSION(:), INTENT(IN) :: PCTSFC
1727 REAL, DIMENSION(:), INTENT(IN) :: PLSTT
1728 !
1729 REAL, DIMENSION(:,:), INTENT(IN) :: PSOILHCAPZ
1730 ! PSOILHCAPZ = soil heat capacity [J/(m3 K)]
1731 REAL, DIMENSION(:,:), INTENT(IN) :: PDZG
1732 ! PDZG = Layer thickness (DIF option)
1733 !
1734 INTEGER, DIMENSION(:), INTENT(IN) :: KWG_LAYER
1735 ! KWG_LAYER = Number of soil moisture layers (DIF option)
1736 !
1737 REAL, DIMENSION(:), INTENT(OUT) :: PPHASEL
1738 ! PPHASEL = Phase changement in litter (W/m2)
1739 REAL, DIMENSION(:), INTENT(OUT) :: PLITCOR
1740 ! PLITCOR = A possible ice mass correction (to be potentially
1741 ! removed from soil) (kg/m2/s)
1742 !
1743 !* 0.2 declarations of local variables
1744 !
1745 INTEGER :: JL ! loop control
1746 !
1747 INTEGER :: INL ! Number of explicit soil layers
1748 !
1749 REAL, DIMENSION(SIZE(PEK%XTG,1)) :: ZEXCESS, ZK, ZHCAPL,ZELITTERI, &
1750  ZDELTAT,ZPHASE,ZPHASEM,ZPHASEF,ZPHASEX, &
1751  ZWRL,ZWRLI,Z0,ZPHASEC
1752 !
1753 REAL :: ZPSI
1754 !
1755 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1756 !
1757 !* 0.3 declaration of local parameters
1758 !
1759 REAL, PARAMETER :: ZERTOL = 1.e-6 ! (-) error tolerance
1760 REAL, PARAMETER :: ZTAUICE = 3300. ! (s) litter phase change characteristic time scale
1761 REAL, PARAMETER :: ZWRLSAT = 0.85 ! (m3/m3) litter porosity
1762 !
1763 !-------------------------------------------------------------------------------
1764 !
1765 IF (lhook) CALL dr_hook('ISBA_MEB:ICE_LITTER',0,zhook_handle)
1766 !
1767 ! Initialization:
1768 ! ---------------
1769 !
1770 IF (SIZE(kwg_layer)>0) THEN
1771  inl = maxval(kwg_layer(:))
1772 ELSE
1773  inl = SIZE(pek%XWG,2)
1774 ENDIF
1775 !
1776 zexcess(:) = 0.0
1777 zphasec(:) = 0.0
1778 plitcor(:) = 0.0
1779 !
1780 zhcapl(:) = 1/(pctsfc(:)*pek%XGNDLITTER(:))
1781 !
1782 !-------------------------------------------------------------------------------
1783 !
1784 ! 1. Surface layer vegetation insulation coefficient (-)
1785 ! ---------------------------------------------------
1786 !
1787 ! 1.1 Convert to m3/m3
1788 ! -----------------
1789 !
1790 zwrl(:)= pek%XWRL(:) /(xrholw*pek%XGNDLITTER(:))
1791 zwrli(:)=pek%XWRLI(:)/(xrholw*pek%XGNDLITTER(:))
1792 !
1793 ! 2. Litter ice evolution computation:
1794 ! --------------------------------
1795 !
1796 zdeltat(:) = pek%XTL(:) - xtt
1797 !
1798 !
1799 ! *Melt* ice if energy and ice available:
1800 zphasem(:) = (ptstep/ztauice)*min((xci*xrholi)*max(0.0,zdeltat),zwrli(:)*(xlmtt*xrholw))
1801 !
1802 ! *Freeze* liquid water if energy and water available and do not exceed porosity:
1803 zphasef(:) = (ptstep/ztauice)*min((xci*xrholi)*max(0.0,-zdeltat),zwrl(:)*(xlmtt*xrholw))
1804 zphasef(:) = min(zphasef(:) , (zwrlsat - zwrli(:)) * (xlmtt*xrholw) ) !!!!! LOOK!!!!!!!!
1805 !
1806 zphase(:) = zphasef(:) - zphasem(:)
1807 
1808 ! Update heat content if melting or freezing
1809 pek%XTL(:) = pek%XTL(:) + zphase(:)/zhcapl(:)
1810 !
1811 ! Get estimate of actual total phase change (J/m3) for equivalent litter water changes:
1812 
1813 zphasex(:) = zphase(:)
1814 !
1815 ! Adjust ice and liquid water conents (m3/m3) accordingly :
1816 !
1817 zwrl(:) = zwrl(:) - zphasex/(xlmtt*xrholw)
1818 zwrli(:) = zwrli(:) + zphasex/(xlmtt*xrholw)
1819 !
1820 ! 2.1 Convert to Kg/m2
1821 ! -----------------
1822 !
1823 pek%XWRL(:) = zwrl(:) * pek%XGNDLITTER(:) * xrholw
1824 pek%XWRLI(:)= zwrli(:) * pek%XGNDLITTER(:) * xrholw
1825 !
1826 ! 3. Adjust litter ice content for sublimation
1827 ! -----------------------------------------
1828 !
1829 zelitteri = plelitteri(:) * (ptstep/plstt(:))
1830 zexcess(:) = max( 0.0 , zelitteri - pek%XWRLI(:) )
1831 plitcor = zexcess / ptstep
1832 pek%XWRLI(:) = pek%XWRLI(:) - ( zelitteri - zexcess )
1833 !
1834 ! 4. Prevent some possible problems
1835 ! ------------------------------
1836 !
1837 pek%XWGI (:,1) = pek%XWGI(:,1)- zexcess / (xrholw * pdzg(:,1))
1838 !
1839 zexcess(:) = max( 0.0, - pek%XWGI(:,1) )
1840 pek%XWGI(:,1) = pek%XWGI(:,1) + zexcess(:)
1841 pek%XWG (:,1) = pek%XWG (:,1) - zexcess(:)
1842 pek%XTG (:,1) = pek%XTG (:,1) + zexcess(:) * (xlmtt*xrholw)/psoilhcapz(:,1)
1843 !
1844 DO jl=1,inl-1
1845  zexcess = max(0.0,-pek%XWG(:,jl))
1846  pek%XWG(:,jl+1) = pek%XWG(:,jl+1) - zexcess*pdzg(:,jl)/pdzg(:,jl+1)
1847  pek%XWG(:,jl) = pek%XWG(:,jl) + zexcess
1848 ENDDO
1849 !
1850 ! 5. Prevent from keeping track of ice in litter
1851 ! -------------------------------------------
1852 !
1853 WHERE (pek%XWRLI(:) < zertol )
1854  pek%XWRL (:) = pek%XWRL(:) + pek%XWRLI(:)
1855  pek%XTL (:) = pek%XTL(:) + pek%XWRLI(:) * xlmtt / pek%XGNDLITTER(:) / zhcapl(:)
1856  zphasec(:) = pek%XWRLI(:) * xlmtt / pek%XGNDLITTER(:)
1857  pek%XWRLI(:) = 0.0
1858 ELSEWHERE
1859  zphasec(:) = 0.0
1860 END WHERE
1861 !
1862 pphasel(:)=(zphase(:) + zphasec(:))/ptstep*pek%XGNDLITTER
1863 !
1864 IF (lhook) CALL dr_hook('ISBA_MEB:ICE_LITTER',1,zhook_handle)
1865 !
1866 END SUBROUTINE ice_litter
1867 
1868 !===============================================================================
1869 
1870 END SUBROUTINE isba_meb
subroutine snowalb_spectral_bands_meb(PVEGTYPE, PEK, PSNOWRHO, PSNOWAGE, PPS
Definition: isba_meb.F90:1287
real, parameter xsw_wght_vis
subroutine isba_meb(IO, KK, PK, PEK, DK, DEK, DMK, G, AG,
Definition: isba_meb.F90:7
real, parameter xsnowdzmin
real, save xcpd
Definition: modd_csts.F90:63
subroutine deallocate_local_vars_prep_grid_soil
Definition: isba_meb.F90:1533
subroutine hydro_veg(HRAIN, PTSTEP, PMUF, PRR, PLEV, PLETR, PVEG, PPSNV, PWR, PWRMAX, PPG, PDRIP, PRRVEG, PLVTT)
Definition: hydro_veg.F90:9
subroutine snow3lradtrans(PSNOWDZMIN, PSPECTRALALBEDO, PSNOWDZ, PS
Definition: isba_meb.F90:1409
real, save xstefan
Definition: modd_csts.F90:59
real, save xlvtt
Definition: modd_csts.F90:70
real, save xlstt
Definition: modd_csts.F90:71
real, parameter xundef
subroutine sum_fluxes_meb_tsplit
Definition: isba_meb.F90:1068
subroutine e_budget_meb(IO, KK, PK, PEK, DK, DEK, DMK, PTSTEP, PLTT, PPS, PCT, PTDEEP_A, PD_G, PS
Definition: e_budget_meb.F90:8
subroutine wet_leaves_frac(PWRM, PVEG, PWRMAX_CF, PZ0, PLAI, PWRMAX, PDELTA)
subroutine allocate_local_vars_prep_grid_soil
Definition: isba_meb.F90:1504
integer, parameter jprb
Definition: parkind1.F90:32
subroutine avg_fluxes_meb_tsplit
Definition: isba_meb.F90:1171
subroutine cotwores(PTSTEP, IO, OSHADE, PK, PEK, PDMAX, PPOI, PCSP, PTG, PF2, PSW_RAD, PQA, PQSAT, PPSNV, PDELTA, PRHOA, PZENITH, PFFV, PIACAN_SUNLIT, PIACAN_SHADE, PFRAC_SUN, PIACAN, PABC, PRS, PGPP, PRESP_LEAF)
Definition: cotwores.F90:10
subroutine init_sum_fluxes_meb_tsplit
Definition: isba_meb.F90:978
subroutine isba_fluxes_meb(KK, PK, PEK, DK, DEK, DMK, PRHOA, PLTT,
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
real, save xci
Definition: modd_csts.F90:65
real, save xday
Definition: modd_csts.F90:45
subroutine ice_litter(PTSTEP, PLELITTERI, PSOILHCAPZ, PEK, KWG_LAYER, PDZG, PPHASEL, PCTSFC, PLSTT, PLITCOR
Definition: isba_meb.F90:1713
real, save xcl
Definition: modd_csts.F90:65
subroutine snow_load_meb(PK, PEK, DEK, PTSTEP, PSR, PWRVNMAX, PKVN, PCHEATV, PMELTVN, PVELC, PSUBVCOR)
logical lhook
Definition: yomhook.F90:15
real, save xrholi
Definition: modd_csts.F90:81
subroutine reshift_meb_soil(OMEB_LITTER, PTGL, PLESFC, PLESFCI, PEK, DEK)
Definition: isba_meb.F90:1557
subroutine preps_for_meb_ebud_rad(PPS, PLAICV, PSNOWRHO, PSNOWSWE, PSNOWHEAT, PSNOWLIQ, PSNOWTEMP, PSNOWDZ, PSCOND, PHEATCAPS, PEMISNOW, PSIGMA_F, PCHIP, PTSTEP, PSR, PTA, PVMOD, PSNOWAGE, PPERMSNOWFRAC)
subroutine veg(PSW_RAD, PTA, PQA, PPS, PRGL, PLAI, PRSMIN, PGAMMA, PF2, PRS)
Definition: veg.F90:8
subroutine drag_meb(IO, PEK, DMK, DK, PTG, PTA, PQA, PVMOD, PWG, PWGI, PWSAT, PWFC, PEXNS, PEXNA, PPS, PRR, PSR, PRHOA, PZ0G_WITHOUT_SNOW, PZ0_MEBV, PZ0H_MEBV, PZ0EFF_MEBV, PZ0_MEBN, PZ0H_MEBN, PZ0EFF_MEBN, PSNOWSWE, PCHIP, PTSTEP, PRS_VG, PRS_VN, PPALPHAN, PZREF, PUREF, PDIRCOSZW, PSNCV, PDELTA, PVELC, PRISNOW, PUSTAR2SNOW, PHUGI, PHVG, PHVN, PLEG_DELTA, PLEGI_DELTA, PHSGL, PHSGF, PFLXC_C_A, PFLXC_N_A, PFLXC_G_C, PFLXC_N_C, PFLXC_VG_C, PFLXC_VN_C, PFLXC_MOM, PQSATG, PQSATV, PQSATC, PQSATN, PDELTAVK)
Definition: drag_meb.F90:17
real, save xrholw
Definition: modd_csts.F90:64
real, save xtt
Definition: modd_csts.F90:66
real, save xlmtt
Definition: modd_csts.F90:72
subroutine isba_lwnet_meb(PLAI, PPSN, PPSNA, PEMIS_N, PEMIS_F, PFF,
subroutine snow_leaves_frac_meb(PPSN, PPALPHAN, PWRVN, PTV, PCHIP, PLAIV, PWRVNMAX, PDELTAVN, PMELTVN)
subroutine prep_meb_soil(OMEB_LITTER, PSOILHCAPZ, PSOILCONDZ, KK, PK, PEK, PD_
Definition: isba_meb.F90:1613
subroutine radiative_transfert(OAGRI_TO_GRASS, PVEGTYPE, PALBVIS_VEG, PALBVIS_SOIL, PALBNIR_VEG, PALBNIR_SOIL, PSW_RAD, PLAI, PZENITH, PABC, PFAPARC, PFAPIRC, PMUS, PLAI_EFFC, OSHADE, PIACAN, PIACAN_SUNLIT, PIACAN_SHADE, PFRAC_SUN, PFAPAR, PFAPIR, PFAPAR_BS, PFAPIR_BS)
real, parameter xsw_wght_nir