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