[7541] | 1 | ! ============================================================================================================================== |
---|
| 2 | ! MODULE : enerbil |
---|
| 3 | ! |
---|
| 4 | ! CONTACT : orchidee-help _at_ listes.ipsl.fr |
---|
| 5 | ! |
---|
| 6 | ! LICENCE : IPSL (2006) |
---|
| 7 | ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
---|
| 8 | ! |
---|
| 9 | !>\BRIEF This module computes the energy balance on |
---|
| 10 | !! continental surfaces. The module contains the following subroutines: enerbil_initialize, enerbil_main, |
---|
| 11 | !! enerbil_clear, enerbil_begin, enerbil_surftemp, enerbil_flux and enerbil_evapveg |
---|
| 12 | !! |
---|
| 13 | !!\n DESCRIPTION : |
---|
| 14 | !! \n |
---|
| 15 | !! \latexonly |
---|
| 16 | !! \input{enerbil_intro2.tex} |
---|
| 17 | !! \endlatexonly |
---|
| 18 | !! |
---|
| 19 | !! IMPORTANT NOTE: The coefficients A and B are defined differently than those in the referenced |
---|
| 20 | !! literature and from those in the code and documentation of the atmospheric model LMDZ. For the |
---|
| 21 | !! avoidance of doubt, the coefficients as described here always refer to the ORCHIDEE coefficients, |
---|
| 22 | !! and are denoted as such (with the marker: ORC). The re-definition of the coefficients takes place |
---|
| 23 | !! within LMDZ before they are passed to ORCHIDEE. The following sequence of expressions is to be |
---|
| 24 | !! found within the LMDZ module 'surf_land_orchidee':\n |
---|
| 25 | !! |
---|
| 26 | !! \latexonly |
---|
| 27 | !! \input{surflandLMDZ1.tex} |
---|
| 28 | !! \input{surflandLMDZ2.tex} |
---|
| 29 | !! \input{surflandLMDZ3.tex} |
---|
| 30 | !! \input{surflandLMDZ4.tex} |
---|
| 31 | !! \endlatexonly |
---|
| 32 | !! \n |
---|
| 33 | !! |
---|
| 34 | !! \latexonly |
---|
| 35 | !! \input{enerbil_symbols.tex} |
---|
| 36 | !! \endlatexonly |
---|
| 37 | !! |
---|
| 38 | !! RECENT CHANGE(S) : None |
---|
| 39 | !! |
---|
| 40 | !! REFERENCE(S) : None |
---|
| 41 | !! |
---|
| 42 | !! SVN : |
---|
| 43 | !! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/enerbil.f90 $ |
---|
| 44 | !! $Date: 2021-07-30 18:36:29 +0200 (Fri, 30 Jul 2021) $ |
---|
| 45 | !! $Revision: 7265 $ |
---|
| 46 | !! \n |
---|
| 47 | !_ ================================================================================================================================ |
---|
| 48 | |
---|
| 49 | MODULE enerbil |
---|
| 50 | |
---|
| 51 | ! routines called : restput, restget |
---|
| 52 | ! |
---|
| 53 | ! modules used |
---|
| 54 | USE ioipsl |
---|
| 55 | USE xios_orchidee |
---|
| 56 | USE ioipsl_para |
---|
| 57 | USE constantes |
---|
| 58 | USE time, ONLY : one_day, dt_sechiba |
---|
| 59 | USE pft_parameters |
---|
| 60 | USE qsat_moisture |
---|
| 61 | USE sechiba_io_p |
---|
| 62 | USE constantes_soil |
---|
| 63 | USE explicitsnow |
---|
| 64 | |
---|
| 65 | IMPLICIT NONE |
---|
| 66 | |
---|
| 67 | PRIVATE |
---|
| 68 | PUBLIC :: enerbil_main, enerbil_initialize, enerbil_finalize, enerbil_clear |
---|
| 69 | |
---|
| 70 | ! variables used inside enerbil module : declaration and initialisation |
---|
| 71 | |
---|
| 72 | ! one dimension array allocated, computed and used in enerbil module exclusively |
---|
| 73 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: psold !! Old surface dry static energy (J kg^{-1}) |
---|
| 74 | !$OMP THREADPRIVATE(psold) |
---|
| 75 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsol_sat !! Saturated specific humudity for old temperature (kg kg^{-1}) |
---|
| 76 | !$OMP THREADPRIVATE(qsol_sat) |
---|
| 77 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: pdqsold !! Deriv. of saturated specific humidity at old temp |
---|
| 78 | !! (kg (kg s)^{-1}) |
---|
| 79 | !$OMP THREADPRIVATE(pdqsold) |
---|
| 80 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: psnew !! New surface static energy (J kg^{-1}) |
---|
| 81 | !$OMP THREADPRIVATE(psnew) |
---|
| 82 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsol_sat_new !! New saturated surface air moisture (kg kg^{-1}) |
---|
| 83 | !$OMP THREADPRIVATE(qsol_sat_new) |
---|
| 84 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: netrad !! Net radiation (W m^{-2}) |
---|
| 85 | !$OMP THREADPRIVATE(netrad) |
---|
| 86 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: lwabs !! LW radiation absorbed by the surface (W m^{-2}) |
---|
| 87 | !$OMP THREADPRIVATE(lwabs) |
---|
| 88 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: lwup !! Long-wave up radiation (W m^{-2}) |
---|
| 89 | !$OMP THREADPRIVATE(lwup) |
---|
| 90 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: lwnet !! Net Long-wave radiation (W m^{-2}) |
---|
| 91 | !$OMP THREADPRIVATE(lwnet) |
---|
| 92 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: fluxsubli !! Energy of sublimation (mm day^{-1}) |
---|
| 93 | !$OMP THREADPRIVATE(fluxsubli) |
---|
| 94 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsat_air !! Air saturated specific humidity (kg kg^{-1}) |
---|
| 95 | !$OMP THREADPRIVATE(qsat_air) |
---|
| 96 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tair !! Air temperature (K) |
---|
| 97 | !$OMP THREADPRIVATE(tair) |
---|
| 98 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: q_sol_pot !! Potential surface humidity |
---|
| 99 | !$OMP THREADPRIVATE(q_sol_pot) |
---|
| 100 | |
---|
| 101 | |
---|
| 102 | CONTAINS |
---|
| 103 | !! ============================================================================================================================= |
---|
| 104 | !! SUBROUTINE: enerbil_initialize |
---|
| 105 | !! |
---|
| 106 | !>\BRIEF Allocate module variables, read from restart file or initialize with default values |
---|
| 107 | !! |
---|
| 108 | !! DESCRIPTION: The purpose of this module is, firstly, to allocate space |
---|
| 109 | !! in memory for key variables within the 'enerbil' module. The second task is to retrieve previous data |
---|
| 110 | !! from the restart file. If the variables are not in restart file, default initialization is done. |
---|
| 111 | !! |
---|
| 112 | !! RECENT CHANGE(S): None |
---|
| 113 | !! |
---|
| 114 | !! MAIN OUTPUT VARIABLE(S): |
---|
| 115 | !! |
---|
| 116 | !! REFERENCE(S): None |
---|
| 117 | !! |
---|
| 118 | !! FLOWCHART: None |
---|
| 119 | !! \n |
---|
| 120 | !_ ============================================================================================================================== |
---|
| 121 | |
---|
| 122 | SUBROUTINE enerbil_initialize (kjit, kjpindex, index, rest_id, & |
---|
| 123 | qair, & |
---|
| 124 | temp_sol, temp_sol_new, tsol_rad, & |
---|
| 125 | evapot, evapot_corr, qsurf, fluxsens, & |
---|
| 126 | fluxlat, vevapp ) |
---|
| 127 | |
---|
| 128 | !! 0 Variable and parameter description |
---|
| 129 | !! 0.1 Input variables |
---|
| 130 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-) |
---|
| 131 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
| 132 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-) |
---|
| 133 | INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier (-) |
---|
| 134 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1}) |
---|
| 135 | |
---|
| 136 | !! 0.2 Output variables |
---|
| 137 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol !! Soil temperature (K) |
---|
| 138 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature (K) |
---|
| 139 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Tsol_rad (W m^{-2}) |
---|
| 140 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: evapot !! Soil Potential Evaporation (mm/tstep) |
---|
| 141 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: evapot_corr !! Soil Potential Evaporation Correction (mm/tstep) |
---|
| 142 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf !! Surface specific humidity (kg kg^{-1}) |
---|
| 143 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux (W m^{-2}) |
---|
| 144 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux (W m^{-2}) |
---|
| 145 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation (mm day^{-1}) |
---|
| 146 | |
---|
| 147 | |
---|
| 148 | !! 0.4 Local variables |
---|
| 149 | INTEGER(i_std) :: ier |
---|
| 150 | |
---|
| 151 | !_ ================================================================================================================================ |
---|
| 152 | |
---|
| 153 | IF (printlev>=3) WRITE (numout,*) 'enerbil_initialize : Start initillization' |
---|
| 154 | |
---|
| 155 | !! 1. Allocate module variables |
---|
| 156 | ALLOCATE (psold(kjpindex),stat=ier) |
---|
| 157 | IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable ','','') |
---|
| 158 | |
---|
| 159 | ALLOCATE (qsol_sat(kjpindex),stat=ier) |
---|
| 160 | IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable ','','') |
---|
| 161 | |
---|
| 162 | ALLOCATE (pdqsold(kjpindex),stat=ier) |
---|
| 163 | IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable ','','') |
---|
| 164 | |
---|
| 165 | ALLOCATE (psnew(kjpindex),stat=ier) |
---|
| 166 | IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable ','','') |
---|
| 167 | |
---|
| 168 | ALLOCATE (qsol_sat_new(kjpindex),stat=ier) |
---|
| 169 | IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable qsol_sat_new','','') |
---|
| 170 | |
---|
| 171 | ALLOCATE (netrad(kjpindex),stat=ier) |
---|
| 172 | IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable netrad','','') |
---|
| 173 | |
---|
| 174 | ALLOCATE (lwabs(kjpindex),stat=ier) |
---|
| 175 | IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable lwabs','','') |
---|
| 176 | |
---|
| 177 | ALLOCATE (lwup(kjpindex),stat=ier) |
---|
| 178 | IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable lwup','','') |
---|
| 179 | |
---|
| 180 | ALLOCATE (lwnet(kjpindex),stat=ier) |
---|
| 181 | IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable lwnet','','') |
---|
| 182 | |
---|
| 183 | ALLOCATE (fluxsubli(kjpindex),stat=ier) |
---|
| 184 | IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable fluxsubli','','') |
---|
| 185 | |
---|
| 186 | ALLOCATE (qsat_air(kjpindex),stat=ier) |
---|
| 187 | IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable qsat_air','','') |
---|
| 188 | |
---|
| 189 | ALLOCATE (tair(kjpindex),stat=ier) |
---|
| 190 | IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable tair','','') |
---|
| 191 | |
---|
| 192 | ALLOCATE (q_sol_pot(kjpindex),stat=ier) |
---|
| 193 | IF (ier /= 0) CALL ipslerr_p(3,'enerbil_initialize','Problem in allocation of variable q_sol_pot','','') |
---|
| 194 | |
---|
| 195 | |
---|
| 196 | !! 2. Initialize variables from restart file or by default values |
---|
| 197 | !! The variables read are: temp_sol (surface temperature), qsurf (near surface specific humidity), |
---|
| 198 | !! evapot (soil potential evaporation), evapot_corr (corrected soil potential evaporation), tsolrad |
---|
| 199 | !! (radiative surface temperature), evapora (evaporation), fluxlat (latent heat flux), fluxsens |
---|
| 200 | !! (sensible heat flux) and temp_sol_new (new surface temperature). |
---|
| 201 | IF (printlev>=3) WRITE (numout,*) 'Read a restart file for ENERBIL variables' |
---|
| 202 | |
---|
| 203 | CALL ioconf_setatt_p('UNITS', 'K') |
---|
| 204 | CALL ioconf_setatt_p('LONG_NAME','Surface temperature') |
---|
| 205 | CALL restget_p (rest_id, 'temp_sol', nbp_glo, 1, 1, kjit, .TRUE., temp_sol, "gather", nbp_glo, index_g) |
---|
| 206 | |
---|
| 207 | !Config Key = ENERBIL_TSURF |
---|
| 208 | !Config Desc = Initial temperature if not found in restart |
---|
| 209 | !Config If = OK_SECHIBA |
---|
| 210 | !Config Def = 280. |
---|
| 211 | !Config Help = The initial value of surface temperature if its value is not found |
---|
| 212 | !Config in the restart file. This should only be used if the model is |
---|
| 213 | !Config started without a restart file. |
---|
| 214 | !Config Units = Kelvin [K] |
---|
| 215 | CALL setvar_p (temp_sol, val_exp,'ENERBIL_TSURF', 280._r_std) |
---|
| 216 | |
---|
| 217 | ! Initialize temp_sol_new with temp_sol. These variables are always equal in the beginning of a new time step. |
---|
| 218 | temp_sol_new(:) = temp_sol(:) |
---|
| 219 | |
---|
| 220 | CALL ioconf_setatt_p('UNITS', 'g/g') |
---|
| 221 | CALL ioconf_setatt_p('LONG_NAME','near surface specific humidity') |
---|
| 222 | CALL restget_p (rest_id, 'qsurf', nbp_glo, 1, 1, kjit, .TRUE., qsurf, "gather", nbp_glo, index_g) |
---|
| 223 | IF ( ALL( qsurf(:) .EQ. val_exp ) ) THEN |
---|
| 224 | qsurf(:) = qair(:) |
---|
| 225 | ENDIF |
---|
| 226 | |
---|
| 227 | CALL ioconf_setatt_p('UNITS', 'mm day^{-1}') |
---|
| 228 | CALL ioconf_setatt_p('LONG_NAME','Soil Potential Evaporation') |
---|
| 229 | CALL restget_p (rest_id, 'evapot', nbp_glo, 1, 1, kjit, .TRUE., evapot, "gather", nbp_glo, index_g) |
---|
| 230 | CALL setvar_p (evapot, val_exp, 'ENERBIL_EVAPOT', zero) |
---|
| 231 | |
---|
| 232 | CALL ioconf_setatt_p('UNITS', 'mm day^{-1}') |
---|
| 233 | CALL ioconf_setatt_p('LONG_NAME','Corrected Soil Potential Evaporation') |
---|
| 234 | CALL restget_p (rest_id, 'evapot_corr', nbp_glo, 1, 1, kjit, .TRUE., evapot_corr, "gather", nbp_glo, index_g) |
---|
| 235 | !Config Key = ENERBIL_EVAPOT |
---|
| 236 | !Config Desc = Initial Soil Potential Evaporation |
---|
| 237 | !Config If = OK_SECHIBA |
---|
| 238 | !Config Def = 0.0 |
---|
| 239 | !Config Help = The initial value of soil potential evaporation if its value |
---|
| 240 | !Config is not found in the restart file. This should only be used if |
---|
| 241 | !Config the model is started without a restart file. |
---|
| 242 | !Config Units = |
---|
| 243 | CALL setvar_p (evapot_corr, val_exp, 'ENERBIL_EVAPOT', zero) |
---|
| 244 | |
---|
| 245 | CALL ioconf_setatt_p('UNITS', 'K') |
---|
| 246 | CALL ioconf_setatt_p('LONG_NAME','Radiative surface temperature') |
---|
| 247 | CALL restget_p (rest_id, 'tsolrad', nbp_glo, 1, 1, kjit, .TRUE., tsol_rad, "gather", nbp_glo, index_g) |
---|
| 248 | IF ( ALL( tsol_rad(:) .EQ. val_exp ) ) THEN |
---|
| 249 | tsol_rad(:) = temp_sol(:) |
---|
| 250 | ENDIF |
---|
| 251 | |
---|
| 252 | !! 1.3 Set the fluxes so that we have something reasonable and not NaN on some machines |
---|
| 253 | CALL ioconf_setatt_p('UNITS', 'Kg/m^2/dt') |
---|
| 254 | CALL ioconf_setatt_p('LONG_NAME','Evaporation') |
---|
| 255 | CALL restget_p (rest_id, 'evapora', nbp_glo, 1, 1, kjit, .TRUE., vevapp, "gather", nbp_glo, index_g) |
---|
| 256 | IF ( ALL( vevapp(:) .EQ. val_exp ) ) THEN |
---|
| 257 | vevapp(:) = zero |
---|
| 258 | ENDIF |
---|
| 259 | |
---|
| 260 | CALL ioconf_setatt_p('UNITS', 'W/m^2') |
---|
| 261 | CALL ioconf_setatt_p('LONG_NAME','Latent heat flux') |
---|
| 262 | CALL restget_p (rest_id, 'fluxlat', nbp_glo, 1, 1, kjit, .TRUE., fluxlat, "gather", nbp_glo, index_g) |
---|
| 263 | IF ( ALL( fluxlat(:) .EQ. val_exp ) ) THEN |
---|
| 264 | fluxlat(:) = zero |
---|
| 265 | ENDIF |
---|
| 266 | |
---|
| 267 | CALL ioconf_setatt_p('UNITS', 'W/m^2') |
---|
| 268 | CALL ioconf_setatt_p('LONG_NAME','Sensible heat flux') |
---|
| 269 | CALL restget_p (rest_id, 'fluxsens', nbp_glo, 1, 1, kjit, .TRUE., fluxsens, "gather", nbp_glo, index_g) |
---|
| 270 | IF ( ALL( fluxsens(:) .EQ. val_exp ) ) THEN |
---|
| 271 | fluxsens(:) = zero |
---|
| 272 | ENDIF |
---|
| 273 | |
---|
| 274 | |
---|
| 275 | CALL ioconf_setatt_p('UNITS', 'kg/m^2') |
---|
| 276 | CALL ioconf_setatt_p('LONG_NAME','Potential saturated surface humidity') |
---|
| 277 | CALL restget_p (rest_id, 'qsolpot', nbp_glo, 1, 1, kjit, .TRUE., q_sol_pot, "gather", nbp_glo, index_g) |
---|
| 278 | IF ( ALL( q_sol_pot(:) .EQ. val_exp ) ) THEN |
---|
| 279 | q_sol_pot = qsurf |
---|
| 280 | ENDIF |
---|
| 281 | |
---|
| 282 | END SUBROUTINE enerbil_initialize |
---|
| 283 | |
---|
| 284 | |
---|
| 285 | |
---|
| 286 | !! =========================================================================================================================== |
---|
| 287 | !! SUBROUTINE : enerbil_main |
---|
| 288 | !! |
---|
| 289 | !! |
---|
| 290 | !>\BRIEF Calls each part of the energy budget calculation in sequence |
---|
| 291 | !! |
---|
| 292 | !! DESCRIPTION : |
---|
| 293 | !! This is the main routine for the 'enerbil' module. It is called |
---|
| 294 | !! once during the initialisation of ORCHIDEE, and then once for each time step. It is called a final |
---|
| 295 | !! time at the culmination of the last time step, for the writing of a restart file.\n |
---|
| 296 | !! |
---|
| 297 | !! The algorithm first calls 'enerbil_begin' for initialisation, followed by 'enerbil_surftemp' to |
---|
| 298 | !! calculate the surface static energy and the saturated surface humidity for the new time step. |
---|
| 299 | !! Next is the module 'enerbil_flux' which calculates the new surface temperature, the net radiation in |
---|
| 300 | !! the surface layer, the total evaporation and the latent and sensible heat fluxes. Finally comes |
---|
| 301 | !! 'enerbil_evapveg', which calculates the evaporation and transpiration from the vegetation.\n |
---|
| 302 | |
---|
| 303 | !! \n |
---|
| 304 | !! |
---|
| 305 | !! RECENT CHANGE(S): None |
---|
| 306 | !! |
---|
| 307 | !! MAIN OUTPUT VARIABLE(S) : evapot, evapot_corr, temp_sol, qsurf, fluxsens, fluxlat, tsol_rad, |
---|
| 308 | !! vevapp, temp_sol_new, vevapnu, vevapsno, transpir, vevapwet |
---|
| 309 | !! |
---|
| 310 | !! REFERENCE(S) : |
---|
| 311 | !! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled |
---|
| 312 | !! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278 |
---|
| 313 | !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations |
---|
| 314 | !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General |
---|
| 315 | !! Circulation Model. Journal of Climate, 6, pp.248-273 |
---|
| 316 | !! - Dufresne, J-L & Ghattas, J, 2009. Description du schéma de la couche limite turbulente et la |
---|
| 317 | !! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11): |
---|
| 318 | !! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf |
---|
| 319 | !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation |
---|
| 320 | !! sur le cycle de l'eau, PhD Thesis, available (22/12/11): |
---|
| 321 | !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf |
---|
| 322 | !! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J |
---|
| 323 | !! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D & |
---|
| 324 | !! Xue, Y, 1998. A proposal for a general interface between land surface schemes and |
---|
| 325 | !! general circulation models. Global and Planetary Change, 19, pp.261-276 |
---|
| 326 | !! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems. |
---|
| 327 | !! Interscience Publishers\n |
---|
| 328 | !! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere |
---|
| 329 | !! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663. |
---|
| 330 | !! |
---|
| 331 | !! FLOWCHART : |
---|
| 332 | !! \latexonly |
---|
| 333 | !! \includegraphics[scale=0.5]{enerbil_main_flowchart.png} |
---|
| 334 | !! \endlatexonly |
---|
| 335 | !! \n |
---|
| 336 | !_ ============================================================================================================================== |
---|
| 337 | |
---|
| 338 | SUBROUTINE enerbil_main (kjit, kjpindex, & |
---|
| 339 | & index, indexveg, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, & |
---|
| 340 | & qair, peqAcoef, peqBcoef, pb, rau, vbeta, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, & |
---|
| 341 | & emis, soilflx, soilcap, q_cdrag, humrel, fluxsens, fluxlat, & |
---|
| 342 | & vevapp, transpir, transpot, vevapnu, vevapwet, vevapsno, vevapflo, temp_sol, tsol_rad, & |
---|
| 343 | & temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id, & |
---|
| 344 | & precip_rain, pgflux, snowdz, temp_sol_add) |
---|
| 345 | |
---|
| 346 | |
---|
| 347 | !! 0 Variable and parameter description |
---|
| 348 | |
---|
| 349 | !! 0.1 Input variables |
---|
| 350 | |
---|
| 351 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-) |
---|
| 352 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
| 353 | INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and _history_ file identifier (-) |
---|
| 354 | INTEGER(i_std),INTENT (in) :: hist2_id !! _history_ file 2 identifier (-) |
---|
| 355 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-) |
---|
| 356 | INTEGER(i_std),DIMENSION(kjpindex*nvm), INTENT(in) :: indexveg !! Indeces of the points on the 3D map |
---|
| 357 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m) |
---|
| 358 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux (W m^{-2}) |
---|
| 359 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2}) |
---|
| 360 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy (?? J) |
---|
| 361 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature (K) |
---|
| 362 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1}) |
---|
| 363 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1}) |
---|
| 364 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! PetAcoef (see note) |
---|
| 365 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! PetBcoef (see note) |
---|
| 366 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1}) |
---|
| 367 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! PeqAcoef (see note) |
---|
| 368 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! PeqBcoef (see note) |
---|
| 369 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa) |
---|
| 370 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air density (kg m^{-3}) |
---|
| 371 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient (-) |
---|
| 372 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-) |
---|
| 373 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta4 !! Bare soil resistance (-) |
---|
| 374 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance |
---|
| 375 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-) |
---|
| 376 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilflx !! Effective ground heat flux including both snow and soil (W m^{-2}) |
---|
| 377 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil calorific capacity including both snow and soil (J K^{-1]) |
---|
| 378 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) |
---|
| 379 | |
---|
| 380 | REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (in) :: humrel !! Soil moisture stress (within range 0 to 1) |
---|
| 381 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta2 !! Interception resistance (-) |
---|
| 382 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3 !! Vegetation resistance (-) |
---|
| 383 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3pot !! Vegetation resistance for potential transpiration |
---|
| 384 | REAL(r_std),DIMENSION (kjpindex),INTENT (in) :: precip_rain !! Rainfall |
---|
| 385 | REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth at each snow layer |
---|
| 386 | |
---|
| 387 | !! 0.2 Output variables |
---|
| 388 | |
---|
| 389 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapnu !! Bare soil evaporation (mm day^{-1}) |
---|
| 390 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapsno !! Snow evaporation (mm day^{-1}) |
---|
| 391 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapflo !! Floodplains evaporation |
---|
| 392 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpir !! Transpiration (mm day^{-1}) |
---|
| 393 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpot !! Potential transpiration |
---|
| 394 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vevapwet !! Interception (mm day^{-1}) |
---|
| 395 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature (K) |
---|
| 396 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K) |
---|
| 397 | |
---|
| 398 | !! 0.3 Modified variables |
---|
| 399 | |
---|
| 400 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot !! Soil Potential Evaporation (mm/tstep) |
---|
| 401 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot_corr !! Soil Potential Evaporation Correction (mm/tstep) |
---|
| 402 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol !! Soil temperature (K) |
---|
| 403 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: qsurf !! Surface specific humidity (kg kg^{-1}) |
---|
| 404 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux (W m^{-2}) |
---|
| 405 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux (W m^{-2}) |
---|
| 406 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: tsol_rad !! Tsol_rad (W m^{-2}) |
---|
| 407 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation (mm day^{-1}) |
---|
| 408 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: pgflux !! Net energy into snowpack(W/m^2) |
---|
| 409 | |
---|
| 410 | !! 0.4 Local variables |
---|
| 411 | |
---|
| 412 | REAL(r_std),DIMENSION (kjpindex) :: epot_air_new, qair_new |
---|
| 413 | REAL(r_std),DIMENSION (kjpindex) :: diffevap !! Difference betwence vevapp and composing fluxes (Kg/m^2/s) |
---|
| 414 | INTEGER(i_std) :: ji,ii,iv |
---|
| 415 | |
---|
| 416 | |
---|
| 417 | !_ ================================================================================================================================ |
---|
| 418 | |
---|
| 419 | !! 1. Computes some initialisation variables |
---|
| 420 | |
---|
| 421 | ! Computes some initialisation variables psold (the old surface static energy), qsol_sat |
---|
| 422 | ! (the saturated surface humidity) and pdqsold (the derivative of the saturated surface humidity, |
---|
| 423 | ! calculated with respect to the surface temperature at the 'old' timestep) |
---|
| 424 | CALL enerbil_begin (kjpindex, temp_sol, lwdown, swnet, pb, psold, qsol_sat, pdqsold, netrad, emis) |
---|
| 425 | |
---|
| 426 | |
---|
| 427 | !! 2. Computes psnew (the surface static energy at the 'new' timestep) |
---|
| 428 | |
---|
| 429 | ! Computes psnew (the surface static energy at the 'new' timestep), qsol_sat_new (the surface |
---|
| 430 | ! saturated humidity at the 'new' timestep), temp_sol_new (the surface temperature at the 'new' |
---|
| 431 | ! timestep), qair_new (the lowest atmospheric humidity at the 'new' timestep) and epot_air_new |
---|
| 432 | ! (the lowest atmospheric evaporation potential at the 'new' timestep) |
---|
| 433 | CALL enerbil_surftemp (kjpindex, zlev, emis, & |
---|
| 434 | & epot_air, petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,& |
---|
| 435 | & vbeta1, vbeta5, soilcap, lwdown, swnet, psnew, qsol_sat_new, temp_sol_new, & |
---|
| 436 | & qair_new, epot_air_new) |
---|
| 437 | |
---|
| 438 | |
---|
| 439 | !! 3. Diagnose components of the energy budget |
---|
| 440 | |
---|
| 441 | ! Diagnoses lwup (upwards long wave radiation), lwnet (net long wave radiation), tsol_rad (radiative |
---|
| 442 | ! ground temperature), netrad (net radiation), qsurf (surface humidity), vevapp (total evaporation), |
---|
| 443 | ! evapot (evaporation potential), evapot_corr (evaporation potential correction factor), |
---|
| 444 | ! fluxlat (latent heat flux), fluxsubli (sublimination heat flux) and fluxsens (sensible heat flux). |
---|
| 445 | |
---|
| 446 | CALL enerbil_flux (kjpindex, emis, temp_sol, rau, u, v, q_cdrag, vbeta, vbeta1, vbeta5, & |
---|
| 447 | & qair_new, epot_air_new, psnew, qsurf, & |
---|
| 448 | & fluxsens , fluxlat , fluxsubli, vevapp, temp_sol_new, lwdown, swnet, & |
---|
| 449 | & lwup, lwnet, pb, tsol_rad, netrad, evapot, evapot_corr, & |
---|
| 450 | & precip_rain,snowdz,temp_air,pgflux, soilcap, temp_sol_add) |
---|
| 451 | |
---|
| 452 | |
---|
| 453 | !! 4. Diagnoses the values for evaporation and transpiration |
---|
| 454 | |
---|
| 455 | ! Diagnoses the values for evaporation and transpiration: vevapsno (snow evaporation), vevapnu |
---|
| 456 | ! (bare soil evaporation), transpir (transpiration) and vevapwet (interception) |
---|
| 457 | CALL enerbil_evapveg (kjpindex, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, & |
---|
| 458 | & rau, u, v, q_cdrag, qair_new, humrel, vevapsno, vevapnu , vevapflo, vevapwet, & |
---|
| 459 | & transpir, transpot, evapot) |
---|
| 460 | |
---|
| 461 | |
---|
| 462 | |
---|
| 463 | !! 5. Write diagnosics |
---|
| 464 | |
---|
| 465 | CALL xios_orchidee_send_field("netrad",netrad) |
---|
| 466 | CALL xios_orchidee_send_field("evapot",evapot/dt_sechiba) |
---|
| 467 | CALL xios_orchidee_send_field("evapot_corr",evapot_corr/dt_sechiba) |
---|
| 468 | CALL xios_orchidee_send_field("lwdown",lwabs) |
---|
| 469 | CALL xios_orchidee_send_field("lwnet",lwnet) |
---|
| 470 | CALL xios_orchidee_send_field("Qv",fluxsubli) |
---|
| 471 | |
---|
| 472 | DO ji=1,kjpindex |
---|
| 473 | diffevap(ji) = vevapp(ji) - ( SUM(vevapwet(ji,:)) + & |
---|
| 474 | SUM(transpir(ji,:)) + vevapnu(ji) + vevapsno(ji) + vevapflo(ji) ) |
---|
| 475 | ENDDO |
---|
| 476 | CALL xios_orchidee_send_field("diffevap",diffevap/dt_sechiba) ! mm/s |
---|
| 477 | |
---|
| 478 | IF ( .NOT. almaoutput ) THEN |
---|
| 479 | CALL histwrite_p(hist_id, 'netrad', kjit, netrad, kjpindex, index) |
---|
| 480 | CALL histwrite_p(hist_id, 'evapot', kjit, evapot, kjpindex, index) |
---|
| 481 | CALL histwrite_p(hist_id, 'evapot_corr', kjit, evapot_corr, kjpindex, index) |
---|
| 482 | CALL histwrite_p(hist_id, 'lwdown', kjit, lwabs, kjpindex, index) |
---|
| 483 | CALL histwrite_p(hist_id, 'lwnet', kjit, lwnet, kjpindex, index) |
---|
| 484 | IF ( hist2_id > 0 ) THEN |
---|
| 485 | CALL histwrite_p(hist2_id, 'netrad', kjit, netrad, kjpindex, index) |
---|
| 486 | CALL histwrite_p(hist2_id, 'evapot', kjit, evapot, kjpindex, index) |
---|
| 487 | CALL histwrite_p(hist2_id, 'evapot_corr', kjit, evapot_corr, kjpindex, index) |
---|
| 488 | CALL histwrite_p(hist2_id, 'lwdown', kjit, lwabs, kjpindex, index) |
---|
| 489 | CALL histwrite_p(hist2_id, 'lwnet', kjit, lwnet, kjpindex, index) |
---|
| 490 | ENDIF |
---|
| 491 | ELSE |
---|
| 492 | CALL histwrite_p(hist_id, 'LWnet', kjit, lwnet, kjpindex, index) |
---|
| 493 | CALL histwrite_p(hist_id, 'Qv', kjit, fluxsubli, kjpindex, index) |
---|
| 494 | CALL histwrite_p(hist_id, 'PotEvap', kjit, evapot_corr, kjpindex, index) |
---|
| 495 | CALL histwrite_p(hist_id, 'PotEvapOld', kjit, evapot, kjpindex, index) |
---|
| 496 | IF ( hist2_id > 0 ) THEN |
---|
| 497 | CALL histwrite_p(hist2_id, 'LWnet', kjit, lwnet, kjpindex, index) |
---|
| 498 | CALL histwrite_p(hist2_id, 'Qv', kjit, fluxsubli, kjpindex, index) |
---|
| 499 | CALL histwrite_p(hist2_id, 'PotEvap', kjit, evapot_corr, kjpindex, index) |
---|
| 500 | ENDIF |
---|
| 501 | ENDIF |
---|
| 502 | |
---|
| 503 | IF (printlev>=3) WRITE (numout,*) ' enerbil_main Done ' |
---|
| 504 | |
---|
| 505 | END SUBROUTINE enerbil_main |
---|
| 506 | |
---|
| 507 | |
---|
| 508 | !! ============================================================================================================================= |
---|
| 509 | !! SUBROUTINE: enerbil_finalize |
---|
| 510 | !! |
---|
| 511 | !>\BRIEF Write to restart file |
---|
| 512 | !! |
---|
| 513 | !! DESCRIPTION: This subroutine writes the module variables and variables calculated in enerbil |
---|
| 514 | !! to restart file |
---|
| 515 | !! |
---|
| 516 | !! RECENT CHANGE(S): None |
---|
| 517 | !! |
---|
| 518 | !! REFERENCE(S): None |
---|
| 519 | !! |
---|
| 520 | !! FLOWCHART: None |
---|
| 521 | !! \n |
---|
| 522 | !_ ============================================================================================================================== |
---|
| 523 | SUBROUTINE enerbil_finalize (kjit, kjpindex, rest_id, & |
---|
| 524 | evapot, evapot_corr, temp_sol, tsol_rad, & |
---|
| 525 | qsurf, fluxsens, fluxlat, vevapp ) |
---|
| 526 | |
---|
| 527 | !! 0 Variable and parameter description |
---|
| 528 | !! 0.1 Input variables |
---|
| 529 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-) |
---|
| 530 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
| 531 | INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier (-) |
---|
| 532 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation (mm/tstep) |
---|
| 533 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr !! Soil Potential Evaporation Correction (mm/tstep) |
---|
| 534 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Soil temperature (K) |
---|
| 535 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tsol_rad !! Tsol_rad (W m^{-2}) |
---|
| 536 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! Surface specific humidity (kg kg^{-1}) |
---|
| 537 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: fluxsens !! Sensible heat flux (W m^{-2}) |
---|
| 538 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: fluxlat !! Latent heat flux (W m^{-2}) |
---|
| 539 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapp !! Total of evaporation (mm day^{-1}) |
---|
| 540 | |
---|
| 541 | !_ ================================================================================================================================ |
---|
| 542 | |
---|
| 543 | !! 1. Write variables to restart file to be used for the next simulation |
---|
| 544 | IF (printlev>=3) WRITE (numout,*) 'Write restart file with ENERBIL variables' |
---|
| 545 | |
---|
| 546 | CALL restput_p(rest_id, 'temp_sol', nbp_glo, 1, 1, kjit, temp_sol, 'scatter', nbp_glo, index_g) |
---|
| 547 | CALL restput_p(rest_id, 'qsurf', nbp_glo, 1, 1, kjit, qsurf, 'scatter', nbp_glo, index_g) |
---|
| 548 | CALL restput_p(rest_id, 'evapot', nbp_glo, 1, 1, kjit, evapot, 'scatter', nbp_glo, index_g) |
---|
| 549 | CALL restput_p(rest_id, 'evapot_corr', nbp_glo, 1, 1, kjit, evapot_corr, 'scatter', nbp_glo, index_g) |
---|
| 550 | CALL restput_p(rest_id, 'tsolrad', nbp_glo, 1, 1, kjit, tsol_rad, 'scatter', nbp_glo, index_g) |
---|
| 551 | CALL restput_p(rest_id, 'evapora', nbp_glo, 1, 1, kjit, vevapp, 'scatter', nbp_glo, index_g) |
---|
| 552 | CALL restput_p(rest_id, 'fluxlat', nbp_glo, 1, 1, kjit, fluxlat, 'scatter', nbp_glo, index_g) |
---|
| 553 | CALL restput_p(rest_id, 'fluxsens', nbp_glo, 1, 1, kjit, fluxsens, 'scatter', nbp_glo, index_g) |
---|
| 554 | CALL restput_p(rest_id, 'qsolpot', nbp_glo, 1, 1, kjit, q_sol_pot, 'scatter', nbp_glo, index_g) |
---|
| 555 | |
---|
| 556 | END SUBROUTINE enerbil_finalize |
---|
| 557 | |
---|
| 558 | |
---|
| 559 | |
---|
| 560 | |
---|
| 561 | !! ============================================================================================================================= |
---|
| 562 | !! SUBROUTINE : enerbil_clear |
---|
| 563 | !! |
---|
| 564 | !>\BRIEF Routine deallocates clear output variables if already allocated. |
---|
| 565 | !! |
---|
| 566 | !! DESCRIPTION : This is a 'housekeeping' routine that deallocates the key output |
---|
| 567 | !! variables, if they have already been allocated. The variables that are deallocated are psold, |
---|
| 568 | !! qsol_sat, pdqsold, psnew, qsol_sat_new, netrad, lwabs, lwup, lwnet, fluxsubli, qsat_air, tair |
---|
| 569 | !! |
---|
| 570 | !! RECENT CHANGE(S) : None |
---|
| 571 | !! |
---|
| 572 | !! MAIN OUTPUT VARIABLE(S) : None |
---|
| 573 | !! |
---|
| 574 | !! REFERENCES : None |
---|
| 575 | !! |
---|
| 576 | !! FLOWCHART : None |
---|
| 577 | !! \n |
---|
| 578 | !_ ============================================================================================================================== |
---|
| 579 | |
---|
| 580 | SUBROUTINE enerbil_clear () |
---|
| 581 | IF ( ALLOCATED (psold)) DEALLOCATE (psold) |
---|
| 582 | IF ( ALLOCATED (qsol_sat)) DEALLOCATE (qsol_sat) |
---|
| 583 | IF ( ALLOCATED (pdqsold)) DEALLOCATE (pdqsold) |
---|
| 584 | IF ( ALLOCATED (psnew)) DEALLOCATE (psnew) |
---|
| 585 | IF ( ALLOCATED (qsol_sat_new)) DEALLOCATE (qsol_sat_new) |
---|
| 586 | IF ( ALLOCATED (netrad)) DEALLOCATE (netrad) |
---|
| 587 | IF ( ALLOCATED (lwabs)) DEALLOCATE (lwabs) |
---|
| 588 | IF ( ALLOCATED (lwup)) DEALLOCATE (lwup) |
---|
| 589 | IF ( ALLOCATED (lwnet)) DEALLOCATE (lwnet) |
---|
| 590 | IF ( ALLOCATED (fluxsubli)) DEALLOCATE (fluxsubli) |
---|
| 591 | IF ( ALLOCATED (qsat_air)) DEALLOCATE (qsat_air) |
---|
| 592 | IF ( ALLOCATED (tair)) DEALLOCATE (tair) |
---|
| 593 | IF ( ALLOCATED (q_sol_pot)) DEALLOCATE (q_sol_pot) |
---|
| 594 | |
---|
| 595 | END SUBROUTINE enerbil_clear |
---|
| 596 | |
---|
| 597 | |
---|
| 598 | !! ============================================================================================================================= |
---|
| 599 | !! SUBROUTINE : enerbil_begin |
---|
| 600 | !! |
---|
| 601 | !>\BRIEF Preliminary variables required for the calculation of |
---|
| 602 | !! the energy budget are derived. |
---|
| 603 | !! |
---|
| 604 | !! DESCRIPTION : This routines computes preliminary variables required |
---|
| 605 | !! for the calculation of the energy budget: the old surface static energy (psold), the surface saturation |
---|
| 606 | !! humidity (qsol_sat), the derivative of satured specific humidity at the old temperature (pdqsold) |
---|
| 607 | !! and the net radiation (netrad). |
---|
| 608 | !! |
---|
| 609 | !! MAIN OUTPUT VARIABLE(S) : psold, qsol_sat, pdqsold, netrad |
---|
| 610 | !! |
---|
| 611 | !! REFERENCE(S) : |
---|
| 612 | !! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled |
---|
| 613 | !! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278 |
---|
| 614 | !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations |
---|
| 615 | !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General |
---|
| 616 | !! Circulation Model. Journal of Climate, 6, pp.248-273 |
---|
| 617 | !! - Dufresne, J-L & Ghattas, J, 2009. Description du schéma de la couche limite turbulente et la |
---|
| 618 | !! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11): |
---|
| 619 | !! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf |
---|
| 620 | !! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J |
---|
| 621 | !! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D & |
---|
| 622 | !! Xue, Y, 1998. A proposal for a general interface between land surface schemes and |
---|
| 623 | !! general circulation models. Global and Planetary Change, 19, pp.261-276 |
---|
| 624 | !! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems. |
---|
| 625 | !! Interscience Publishers\n |
---|
| 626 | !! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere |
---|
| 627 | !! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663. |
---|
| 628 | !! |
---|
| 629 | !! FLOWCHART : None |
---|
| 630 | !! \n |
---|
| 631 | !_ ============================================================================================================================== |
---|
| 632 | |
---|
| 633 | SUBROUTINE enerbil_begin (kjpindex, temp_sol, lwdown, swnet, pb, psold, qsol_sat, pdqsold, netrad, emis) |
---|
| 634 | |
---|
| 635 | !! 0. Variable and parameter declaration |
---|
| 636 | |
---|
| 637 | !! 0.1 Input variables |
---|
| 638 | |
---|
| 639 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
| 640 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Soil temperature (K) |
---|
| 641 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux (W m^{-2}) |
---|
| 642 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2}) |
---|
| 643 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa) |
---|
| 644 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-) |
---|
| 645 | |
---|
| 646 | !! 0.2 Output variables |
---|
| 647 | |
---|
| 648 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: psold !! Old surface dry static energy (J kg^{-1}) |
---|
| 649 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsol_sat !! Saturated specific humudity for old temperature |
---|
| 650 | !! (kg kg^{-1}) |
---|
| 651 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: pdqsold !! Derivative of satured specific humidity at the old |
---|
| 652 | !! temperature (kg (kg s)^{-1}) |
---|
| 653 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: netrad !! Net radiation (W m^{-2}) |
---|
| 654 | |
---|
| 655 | !! 0.3 Modified variables |
---|
| 656 | |
---|
| 657 | !! 0.4 Local variables |
---|
| 658 | |
---|
| 659 | INTEGER(i_std) :: ji |
---|
| 660 | REAL(r_std), DIMENSION(kjpindex) :: dev_qsol |
---|
| 661 | REAL(r_std), PARAMETER :: missing = 999998. |
---|
| 662 | |
---|
| 663 | !_ ================================================================================================================================ |
---|
| 664 | |
---|
| 665 | !! 1. Computes psold (the surface static energy for the old timestep) |
---|
| 666 | |
---|
| 667 | !! We here define the surface static energy for the 'old' timestep, in terms of the surface |
---|
| 668 | !! temperature and heat capacity. |
---|
| 669 | !! \latexonly |
---|
| 670 | !! \input{enerbilbegin1.tex} |
---|
| 671 | !! \endlatexonly |
---|
| 672 | psold(:) = temp_sol(:)*cp_air |
---|
| 673 | !! 2. Computes qsol_sat (the surface saturated humidity). |
---|
| 674 | |
---|
| 675 | !! We call the routine 'qsatcalc' from within the module 'src_parameters/constantes_veg'. |
---|
| 676 | CALL qsatcalc (kjpindex, temp_sol, pb, qsol_sat) |
---|
| 677 | IF ( diag_qsat ) THEN |
---|
| 678 | IF ( ANY(ABS(qsol_sat(:)) .GT. missing) ) THEN |
---|
| 679 | DO ji = 1, kjpindex |
---|
| 680 | IF ( ABS(qsol_sat(ji)) .GT. missing) THEN |
---|
| 681 | WRITE(numout,*) 'ERROR on ji = ', ji |
---|
| 682 | WRITE(numout,*) 'temp_sol(ji), pb(ji) :', temp_sol(ji), pb(ji) |
---|
| 683 | CALL ipslerr_p (3,'enerbil_begin', & |
---|
| 684 | & 'qsol too high ','','') |
---|
| 685 | ENDIF |
---|
| 686 | ENDDO |
---|
| 687 | ENDIF |
---|
| 688 | ENDIF |
---|
| 689 | |
---|
| 690 | !! 3. Computes pdqsold |
---|
| 691 | |
---|
| 692 | !! Computes pdqsold (the derivative of the saturated humidity with respect to temperature |
---|
| 693 | !! analysed at the surface temperature at the 'old' timestep. |
---|
| 694 | !! We call the routine 'dev_qsatcalc' from qsat_moisture module. |
---|
| 695 | CALL dev_qsatcalc (kjpindex, temp_sol, pb, dev_qsol) |
---|
| 696 | |
---|
| 697 | !! \latexonly |
---|
| 698 | !! \input{enerbilbegin2.tex} |
---|
| 699 | !! \endlatexonly |
---|
| 700 | pdqsold(:) = dev_qsol(:) |
---|
| 701 | |
---|
| 702 | |
---|
| 703 | IF ( diag_qsat ) THEN |
---|
| 704 | IF ( ANY(ABS( pdqsold(:)) .GT. missing) ) THEN |
---|
| 705 | DO ji = 1, kjpindex |
---|
| 706 | IF ( ABS( pdqsold(ji)) .GT. missing ) THEN |
---|
| 707 | WRITE(numout,*) 'ERROR on ji = ', ji |
---|
| 708 | WRITE(numout,*) 'temp_sol(ji), pb(ji) :', temp_sol(ji), pb(ji) |
---|
| 709 | CALL ipslerr_p (3,'enerbil_begin', & |
---|
| 710 | & 'pdqsold too high ','','') |
---|
| 711 | ENDIF |
---|
| 712 | ENDDO |
---|
| 713 | ENDIF |
---|
| 714 | ENDIF |
---|
| 715 | |
---|
| 716 | !! 4. Computes the net radiation and the absorbed LW radiation absorbed at the surface. |
---|
| 717 | |
---|
| 718 | !! Long wave radiation absorbed by the surface is the product of emissivity and downwelling LW radiation |
---|
| 719 | !! \latexonly |
---|
| 720 | !! \input{enerbilbegin3.tex} |
---|
| 721 | !! \endlatexonly |
---|
| 722 | lwabs(:) = emis(:) * lwdown(:) |
---|
| 723 | |
---|
| 724 | !! Net radiation is calculated as: |
---|
| 725 | !! \latexonly |
---|
| 726 | !! \input{enerbilbegin4.tex} |
---|
| 727 | !! \endlatexonly |
---|
| 728 | netrad(:) = lwdown(:) + swnet (:) - (emis(:) * c_stefan * temp_sol(:)**4 + (un - emis(:)) * lwdown(:)) |
---|
| 729 | |
---|
| 730 | IF (printlev>=3) WRITE (numout,*) ' enerbil_begin done ' |
---|
| 731 | |
---|
| 732 | END SUBROUTINE enerbil_begin |
---|
| 733 | |
---|
| 734 | |
---|
| 735 | !! ============================================================================================================================= |
---|
| 736 | !! SUBROUTINE : enerbil_surftemp |
---|
| 737 | !! |
---|
| 738 | !>\BRIEF This routine computes the new surface static energy |
---|
| 739 | !! (psnew) and the saturated humidity at the surface (qsol_sat_new). |
---|
| 740 | !! |
---|
| 741 | !! DESCRIPTION : This is the key part of the enerbil module, for which |
---|
| 742 | !! the energy budget in the surface layer is solved and changes over the model timestep 'dt_sechiba' are |
---|
| 743 | !! quantified for the surface static energy, surface temperature, surface humidity, the 'air' (or lowest |
---|
| 744 | !! level atmospheric model) temperature, the 'air' (or lowest level atmospheric model) humidity, |
---|
| 745 | !! according to the method that is laid out by Dufresne \& Ghattas (2009) and Shultz et al. (2001). |
---|
| 746 | !! |
---|
| 747 | !! It computes the energy balance at the surface with an implicit scheme, that is connected to the |
---|
| 748 | !! Richtmyer and Morton algorithm of the PBL. By computing the surface temperature and surface humidity |
---|
| 749 | !! the routine also implicitly estimates the various fluxes which balance in order to give the new |
---|
| 750 | !! temperature. Thus once the surface temperature has been obtained all the fluxes need to be diagnosed. |
---|
| 751 | !! |
---|
| 752 | !! Since the explicit snow scheme is adapted, the notion of skin layer in ORCHIDEE is abandoned and |
---|
| 753 | !! the first thin snow layer is adopted to solve the surface energy budget. |
---|
| 754 | !! Implicit method is still used for coupling to the atmosphere. |
---|
| 755 | !! In this new scheme, the snow temperature profile is simulatenously solved |
---|
| 756 | !! along with the surface snow temperature, and the details are referenced to Boone et al. (2010) |
---|
| 757 | !! |
---|
| 758 | !! MAIN OUTPUT VARIABLE(S) : psnew, qsol_sat_new, temp_sol_new, qair_new, epot_air_new |
---|
| 759 | !! |
---|
| 760 | !! REFERENCE(S) : |
---|
| 761 | !! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled |
---|
| 762 | !! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278 |
---|
| 763 | !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations |
---|
| 764 | !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General |
---|
| 765 | !! Circulation Model. Journal of Climate, 6, pp.248-273 |
---|
| 766 | !! - Dufresne, J-L & Ghattas, J, 2009. Description du schéma de la couche limite turbulente et la |
---|
| 767 | !! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11): |
---|
| 768 | !! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf |
---|
| 769 | !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation |
---|
| 770 | !! sur le cycle de l'eau, PhD Thesis, available (22/12/11): |
---|
| 771 | !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pd |
---|
| 772 | !! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J |
---|
| 773 | !! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D & |
---|
| 774 | !! Xue, Y, 1998. A proposal for a general interface between land surface schemes and |
---|
| 775 | !! general circulation models. Global and Planetary Change, 19, pp.261-276 |
---|
| 776 | !! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems. |
---|
| 777 | !! Interscience Publishers |
---|
| 778 | !! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere |
---|
| 779 | !! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663. |
---|
| 780 | !! |
---|
| 781 | !! FLOWCHART : None |
---|
| 782 | !! |
---|
| 783 | !_ ============================================================================================================================== |
---|
| 784 | |
---|
| 785 | |
---|
| 786 | SUBROUTINE enerbil_surftemp (kjpindex, zlev, emis, epot_air, & |
---|
| 787 | & petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,& |
---|
| 788 | & vbeta1, vbeta5, soilcap, lwdown, swnet, psnew, qsol_sat_new, temp_sol_new, & |
---|
| 789 | & qair_new, epot_air_new) |
---|
| 790 | |
---|
| 791 | |
---|
| 792 | |
---|
| 793 | !! 0. Variable and parameter declaration |
---|
| 794 | |
---|
| 795 | !! 0.1 Input variables |
---|
| 796 | |
---|
| 797 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
| 798 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m) |
---|
| 799 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-) |
---|
| 800 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy (?? J) |
---|
| 801 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! PetAcoef (see note) |
---|
| 802 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! PetBcoef (see note) |
---|
| 803 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1}) |
---|
| 804 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! PeqAcoef (see note) |
---|
| 805 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! PeqBcoef (see note) |
---|
| 806 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilflx !! Soil flux (W m^{-2}) |
---|
| 807 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air density (kg m^{-3}) |
---|
| 808 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u, v !! Wind velocity by directional components |
---|
| 809 | !! u (Eastwards) and v (Northwards) (m s^{-1}) |
---|
| 810 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) |
---|
| 811 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient (-) |
---|
| 812 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-) |
---|
| 813 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance |
---|
| 814 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil calorific capacity (J K^{-1]) |
---|
| 815 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux (W m^{-2}) |
---|
| 816 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2}) |
---|
| 817 | |
---|
| 818 | !! 0.2 Output variables |
---|
| 819 | |
---|
| 820 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: psnew !! New surface static energy (J kg^{-1}) |
---|
| 821 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsol_sat_new !! New saturated surface air moisture (kg kg^{-1}) |
---|
| 822 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature (K) |
---|
| 823 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qair_new !! New air moisture (kg kg^{-1}) |
---|
| 824 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: epot_air_new !! New air temperature (K) |
---|
| 825 | |
---|
| 826 | |
---|
| 827 | !! 0.4 Local variables |
---|
| 828 | |
---|
| 829 | INTEGER(i_std) :: ji |
---|
| 830 | REAL(r_std),DIMENSION (kjpindex) :: zicp |
---|
| 831 | REAL(r_std) :: fevap |
---|
| 832 | REAL(r_std) :: sensfl_old, larsub_old, lareva_old, dtheta, sum_old, sum_sns |
---|
| 833 | REAL(r_std) :: zikt, zikq, netrad_sns, sensfl_sns, larsub_sns, lareva_sns |
---|
| 834 | REAL(r_std) :: speed |
---|
| 835 | REAL(r_std),DIMENSION (kjpindex) :: dtes !! Diagnostic output variable for the change of the thermal energy content of surface for which the energy balance is calculated (J m-2) |
---|
| 836 | |
---|
| 837 | !_ ================================================================================================================================ |
---|
| 838 | |
---|
| 839 | zicp = un / cp_air |
---|
| 840 | |
---|
| 841 | DO ji=1,kjpindex |
---|
| 842 | !! 1. Derivation of auxiliary variables |
---|
| 843 | |
---|
| 844 | !! \latexonly |
---|
| 845 | !! \input{enerbilsurftemp1.tex} |
---|
| 846 | !! \endlatexonly |
---|
| 847 | speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) |
---|
| 848 | ! |
---|
| 849 | !! \latexonly |
---|
| 850 | !! \input{enerbilsurftemp2.tex} |
---|
| 851 | !! \endlatexonly |
---|
| 852 | zikt = 1/(rau(ji) * speed * q_cdrag(ji)) |
---|
| 853 | zikq = 1/(rau(ji) * speed * q_cdrag(ji)) |
---|
| 854 | |
---|
| 855 | !! 2. Computation of fluxes for the old surface conditions |
---|
| 856 | |
---|
| 857 | !! As a reference, we first calculate the sensible and latent heat for the 'old' timestep. |
---|
| 858 | |
---|
| 859 | !! 2.1 Sensible heat for the 'old' timestep |
---|
| 860 | !! This is equation (64) of (Dufresne & Ghattas, 2009), rewritten in terms of H_old. We make the |
---|
| 861 | !! approximation that (P_r/P)^{kappa}=1 and convert the A and B coefficients to the ORCHIDEE |
---|
| 862 | !! format (see introductory note, above). Also, the equation here is in terms of surface static |
---|
| 863 | !! energy, which is defined as per section 3.1, above. |
---|
| 864 | !! \latexonly |
---|
| 865 | !! \input{enerbilsurftemp3.tex} |
---|
| 866 | !! \endlatexonly |
---|
| 867 | sensfl_old = (petBcoef(ji) - psold(ji)) / (zikt - petAcoef(ji)) |
---|
| 868 | !! 2.2 Latent heat for the 'old' timestep |
---|
| 869 | !! This is equation (70) of (Dufresne & Ghattas, 2009), rewritten in terms of {\lambda E}_{old}. |
---|
| 870 | !! Again we convert the A and B coefficients from the LMDZ format to the ORCHIDEE format.\n |
---|
| 871 | !! There are two forms of the equation - first for the latent heat as a result of sublimination |
---|
| 872 | !! processes and then as a result of evaporative processes: |
---|
| 873 | !! \latexonly |
---|
| 874 | !! \input{enerbilsurftemp4.tex} |
---|
| 875 | !! \endlatexonly |
---|
| 876 | larsub_old = chalsu0 * vbeta1(ji) * (un - vbeta5(ji)) * & |
---|
| 877 | (peqBcoef(ji) - qsol_sat(ji)) / (zikq - vbeta1(ji) * (un - vbeta5(ji))* peqAcoef(ji)) |
---|
| 878 | !! \latexonly |
---|
| 879 | !! \input{enerbilsurftemp5.tex} |
---|
| 880 | !! \endlatexonly |
---|
| 881 | lareva_old = chalev0 * (un - vbeta1(ji)) * (un - vbeta5(ji)) * vbeta(ji) * & |
---|
| 882 | (peqBcoef(ji) - qsol_sat(ji)) / & |
---|
| 883 | (zikq - (un - vbeta1(ji)) * (un - vbeta5(ji)) * vbeta(ji) * peqAcoef(ji)) & |
---|
| 884 | + chalev0 * vbeta5(ji) * (peqBcoef(ji) - qsol_sat(ji)) / (zikq - vbeta5(ji) * peqAcoef(ji)) |
---|
| 885 | |
---|
| 886 | !! 3. Calculation of sensitivity terms |
---|
| 887 | |
---|
| 888 | !! We here calculate the sensitivity for the different fluxes to changes in dtheta, which is the |
---|
| 889 | !! change in the surface static energy over the model time step (dt_sechiba). |
---|
| 890 | |
---|
| 891 | !! 3.1 Net radiation sensitivity |
---|
| 892 | !! This is an explicit-implicit representation of the Stefan-Boltzmann law - the explicit terms |
---|
| 893 | !! are ${ps}_{old}^3$, and it is completed when multiplied by dtheta (defined above). |
---|
| 894 | !! \latexonly |
---|
| 895 | !! \input{enerbilsurftemp6.tex} |
---|
| 896 | !! \endlatexonly |
---|
| 897 | netrad_sns = zicp(ji) * quatre * emis(ji) * c_stefan * ((zicp(ji) * psold(ji))**3) |
---|
| 898 | |
---|
| 899 | !! 3.2 Sensible heat flux sensitivity |
---|
| 900 | !! This is the temperature sensitivity term N_1^h derived in equation (66) of (Dufresne & Ghattas, 2009), |
---|
| 901 | !! where we again assume that (P_r/P)^{kappa}=1, convert the A and B coefficients to the ORCHIDEE format, |
---|
| 902 | !! and rewrite in terms of surface static energy, rather than temperature (see section 3.1, above). |
---|
| 903 | !! \latexonly |
---|
| 904 | !! \input{enerbilsurftemp7.tex} |
---|
| 905 | !! \endlatexonly |
---|
| 906 | sensfl_sns = un / (zikt - petAcoef(ji)) |
---|
| 907 | !! 3.3 Latent heat flux sensitivity |
---|
| 908 | !! This is the humidity sensitivity term N_1^q derived in equation (72) of (Dufresne & Ghattas, 2009), where |
---|
| 909 | !! the A and B coefficients are written in the ORCHIDEE format. |
---|
| 910 | !! larsub_sns is the latent heat sensitivity for sublimination. The coefficient vbeta1 is the evaporation |
---|
| 911 | !! coefficient. It represents the relationship between the real and potential evaporation. It is derived |
---|
| 912 | !! in the module 'src_sechiba/diffuco_snow'. |
---|
| 913 | !! \latexonly |
---|
| 914 | !! \input{enerbilsurftemp8.tex} |
---|
| 915 | !! \endlatexonly |
---|
| 916 | larsub_sns = chalsu0 * vbeta1(ji) * (un - vbeta5(ji)) * zicp(ji) * & |
---|
| 917 | pdqsold(ji) / (zikq - vbeta1(ji) * (un - vbeta5(ji)) * peqAcoef(ji)) |
---|
| 918 | |
---|
| 919 | !! lareva_sns is the latent heat sensitivity for evaporation. vbeta1 (src_sechiba/diffuco_snow), |
---|
| 920 | !! and vbeta (src_sechiba/diffuco_comb) are the evaporation |
---|
| 921 | !! coefficients. |
---|
| 922 | !! \latexonly |
---|
| 923 | !! \input{enerbilsurftemp9.tex} |
---|
| 924 | !! \endlatexonly |
---|
| 925 | lareva_sns = chalev0 * ((un - vbeta1(ji))*(un - vbeta5(ji)) * vbeta(ji) + vbeta5(ji)) * & |
---|
| 926 | & zicp(ji) * pdqsold(ji) / & |
---|
| 927 | (zikq - ((un - vbeta1(ji))*(un - vbeta5(ji)) * vbeta(ji) + vbeta5(ji))* peqAcoef(ji)) |
---|
| 928 | |
---|
| 929 | !! 4. Solution of the energy balance |
---|
| 930 | |
---|
| 931 | !! 4.1 We calculate the total flux for the 'old' timestep. |
---|
| 932 | !! \latexonly |
---|
| 933 | !! \input{enerbilsurftemp10.tex} |
---|
| 934 | !! \endlatexonly |
---|
| 935 | sum_old = netrad(ji) + sensfl_old + larsub_old + lareva_old + soilflx(ji) |
---|
| 936 | |
---|
| 937 | !! 4.2 We calculate the total sensitivity dtheta (the change in the |
---|
| 938 | !! surface static energy over the timestep. |
---|
| 939 | !! \latexonly |
---|
| 940 | !! \input{enerbilsurftemp11.tex} |
---|
| 941 | !! \endlatexonly |
---|
| 942 | sum_sns = netrad_sns + sensfl_sns + larsub_sns + lareva_sns |
---|
| 943 | !! 4.3 Calculation of dtheta (change in surface static energy over the |
---|
| 944 | !! timestep. |
---|
| 945 | !! \latexonly |
---|
| 946 | !! \input{enerbilsurftemp12.tex} |
---|
| 947 | !! \endlatexonly |
---|
| 948 | |
---|
| 949 | dtheta = dt_sechiba * sum_old / (zicp(ji) * soilcap(ji) + dt_sechiba * sum_sns) |
---|
| 950 | |
---|
| 951 | !! 4.4 Determination of state variables for the 'new' timestep |
---|
| 952 | !! No we have dtheta, we can determine the surface static energy that |
---|
| 953 | !! corresponds to the 'new' timestep. |
---|
| 954 | !! \latexonly |
---|
| 955 | !! \input{enerbilsurftemp13.tex} |
---|
| 956 | !! \endlatexonly |
---|
| 957 | psnew(ji) = psold(ji) + dtheta |
---|
| 958 | |
---|
| 959 | !! The new surface saturated humidity can be calculated by equation (69) |
---|
| 960 | !! of (Dufresne & Ghattas, 2009), in which we substitute dtheta for the |
---|
| 961 | !! change between old and new temperature using the relationship from 3.1, |
---|
| 962 | !! above. |
---|
| 963 | !! \latexonly |
---|
| 964 | !! \input{enerbilsurftemp14.tex} |
---|
| 965 | !! \endlatexonly |
---|
| 966 | qsol_sat_new(ji) = qsol_sat(ji) + zicp(ji) * pdqsold(ji) * dtheta |
---|
| 967 | !! The new surface temperature is determined from the new surface static temperature, |
---|
| 968 | !! again by using the relationship in section 3.1. |
---|
| 969 | !! \latexonly |
---|
| 970 | !! \input{enerbilsurftemp15.tex} |
---|
| 971 | !! \endlatexonly |
---|
| 972 | temp_sol_new(ji) = psnew(ji) / cp_air |
---|
| 973 | |
---|
| 974 | !! 4.5 Calculation of new evaporation potential and new evaporation latent heat |
---|
| 975 | !! flux (???) |
---|
| 976 | !! \latexonly |
---|
| 977 | !! \input{enerbilsurftemp16.tex} |
---|
| 978 | !! \endlatexonly |
---|
| 979 | epot_air_new(ji) = zikt * (sensfl_old - sensfl_sns * dtheta) + psnew(ji) |
---|
| 980 | |
---|
| 981 | !! \latexonly |
---|
| 982 | !! \input{enerbilsurftemp17.tex} |
---|
| 983 | !! \endlatexonly |
---|
| 984 | fevap = (lareva_old - lareva_sns * dtheta) + (larsub_old - larsub_sns * dtheta) |
---|
| 985 | |
---|
| 986 | IF ( ABS(fevap) < EPSILON(un) ) THEN |
---|
| 987 | |
---|
| 988 | !! \latexonly |
---|
| 989 | !! \input{enerbilsurftemp18.tex} |
---|
| 990 | !! \endlatexonly |
---|
| 991 | qair_new(ji) = qair(ji) |
---|
| 992 | ELSE |
---|
| 993 | !! \latexonly |
---|
| 994 | !! \input{enerbilsurftemp19.tex} |
---|
| 995 | !! \endlatexonly |
---|
| 996 | qair_new(ji) = zikq * un / ( chalsu0 * vbeta1(ji) * (un - vbeta5(ji)) + & |
---|
| 997 | & chalev0 * ((un - vbeta1(ji))*(un - vbeta5(ji)) * vbeta(ji) + vbeta5(ji)) ) & |
---|
| 998 | & * fevap + qsol_sat_new(ji) |
---|
| 999 | ENDIF |
---|
| 1000 | |
---|
| 1001 | ! Calculate dtes only for diagnostic output and send it to xios |
---|
| 1002 | dtes(ji) = (soilcap(ji)/cp_air)*dtheta/dt_sechiba |
---|
| 1003 | |
---|
| 1004 | ENDDO |
---|
| 1005 | |
---|
| 1006 | ! Send diagnostic variable to XIOS |
---|
| 1007 | CALL xios_orchidee_send_field("dtes", dtes) |
---|
| 1008 | |
---|
| 1009 | IF (printlev>=3) WRITE (numout,*) ' enerbil_surftemp done ' |
---|
| 1010 | |
---|
| 1011 | END SUBROUTINE enerbil_surftemp |
---|
| 1012 | |
---|
| 1013 | |
---|
| 1014 | !! ============================================================================================================================= |
---|
| 1015 | !! SUBROUTINE : enerbil_pottemp |
---|
| 1016 | !! |
---|
| 1017 | !>\BRIEF This subroutine computes the surface temperature and humidity should the surface been unstressed |
---|
| 1018 | !! |
---|
| 1019 | !! DESCRIPTION : |
---|
| 1020 | !! |
---|
| 1021 | !! MAIN OUTPUT VARIABLE(S) : |
---|
| 1022 | !! |
---|
| 1023 | !! REFERENCE(S) : |
---|
| 1024 | !! |
---|
| 1025 | !! FLOWCHART : None |
---|
| 1026 | !! |
---|
| 1027 | !_ ============================================================================================================================== |
---|
| 1028 | |
---|
| 1029 | SUBROUTINE enerbil_pottemp (kjpindex, zlev, emis, epot_air, & |
---|
| 1030 | & petAcoef, petBcoef, qair, peqAcoef, peqBcoef, soilflx, rau, u, v, q_cdrag, vbeta,& |
---|
| 1031 | & vbeta1, vbeta5, soilcap, lwdown, swnet, q_sol_pot, temp_sol_pot) |
---|
| 1032 | |
---|
| 1033 | ! interface |
---|
| 1034 | ! input scalar |
---|
| 1035 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size |
---|
| 1036 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer |
---|
| 1037 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity |
---|
| 1038 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy |
---|
| 1039 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! PetAcoef |
---|
| 1040 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! PetBcoef |
---|
| 1041 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity |
---|
| 1042 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! PeqAcoef |
---|
| 1043 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! PeqBcoef |
---|
| 1044 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilflx !! Soil flux |
---|
| 1045 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density |
---|
| 1046 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u, v !! Wind |
---|
| 1047 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) |
---|
| 1048 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient |
---|
| 1049 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance |
---|
| 1050 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance |
---|
| 1051 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil calorific capacity |
---|
| 1052 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux |
---|
| 1053 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux |
---|
| 1054 | ! output fields |
---|
| 1055 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: q_sol_pot !! Potential surface air moisture |
---|
| 1056 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_pot !! Potential soil temperature |
---|
| 1057 | |
---|
| 1058 | |
---|
| 1059 | ! local declaration |
---|
| 1060 | INTEGER(i_std) :: ji |
---|
| 1061 | REAL(r_std),DIMENSION (kjpindex) :: zicp |
---|
| 1062 | REAL(r_std) :: fevap |
---|
| 1063 | REAL(r_std) :: sensfl_old, larsub_old, lareva_old, dtheta, sum_old, sum_sns |
---|
| 1064 | REAL(r_std) :: zikt, zikq, netrad_sns, sensfl_sns, larsub_sns, lareva_sns |
---|
| 1065 | REAL(r_std) :: speed |
---|
| 1066 | |
---|
| 1067 | !_ ============================================================================================================================== |
---|
| 1068 | |
---|
| 1069 | zicp = un / cp_air |
---|
| 1070 | ! |
---|
| 1071 | DO ji=1,kjpindex |
---|
| 1072 | |
---|
| 1073 | dtheta = zero |
---|
| 1074 | fevap = zero |
---|
| 1075 | |
---|
| 1076 | temp_sol_pot(ji) = temp_sol_pot(ji) + dtheta |
---|
| 1077 | |
---|
| 1078 | q_sol_pot(ji) = q_sol_pot(ji) + fevap |
---|
| 1079 | |
---|
| 1080 | ENDDO |
---|
| 1081 | |
---|
| 1082 | END SUBROUTINE enerbil_pottemp |
---|
| 1083 | |
---|
| 1084 | |
---|
| 1085 | !! ============================================================================================================================= |
---|
| 1086 | !! SUBROUTINE : enerbil_flux |
---|
| 1087 | !! |
---|
| 1088 | !>\BRIEF Computes the new soil temperature, net radiation and |
---|
| 1089 | !! latent and sensible heat flux for the new time step. |
---|
| 1090 | !! |
---|
| 1091 | !! DESCRIPTION : This routine diagnoses, based on the new soil temperature, |
---|
| 1092 | !! the net radiation, the total evaporation, the latent heat flux, the sensible heat flux and the |
---|
| 1093 | !! sublimination flux. It also diagnoses the potential evaporation used for the fluxes (evapot) and the potential |
---|
| 1094 | !! as defined by Penman & Monteith (Monteith, 1965) based on the correction term developed by Chris |
---|
| 1095 | !! Milly (1992). This Penman-Monteith formulation is required for the estimation of bare soil evaporation |
---|
| 1096 | !! when the 11 layer CWRR moisture scheme is used for the soil hydrology. |
---|
| 1097 | !! |
---|
| 1098 | !! MAIN OUTPUT VARIABLE(S) : qsurf, fluxsens, fluxlat, fluxsubli, vevapp, lwup, lwnet, |
---|
| 1099 | !! tsol_rad, netrad, evapot, evapot_corr |
---|
| 1100 | !! |
---|
| 1101 | !! REFERENCE(S) : |
---|
| 1102 | !! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled |
---|
| 1103 | !! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278 |
---|
| 1104 | !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations |
---|
| 1105 | !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General |
---|
| 1106 | !! Circulation Model. Journal of Climate, 6, pp.248-273 |
---|
| 1107 | !! - Dufresne, J-L & Ghattas, J, 2009. Description du schéma de la couche limite turbulente et la |
---|
| 1108 | !! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11): |
---|
| 1109 | !! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf |
---|
| 1110 | !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation |
---|
| 1111 | !! sur le cycle de l'eau, PhD Thesis, available (22/12/11): |
---|
| 1112 | !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf |
---|
| 1113 | !! - Monteith, JL, 1965. Evaporation and Environment, paper presented at Symposium of the Society |
---|
| 1114 | !! for Experimental Biology |
---|
| 1115 | !! - Monteith & Unsworth, 2008. Principles of Environmental Physics (third edition), published Elsevier |
---|
| 1116 | !! ISBN 978-0-12-505103-3 |
---|
| 1117 | !! - Milly, P. C. D., 1992: Potential Evaporation and Soil Moisture in General Circulation Models. |
---|
| 1118 | !! Journal of Climate, 5, pp. 209â226. |
---|
| 1119 | !! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J |
---|
| 1120 | !! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D & |
---|
| 1121 | !! Xue, Y, 1998. A proposal for a general interface between land surface schemes and |
---|
| 1122 | !! general circulation models. Global and Planetary Change, 19, pp.261-276 |
---|
| 1123 | !! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems. |
---|
| 1124 | !! Interscience Publishers |
---|
| 1125 | !! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere |
---|
| 1126 | !! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663. |
---|
| 1127 | !! |
---|
| 1128 | !! FLOWCHART : None |
---|
| 1129 | !! \n |
---|
| 1130 | !_ ============================================================================================================================== |
---|
| 1131 | |
---|
| 1132 | SUBROUTINE enerbil_flux (kjpindex, emis, temp_sol, rau, u, v, q_cdrag, vbeta, vbeta1, vbeta5, & |
---|
| 1133 | & qair, epot_air, psnew, qsurf, fluxsens, fluxlat, fluxsubli, vevapp, temp_sol_new, & |
---|
| 1134 | & lwdown, swnet, lwup, lwnet, pb, tsol_rad, netrad, evapot, evapot_corr,& |
---|
| 1135 | & precip_rain,snowdz,temp_air,pgflux, soilcap, temp_sol_add) |
---|
| 1136 | |
---|
| 1137 | !! 0. Variable and parameter declaration |
---|
| 1138 | |
---|
| 1139 | !! 0.1 Input variables |
---|
| 1140 | |
---|
| 1141 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
| 1142 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: emis !! Emissivity (-) |
---|
| 1143 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Surface temperature (K) |
---|
| 1144 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density (kg m^{-3}) |
---|
| 1145 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u,v !! Wind velocity in components u and v (m s^{-1}) |
---|
| 1146 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) |
---|
| 1147 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta !! Resistance coefficient (-) |
---|
| 1148 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-) |
---|
| 1149 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Flood resistance |
---|
| 1150 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1}) |
---|
| 1151 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy (J) |
---|
| 1152 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: psnew !! New surface static energy (J kg^{-1}) |
---|
| 1153 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature (K) |
---|
| 1154 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa) |
---|
| 1155 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth |
---|
| 1156 | REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: precip_rain !! Rainfall |
---|
| 1157 | REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: temp_air !! Air temperature |
---|
| 1158 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Downward long wave radiation (W m^{-2}) |
---|
| 1159 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net short wave radiation (W m^{-2}) |
---|
| 1160 | REAL(r_std), DIMENSION (kjpindex),INTENT(in) :: soilcap !! Soil calorific capacity including snow and soil (J K^{-1}) |
---|
| 1161 | |
---|
| 1162 | !! 0.2 Output variables |
---|
| 1163 | |
---|
| 1164 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf !! Surface specific humidity (kg kg^{-1}) |
---|
| 1165 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible heat flux (W m^{-2}) |
---|
| 1166 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent heat flux (W m^{-2}) |
---|
| 1167 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsubli !! Energy of sublimation (W m^{-2}) |
---|
| 1168 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation (mm day^{-1}) |
---|
| 1169 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: lwup !! Long-wave up radiation (W m^{-2}) |
---|
| 1170 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: lwnet !! Long-wave net radiation (W m^{-2}) |
---|
| 1171 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Radiative surface temperature (W m^{-2}) |
---|
| 1172 | REAL(r_std), DIMENSION (kjpindex),INTENT(out) :: temp_sol_add !! Additional energy to melt snow for snow ablation case (K) |
---|
| 1173 | |
---|
| 1174 | |
---|
| 1175 | !! 0.3 Modified variables |
---|
| 1176 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: netrad !! Net radiation (W m^{-2}) |
---|
| 1177 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot !! Soil Potential Evaporation (mm/tstep) |
---|
| 1178 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evapot_corr !! Soil Potential Evaporation Correction (mm/tstep) |
---|
| 1179 | REAL(r_std), DIMENSION (kjpindex),INTENT(inout) :: pgflux !! Net energy into snowpack(W/m^2) |
---|
| 1180 | |
---|
| 1181 | !! 0.4 Local variables |
---|
| 1182 | |
---|
| 1183 | INTEGER(i_std) :: ji |
---|
| 1184 | REAL(r_std),DIMENSION (kjpindex) :: grad_qsat |
---|
| 1185 | REAL(r_std) :: correction |
---|
| 1186 | REAL(r_std) :: speed !! Speed (m s^{-1}), |
---|
| 1187 | REAL(r_std) :: qc !! Surface drag coefficient (??) |
---|
| 1188 | LOGICAL,DIMENSION (kjpindex) :: warning_correction |
---|
| 1189 | REAL(r_std), DIMENSION (kjpindex) :: zgflux |
---|
| 1190 | REAL(r_std), DIMENSION (kjpindex) :: qsol_sat_tmp |
---|
| 1191 | REAL(r_std), DIMENSION (kjpindex) :: zerocelcius |
---|
| 1192 | REAL(r_std), DIMENSION (kjpindex) :: PHPSNOW |
---|
| 1193 | REAL(r_std), DIMENSION (kjpindex) :: lwup_tmp |
---|
| 1194 | REAL(r_std), DIMENSION (kjpindex) :: netrad_tmp |
---|
| 1195 | REAL(r_std), DIMENSION (kjpindex) :: fluxsens_tmp |
---|
| 1196 | REAL(r_std), DIMENSION (kjpindex) :: fluxlat_tmp |
---|
| 1197 | REAL(r_std), DIMENSION (kjpindex) :: betatot !! beta used in the computation of qsurf (-) |
---|
| 1198 | !_ ================================================================================================================================ |
---|
| 1199 | |
---|
| 1200 | zerocelcius(:) = tp_00 |
---|
| 1201 | CALL qsatcalc (kjpindex, zerocelcius, pb, qsol_sat_tmp) |
---|
| 1202 | |
---|
| 1203 | DO ji=1,kjpindex |
---|
| 1204 | !! 1. Determination of 'housekeeping' variables |
---|
| 1205 | |
---|
| 1206 | !! The horizontal wind speed is calculated from the velocities along each axis using Pythagorus' theorem. |
---|
| 1207 | !! \latexonly |
---|
| 1208 | !! \input{enerbilflux1.tex} |
---|
| 1209 | !! \endlatexonly |
---|
| 1210 | speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) |
---|
| 1211 | |
---|
| 1212 | !! From definition, the surface drag coefficient is defined: |
---|
| 1213 | !! \latexonly |
---|
| 1214 | !! \input{enerbilflux2.tex} |
---|
| 1215 | !! \endlatexonly |
---|
| 1216 | qc = speed * q_cdrag(ji) |
---|
| 1217 | |
---|
| 1218 | !! 2. Calculation of the net upward flux of longwave radiation |
---|
| 1219 | |
---|
| 1220 | !! We first of all calculate the radiation as a result of the Stefan-Boltzmann equation, |
---|
| 1221 | !! which is the sum of the calculated values at the surface temperature at the 'old' |
---|
| 1222 | !! temperature and the value that corresponds to the difference between the old temperature |
---|
| 1223 | !! and the temperature at the 'new' timestep. |
---|
| 1224 | !! \latexonly |
---|
| 1225 | !! \input{enerbilflux3.tex} |
---|
| 1226 | !! \endlatexonly |
---|
| 1227 | lwup(ji) = emis(ji) * c_stefan * temp_sol(ji)**4 + & |
---|
| 1228 | & quatre * emis(ji) * c_stefan * temp_sol(ji)**3 * & |
---|
| 1229 | & (temp_sol_new(ji) - temp_sol(ji)) |
---|
| 1230 | |
---|
| 1231 | !! We then add to this value that from the reflected longwave radiation: |
---|
| 1232 | !! \latexonly |
---|
| 1233 | !! \input{enerbilflux4.tex} |
---|
| 1234 | !! \endlatexonly |
---|
| 1235 | lwup(ji) = lwup(ji) + (un - emis(ji)) * lwdown(ji) |
---|
| 1236 | |
---|
| 1237 | !! The radiative surface temperature is calculated by inverting the Stefan-Boltzmann relation: |
---|
| 1238 | !! \latexonly |
---|
| 1239 | !! \input{enerbilflux5.tex} |
---|
| 1240 | !! \endlatexonly |
---|
| 1241 | !! The implicit solution computes an emitted long-wave flux which is a limited Taylor expansion |
---|
| 1242 | !! of the future surface temperature around the current values. Thus the long-wave flux does not |
---|
| 1243 | !! correspond to the new surface temperature but some intermediate value. So we need to deduce the |
---|
| 1244 | !! radiative surface temperature to which this flux corresponds. |
---|
| 1245 | !! |
---|
| 1246 | tsol_rad(ji) = (lwup(ji)/ (emis(ji) * c_stefan)) **(1./quatre) |
---|
| 1247 | |
---|
| 1248 | !! qsurf (the surface specific humidity) is a simple diagnostic which will be used by the |
---|
| 1249 | !! GCM to compute the dependence of of the surface layer stability on moisture. |
---|
| 1250 | !! \latexonly |
---|
| 1251 | !! \input{enerbilflux6.tex} |
---|
| 1252 | !! \endlatexonly |
---|
| 1253 | |
---|
| 1254 | betatot(ji) = vbeta1(ji) * (un - vbeta5(ji)) + vbeta5(ji) + & |
---|
| 1255 | & (un - vbeta1(ji))*(un - vbeta5(ji)) * vbeta(ji) |
---|
| 1256 | |
---|
| 1257 | qsurf(ji) = betatot(ji) * qsol_sat_new(ji) + (un - betatot(ji)) * qair(ji) |
---|
| 1258 | |
---|
| 1259 | ! When no turbulence and no wind (qc small) the near surface air moisture is equal to qair. |
---|
| 1260 | ! This is related to problem with interpolation in stdlevvar in LMDZ. |
---|
| 1261 | IF ( qc .LE. min_qc ) qsurf(ji) = qair(ji) |
---|
| 1262 | |
---|
| 1263 | |
---|
| 1264 | !! \latexonly |
---|
| 1265 | !! \input{enerbilflux7.tex} |
---|
| 1266 | !! \endlatexonly |
---|
| 1267 | qsurf(ji) = MAX(qsurf(ji), qair(ji)) |
---|
| 1268 | |
---|
| 1269 | !! Net downward radiation is the sum of the down-welling less the up-welling long wave flux, plus |
---|
| 1270 | !! the short wave radiation. |
---|
| 1271 | !! \latexonly |
---|
| 1272 | !! \emissivity absorbed lw radiationinput{enerbilflux8.tex} |
---|
| 1273 | !! \endlatexonly |
---|
| 1274 | netrad(ji) = lwdown(ji) + swnet(ji) - lwup(ji) |
---|
| 1275 | |
---|
| 1276 | !! 'vevapp' is the sum of the total evaporative processes (snow plus non-snow processes). |
---|
| 1277 | !! \latexonly |
---|
| 1278 | !! \input{enerbilflux9.tex} |
---|
| 1279 | !! \endlatexonly |
---|
| 1280 | vevapp(ji) = dt_sechiba * rau(ji) * qc * (vbeta1(ji) * (un - vbeta5(ji)) + vbeta5(ji)) * & |
---|
| 1281 | & (qsol_sat_new(ji) - qair(ji)) + & |
---|
| 1282 | & dt_sechiba * rau(ji) * qc * (un - vbeta1(ji))*(un-vbeta5(ji)) * vbeta(ji) * & |
---|
| 1283 | & (qsol_sat_new(ji) - qair(ji)) |
---|
| 1284 | |
---|
| 1285 | !! The total latent heat flux is the sum of the snow plus non-snow processes. |
---|
| 1286 | !! \latexonly |
---|
| 1287 | !! \input{enerbilflux10.tex} |
---|
| 1288 | !! \endlatexonly |
---|
| 1289 | fluxlat(ji) = chalsu0 * rau(ji) * qc * vbeta1(ji) * (un - vbeta5(ji)) * & |
---|
| 1290 | & (qsol_sat_new(ji) - qair(ji)) + & |
---|
| 1291 | & chalev0 * rau(ji) * qc * vbeta5(ji) *& |
---|
| 1292 | & (qsol_sat_new(ji) - qair(ji)) + & |
---|
| 1293 | & chalev0 * rau(ji) * qc * (un - vbeta1(ji)) * (un - vbeta5(ji)) * vbeta(ji) * & |
---|
| 1294 | & (qsol_sat_new(ji) - qair(ji)) |
---|
| 1295 | |
---|
| 1296 | !! The sublimination flux concerns is calculated using vbeta1, the snow resistance. |
---|
| 1297 | !! \latexonly |
---|
| 1298 | !! \input{enerbilflux11.tex} |
---|
| 1299 | !! \endlatexonly |
---|
| 1300 | fluxsubli(ji) = chalsu0 * rau(ji) * qc * vbeta1(ji) * (un - vbeta5(ji)) * & |
---|
| 1301 | & (qsol_sat_new(ji) - qair(ji)) |
---|
| 1302 | |
---|
| 1303 | !! The sensible heat flux is a factor of the difference between the new surface static energy |
---|
| 1304 | !! and the potential energy of air. |
---|
| 1305 | !! \latexonly |
---|
| 1306 | !! \input{enerbilflux12.tex} |
---|
| 1307 | !! \endlatexonly |
---|
| 1308 | fluxsens(ji) = rau(ji) * qc * (psnew(ji) - epot_air(ji)) |
---|
| 1309 | |
---|
| 1310 | !! This is the net longwave downwards radiation. |
---|
| 1311 | !! \latexonly |
---|
| 1312 | !! \input{enerbilflux13.tex} |
---|
| 1313 | !! \endlatexonly |
---|
| 1314 | lwnet(ji) = lwdown(ji) - lwup(ji) |
---|
| 1315 | |
---|
| 1316 | !! Diagnoses the potential evaporation used for the fluxes (evapot) |
---|
| 1317 | !! \latexonly |
---|
| 1318 | !! \input{enerbilflux14.tex} |
---|
| 1319 | !! \endlatexonly |
---|
| 1320 | evapot(ji) = MAX(zero, dt_sechiba * rau(ji) * qc * (qsol_sat_new(ji) - qair(ji))) |
---|
| 1321 | |
---|
| 1322 | !! From definition we can say: |
---|
| 1323 | !! \latexonly |
---|
| 1324 | !! \input{enerbilflux15.tex} |
---|
| 1325 | !! \endlatexonly |
---|
| 1326 | tair(ji) = epot_air(ji) / cp_air |
---|
| 1327 | |
---|
| 1328 | |
---|
| 1329 | !! To calculate net energy flux into the snowpack |
---|
| 1330 | PHPSNOW(ji) = precip_rain(ji)*(4.218E+3)*(MAX(tp_00,temp_air(ji))-tp_00)/dt_sechiba ! (w/m2) |
---|
| 1331 | pgflux(ji) = netrad(ji) - fluxsens(ji) - fluxlat(ji) + PHPSNOW(ji) |
---|
| 1332 | |
---|
| 1333 | |
---|
| 1334 | !! To get the extra energy used to melt the snowpack |
---|
| 1335 | IF (temp_sol_new (ji) > tp_00 .AND. & |
---|
| 1336 | SUM(snowdz(ji,:)) .GT. zero .AND. soilcap(ji) .GT. min_sechiba) THEN |
---|
| 1337 | |
---|
| 1338 | lwup_tmp(ji) = emis(ji) * c_stefan * temp_sol(ji)**4 + & |
---|
| 1339 | quatre * emis(ji) * c_stefan * temp_sol(ji)**3 * (tp_00 - temp_sol(ji)) |
---|
| 1340 | |
---|
| 1341 | lwup_tmp(ji) = lwup_tmp(ji) + (un - emis(ji)) * lwdown(ji) |
---|
| 1342 | netrad_tmp(ji) = lwdown(ji) + swnet(ji) - lwup_tmp(ji) |
---|
| 1343 | fluxsens_tmp(ji) = rau(ji) * qc * cp_air * (tp_00 - epot_air(ji)/cp_air) |
---|
| 1344 | fluxlat_tmp(ji) = chalsu0 * rau(ji) * qc * vbeta1(ji) * (un-vbeta5(ji)) * & |
---|
| 1345 | & (qsol_sat_tmp(ji) - qair(ji)) + & |
---|
| 1346 | & chalev0 * rau(ji) * qc * vbeta5(ji) *& |
---|
| 1347 | & (qsol_sat_tmp(ji) - qair(ji)) + & |
---|
| 1348 | & chalev0 * rau(ji) * qc * (un - vbeta1(ji)) * (un-vbeta5(ji))* vbeta(ji) * & |
---|
| 1349 | & (qsol_sat_tmp(ji) - qair(ji)) |
---|
| 1350 | |
---|
| 1351 | zgflux(ji) = netrad_tmp(ji) - fluxsens_tmp(ji) - fluxlat_tmp(ji)+PHPSNOW(ji) |
---|
| 1352 | |
---|
| 1353 | temp_sol_add(ji) = -(pgflux(ji) - zgflux(ji))*dt_sechiba/soilcap(ji) |
---|
| 1354 | |
---|
| 1355 | pgflux(ji) = zgflux(ji) |
---|
| 1356 | |
---|
| 1357 | ELSE |
---|
| 1358 | |
---|
| 1359 | temp_sol_add(ji) = zero |
---|
| 1360 | |
---|
| 1361 | ENDIF |
---|
| 1362 | |
---|
| 1363 | ENDDO |
---|
| 1364 | |
---|
| 1365 | !! 3. Define qsat_air with the subroutine src_parameter: |
---|
| 1366 | |
---|
| 1367 | CALL qsatcalc(kjpindex, tair, pb, qsat_air) |
---|
| 1368 | |
---|
| 1369 | CALL dev_qsatcalc(kjpindex, tair, pb, grad_qsat) |
---|
| 1370 | |
---|
| 1371 | ! grad_qsat(:)= (qsol_sat_new(:)- qsat_air(:)) / ((psnew(:) - epot_air(:)) / cp_air) ! * dt_sechiba |
---|
| 1372 | |
---|
| 1373 | warning_correction(:)=.FALSE. |
---|
| 1374 | DO ji=1,kjpindex |
---|
| 1375 | |
---|
| 1376 | !! \latexonly |
---|
| 1377 | !! \input{enerbilflux16.tex} |
---|
| 1378 | !! \endlatexonly |
---|
| 1379 | speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) |
---|
| 1380 | |
---|
| 1381 | !! \latexonly |
---|
| 1382 | !! \input{enerbilflux17.tex} |
---|
| 1383 | !! \endlatexonly |
---|
| 1384 | qc = speed * q_cdrag(ji) |
---|
| 1385 | |
---|
| 1386 | !! Derive the potential as defined by Penman & Monteith (Monteith, 1965) based on the correction |
---|
| 1387 | !! term developed by Chris Milly (1992). |
---|
| 1388 | IF ((evapot(ji) .GT. min_sechiba) .AND. ((psnew(ji) - epot_air(ji)) .NE. zero )) THEN |
---|
| 1389 | ! |
---|
| 1390 | !! \latexonly |
---|
| 1391 | !! \input{enerbilflux18.tex} |
---|
| 1392 | !! \endlatexonly |
---|
| 1393 | correction = (quatre * emis(ji) * c_stefan * tair(ji)**3 + rau(ji) * qc * cp_air + & |
---|
| 1394 | & chalev0 * rau(ji) * qc * grad_qsat(ji) * vevapp(ji) / evapot(ji) ) |
---|
| 1395 | |
---|
| 1396 | !! \latexonly |
---|
| 1397 | !! \input{enerbilflux19.tex} |
---|
| 1398 | !! \endlatexonly |
---|
| 1399 | IF (ABS(correction) .GT. min_sechiba) THEN |
---|
| 1400 | correction = chalev0 * rau(ji) * qc * grad_qsat(ji) * (un - vevapp(ji)/evapot(ji)) / correction |
---|
| 1401 | ELSE |
---|
| 1402 | warning_correction(ji)=.TRUE. |
---|
| 1403 | ENDIF |
---|
| 1404 | ELSE |
---|
| 1405 | correction = zero |
---|
| 1406 | ENDIF |
---|
| 1407 | correction = MAX (zero, correction) |
---|
| 1408 | |
---|
| 1409 | !! \latexonly |
---|
| 1410 | !! \input{enerbilflux20.tex} |
---|
| 1411 | !! \endlatexonly |
---|
| 1412 | evapot_corr(ji) = evapot(ji) / (un + correction) |
---|
| 1413 | |
---|
| 1414 | ENDDO |
---|
| 1415 | IF ( ANY(warning_correction) ) THEN |
---|
| 1416 | DO ji=1,kjpindex |
---|
| 1417 | IF ( warning_correction(ji) ) THEN |
---|
| 1418 | WRITE(numout,*) ji,"Denominateur de la correction de milly nul! Aucune correction appliquee" |
---|
| 1419 | ENDIF |
---|
| 1420 | ENDDO |
---|
| 1421 | ENDIF |
---|
| 1422 | IF (printlev>=3) WRITE (numout,*) ' enerbil_flux done ' |
---|
| 1423 | |
---|
| 1424 | END SUBROUTINE enerbil_flux |
---|
| 1425 | |
---|
| 1426 | |
---|
| 1427 | !! =========================================================================================================================== |
---|
| 1428 | !! SUBROUTINE : enerbil_evapveg |
---|
| 1429 | !! |
---|
| 1430 | !>\BRIEF Splits total evaporation loss into individual process |
---|
| 1431 | !! components and calculates the assimilation. |
---|
| 1432 | !! |
---|
| 1433 | !! DESCRIPTION : Based on the estimation of the fluxes in enerbil_flux, |
---|
| 1434 | !! this routine splits the total evaporation into transpiration interception loss from vegetation, |
---|
| 1435 | !! bare soil evaporation and snow sublimation. It then calculates the assimilation. |
---|
| 1436 | !! |
---|
| 1437 | !! MAIN OUTPUT VARIABLE(S) : vevapsno, vevapnu, transpir, vevapwet |
---|
| 1438 | !! |
---|
| 1439 | !! REFERENCE(S) : |
---|
| 1440 | !! - Best, MJ, Beljaars, A, Polcher, J & Viterbo, P, 2004. A proposed structure for coupling tiled |
---|
| 1441 | !! surfaces with the planetary boundary layer. Journal of Hydrometeorology, 5, pp.1271-1278 |
---|
| 1442 | !! - de Noblet-Ducoudré, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations |
---|
| 1443 | !! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General |
---|
| 1444 | !! Circulation Model. Journal of Climate, 6, pp.248-273 |
---|
| 1445 | !! - Dufresne, J-L & Ghattas, J, 2009. Description du schéma de la couche limite turbulente et la |
---|
| 1446 | !! interface avec la surface planetaire dans LMDZ, Technical note, available (22/12/11): |
---|
| 1447 | !! http://lmdz.lmd.jussieu.fr/developpeurs/notes-techniques/ressources/pbl_surface.pdf |
---|
| 1448 | !! - Guimberteau, M, 2010. Modélisation de l'hydrologie continentale et influences de l'irrigation |
---|
| 1449 | !! sur le cycle de l'eau, PhD Thesis, available (22/12/11): |
---|
| 1450 | !! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf |
---|
| 1451 | !! - Polcher, J. McAvaney, B, Viterbo, P, Gaertner, MA, Hahmann, A, Mahfouf, J-F, Noilhan, J |
---|
| 1452 | !! Phillips, TJ, Pitman, AJ, Schlosser, CA, Schulz, J-P, Timbal, B, Verseghy, D & |
---|
| 1453 | !! Xue, Y, 1998. A proposal for a general interface between land surface schemes and |
---|
| 1454 | !! general circulation models. Global and Planetary Change, 19, pp.261-276 |
---|
| 1455 | !! - Richtmeyer, RD, Morton, KW, 1967. Difference Methods for Initial-Value Problems. |
---|
| 1456 | !! Interscience Publishers |
---|
| 1457 | !! - Schulz, Jan-Peter, Lydia DÃŒmenil, Jan Polcher, 2001: On the Land SurfaceâAtmosphere |
---|
| 1458 | !! Coupling and Its Impact in a Single-Column Atmospheric Model. J. Appl. Meteor., 40, 642â663. |
---|
| 1459 | !! |
---|
| 1460 | !! FLOWCHART : None |
---|
| 1461 | !! \n |
---|
| 1462 | !_ ============================================================================================================================== |
---|
| 1463 | |
---|
| 1464 | SUBROUTINE enerbil_evapveg (kjpindex, vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, & |
---|
| 1465 | & rau, u, v, q_cdrag, qair, humrel, vevapsno, vevapnu , vevapflo, vevapwet, & |
---|
| 1466 | & transpir, transpot, evapot) |
---|
| 1467 | |
---|
| 1468 | !! 0. Variable and parameter declaration |
---|
| 1469 | |
---|
| 1470 | !! 0.1 Input variables |
---|
| 1471 | |
---|
| 1472 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-) |
---|
| 1473 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta1 !! Snow resistance (-) |
---|
| 1474 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta4 !! Bare soil resistance (-) |
---|
| 1475 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vbeta5 !! Floodplains resistance |
---|
| 1476 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density (kg m^{-3}) |
---|
| 1477 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u, v !! Wind velocity in directions u and v (m s^{-1}) |
---|
| 1478 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-) |
---|
| 1479 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity (kg kg^{-1}) |
---|
| 1480 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation (mm/tstep) |
---|
| 1481 | REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (in) :: humrel !! Soil moisture stress (within range 0 to 1) |
---|
| 1482 | !!$ DS 15022011 humrel was used in a previous version of Orchidee, developped by Nathalie. Need to be discussed if it should be introduce again |
---|
| 1483 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta2 !! Interception resistance (-) |
---|
| 1484 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3 !! Vegetation resistance (-) |
---|
| 1485 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3pot !! Vegetation resistance for potential transpiration |
---|
| 1486 | |
---|
| 1487 | !! 0.2 Output variables |
---|
| 1488 | |
---|
| 1489 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapsno !! Snow evaporation (mm day^{-1}) |
---|
| 1490 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapnu !! Bare soil evaporation (mm day^{-1}) |
---|
| 1491 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapflo !! Floodplains evaporation |
---|
| 1492 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpir !! Transpiration (mm day^{-1}) |
---|
| 1493 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: transpot !! Potential Transpiration |
---|
| 1494 | REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vevapwet !! Interception (mm day^{-1}) |
---|
| 1495 | |
---|
| 1496 | !! 0.3 Modified variables |
---|
| 1497 | |
---|
| 1498 | !! 0.4 Local variables |
---|
| 1499 | |
---|
| 1500 | INTEGER(i_std) :: ji, jv |
---|
| 1501 | REAL(r_std), DIMENSION(kjpindex) :: xx |
---|
| 1502 | REAL(r_std), DIMENSION(kjpindex) :: vbeta2sum, vbeta3sum |
---|
| 1503 | REAL(r_std) :: speed !! Speed (m s^{-1}) |
---|
| 1504 | |
---|
| 1505 | !_ ============================================================================================================================== |
---|
| 1506 | ! initialisation: utile pour calculer l'evaporation des floodplains dans lesquelles il y a de la vegetation |
---|
| 1507 | vbeta2sum(:) = 0. |
---|
| 1508 | vbeta3sum(:) = 0. |
---|
| 1509 | DO jv = 1, nvm |
---|
| 1510 | vbeta2sum(:) = vbeta2sum(:) + vbeta2(:,jv) |
---|
| 1511 | vbeta3sum(:) = vbeta3sum(:) + vbeta3(:,jv) |
---|
| 1512 | ENDDO |
---|
| 1513 | |
---|
| 1514 | !! 1. Compute vevapsno (evaporation from snow) and vevapnu (bare soil evaporation) |
---|
| 1515 | |
---|
| 1516 | DO ji=1,kjpindex |
---|
| 1517 | |
---|
| 1518 | !! \latexonly |
---|
| 1519 | !! \input{enerbilevapveg1.tex} |
---|
| 1520 | !! \endlatexonly |
---|
| 1521 | speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) |
---|
| 1522 | |
---|
| 1523 | |
---|
| 1524 | !! 1.1 Snow sublimation |
---|
| 1525 | !! \latexonly |
---|
| 1526 | !! \input{enerbilevapveg2.tex} |
---|
| 1527 | !! \endlatexonly |
---|
| 1528 | vevapsno(ji) = (un - vbeta5(ji)) * vbeta1(ji) * dt_sechiba * rau(ji) * speed * q_cdrag(ji) * (qsol_sat_new(ji) - qair(ji)) |
---|
| 1529 | |
---|
| 1530 | !! 1.2 Bare soil evaporation |
---|
| 1531 | !! \latexonly |
---|
| 1532 | !! \input{enerbilevapveg3.tex} |
---|
| 1533 | !! \endlatexonly |
---|
| 1534 | vevapnu(ji) = (un - vbeta1(ji)) * (un-vbeta5(ji)) * vbeta4(ji) * dt_sechiba * rau(ji) * speed * q_cdrag(ji) & |
---|
| 1535 | & * (qsol_sat_new(ji) - qair(ji)) |
---|
| 1536 | ! |
---|
| 1537 | ! 1.3 floodplains evaporation - transpiration et interception prioritaires dans les floodplains |
---|
| 1538 | ! |
---|
| 1539 | vevapflo(ji) = vbeta5(ji) & |
---|
| 1540 | & * dt_sechiba * rau(ji) * speed * q_cdrag(ji) * (qsol_sat_new(ji) - qair(ji)) |
---|
| 1541 | |
---|
| 1542 | END DO |
---|
| 1543 | |
---|
| 1544 | !! 2. Compute transpir (transpiration) and vevapwet (interception) |
---|
| 1545 | |
---|
| 1546 | !! Preliminaries |
---|
| 1547 | DO ji = 1, kjpindex |
---|
| 1548 | |
---|
| 1549 | !! \latexonly |
---|
| 1550 | !! \input{enerbilevapveg4.tex} |
---|
| 1551 | !! \endlatexonly |
---|
| 1552 | speed = MAX(min_wind, SQRT (u(ji)*u(ji) + v(ji)*v(ji))) |
---|
| 1553 | |
---|
| 1554 | !! \latexonly |
---|
| 1555 | !! \input{enerbilevapveg5.tex} |
---|
| 1556 | !! \endlatexonly |
---|
| 1557 | xx(ji) = dt_sechiba * (un-vbeta1(ji)) * (un-vbeta5(ji)) * (qsol_sat_new(ji)-qair(ji)) * rau(ji) * speed * q_cdrag(ji) |
---|
| 1558 | |
---|
| 1559 | ENDDO |
---|
| 1560 | |
---|
| 1561 | DO jv=1,nvm |
---|
| 1562 | DO ji=1,kjpindex |
---|
| 1563 | |
---|
| 1564 | !! 2.1 Calculate interception loss |
---|
| 1565 | !! \latexonly |
---|
| 1566 | !! \input{enerbilevapveg6.tex} |
---|
| 1567 | !! \endlatexonly |
---|
| 1568 | vevapwet(ji,jv) = xx(ji) * vbeta2(ji,jv) |
---|
| 1569 | ! |
---|
| 1570 | !! 2.2 Calculate transpiration |
---|
| 1571 | !! \latexonly |
---|
| 1572 | !! \input{enerbilevapveg7.tex} |
---|
| 1573 | !! \endlatexonly |
---|
| 1574 | transpir (ji,jv) = xx(ji) * vbeta3(ji,jv) |
---|
| 1575 | |
---|
| 1576 | transpot(ji,jv) = xx(ji) * vbeta3pot(ji,jv) |
---|
| 1577 | |
---|
| 1578 | END DO |
---|
| 1579 | END DO |
---|
| 1580 | |
---|
| 1581 | |
---|
| 1582 | |
---|
| 1583 | IF (printlev>=3) WRITE (numout,*) ' enerbil_evapveg done ' |
---|
| 1584 | |
---|
| 1585 | END SUBROUTINE enerbil_evapveg |
---|
| 1586 | |
---|
| 1587 | END MODULE enerbil |
---|