[7541] | 1 | ! ================================================================================================================================= |
---|
| 2 | ! MODULE : thermosoil |
---|
| 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 Calculates the soil temperatures by solving the heat |
---|
| 10 | !! diffusion equation within the soil. This module is only used with CWRR hydrology. |
---|
| 11 | !! |
---|
| 12 | !!\n DESCRIPTION : General important informations about the numerical scheme and |
---|
| 13 | !! the soil vertical discretization:\n |
---|
| 14 | !! - the soil is zmaxt deep (by default 10m) and divided into "ngrnd" layers. |
---|
| 15 | !! From 0-zmaxh(default 2m), the discretization is the same as for hydrology. |
---|
| 16 | !! From zmaxh(2m) and below, the depth increase linearly (by default) or geometrically. \n |
---|
| 17 | !! - "jg" is usually used as the index going from 1 to ngrnd to describe the |
---|
| 18 | !! layers, from top (jg=1) to bottom (jg=ngrnd)\n |
---|
| 19 | !! - the thermal numerical scheme is implicit finite differences.\n |
---|
| 20 | !! -- When it is resolved in thermosoil_profile at the present timestep t, the |
---|
| 21 | !! dependancy from the previous timestep (t-1) is hidden in the |
---|
| 22 | !! integration coefficients cgrnd and dgrnd, which are therefore |
---|
| 23 | !! calculated at the very end of thermosoil_main (call to |
---|
| 24 | !! thermosoil_coef) for use in the next timestep.\n |
---|
| 25 | !! -- At timestep t, the system becomes :\n |
---|
| 26 | !! |
---|
| 27 | !! T(k+1)=cgrnd(k)+dgrnd(k)*T(k) \n |
---|
| 28 | !! -- EQ1 -- \n |
---|
| 29 | !! |
---|
| 30 | !! (the bottom boundary condition has been used to obtained this equation).\n |
---|
| 31 | !! To solve it, the uppermost soil temperature T(1) is required. |
---|
| 32 | !! It is obtained from the surface temperature Ts, which is |
---|
| 33 | !! considered a linear extrapolation of T(1) and T(2)\n |
---|
| 34 | !! |
---|
| 35 | !! Ts=(1+lambda)*T(1) -lambda*T(2) \n |
---|
| 36 | !! -- EQ2--\n |
---|
| 37 | !! |
---|
| 38 | !! -- caveat 1 : Ts is called 'temp_soil_new' in this routine, |
---|
| 39 | !! don' t act.\n |
---|
| 40 | !! -- caveat 2 : actually, the surface temperature at time t Ts |
---|
| 41 | !! depends on the soil temperature at time t through the |
---|
| 42 | !! ground heat flux. This is again implicitly solved, with Ts(t) |
---|
| 43 | !! expressed as :\n |
---|
| 44 | !! |
---|
| 45 | !! soilcap*(Ts(t)-Ts(t-1))/dt=soilflx+otherfluxes(Ts(t))\n |
---|
| 46 | !! -- EQ3 --\n |
---|
| 47 | !! |
---|
| 48 | !! and the dependency from the previous timestep is hidden in |
---|
| 49 | !! soilcap and soilflx (apparent surface heat capacity and heat |
---|
| 50 | !! flux respectively). Soilcap and soilflx are therefore |
---|
| 51 | !! calculated at the previous timestep, at the very end of thermosoil |
---|
| 52 | !! (final call to thermosoil_coef) and stored to be used at the next time step. |
---|
| 53 | !! At timestep t, EQ3 is solved for Ts in enerbil, and Ts |
---|
| 54 | !! is used in thermosoil to get T(1) and solve EQ1.\n |
---|
| 55 | !! |
---|
| 56 | !! - lambda is the @tex $\mu$ @endtex of F. Hourdin' s PhD thesis, equation (A28); ie the |
---|
| 57 | !! coefficient of the linear extrapolation of Ts (surface temperature) from T1 and T2 (ptn(jg=1) and ptn(jg=2)), so that:\n |
---|
| 58 | !! Ts= (1+lambda)*T(1)-lambda*T(2) --EQ2-- \n |
---|
| 59 | !! lambda = (zlt(1))/((zlt(2)-zlt(1))) \n |
---|
| 60 | !! |
---|
| 61 | !! RECENT CHANGE(S) : - Change soil thermal properties to consider also soil texture, rev 2922. |
---|
| 62 | !! - Change vertical discretization, rev 2917. Note: In the revised thermosoil, |
---|
| 63 | !! cstgrnd and lskin are not needed any more. The depth znt, zlt and dlt |
---|
| 64 | !! are computed in vertical_soil and are in meter |
---|
| 65 | !! |
---|
| 66 | !! REFERENCE(S) : None |
---|
| 67 | !! |
---|
| 68 | !! SVN : |
---|
| 69 | !! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90 $ |
---|
| 70 | !! $Date: 2021-05-31 17:59:46 +0200 (Mon, 31 May 2021) $ |
---|
| 71 | !! $Revision: 7207 $ |
---|
| 72 | !! \n |
---|
| 73 | !_ ================================================================================================================================ |
---|
| 74 | |
---|
| 75 | MODULE thermosoil |
---|
| 76 | |
---|
| 77 | ! modules used : |
---|
| 78 | USE ioipsl |
---|
| 79 | USE ioipsl_para |
---|
| 80 | USE xios_orchidee |
---|
| 81 | USE constantes |
---|
| 82 | USE time, ONLY : one_day, dt_sechiba |
---|
| 83 | USE constantes_soil |
---|
| 84 | USE sechiba_io_p |
---|
| 85 | USE grid |
---|
| 86 | |
---|
| 87 | IMPLICIT NONE |
---|
| 88 | |
---|
| 89 | !private and public routines : |
---|
| 90 | PRIVATE |
---|
| 91 | PUBLIC :: thermosoil_main, thermosoil_clear, thermosoil_initialize, thermosoil_finalize, thermosoil_xios_initialize |
---|
| 92 | |
---|
| 93 | REAL(r_std), SAVE :: lambda !! See Module description |
---|
| 94 | !$OMP THREADPRIVATE(lambda) |
---|
| 95 | REAL(r_std), SAVE :: fz1, zalph !! usefull constants for diverse use |
---|
| 96 | !$OMP THREADPRIVATE(fz1, zalph) |
---|
| 97 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn !! vertically discretized |
---|
| 98 | !! soil temperatures @tex ($K$) @endtex. |
---|
| 99 | !$OMP THREADPRIVATE(ptn) |
---|
| 100 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: dz1 !! numerical constant used in the thermal numerical |
---|
| 101 | !! scheme @tex ($m^{-1}$) @endtex. ; it corresponds |
---|
| 102 | !! to the coefficient @tex $d_k$ @endtex of equation |
---|
| 103 | !! (A.12) in F. Hourdin PhD thesis. |
---|
| 104 | !$OMP THREADPRIVATE(dz1) |
---|
| 105 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: cgrnd !! integration coefficient for the numerical scheme, |
---|
| 106 | !! see eq.1 |
---|
| 107 | !$OMP THREADPRIVATE(cgrnd) |
---|
| 108 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dgrnd !! integration coefficient for the numerical scheme, |
---|
| 109 | !! see eq.1 |
---|
| 110 | !$OMP THREADPRIVATE(dgrnd) |
---|
| 111 | |
---|
| 112 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa !! volumetric vertically discretized soil heat |
---|
| 113 | !! capacity @tex ($J K^{-1} m^{-3}$) @endtex. |
---|
| 114 | !! It depends on the soil |
---|
| 115 | !! moisture content (shum_ngrnd_perma) and is calculated at |
---|
| 116 | !! each time step in thermosoil_coef. |
---|
| 117 | !$OMP THREADPRIVATE(pcapa) |
---|
| 118 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa !! vertically discretized soil thermal conductivity |
---|
| 119 | !! @tex ($W K^{-1} m^{-1}$) @endtex. Same as pcapa. |
---|
| 120 | !$OMP THREADPRIVATE(pkappa) |
---|
| 121 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_snow !! volumetric vertically discretized snow heat |
---|
| 122 | !! capacity @tex ($J K^{-1} m^{-3}$) @endtex. |
---|
| 123 | !$OMP THREADPRIVATE(pcapa_snow) |
---|
| 124 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pkappa_snow !! vertically discretized snow thermal conductivity |
---|
| 125 | !! @tex ($W K^{-1} m^{-1}$) @endtex. |
---|
| 126 | !$OMP THREADPRIVATE(pkappa_snow) |
---|
| 127 | |
---|
| 128 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcapa_en !! heat capacity used for surfheat_incr and |
---|
| 129 | !! coldcont_incr |
---|
| 130 | !$OMP THREADPRIVATE(pcapa_en) |
---|
| 131 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ptn_beg !! ptn as it is after thermosoil_profile but before thermosoil_coef, |
---|
| 132 | !! used in thermosoil_readjust |
---|
| 133 | !$OMP THREADPRIVATE(ptn_beg) |
---|
| 134 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: temp_sol_beg !! Surface temperature at previous timestep (K) |
---|
| 135 | !$OMP THREADPRIVATE(temp_sol_beg) |
---|
| 136 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: surfheat_incr !! Change in soil heat content during the timestep |
---|
| 137 | !! @tex ($J$) @endtex. |
---|
| 138 | !$OMP THREADPRIVATE(surfheat_incr) |
---|
| 139 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: coldcont_incr !! Change in snow heat content @tex ($J$) @endtex. |
---|
| 140 | !$OMP THREADPRIVATE(coldcont_incr) |
---|
| 141 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: shum_ngrnd_perma !! Saturation degree on the thermal axes (0-1, dimensionless) |
---|
| 142 | !$OMP THREADPRIVATE(shum_ngrnd_perma) |
---|
| 143 | |
---|
| 144 | ! Variables related to soil freezing |
---|
| 145 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: profil_froz !! Frozen fraction of the soil on hydrological levels (-) |
---|
| 146 | !$OMP THREADPRIVATE(profil_froz) |
---|
| 147 | REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:) :: e_soil_lat !! Latent heat released or consumed in the freezing/thawing processes summed vertically |
---|
| 148 | !! for the whole soil (J/m2) and on the whole simulation to check/correct energy conservation |
---|
| 149 | !$OMP THREADPRIVATE(e_soil_lat) |
---|
| 150 | REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:,:) :: pcappa_supp !! Additional surfacic heat capacity due to soil freezing for each soil layer (J/K/m2) |
---|
| 151 | !$OMP THREADPRIVATE(pcappa_supp) |
---|
| 152 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: dz5 !! Used for numerical calculation [-] |
---|
| 153 | !$OMP THREADPRIVATE(dz5) |
---|
| 154 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: QZ !! quartz content [-] |
---|
| 155 | !$OMP THREADPRIVATE(QZ) |
---|
| 156 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: so_capa_dry_ns !! Dry soil Heat capacity of soils,J.m^{-3}.K^{-1} |
---|
| 157 | !$OMP THREADPRIVATE(so_capa_dry_ns) |
---|
| 158 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: so_capa_ice !! Heat capacity of saturated frozen soil (J/K/m3) |
---|
| 159 | !$OMP THREADPRIVATE(so_capa_ice) |
---|
| 160 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mc_layt !! Volumetric soil moisture (liquid+ice) (m3/m3) on the thermodynamical levels at interface |
---|
| 161 | !$OMP THREADPRIVATE(mc_layt) |
---|
| 162 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mcl_layt !! Volumetric soil moisture (liquid) (m3/m3) on the thermodynamical levels at interface |
---|
| 163 | !$OMP THREADPRIVATE(mcl_layt) |
---|
| 164 | REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_layt !! Total soil moisture content for each layer (liquid+ice) (mm) on the thermodynamical levels |
---|
| 165 | !$OMP THREADPRIVATE(tmc_layt) |
---|
| 166 | INTEGER(i_std), SAVE :: brk_flag = 0 !! Flag to consider bedrock: 0.no; 1.yes |
---|
| 167 | !$OMP THREADPRIVATE(brk_flag) |
---|
| 168 | |
---|
| 169 | |
---|
| 170 | CONTAINS |
---|
| 171 | |
---|
| 172 | |
---|
| 173 | !! ============================================================================================================================= |
---|
| 174 | !! SUBROUTINE: thermosoil_xios_initialize |
---|
| 175 | !! |
---|
| 176 | !>\BRIEF Initialize xios dependant defintion before closing context defintion |
---|
| 177 | !! |
---|
| 178 | !! DESCRIPTION: Initialize xios dependant defintion before closing context defintion |
---|
| 179 | !! Reading is deactivated if the sechiba restart file exists because the |
---|
| 180 | !! variable should be in the restart file already. |
---|
| 181 | !! |
---|
| 182 | !! \n |
---|
| 183 | !_ ============================================================================================================================== |
---|
| 184 | |
---|
| 185 | SUBROUTINE thermosoil_xios_initialize |
---|
| 186 | |
---|
| 187 | CHARACTER(LEN=255) :: filename, name |
---|
| 188 | |
---|
| 189 | filename = 'reftemp.nc' |
---|
| 190 | CALL getin_p('REFTEMP_FILE',filename) |
---|
| 191 | |
---|
| 192 | name = filename(1:LEN_TRIM(FILENAME)-3) |
---|
| 193 | CALL xios_orchidee_set_file_attr("reftemp_file",name=name) |
---|
| 194 | |
---|
| 195 | ! Check if the reftemp file will be read by XIOS, by IOIPSL or not at all |
---|
| 196 | IF (xios_interpolation .AND. read_reftemp .AND. restname_in=='NONE') THEN |
---|
| 197 | ! The reftemp file will be read using XIOS |
---|
| 198 | IF (printlev>=2) WRITE(numout,*) 'Reading of reftemp file will be done later using XIOS. The filename is ', filename |
---|
| 199 | ELSE |
---|
| 200 | IF (.NOT. read_reftemp) THEN |
---|
| 201 | IF (printlev>=2) WRITE (numout,*) 'No reading of reftemp will be done because read_reftemp=FALSE' |
---|
| 202 | ELSE IF (restname_in=='NONE') THEN |
---|
| 203 | IF (printlev>=2) WRITE (numout,*) 'The reftemp file will be read later by IOIPSL' |
---|
| 204 | ELSE |
---|
| 205 | IF (printlev>=2) WRITE (numout,*) 'The reftemp file will not be read because the restart file exists.' |
---|
| 206 | END IF |
---|
| 207 | |
---|
| 208 | ! The reftemp file will not be read by XIOS. Now deactivate albedo for XIOS. |
---|
| 209 | CALL xios_orchidee_set_file_attr("reftemp_file",enabled=.FALSE.) |
---|
| 210 | CALL xios_orchidee_set_field_attr("reftemp_interp",enabled=.FALSE.) |
---|
| 211 | ENDIF |
---|
| 212 | |
---|
| 213 | END SUBROUTINE thermosoil_xios_initialize |
---|
| 214 | |
---|
| 215 | !! ============================================================================================================================= |
---|
| 216 | !! SUBROUTINE : thermosoil_initialize |
---|
| 217 | !! |
---|
| 218 | !>\BRIEF Allocate module variables, read from restart file or initialize with default values |
---|
| 219 | !! |
---|
| 220 | !! DESCRIPTION : Allocate module variables, read from restart file or initialize with default values. |
---|
| 221 | !! Call thermosoil_var_init to calculate physical constants. |
---|
| 222 | !! Call thermosoil_coef to calculate thermal soil properties. |
---|
| 223 | !! |
---|
| 224 | !! RECENT CHANGE(S) : None |
---|
| 225 | !! |
---|
| 226 | !! REFERENCE(S) : None |
---|
| 227 | !! |
---|
| 228 | !! FLOWCHART : None |
---|
| 229 | !! \n |
---|
| 230 | !_ ============================================================================================================================== |
---|
| 231 | SUBROUTINE thermosoil_initialize (kjit, kjpindex, rest_id, mcs, & |
---|
| 232 | temp_sol_new, snow, shumdiag_perma, & |
---|
| 233 | soilcap, soilflx, stempdiag, & |
---|
| 234 | gtemp, & |
---|
| 235 | mc_layh, mcl_layh, tmc_layh, njsc, & |
---|
| 236 | frac_snow_veg,frac_snow_nobio,totfrac_nobio, & |
---|
| 237 | snowdz, snowrho, snowtemp, lambda_snow, cgrnd_snow, dgrnd_snow, pb) |
---|
| 238 | |
---|
| 239 | !! 0. Variable and parameter declaration |
---|
| 240 | !! 0.1 Input variables |
---|
| 241 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless) |
---|
| 242 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
| 243 | INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier (unitless) |
---|
| 244 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3) |
---|
| 245 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step, |
---|
| 246 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg) |
---|
| 247 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree (0-1, unitless) |
---|
| 248 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mc_layh !! Volumetric soil moisture content (liquid+ice) for hydrological layers, at node (m3/m3) |
---|
| 249 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content (liquid) for hydrological layers, at node (m3/m3) |
---|
| 250 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content(liquid+ice) for hydrological layers (mm) |
---|
| 251 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) |
---|
| 252 | REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area |
---|
| 253 | REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area |
---|
| 254 | REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+... |
---|
| 255 | !! (unitless,0-1) |
---|
| 256 | REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowdz !! Snow depth |
---|
| 257 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density |
---|
| 258 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K) |
---|
| 259 | REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa) |
---|
| 260 | |
---|
| 261 | !! 0.2 Output variables |
---|
| 262 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilcap !! apparent surface heat capacity considering snow and soil surface (J m-2 K-1) |
---|
| 263 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilflx !! apparent soil heat flux considering snow and soil surface (W m-2) |
---|
| 264 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! temperature profile on the levels in hydrol(K) |
---|
| 265 | REAL(r_std),DIMENSION (kjpindex),INTENT(out) :: gtemp !! First soil layer temperature |
---|
| 266 | |
---|
| 267 | !! 0.3 Modified variables |
---|
| 268 | REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature |
---|
| 269 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme |
---|
| 270 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme |
---|
| 271 | |
---|
| 272 | !! 0.4 Local variables |
---|
| 273 | INTEGER(i_std) :: ier, i, jg, jsc |
---|
| 274 | LOGICAL :: calculate_coef !! Local flag to initialize variables by call to thermosoil_coef |
---|
| 275 | !_ ================================================================================================================================ |
---|
| 276 | |
---|
| 277 | |
---|
| 278 | ! |
---|
| 279 | ! !! Flag to consider bedrock at deeper layers |
---|
| 280 | ! !! It affects heat capacity and thermal conductivity (energy balance). |
---|
| 281 | ! |
---|
| 282 | !Config Key = BEDROCK_FLAG |
---|
| 283 | !Config Desc = Flag to consider bedrock at deeper layers. |
---|
| 284 | !Config If = |
---|
| 285 | !Config Def = 0 |
---|
| 286 | !Config Help = 0, no, 1, yes. |
---|
| 287 | !Config Units = [FLAG] |
---|
| 288 | brk_flag = 0 |
---|
| 289 | CALL getin_p('BEDROCK_FLAG', brk_flag) |
---|
| 290 | |
---|
| 291 | IF (printlev >= 3) WRITE (numout,*) 'Start thermosoil_initialize ' |
---|
| 292 | |
---|
| 293 | !! 1. Allocate soil temperatures variables |
---|
| 294 | ALLOCATE (ptn(kjpindex,ngrnd),stat=ier) |
---|
| 295 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ptn','','') |
---|
| 296 | |
---|
| 297 | ALLOCATE (dz1(ngrnd),stat=ier) |
---|
| 298 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dz1','','') |
---|
| 299 | |
---|
| 300 | ALLOCATE (cgrnd(kjpindex,ngrnd-1),stat=ier) |
---|
| 301 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of cgrnd','','') |
---|
| 302 | |
---|
| 303 | ALLOCATE (dgrnd(kjpindex,ngrnd-1),stat=ier) |
---|
| 304 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dgrnd','','') |
---|
| 305 | |
---|
| 306 | ALLOCATE (pcapa(kjpindex,ngrnd),stat=ier) |
---|
| 307 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa','','') |
---|
| 308 | |
---|
| 309 | ALLOCATE (pkappa(kjpindex,ngrnd),stat=ier) |
---|
| 310 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pkappa','','') |
---|
| 311 | |
---|
| 312 | ALLOCATE (pcapa_snow(kjpindex,nsnow),stat=ier) |
---|
| 313 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa_snow','','') |
---|
| 314 | |
---|
| 315 | ALLOCATE (pkappa_snow(kjpindex,nsnow),stat=ier) |
---|
| 316 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pkappa_snow','','') |
---|
| 317 | |
---|
| 318 | ! Temporary fix: Initialize following variable because they are output to xios before the first calculation |
---|
| 319 | pcapa = 0 |
---|
| 320 | pkappa = 0 |
---|
| 321 | pcapa_snow = 0 |
---|
| 322 | pkappa_snow = 0 |
---|
| 323 | |
---|
| 324 | ALLOCATE (surfheat_incr(kjpindex),stat=ier) |
---|
| 325 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of surfheat_incr','','') |
---|
| 326 | |
---|
| 327 | ALLOCATE (coldcont_incr(kjpindex),stat=ier) |
---|
| 328 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of coldcont_incr','','') |
---|
| 329 | |
---|
| 330 | ALLOCATE (pcapa_en(kjpindex,ngrnd),stat=ier) |
---|
| 331 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of pcapa_en','','') |
---|
| 332 | ! Initialization to zero used at first time step in thermosoil_energy_diag, only for diagnostic variables coldcont_incr and surfheat_incr |
---|
| 333 | pcapa_en(:,:) = 0. |
---|
| 334 | |
---|
| 335 | ALLOCATE (ptn_beg(kjpindex,ngrnd),stat=ier) |
---|
| 336 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ptn_beg','','') |
---|
| 337 | |
---|
| 338 | ALLOCATE (temp_sol_beg(kjpindex),stat=ier) |
---|
| 339 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of temp_sol_beg','','') |
---|
| 340 | |
---|
| 341 | ALLOCATE (shum_ngrnd_perma(kjpindex,ngrnd),stat=ier) |
---|
| 342 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of shum_ngrnd_perma','','') |
---|
| 343 | |
---|
| 344 | ALLOCATE (profil_froz(kjpindex,ngrnd),stat=ier) |
---|
| 345 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of profil_froz','','') |
---|
| 346 | |
---|
| 347 | IF (ok_freeze_thermix) THEN |
---|
| 348 | ALLOCATE (pcappa_supp(kjpindex,ngrnd),stat=ier) |
---|
| 349 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of ok_freeze_termix','','') |
---|
| 350 | ! Initialization to zero used at first time step only for diagnostic output. |
---|
| 351 | ! This variable is only used in thermosoil_readajust and always calculated before in thermosoil_getdiff. |
---|
| 352 | pcappa_supp(:,:) = 0. |
---|
| 353 | END IF |
---|
| 354 | IF (ok_Ecorr) THEN |
---|
| 355 | ALLOCATE (e_soil_lat(kjpindex),stat=ier) |
---|
| 356 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of e_soil_lat','','') |
---|
| 357 | END IF |
---|
| 358 | |
---|
| 359 | ALLOCATE (dz5(ngrnd),stat=ier) |
---|
| 360 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of dz5','','') |
---|
| 361 | |
---|
| 362 | ALLOCATE (mc_layt(kjpindex,ngrnd),stat=ier) |
---|
| 363 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mc_layt','','') |
---|
| 364 | |
---|
| 365 | ALLOCATE (mcl_layt(kjpindex,ngrnd),stat=ier) |
---|
| 366 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of mcl_layt','','') |
---|
| 367 | |
---|
| 368 | ALLOCATE (tmc_layt(kjpindex,ngrnd),stat=ier) |
---|
| 369 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of tmc_layt','','') |
---|
| 370 | |
---|
| 371 | ALLOCATE (QZ(nscm),stat=ier) |
---|
| 372 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of QZ','','') |
---|
| 373 | |
---|
| 374 | ALLOCATE (so_capa_dry_ns(nscm),stat=ier) |
---|
| 375 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_dry_ns','','') |
---|
| 376 | |
---|
| 377 | ALLOCATE (so_capa_ice(kjpindex),stat=ier) |
---|
| 378 | IF (ier /= 0) CALL ipslerr_p(3,'thermosoil_initialize', 'Error in allocation of so_capa_ice','','') |
---|
| 379 | |
---|
| 380 | |
---|
| 381 | !! Soil texture choose |
---|
| 382 | SELECTCASE (nscm) |
---|
| 383 | CASE (3) |
---|
| 384 | QZ(:) = QZ_fao(:) |
---|
| 385 | so_capa_dry_ns(:) = so_capa_dry_ns_fao(:) |
---|
| 386 | CASE (13) !Salma changed from 12 to 13 for the new class Oxisols |
---|
| 387 | QZ(:) = QZ_usda(:) |
---|
| 388 | so_capa_dry_ns(:) = so_capa_dry_ns_usda(:) |
---|
| 389 | CASE DEFAULT |
---|
| 390 | WRITE (numout,*) 'Unsupported soil type classification. Choose between zobler, fao and usda according to the map' |
---|
| 391 | STOP 'thermosoil_initialize' |
---|
| 392 | ENDSELECT |
---|
| 393 | |
---|
| 394 | |
---|
| 395 | !! 2. Initialize variable from restart file or with default values |
---|
| 396 | |
---|
| 397 | !! Reads restart files for soil temperatures only. If no restart file is |
---|
| 398 | !! found, the initial soil temperature is by default set to 280K at all depths. The user |
---|
| 399 | !! can decide to initialize soil temperatures at an other value, in which case he should set the flag THERMOSOIL_TPRO |
---|
| 400 | !! to this specific value in the run.def. |
---|
| 401 | IF (printlev>=3) WRITE (numout,*) 'Read restart file for THERMOSOIL variables' |
---|
| 402 | |
---|
| 403 | CALL ioconf_setatt_p('UNITS', 'K') |
---|
| 404 | CALL ioconf_setatt_p('LONG_NAME','Soil Temperature profile') |
---|
| 405 | CALL restget_p (rest_id, 'ptn', nbp_glo, ngrnd, 1, kjit, .TRUE., ptn, "gather", nbp_glo, index_g) |
---|
| 406 | |
---|
| 407 | ! Initialize ptn if it was not found in restart file |
---|
| 408 | IF (ALL(ptn(:,:)==val_exp)) THEN |
---|
| 409 | ! ptn was not found in restart file |
---|
| 410 | |
---|
| 411 | IF (read_reftemp) THEN |
---|
| 412 | ! Read variable ptn from file |
---|
| 413 | CALL thermosoil_read_reftempfile(kjpindex,lalo,ptn) |
---|
| 414 | ELSE |
---|
| 415 | ! Initialize ptn with a constant value which can be set in run.def |
---|
| 416 | |
---|
| 417 | !Config Key = THERMOSOIL_TPRO |
---|
| 418 | !Config Desc = Initial soil temperature profile if not found in restart |
---|
| 419 | !Config Def = 280. |
---|
| 420 | !Config If = OK_SECHIBA |
---|
| 421 | !Config Help = The initial value of the temperature profile in the soil if |
---|
| 422 | !Config its value is not found in the restart file. Here |
---|
| 423 | !Config we only require one value as we will assume a constant |
---|
| 424 | !Config throughout the column. |
---|
| 425 | !Config Units = Kelvin [K] |
---|
| 426 | CALL setvar_p (ptn, val_exp,'THERMOSOIL_TPRO',280._r_std) |
---|
| 427 | END IF |
---|
| 428 | END IF |
---|
| 429 | |
---|
| 430 | ! Initialize ptn_beg (variable needed in thermosoil_readadjust called from thermosoil_coef) |
---|
| 431 | ptn_beg(:,:) = ptn(:,:) |
---|
| 432 | |
---|
| 433 | ! Initialize temp_sol_beg with values from previous time-step |
---|
| 434 | temp_sol_beg(:) = temp_sol_new(:) |
---|
| 435 | |
---|
| 436 | ! Read e_soil_lat from restart file or initialize |
---|
| 437 | IF (ok_Ecorr) THEN |
---|
| 438 | CALL restget_p (rest_id, 'e_soil_lat', nbp_glo, 1, 1, kjit, .TRUE., & |
---|
| 439 | e_soil_lat, "gather", nbp_glo, index_g) |
---|
| 440 | CALL setvar_p (e_soil_lat, val_exp,'NO_KEYWORD',zero) |
---|
| 441 | END IF |
---|
| 442 | |
---|
| 443 | ! Read gtemp from restart file |
---|
| 444 | CALL restget_p (rest_id, 'gtemp', nbp_glo, 1, 1, kjit, .TRUE., & |
---|
| 445 | gtemp, "gather", nbp_glo, index_g) |
---|
| 446 | CALL setvar_p (gtemp, val_exp,'NO_KEYWORD',zero) |
---|
| 447 | |
---|
| 448 | |
---|
| 449 | ! Read variables calculated in thermosoil_coef from restart file |
---|
| 450 | ! If the variables were not found in the restart file, the logical |
---|
| 451 | ! calculate_coef will be true and thermosoil_coef will be called further below. |
---|
| 452 | ! These variables need to be in the restart file to avoid a time shift that |
---|
| 453 | ! would be done using thermosoil_coef at this stage. |
---|
| 454 | calculate_coef=.FALSE. |
---|
| 455 | CALL ioconf_setatt_p('UNITS', 'J m-2 K-1') |
---|
| 456 | CALL ioconf_setatt_p('LONG_NAME','Apparent surface heat capacity') |
---|
| 457 | CALL restget_p (rest_id, 'soilcap', nbp_glo, 1, 1, kjit, .TRUE., & |
---|
| 458 | soilcap, "gather", nbp_glo, index_g) |
---|
| 459 | IF (ALL(soilcap(:)==val_exp)) calculate_coef=.TRUE. |
---|
| 460 | |
---|
| 461 | CALL ioconf_setatt_p('UNITS', 'W m-2') |
---|
| 462 | CALL ioconf_setatt_p('LONG_NAME','Apparent soil heat flux') |
---|
| 463 | CALL restget_p (rest_id, 'soilflx', nbp_glo, 1, 1, kjit, .TRUE., & |
---|
| 464 | soilflx, "gather", nbp_glo, index_g) |
---|
| 465 | IF (ALL(soilflx(:)==val_exp)) calculate_coef=.TRUE. |
---|
| 466 | |
---|
| 467 | CALL ioconf_setatt_p('UNITS', '') |
---|
| 468 | CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme') |
---|
| 469 | CALL restget_p (rest_id, 'cgrnd', nbp_glo, ngrnd-1, 1, kjit, .TRUE., & |
---|
| 470 | cgrnd, "gather", nbp_glo, index_g) |
---|
| 471 | IF (ALL(cgrnd(:,:)==val_exp)) calculate_coef=.TRUE. |
---|
| 472 | |
---|
| 473 | CALL ioconf_setatt_p('UNITS', '') |
---|
| 474 | CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme') |
---|
| 475 | CALL restget_p (rest_id, 'dgrnd', nbp_glo, ngrnd-1, 1, kjit, .TRUE., & |
---|
| 476 | dgrnd, "gather", nbp_glo, index_g) |
---|
| 477 | IF (ALL(dgrnd(:,:)==val_exp)) calculate_coef=.TRUE. |
---|
| 478 | |
---|
| 479 | CALL ioconf_setatt_p('UNITS', '') |
---|
| 480 | CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme') |
---|
| 481 | CALL restget_p (rest_id, 'cgrnd_snow', nbp_glo, nsnow, 1, kjit, .TRUE., & |
---|
| 482 | cgrnd_snow, "gather", nbp_glo, index_g) |
---|
| 483 | IF (ALL(cgrnd_snow(:,:)==val_exp)) calculate_coef=.TRUE. |
---|
| 484 | |
---|
| 485 | CALL ioconf_setatt_p('UNITS', '') |
---|
| 486 | CALL ioconf_setatt_p('LONG_NAME','Integration coefficient for the numerical scheme') |
---|
| 487 | CALL restget_p (rest_id, 'dgrnd_snow', nbp_glo, nsnow, 1, kjit, .TRUE., & |
---|
| 488 | dgrnd_snow, "gather", nbp_glo, index_g) |
---|
| 489 | IF (ALL(dgrnd_snow(:,:)==val_exp)) calculate_coef=.TRUE. |
---|
| 490 | |
---|
| 491 | CALL ioconf_setatt_p('UNITS', '') |
---|
| 492 | CALL ioconf_setatt_p('LONG_NAME','Coefficient of the linear extrapolation of surface temperature') |
---|
| 493 | CALL restget_p (rest_id, 'lambda_snow', nbp_glo, 1, 1, kjit, .TRUE., & |
---|
| 494 | lambda_snow, "gather", nbp_glo, index_g) |
---|
| 495 | IF (ALL(lambda_snow(:)==val_exp)) calculate_coef=.TRUE. |
---|
| 496 | |
---|
| 497 | !! 2.2 Computes some physical constants and arrays depending on the soil vertical discretization |
---|
| 498 | |
---|
| 499 | ! Calculate so_capa_ice |
---|
| 500 | so_capa_ice(:) = so_capa_dry + mcs(:)*capa_ice*rho_ice |
---|
| 501 | IF (printlev>=2) WRITE(numout,*) 'Calculation of so_capa_ice(:)=', so_capa_ice(:),' and capa_ice=',capa_ice |
---|
| 502 | |
---|
| 503 | ! Computing some usefull constants for the numerical scheme |
---|
| 504 | ! Use znt(depth of nodes) and zlt(depth of deeper layer interface) from vertical_soil module. |
---|
| 505 | DO jg=1,ngrnd-1 |
---|
| 506 | dz1(jg) = un / (znt(jg+1) - znt(jg)) |
---|
| 507 | dz5(jg) = (zlt(jg) - znt(jg)) * dz1(jg) |
---|
| 508 | ENDDO |
---|
| 509 | dz5(ngrnd) = 0.0 |
---|
| 510 | lambda = znt(1) * dz1(1) |
---|
| 511 | |
---|
| 512 | ! Send out the temperature profile on the first nslm levels(the levels treated in hydrol) |
---|
| 513 | stempdiag(:,:) = ptn(:,1:nslm) |
---|
| 514 | |
---|
| 515 | |
---|
| 516 | !! 2.3. Computes cgrnd, dgrnd, soilflx and soilcap coefficients only if they were not found in restart file. |
---|
| 517 | IF (calculate_coef) THEN |
---|
| 518 | ! Interpolate variables needed by thermosoil_coef to the thermal levels |
---|
| 519 | CALL thermosoil_humlev(kjpindex, shumdiag_perma, mc_layh, mcl_layh, tmc_layh) |
---|
| 520 | |
---|
| 521 | IF (printlev>=3) WRITE (numout,*) 'thermosoil_coef will be called in the intialization phase' |
---|
| 522 | CALL thermosoil_coef (& |
---|
| 523 | kjpindex, temp_sol_new, snow, njsc, & |
---|
| 524 | mcs, frac_snow_veg, frac_snow_nobio, totfrac_nobio, & |
---|
| 525 | snowdz, snowrho, snowtemp, pb, & |
---|
| 526 | ptn, & |
---|
| 527 | soilcap, soilflx, cgrnd, dgrnd,& |
---|
| 528 | lambda_snow, cgrnd_snow, dgrnd_snow) |
---|
| 529 | END IF |
---|
| 530 | |
---|
| 531 | END SUBROUTINE thermosoil_initialize |
---|
| 532 | |
---|
| 533 | |
---|
| 534 | !! ================================================================================================================================ |
---|
| 535 | !! SUBROUTINE : thermosoil_main |
---|
| 536 | !! |
---|
| 537 | !>\BRIEF Thermosoil_main computes the soil thermal properties and dynamics, ie solves |
---|
| 538 | !! the heat diffusion equation within the soil. |
---|
| 539 | !! |
---|
| 540 | !! DESCRIPTION : The resolution of the soil heat diffusion equation |
---|
| 541 | !! relies on a numerical finite-difference implicit scheme |
---|
| 542 | !! fully described in the reference and in the header of the thermosoil module. |
---|
| 543 | !! - The dependency of the previous timestep hidden in the |
---|
| 544 | !! integration coefficients cgrnd and dgrnd (EQ1), calculated in thermosoil_coef, and |
---|
| 545 | !! called at the end of the routine to prepare for the next timestep. |
---|
| 546 | !! - The effective computation of the new soil temperatures is performed in thermosoil_profile. |
---|
| 547 | !! |
---|
| 548 | !! - thermosoil_coef calculates the coefficients for the numerical scheme for the very first iteration of thermosoil; |
---|
| 549 | !! after that, thermosoil_coef is called only at the end of the module to calculate the coefficients for the next timestep. |
---|
| 550 | !! - thermosoil_profile solves the numerical scheme.\n |
---|
| 551 | !! |
---|
| 552 | !! - Flags : one unique flag : THERMOSOIL_TPRO (to be set to the desired initial soil in-depth temperature in K; by default 280K) |
---|
| 553 | !! |
---|
| 554 | !! RECENT CHANGE(S) : Change vertical discretization (consistent with hydrology layers) and soil thermal properties (taking into account soil texture effects). |
---|
| 555 | !! |
---|
| 556 | !! MAIN OUTPUT VARIABLE(S): vertically discretized soil temperatures ptn, soil |
---|
| 557 | !! thermal properties (pcapa, pkappa), apparent surface heat capacity (soilcap) |
---|
| 558 | !! and heat flux (soilflx) to be used in enerbil at the next timestep to solve |
---|
| 559 | !! the surface energy balance. |
---|
| 560 | !! |
---|
| 561 | !! REFERENCE(S) : |
---|
| 562 | !! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres, |
---|
| 563 | !! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin' s PhD thesis relative to the thermal |
---|
| 564 | !! integration scheme has been scanned and is provided along with the documentation, with name : |
---|
| 565 | !! Hourdin_1992_PhD_thermal_scheme.pdf |
---|
| 566 | !! |
---|
| 567 | !! FLOWCHART : |
---|
| 568 | !! \latexonly |
---|
| 569 | !! \includegraphics[scale = 1]{thermosoil_flowchart.png} |
---|
| 570 | !! \endlatexonly |
---|
| 571 | !! |
---|
| 572 | !! \n |
---|
| 573 | !_ ================================================================================================================================ |
---|
| 574 | |
---|
| 575 | SUBROUTINE thermosoil_main (kjit, kjpindex, & |
---|
| 576 | index, indexgrnd, mcs, & |
---|
| 577 | temp_sol_new, snow, soilcap, soilflx, & |
---|
| 578 | shumdiag_perma, stempdiag, ptnlev1, rest_id, hist_id, hist2_id, & |
---|
| 579 | snowdz,snowrho,snowtemp,gtemp,pb,& |
---|
| 580 | mc_layh, mcl_layh, tmc_layh, njsc, frac_snow_veg,frac_snow_nobio,totfrac_nobio,temp_sol_add, & |
---|
| 581 | lambda_snow, cgrnd_snow, dgrnd_snow) |
---|
| 582 | |
---|
| 583 | !! 0. Variable and parameter declaration |
---|
| 584 | |
---|
| 585 | !! 0.1 Input variables |
---|
| 586 | |
---|
| 587 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless) |
---|
| 588 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
| 589 | INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! Restart_ file and history file identifier |
---|
| 590 | !! (unitless) |
---|
| 591 | INTEGER(i_std),INTENT (in) :: hist2_id !! history file 2 identifier (unitless) |
---|
| 592 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (unitless) |
---|
| 593 | INTEGER(i_std),DIMENSION (kjpindex*ngrnd), INTENT (in):: indexgrnd !! Indeces of the points on the 3D map (vertical |
---|
| 594 | !! dimension towards the ground) (unitless) |
---|
| 595 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3) |
---|
| 596 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_new !! Surface temperature at the present time-step, |
---|
| 597 | !! Ts @tex ($K$) @endtex |
---|
| 598 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass @tex ($kg$) @endtex. |
---|
| 599 | !! Caveat: when there is snow on the |
---|
| 600 | !! ground, the snow is integrated into the soil for |
---|
| 601 | !! the calculation of the thermal dynamics. It means |
---|
| 602 | !! that the uppermost soil layers can completely or |
---|
| 603 | !! partially consist in snow. In the second case, zx1 |
---|
| 604 | !! and zx2 are the fraction of the soil layer |
---|
| 605 | !! consisting in snow and 'normal' soil, respectively |
---|
| 606 | !! This is calculated in thermosoil_coef. |
---|
| 607 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree (0-1, unitless) |
---|
| 608 | REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowdz !! Snow depth |
---|
| 609 | REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (in) :: snowrho !! Snow density |
---|
| 610 | REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT (inout) :: snowtemp !! Snow temperature (K) |
---|
| 611 | REAL(r_std), DIMENSION (kjpindex),INTENT (in) :: pb !! Surface presure (hPa) |
---|
| 612 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mc_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid + ice) (m3/m3) |
---|
| 613 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid) (m3/m3) |
---|
| 614 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content for each layer in hydrol(liquid + ice) (mm) |
---|
| 615 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) |
---|
| 616 | REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area |
---|
| 617 | REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area |
---|
| 618 | REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+... |
---|
| 619 | !!(unitless,0-1) |
---|
| 620 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: temp_sol_add !! additional surface temperature due to the melt of first layer |
---|
| 621 | !! at the present time-step @tex ($K$) @endtex |
---|
| 622 | |
---|
| 623 | !! 0.2 Output variables |
---|
| 624 | |
---|
| 625 | REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: ptnlev1 !! 1st level soil temperature |
---|
| 626 | REAL(r_std),DIMENSION (kjpindex),INTENT(out) :: gtemp !! First soil layer temperature |
---|
| 627 | |
---|
| 628 | |
---|
| 629 | !! 0.3 Modified variables |
---|
| 630 | |
---|
| 631 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilcap !! apparent surface heat capacity considering snow and soil surface |
---|
| 632 | !! @tex ($J m^{-2} K^{-1}$) @endtex |
---|
| 633 | REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: soilflx !! apparent soil heat flux considering snow and soil surface |
---|
| 634 | !! @tex ($W m^{-2}$) @endtex |
---|
| 635 | !! , positive |
---|
| 636 | !! towards the soil, writen as Qg (ground heat flux) |
---|
| 637 | !! in the history files, and computed at the end of |
---|
| 638 | !! thermosoil for the calculation of Ts in enerbil, |
---|
| 639 | !! see EQ3. |
---|
| 640 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! temperature profile @tex ($K$) @endtex |
---|
| 641 | REAL(r_std),DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature |
---|
| 642 | REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme |
---|
| 643 | REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme |
---|
| 644 | |
---|
| 645 | !! 0.4 Local variables |
---|
| 646 | |
---|
| 647 | INTEGER(i_std) :: jv,ji,ii |
---|
| 648 | REAL(r_std),DIMENSION (kjpindex) :: snowtemp_weighted!! Snow temperature weighted by snow density, only for diag (K) |
---|
| 649 | REAL(r_std),DIMENSION (kjpindex, nsnow) :: pkappa_snow_diag !! Only for diag, containing xios_default_val |
---|
| 650 | REAL(r_std),DIMENSION (kjpindex, nsnow) :: pcapa_snow_diag !! Only for diag, containing xios_default_val |
---|
| 651 | REAL(r_std),DIMENSION (kjpindex, nsnow) :: snowtemp_diag !! Only for diag, containing xios_default_val |
---|
| 652 | |
---|
| 653 | !_ ================================================================================================================================ |
---|
| 654 | |
---|
| 655 | !! 3. Put the soil wetness diagnostic on the levels of the soil temperature |
---|
| 656 | |
---|
| 657 | !!?? this could logically be put just before the last call to |
---|
| 658 | !!thermosoil_coef, as the results are used there... |
---|
| 659 | CALL thermosoil_humlev(kjpindex, shumdiag_perma, mc_layh, mcl_layh, tmc_layh) |
---|
| 660 | |
---|
| 661 | |
---|
| 662 | !! 4. Effective computation of the soil temperatures profile. |
---|
| 663 | !! cgrnd and dgrnd have been calculated in thermosoil_coef at the previous time step |
---|
| 664 | !! but they are correct for the actual time-step. |
---|
| 665 | CALL thermosoil_profile (kjpindex, temp_sol_new, & |
---|
| 666 | frac_snow_veg, frac_snow_nobio, totfrac_nobio, & |
---|
| 667 | ptn, stempdiag, snowtemp, & |
---|
| 668 | cgrnd_snow, dgrnd_snow) |
---|
| 669 | |
---|
| 670 | |
---|
| 671 | !! 5. Call to thermosoil_energy_diag for calculation of diagnostic variables |
---|
| 672 | CALL thermosoil_energy_diag(kjpindex, temp_sol_new, soilcap) |
---|
| 673 | |
---|
| 674 | !! Save ptn at current stage, to be used in thermosoil_readjust |
---|
| 675 | ptn_beg(:,:) = ptn(:,:) |
---|
| 676 | |
---|
| 677 | !! 6. Writing the history files according to the ALMA standards (or not..) |
---|
| 678 | |
---|
| 679 | ! Add XIOS default value where no snow |
---|
| 680 | DO ji=1,kjpindex |
---|
| 681 | IF (snow(ji) .GT. zero) THEN |
---|
| 682 | pkappa_snow_diag(ji,:) = pkappa_snow(ji,:) |
---|
| 683 | pcapa_snow_diag(ji,:) = pcapa_snow(ji,:) |
---|
| 684 | snowtemp_diag(ji,:) = snowtemp(ji,:) |
---|
| 685 | ELSE |
---|
| 686 | pkappa_snow_diag(ji,:) = xios_default_val |
---|
| 687 | pcapa_snow_diag(ji,:) = xios_default_val |
---|
| 688 | snowtemp_diag(ji,:) = xios_default_val |
---|
| 689 | END IF |
---|
| 690 | END DO |
---|
| 691 | |
---|
| 692 | DO ji=1,kjpindex |
---|
| 693 | ! Use min_sechiba instead of zero to avoid problem with division by zero |
---|
| 694 | IF (snow(ji) .GT. min_sechiba) THEN |
---|
| 695 | snowtemp_weighted(ji) = SUM(snowtemp(ji,:)*snowrho(ji,:))/SUM(snowrho(ji,:)) |
---|
| 696 | ELSE |
---|
| 697 | snowtemp_weighted(ji) = xios_default_val |
---|
| 698 | END IF |
---|
| 699 | END DO |
---|
| 700 | CALL xios_orchidee_send_field("snowtemp_weighted",snowtemp_weighted) |
---|
| 701 | |
---|
| 702 | CALL xios_orchidee_send_field("ptn",ptn) |
---|
| 703 | CALL xios_orchidee_send_field("soilflx",soilflx) |
---|
| 704 | CALL xios_orchidee_send_field("surfheat_incr",surfheat_incr) |
---|
| 705 | CALL xios_orchidee_send_field("coldcont_incr",coldcont_incr) |
---|
| 706 | CALL xios_orchidee_send_field("pkappa",pkappa) |
---|
| 707 | CALL xios_orchidee_send_field("pkappa_snow",pkappa_snow_diag) |
---|
| 708 | CALL xios_orchidee_send_field("pcapa",pcapa) |
---|
| 709 | CALL xios_orchidee_send_field("pcapa_snow",pcapa_snow_diag) |
---|
| 710 | CALL xios_orchidee_send_field("snowtemp",snowtemp_diag) |
---|
| 711 | |
---|
| 712 | |
---|
| 713 | IF ( .NOT. almaoutput ) THEN |
---|
| 714 | CALL histwrite_p(hist_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd) |
---|
| 715 | CALL histwrite_p(hist_id, 'Qg', kjit, soilflx, kjpindex, index) |
---|
| 716 | CALL histwrite_p(hist_id, 'ptn_beg', kjit, ptn_beg, kjpindex*ngrnd, indexgrnd) |
---|
| 717 | CALL histwrite_p(hist_id, 'pkappa', kjit, pkappa, kjpindex*ngrnd, indexgrnd) |
---|
| 718 | CALL histwrite_p(hist_id, 'pcapa', kjit, pcapa, kjpindex*ngrnd, indexgrnd) |
---|
| 719 | |
---|
| 720 | IF (ok_freeze_thermix) THEN |
---|
| 721 | CALL histwrite_p(hist_id, 'profil_froz', kjit, profil_froz, kjpindex*ngrnd, indexgrnd) |
---|
| 722 | CALL histwrite_p(hist_id, 'pcappa_supp', kjit, pcappa_supp, kjpindex*ngrnd, indexgrnd) |
---|
| 723 | END IF |
---|
| 724 | CALL histwrite_p(hist_id, 'shum_ngrnd_perma', kjit, shum_ngrnd_perma(:,:), kjpindex*ngrnd, indexgrnd) |
---|
| 725 | |
---|
| 726 | ELSE |
---|
| 727 | CALL histwrite_p(hist_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd) |
---|
| 728 | CALL histwrite_p(hist_id, 'Qg', kjit, soilflx, kjpindex, index) |
---|
| 729 | CALL histwrite_p(hist_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index) |
---|
| 730 | CALL histwrite_p(hist_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index) |
---|
| 731 | ENDIF |
---|
| 732 | IF ( hist2_id > 0 ) THEN |
---|
| 733 | IF ( .NOT. almaoutput ) THEN |
---|
| 734 | CALL histwrite_p(hist2_id, 'ptn', kjit, ptn, kjpindex*ngrnd, indexgrnd) |
---|
| 735 | ELSE |
---|
| 736 | CALL histwrite_p(hist2_id, 'SoilTemp', kjit, ptn, kjpindex*ngrnd, indexgrnd) |
---|
| 737 | CALL histwrite_p(hist2_id, 'Qg', kjit, soilflx, kjpindex, index) |
---|
| 738 | CALL histwrite_p(hist2_id, 'DelSurfHeat', kjit, surfheat_incr, kjpindex, index) |
---|
| 739 | CALL histwrite_p(hist2_id, 'DelColdCont', kjit, coldcont_incr, kjpindex, index) |
---|
| 740 | ENDIF |
---|
| 741 | ENDIF |
---|
| 742 | |
---|
| 743 | !! 7. A last final call to thermosoil_coef |
---|
| 744 | |
---|
| 745 | !! A last final call to thermosoil_coef, which calculates the different |
---|
| 746 | !!coefficients (cgrnd, dgrnd, soilcap, soilflx) from this time step to be |
---|
| 747 | !!used at the next time step, either in the surface temperature calculation |
---|
| 748 | !!(soilcap, soilflx) or in the soil thermal numerical scheme. |
---|
| 749 | CALL thermosoil_coef (& |
---|
| 750 | kjpindex, temp_sol_new, snow, njsc, & |
---|
| 751 | mcs, frac_snow_veg, frac_snow_nobio, totfrac_nobio, & |
---|
| 752 | snowdz, snowrho, snowtemp, pb, & |
---|
| 753 | ptn, & |
---|
| 754 | soilcap, soilflx, cgrnd, dgrnd,& |
---|
| 755 | lambda_snow, cgrnd_snow, dgrnd_snow) |
---|
| 756 | |
---|
| 757 | |
---|
| 758 | ! Save variables for explicit snow model |
---|
| 759 | gtemp(:) = ptn(:,1) |
---|
| 760 | |
---|
| 761 | !! Initialize output arguments to be used in sechiba |
---|
| 762 | ptnlev1(:) = ptn(:,1) |
---|
| 763 | |
---|
| 764 | !! Surface temperature is forced to zero celcius if its value is larger than melting point |
---|
| 765 | DO ji=1,kjpindex |
---|
| 766 | IF (SUM(snowdz(ji,:)) .GT. 0.0) THEN |
---|
| 767 | IF (temp_sol_new(ji) .GE. tp_00) THEN |
---|
| 768 | temp_sol_new(ji) = tp_00 |
---|
| 769 | ENDIF |
---|
| 770 | END IF |
---|
| 771 | END DO |
---|
| 772 | |
---|
| 773 | IF (printlev>=3) WRITE (numout,*) ' thermosoil_main done ' |
---|
| 774 | |
---|
| 775 | END SUBROUTINE thermosoil_main |
---|
| 776 | |
---|
| 777 | !! ============================================================================================================================= |
---|
| 778 | !! SUBROUTINE : thermosoil_finalize |
---|
| 779 | !! |
---|
| 780 | !>\BRIEF Write to restart file |
---|
| 781 | !! |
---|
| 782 | !! DESCRIPTION : This subroutine writes the module variables and variables calculated in thermosoil |
---|
| 783 | !! to restart file |
---|
| 784 | !! \n |
---|
| 785 | !_ ============================================================================================================================== |
---|
| 786 | SUBROUTINE thermosoil_finalize (kjit, kjpindex, rest_id, gtemp, & |
---|
| 787 | soilcap, soilflx, lambda_snow, cgrnd_snow, dgrnd_snow) |
---|
| 788 | |
---|
| 789 | !! 0. Variable and parameter declaration |
---|
| 790 | !! 0.1 Input variables |
---|
| 791 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless) |
---|
| 792 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
| 793 | INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier(unitless) |
---|
| 794 | REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: gtemp !! First soil layer temperature |
---|
| 795 | REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilcap |
---|
| 796 | REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: soilflx |
---|
| 797 | REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature |
---|
| 798 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme |
---|
| 799 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme |
---|
| 800 | |
---|
| 801 | !_ ================================================================================================================================ |
---|
| 802 | |
---|
| 803 | !! 1. Write variables to restart file to be used for the next simulation |
---|
| 804 | IF (printlev>=3) WRITE (numout,*) 'Write restart file with THERMOSOIL variables' |
---|
| 805 | |
---|
| 806 | CALL restput_p(rest_id, 'ptn', nbp_glo, ngrnd, 1, kjit, ptn, 'scatter', nbp_glo, index_g) |
---|
| 807 | |
---|
| 808 | IF (ok_Ecorr) THEN |
---|
| 809 | CALL restput_p(rest_id, 'e_soil_lat', nbp_glo, 1 , 1, kjit, e_soil_lat, 'scatter', nbp_glo, index_g) |
---|
| 810 | END IF |
---|
| 811 | |
---|
| 812 | CALL restput_p(rest_id, 'gtemp', nbp_glo, 1, 1, kjit, gtemp, 'scatter', nbp_glo, index_g) |
---|
| 813 | CALL restput_p(rest_id, 'soilcap', nbp_glo, 1, 1, kjit, soilcap, 'scatter', nbp_glo, index_g) |
---|
| 814 | CALL restput_p(rest_id, 'soilflx', nbp_glo, 1, 1, kjit, soilflx, 'scatter', nbp_glo, index_g) |
---|
| 815 | CALL restput_p(rest_id, 'cgrnd', nbp_glo, ngrnd-1, 1, kjit, cgrnd, 'scatter', nbp_glo, index_g) |
---|
| 816 | CALL restput_p(rest_id, 'dgrnd', nbp_glo, ngrnd-1, 1, kjit, dgrnd, 'scatter', nbp_glo, index_g) |
---|
| 817 | CALL restput_p(rest_id, 'cgrnd_snow', nbp_glo, nsnow, 1, kjit, cgrnd_snow, 'scatter', nbp_glo, index_g) |
---|
| 818 | CALL restput_p(rest_id, 'dgrnd_snow', nbp_glo, nsnow, 1, kjit, dgrnd_snow, 'scatter', nbp_glo, index_g) |
---|
| 819 | CALL restput_p(rest_id, 'lambda_snow', nbp_glo, 1, 1, kjit, lambda_snow, 'scatter', nbp_glo, index_g) |
---|
| 820 | |
---|
| 821 | END SUBROUTINE thermosoil_finalize |
---|
| 822 | |
---|
| 823 | |
---|
| 824 | !! ================================================================================================================================ |
---|
| 825 | !! SUBROUTINE : thermosoil_clear |
---|
| 826 | !! |
---|
| 827 | !>\BRIEF Deallocates the allocated arrays. |
---|
| 828 | !! The call of thermosoil_clear originates from sechiba_clear but the calling sequence and |
---|
| 829 | !! its purpose require further investigation. |
---|
| 830 | !! |
---|
| 831 | !! DESCRIPTION : None |
---|
| 832 | !! |
---|
| 833 | !! RECENT CHANGE(S) : None |
---|
| 834 | !! |
---|
| 835 | !! MAIN OUTPUT VARIABLE(S): None |
---|
| 836 | !! |
---|
| 837 | !! REFERENCE(S) : None |
---|
| 838 | !! |
---|
| 839 | !! FLOWCHART : None |
---|
| 840 | !! \n |
---|
| 841 | !_ ================================================================================================================================ |
---|
| 842 | |
---|
| 843 | SUBROUTINE thermosoil_clear() |
---|
| 844 | |
---|
| 845 | IF ( ALLOCATED (ptn)) DEALLOCATE (ptn) |
---|
| 846 | IF ( ALLOCATED (cgrnd)) DEALLOCATE (cgrnd) |
---|
| 847 | IF ( ALLOCATED (dgrnd)) DEALLOCATE (dgrnd) |
---|
| 848 | IF ( ALLOCATED (pcapa)) DEALLOCATE (pcapa) |
---|
| 849 | IF ( ALLOCATED (pkappa)) DEALLOCATE (pkappa) |
---|
| 850 | IF ( ALLOCATED (pcapa_snow)) DEALLOCATE (pcapa_snow) |
---|
| 851 | IF ( ALLOCATED (pkappa_snow)) DEALLOCATE (pkappa_snow) |
---|
| 852 | IF ( ALLOCATED (pcapa_en)) DEALLOCATE (pcapa_en) |
---|
| 853 | IF ( ALLOCATED (ptn_beg)) DEALLOCATE (ptn_beg) |
---|
| 854 | IF ( ALLOCATED (temp_sol_beg)) DEALLOCATE (temp_sol_beg) |
---|
| 855 | IF ( ALLOCATED (surfheat_incr)) DEALLOCATE (surfheat_incr) |
---|
| 856 | IF ( ALLOCATED (coldcont_incr)) DEALLOCATE (coldcont_incr) |
---|
| 857 | IF ( ALLOCATED (shum_ngrnd_perma)) DEALLOCATE (shum_ngrnd_perma) |
---|
| 858 | IF ( ALLOCATED (profil_froz)) DEALLOCATE (profil_froz) |
---|
| 859 | IF ( ALLOCATED (mc_layt)) DEALLOCATE (mc_layt) |
---|
| 860 | IF ( ALLOCATED (mcl_layt)) DEALLOCATE (mcl_layt) |
---|
| 861 | IF ( ALLOCATED (tmc_layt)) DEALLOCATE (tmc_layt) |
---|
| 862 | END SUBROUTINE thermosoil_clear |
---|
| 863 | |
---|
| 864 | |
---|
| 865 | !! ================================================================================================================================ |
---|
| 866 | !! SUBROUTINE : thermosoil_coef |
---|
| 867 | !! |
---|
| 868 | !>\BRIEF Calculate soil thermal properties, integration coefficients, apparent heat flux, |
---|
| 869 | !! surface heat capacity, |
---|
| 870 | !! |
---|
| 871 | !! DESCRIPTION : This routine computes : \n |
---|
| 872 | !! 1. the soil thermal properties. \n |
---|
| 873 | !! 2. the integration coefficients of the thermal numerical scheme, cgrnd and dgrnd, |
---|
| 874 | !! which depend on the vertical grid and on soil properties, and are used at the next |
---|
| 875 | !! timestep.\n |
---|
| 876 | !! 3. the soil apparent heat flux and surface heat capacity (soilflx |
---|
| 877 | !! and soilcap), used by enerbil to compute the surface temperature at the next |
---|
| 878 | !! timestep.\n |
---|
| 879 | !! - The soil thermal properties depend on water content (shum_ngrnd_perma, shumdiag_perma, |
---|
| 880 | !! mc_layt, mcl_layt, tmc_layt), dominant soil texture(njsc), and on the presence |
---|
| 881 | !! of snow : snow is integrated into the soil for the thermal calculations, ie if there |
---|
| 882 | !! is snow on the ground, the first thermal layer(s) consist in snow, depending on the |
---|
| 883 | !! snow-depth. If a layer consists out of snow and soil, wheighed soil properties are |
---|
| 884 | !! calculated\n |
---|
| 885 | !! - The coefficients cgrnd and dgrnd are the integration |
---|
| 886 | !! coefficients for the thermal scheme \n |
---|
| 887 | !! T(k+1)=cgrnd(k)+dgrnd(k)*T(k) \n |
---|
| 888 | !! -- EQ1 -- \n |
---|
| 889 | !! They correspond respectively to $\beta$ and $\alpha$ from F. Hourdin\'s thesis and |
---|
| 890 | !! their expression can be found in this document (eq A19 and A20) |
---|
| 891 | !! - soilcap and soilflx are the apparent surface heat capacity and flux |
---|
| 892 | !! used in enerbil at the next timestep to solve the surface |
---|
| 893 | !! balance for Ts (EQ3); they correspond to $C_s$ and $F_s$ in F. |
---|
| 894 | !! Hourdin\'s PhD thesis and are expressed in eq. A30 and A31. \n |
---|
| 895 | !! soilcap*(Ts(t)-Ts(t-1))/dt=soilflx+otherfluxes(Ts(t)) \n |
---|
| 896 | !! -- EQ3 --\n |
---|
| 897 | !! |
---|
| 898 | !! RECENT CHANGE(S) : None |
---|
| 899 | !! |
---|
| 900 | !! MAIN OUTPUT VARIABLE(S): cgrnd, dgrnd, pcapa, pkappa, soilcap, soilflx |
---|
| 901 | !! |
---|
| 902 | !! REFERENCE(S) : |
---|
| 903 | !! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres, |
---|
| 904 | !! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal |
---|
| 905 | !! integration scheme has been scanned and is provided along with the documentation, with name : |
---|
| 906 | !! Hourdin_1992_PhD_thermal_scheme.pdf |
---|
| 907 | !! |
---|
| 908 | !! FLOWCHART : None |
---|
| 909 | !! \n |
---|
| 910 | !_ ================================================================================================================================ |
---|
| 911 | |
---|
| 912 | SUBROUTINE thermosoil_coef (kjpindex, temp_sol_new, snow, njsc, & |
---|
| 913 | mcs, frac_snow_veg, frac_snow_nobio,totfrac_nobio, & |
---|
| 914 | snowdz, snowrho, snowtemp, pb, & |
---|
| 915 | ptn, & |
---|
| 916 | soilcap, soilflx, cgrnd, dgrnd,& |
---|
| 917 | lambda_snow, cgrnd_snow, dgrnd_snow) |
---|
| 918 | |
---|
| 919 | !! 0. Variables and parameter declaration |
---|
| 920 | |
---|
| 921 | !! 0.1 Input variables |
---|
| 922 | |
---|
| 923 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
| 924 | REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! soil surface temperature @tex ($K$) @endtex |
---|
| 925 | REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: snow !! snow mass @tex ($Kg$) @endtex |
---|
| 926 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class |
---|
| 927 | !! in the grid cell (1-nscm, unitless) |
---|
| 928 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3) |
---|
| 929 | REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area |
---|
| 930 | REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Snow cover fraction on non-vegeted area |
---|
| 931 | REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+... |
---|
| 932 | !!(unitless,0-1) |
---|
| 933 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow depth (m) |
---|
| 934 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density |
---|
| 935 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K) |
---|
| 936 | REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Surface presure (hPa) |
---|
| 937 | |
---|
| 938 | !! 0.2 Output variables |
---|
| 939 | |
---|
| 940 | REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilcap !! surface heat capacity considering snow and soil surface |
---|
| 941 | !! @tex ($J m^{-2} K^{-1}$) @endtex |
---|
| 942 | REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: soilflx !! surface heat flux considering snow and soil surface @tex ($W m^{-2}$) @endtex, |
---|
| 943 | !! positive towards the |
---|
| 944 | !! soil, writen as Qg (ground heat flux) in the history |
---|
| 945 | !! files. |
---|
| 946 | REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: cgrnd !! matrix coefficient for the computation of soil |
---|
| 947 | !! temperatures (beta in F. Hourdin thesis) |
---|
| 948 | REAL(r_std), DIMENSION (kjpindex,ngrnd-1), INTENT(out) :: dgrnd !! matrix coefficient for the computation of soil |
---|
| 949 | !! temperatures (alpha in F. Hourdin thesis) |
---|
| 950 | |
---|
| 951 | |
---|
| 952 | !! 0.3 Modified variable |
---|
| 953 | |
---|
| 954 | REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT (inout):: ptn !! vertically discretized soil temperatures. ptn is only modified if ok_Ecorr. |
---|
| 955 | REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature |
---|
| 956 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: cgrnd_snow !! Integration coefficient for snow numerical scheme |
---|
| 957 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT (inout):: dgrnd_snow !! Integration coefficient for snow numerical scheme |
---|
| 958 | |
---|
| 959 | !! 0.4 Local variables |
---|
| 960 | |
---|
| 961 | INTEGER(i_std) :: ji, jg |
---|
| 962 | REAL(r_std), DIMENSION (kjpindex,ngrnd-1) :: zdz1 !! numerical (buffer) constant |
---|
| 963 | !! @tex ($W m^{-1} K^{-1}$) @endtex |
---|
| 964 | REAL(r_std), DIMENSION (kjpindex,ngrnd) :: zdz2 !! numerical (buffer) constant |
---|
| 965 | !! @tex ($W m^{-1} K^{-1}$) @endtex |
---|
| 966 | REAL(r_std), DIMENSION (kjpindex) :: z1 !! numerical constant @tex ($W m^{-1} K^{-1}$) @endtex |
---|
| 967 | REAL(r_std), DIMENSION (kjpindex) :: soilcap_nosnow !! surface heat capacity |
---|
| 968 | !! @tex ($J m^{-2} K^{-1}$) |
---|
| 969 | !! @endtex |
---|
| 970 | REAL(r_std), DIMENSION (kjpindex) :: soilflx_nosnow !! surface heat flux @tex ($W m^{-2}$) @endtex, |
---|
| 971 | !! positive towards the soil, written as Qg |
---|
| 972 | !!(ground heat flux in the history files). |
---|
| 973 | REAL(r_std), DIMENSION (kjpindex) :: snowcap !! apparent snow heat capacity @tex ($J m^{-2} K^{-1}$) |
---|
| 974 | REAL(r_std), DIMENSION (kjpindex) :: snowflx !! apparent snow-atmosphere heat flux @tex ($W m^{-2}$) @endtex |
---|
| 975 | REAL(r_std), DIMENSION (kjpindex,nsnow) :: dz1_snow |
---|
| 976 | REAL(r_std), DIMENSION (kjpindex,nsnow) :: ZSNOWDZM |
---|
| 977 | REAL(r_std), DIMENSION (kjpindex,nsnow) :: dz2_snow |
---|
| 978 | REAL(r_std), DIMENSION (kjpindex,nsnow) :: zdz1_snow |
---|
| 979 | REAL(r_std), DIMENSION (kjpindex,nsnow) :: zdz2_snow |
---|
| 980 | REAL(r_std), DIMENSION (kjpindex) :: z1_snow |
---|
| 981 | REAL(r_std), DIMENSION (kjpindex) :: snowflxtot !! Total snow flux (including snow on vegetated and bare soil and nobio areas) |
---|
| 982 | !! @tex ($W m^{-2}$) @endtex |
---|
| 983 | !! positive towards the soil |
---|
| 984 | |
---|
| 985 | !_ ================================================================================================================================ |
---|
| 986 | |
---|
| 987 | !! 1. Computation of the soil thermal properties |
---|
| 988 | |
---|
| 989 | ! Computation of the soil thermal properties; snow properties are also accounted for |
---|
| 990 | IF (ok_freeze_thermix) THEN |
---|
| 991 | CALL thermosoil_getdiff( kjpindex, snow, ptn, mcs, njsc, snowrho, snowtemp, pb ) |
---|
| 992 | ELSE |
---|
| 993 | ! Special case without soil freezing |
---|
| 994 | CALL thermosoil_getdiff_old_thermix_without_snow( kjpindex, mcs, njsc, snowrho, snowtemp, pb ) |
---|
| 995 | ENDIF |
---|
| 996 | |
---|
| 997 | ! Energy conservation : Correction to make sure that the same latent heat is released and |
---|
| 998 | ! consumed during freezing and thawing |
---|
| 999 | IF (ok_Ecorr) THEN |
---|
| 1000 | CALL thermosoil_readjust(kjpindex, ptn) |
---|
| 1001 | ENDIF |
---|
| 1002 | |
---|
| 1003 | |
---|
| 1004 | !! 2. Computation of the coefficients of the numerical integration scheme for the soil layers |
---|
| 1005 | |
---|
| 1006 | !! 2.1 Calculate numerical coefficients zdz1 and zdz2 |
---|
| 1007 | DO jg=1,ngrnd |
---|
| 1008 | DO ji=1,kjpindex |
---|
| 1009 | zdz2(ji,jg)=pcapa(ji,jg) * dlt(jg)/dt_sechiba |
---|
| 1010 | ENDDO |
---|
| 1011 | ENDDO |
---|
| 1012 | |
---|
| 1013 | DO jg=1,ngrnd-1 |
---|
| 1014 | DO ji=1,kjpindex |
---|
| 1015 | zdz1(ji,jg) = dz1(jg) * pkappa(ji,jg) |
---|
| 1016 | ENDDO |
---|
| 1017 | ENDDO |
---|
| 1018 | |
---|
| 1019 | !! 2.2 Calculate coefficients cgrnd and dgrnd for soil |
---|
| 1020 | DO ji = 1,kjpindex |
---|
| 1021 | z1(ji) = zdz2(ji,ngrnd) + zdz1(ji,ngrnd-1) |
---|
| 1022 | cgrnd(ji,ngrnd-1) = zdz2(ji,ngrnd) * ptn(ji,ngrnd) / z1(ji) |
---|
| 1023 | dgrnd(ji,ngrnd-1) = zdz1(ji,ngrnd-1) / z1(ji) |
---|
| 1024 | ENDDO |
---|
| 1025 | |
---|
| 1026 | DO jg = ngrnd-1,2,-1 |
---|
| 1027 | DO ji = 1,kjpindex |
---|
| 1028 | z1(ji) = un / (zdz2(ji,jg) + zdz1(ji,jg-1) + zdz1(ji,jg) * (un - dgrnd(ji,jg))) |
---|
| 1029 | cgrnd(ji,jg-1) = (ptn(ji,jg) * zdz2(ji,jg) + zdz1(ji,jg) * cgrnd(ji,jg)) * z1(ji) |
---|
| 1030 | dgrnd(ji,jg-1) = zdz1(ji,jg-1) * z1(ji) |
---|
| 1031 | ENDDO |
---|
| 1032 | ENDDO |
---|
| 1033 | |
---|
| 1034 | |
---|
| 1035 | !! 3. Computation of the coefficients of the numerical integration scheme for the snow layers |
---|
| 1036 | |
---|
| 1037 | !! 3.1 Calculate numerical coefficients zdz1_snow, zdz2_snow and lambda_snow |
---|
| 1038 | DO ji = 1, kjpindex |
---|
| 1039 | |
---|
| 1040 | ! Calculate internal values |
---|
| 1041 | DO jg = 1, nsnow |
---|
| 1042 | ZSNOWDZM(ji,jg) = MAX(snowdz(ji,jg),psnowdzmin) |
---|
| 1043 | ENDDO |
---|
| 1044 | dz2_snow(ji,:)=ZSNOWDZM(ji,:) |
---|
| 1045 | |
---|
| 1046 | DO jg = 1, nsnow-1 |
---|
| 1047 | dz1_snow(ji,jg) = 2.0 / (dz2_snow(ji,jg+1)+dz2_snow(ji,jg)) |
---|
| 1048 | ENDDO |
---|
| 1049 | |
---|
| 1050 | lambda_snow(ji) = dz2_snow(ji,1)/2.0 * dz1_snow(ji,1) |
---|
| 1051 | |
---|
| 1052 | DO jg=1,nsnow |
---|
| 1053 | zdz2_snow(ji,jg)=pcapa_snow(ji,jg) * dz2_snow(ji,jg)/dt_sechiba |
---|
| 1054 | ENDDO |
---|
| 1055 | |
---|
| 1056 | DO jg=1,nsnow-1 |
---|
| 1057 | zdz1_snow(ji,jg) = dz1_snow(ji,jg) * pkappa_snow(ji,jg) |
---|
| 1058 | ENDDO |
---|
| 1059 | |
---|
| 1060 | ! the bottom snow |
---|
| 1061 | zdz1_snow(ji,nsnow) = pkappa_snow(ji,nsnow) / ( zlt(1) + dz2_snow(ji,nsnow)/2 ) |
---|
| 1062 | |
---|
| 1063 | ENDDO |
---|
| 1064 | |
---|
| 1065 | !! 3.2 Calculate coefficients cgrnd_snow and dgrnd_snow for snow |
---|
| 1066 | DO ji = 1,kjpindex |
---|
| 1067 | ! bottom level |
---|
| 1068 | z1_snow(ji) = zdz2(ji,1)+(un-dgrnd(ji,1))*zdz1(ji,1)+zdz1_snow(ji,nsnow) |
---|
| 1069 | cgrnd_snow(ji,nsnow) = (zdz2(ji,1) * ptn(ji,1) + zdz1(ji,1) * cgrnd(ji,1) ) / z1_snow(ji) |
---|
| 1070 | dgrnd_snow(ji,nsnow) = zdz1_snow(ji,nsnow) / z1_snow(ji) |
---|
| 1071 | |
---|
| 1072 | ! next-to-bottom level |
---|
| 1073 | z1_snow(ji) = zdz2_snow(ji,nsnow)+(un-dgrnd_snow(ji,nsnow))*zdz1_snow(ji,nsnow)+zdz1_snow(ji,nsnow-1) |
---|
| 1074 | cgrnd_snow(ji,nsnow-1) = (zdz2_snow(ji,nsnow)*snowtemp(ji,nsnow)+& |
---|
| 1075 | zdz1_snow(ji,nsnow)*cgrnd_snow(ji,nsnow))/z1_snow(ji) |
---|
| 1076 | dgrnd_snow(ji,nsnow-1) = zdz1_snow(ji,nsnow-1) / z1_snow(ji) |
---|
| 1077 | |
---|
| 1078 | DO jg = nsnow-1,2,-1 |
---|
| 1079 | z1_snow(ji) = un / (zdz2_snow(ji,jg) + zdz1_snow(ji,jg-1) + zdz1_snow(ji,jg) * (un - dgrnd_snow(ji,jg))) |
---|
| 1080 | cgrnd_snow(ji,jg-1) = (snowtemp(ji,jg) * zdz2_snow(ji,jg) + zdz1_snow(ji,jg) * cgrnd_snow(ji,jg)) * z1_snow(ji) |
---|
| 1081 | dgrnd_snow(ji,jg-1) = zdz1_snow(ji,jg-1) * z1_snow(ji) |
---|
| 1082 | ENDDO |
---|
| 1083 | ENDDO |
---|
| 1084 | |
---|
| 1085 | |
---|
| 1086 | |
---|
| 1087 | !! 4. Computation of the apparent ground heat flux |
---|
| 1088 | !! Computation of apparent snow-atmosphere flux |
---|
| 1089 | DO ji = 1,kjpindex |
---|
| 1090 | snowflx(ji) = zdz1_snow(ji,1) * (cgrnd_snow(ji,1) + (dgrnd_snow(ji,1)-1.) * snowtemp(ji,1)) |
---|
| 1091 | snowcap(ji) = (zdz2_snow(ji,1) * dt_sechiba + dt_sechiba * (un - dgrnd_snow(ji,1)) * zdz1_snow(ji,1)) |
---|
| 1092 | z1_snow(ji) = lambda_snow(ji) * (un - dgrnd_snow(ji,1)) + un |
---|
| 1093 | snowcap(ji) = snowcap(ji) / z1_snow(ji) |
---|
| 1094 | snowflx(ji) = snowflx(ji) + & |
---|
| 1095 | snowcap(ji) * (snowtemp(ji,1) * z1_snow(ji) - lambda_snow(ji) * cgrnd_snow(ji,1) - temp_sol_new(ji)) / dt_sechiba |
---|
| 1096 | ENDDO |
---|
| 1097 | |
---|
| 1098 | |
---|
| 1099 | !! Computation of the apparent ground heat flux (> towards the soil) and |
---|
| 1100 | !! apparent surface heat capacity, used at the next timestep by enerbil to |
---|
| 1101 | !! compute the surface temperature. |
---|
| 1102 | DO ji = 1,kjpindex |
---|
| 1103 | soilflx_nosnow(ji) = zdz1(ji,1) * (cgrnd(ji,1) + (dgrnd(ji,1)-1.) * ptn(ji,1)) |
---|
| 1104 | soilcap_nosnow(ji) = (zdz2(ji,1) * dt_sechiba + dt_sechiba * (un - dgrnd(ji,1)) * zdz1(ji,1)) |
---|
| 1105 | z1(ji) = lambda * (un - dgrnd(ji,1)) + un |
---|
| 1106 | soilcap_nosnow(ji) = soilcap_nosnow(ji) / z1(ji) |
---|
| 1107 | soilflx_nosnow(ji) = soilflx_nosnow(ji) + & |
---|
| 1108 | & soilcap_nosnow(ji) * (ptn(ji,1) * z1(ji) - lambda * cgrnd(ji,1) - temp_sol_new(ji)) / dt_sechiba |
---|
| 1109 | ENDDO |
---|
| 1110 | |
---|
| 1111 | !! Add snow fraction |
---|
| 1112 | ! Using an effective heat capacity and heat flux by a simple pondering of snow and soil fraction |
---|
| 1113 | DO ji = 1, kjpindex |
---|
| 1114 | soilcap(ji) = snowcap(ji)*frac_snow_veg(ji)*(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation |
---|
| 1115 | soilcap_nosnow(ji)*SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji)+ & ! weights related to SCF on nobio |
---|
| 1116 | soilcap_nosnow(ji)*(1-(frac_snow_veg(ji)*(1-totfrac_nobio(ji))+SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji))) ! weights related to non snow fraction |
---|
| 1117 | soilflx(ji) = snowflx(ji)*frac_snow_veg(ji)*(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation |
---|
| 1118 | soilflx_nosnow(ji)*SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji)+ & ! weights related to SCF on nobio |
---|
| 1119 | soilflx_nosnow(ji)*(1-(frac_snow_veg(ji)*(1-totfrac_nobio(ji))+SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji))) ! weights related to non snow fraction |
---|
| 1120 | ENDDO |
---|
| 1121 | |
---|
| 1122 | ! Total snow flux (including snow on vegetated and bare soil and nobio areas) |
---|
| 1123 | DO ji = 1, kjpindex |
---|
| 1124 | snowflxtot(ji) = snowflx(ji)*frac_snow_veg(ji)*(1-totfrac_nobio(ji)) + & |
---|
| 1125 | soilflx_nosnow(ji)*SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji) |
---|
| 1126 | ENDDO |
---|
| 1127 | CALL xios_orchidee_send_field("snowflxtot",snowflxtot(:)) |
---|
| 1128 | |
---|
| 1129 | IF (printlev>=3) WRITE (numout,*) ' thermosoil_coef done ' |
---|
| 1130 | |
---|
| 1131 | END SUBROUTINE thermosoil_coef |
---|
| 1132 | |
---|
| 1133 | |
---|
| 1134 | !! ================================================================================================================================ |
---|
| 1135 | !! SUBROUTINE : thermosoil_profile |
---|
| 1136 | !! |
---|
| 1137 | !>\BRIEF In this routine solves the numerical soil thermal scheme, ie calculates the new soil temperature profile; |
---|
| 1138 | !! |
---|
| 1139 | !! |
---|
| 1140 | !! DESCRIPTION : The calculation of the new soil temperature profile is based on |
---|
| 1141 | !! the cgrnd and dgrnd values from the previous timestep and the surface temperature Ts aka temp_sol_new. (see detailed |
---|
| 1142 | !! explanation in the header of the thermosoil module or in the reference).\n |
---|
| 1143 | !! T(k+1)=cgrnd(k)+dgrnd(k)*T(k)\n |
---|
| 1144 | !! -- EQ1 --\n |
---|
| 1145 | !! Ts=(1+lambda)*T(1) -lambda*T(2)\n |
---|
| 1146 | !! -- EQ2--\n |
---|
| 1147 | !! |
---|
| 1148 | !! RECENT CHANGE(S) : None |
---|
| 1149 | !! |
---|
| 1150 | !! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis), |
---|
| 1151 | !! stempdiag (soil temperature profile on the diagnostic axis) |
---|
| 1152 | !! |
---|
| 1153 | !! REFERENCE(S) : |
---|
| 1154 | !! - Hourdin, F. (1992). Study and numerical simulation of the general circulation of planetary atmospheres, |
---|
| 1155 | !! Ph.D. thesis, Paris VII University. Remark: the part of F. Hourdin's PhD thesis relative to the thermal |
---|
| 1156 | !! integration scheme has been scanned and is provided along with the documentation, with name : |
---|
| 1157 | !! Hourdin_1992_PhD_thermal_scheme.pdf |
---|
| 1158 | !! |
---|
| 1159 | !! FLOWCHART : None |
---|
| 1160 | !! \n |
---|
| 1161 | !_ ================================================================================================================================ |
---|
| 1162 | |
---|
| 1163 | SUBROUTINE thermosoil_profile (kjpindex, temp_sol_new, & |
---|
| 1164 | frac_snow_veg, frac_snow_nobio, totfrac_nobio, & |
---|
| 1165 | ptn, stempdiag, snowtemp, & |
---|
| 1166 | cgrnd_snow, dgrnd_snow) |
---|
| 1167 | |
---|
| 1168 | !! 0. Variables and parameter declaration |
---|
| 1169 | |
---|
| 1170 | !! 0.1 Input variables |
---|
| 1171 | |
---|
| 1172 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
| 1173 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step |
---|
| 1174 | !! @tex ($K$) @endtex |
---|
| 1175 | REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area |
---|
| 1176 | REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio!! Snow cover fraction on non-vegeted area |
---|
| 1177 | REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+... |
---|
| 1178 | !! (unitless,0-1) |
---|
| 1179 | REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K) |
---|
| 1180 | REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme |
---|
| 1181 | REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT(in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme |
---|
| 1182 | |
---|
| 1183 | !! 0.3 Modified variables |
---|
| 1184 | |
---|
| 1185 | |
---|
| 1186 | !! 0.2 Output variables |
---|
| 1187 | REAL(r_std),DIMENSION (kjpindex,ngrnd), INTENT (out) :: ptn !! vertically discretized soil temperatures |
---|
| 1188 | !! @tex ($K$) @endtex |
---|
| 1189 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: stempdiag !! diagnostic temperature profile |
---|
| 1190 | !! @tex ($K$) @endtex |
---|
| 1191 | |
---|
| 1192 | !! 0.4 Local variables |
---|
| 1193 | |
---|
| 1194 | INTEGER(i_std) :: ji, jg |
---|
| 1195 | REAL(r_std) :: temp_sol_eff !! effective surface temperature including snow and soil |
---|
| 1196 | |
---|
| 1197 | !_ ================================================================================================================================ |
---|
| 1198 | |
---|
| 1199 | !! 1. Computes the soil temperatures ptn. |
---|
| 1200 | |
---|
| 1201 | !! 1.1. ptn(jg=1) using EQ1 and EQ2 |
---|
| 1202 | DO ji = 1,kjpindex |
---|
| 1203 | |
---|
| 1204 | ! Using an effective surface temperature by a simple pondering |
---|
| 1205 | temp_sol_eff=snowtemp(ji,nsnow)*frac_snow_veg(ji)*(1-totfrac_nobio(ji))+ & ! weights related to snow cover fraction on vegetation |
---|
| 1206 | temp_sol_new(ji)*SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji)+ & ! weights related to SCF on nobio |
---|
| 1207 | temp_sol_new(ji)*(1-(frac_snow_veg(ji)*(1-totfrac_nobio(ji))+SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji))) ! weights related to non snow fraction |
---|
| 1208 | ! Soil temperature calculation with explicit snow if there is snow on the ground |
---|
| 1209 | ptn(ji,1) = cgrnd_snow(ji,nsnow) + dgrnd_snow(ji,nsnow) * temp_sol_eff |
---|
| 1210 | ENDDO |
---|
| 1211 | |
---|
| 1212 | !! 1.2. ptn(jg=2:ngrnd) using EQ1. |
---|
| 1213 | DO jg = 1,ngrnd-1 |
---|
| 1214 | DO ji = 1,kjpindex |
---|
| 1215 | ptn(ji,jg+1) = cgrnd(ji,jg) + dgrnd(ji,jg) * ptn(ji,jg) |
---|
| 1216 | ENDDO |
---|
| 1217 | ENDDO |
---|
| 1218 | |
---|
| 1219 | !! 2. Assigne the soil temperature to the output variable. It is already on the right axis. |
---|
| 1220 | stempdiag(:,:) = ptn(:,1:nslm) |
---|
| 1221 | |
---|
| 1222 | IF (printlev>=3) WRITE (numout,*) ' thermosoil_profile done ' |
---|
| 1223 | |
---|
| 1224 | END SUBROUTINE thermosoil_profile |
---|
| 1225 | |
---|
| 1226 | !================================================================================================================================ |
---|
| 1227 | !! SUBROUTINE : thermosoil_cond |
---|
| 1228 | !! |
---|
| 1229 | !>\BRIEF Calculate soil thermal conductivity. |
---|
| 1230 | !! |
---|
| 1231 | !! DESCRIPTION : This routine computes soil thermal conductivity |
---|
| 1232 | !! Code introduced from NOAH LSM. |
---|
| 1233 | !! |
---|
| 1234 | !! RECENT CHANGE(S) : None |
---|
| 1235 | !! |
---|
| 1236 | !! MAIN OUTPUT VARIABLE(S): cnd |
---|
| 1237 | !! |
---|
| 1238 | !! REFERENCE(S) : |
---|
| 1239 | !! Farouki, O.T.,1986: Thermal Properties of Soils. Series on Rock |
---|
| 1240 | !! and Soil Mechanics, Vol. 11, Trans Tech, 136 PP. |
---|
| 1241 | !! Johansen, O., 1975: Thermal Conductivity of Soils. Ph.D. Thesis, |
---|
| 1242 | !! University of Trondheim, |
---|
| 1243 | !! Peters-Lidard, C. D., Blackburn, E., Liang, X., & Wood, E. F., |
---|
| 1244 | !! 1998: The effect of soil thermal conductivity |
---|
| 1245 | !! Parameterization on Surface Energy fluxes |
---|
| 1246 | !! and Temperatures. J. of The Atmospheric Sciences, |
---|
| 1247 | !! Vol. 55, pp. 1209-1224. |
---|
| 1248 | !! Modify histroy: |
---|
| 1249 | !! |
---|
| 1250 | !! FLOWCHART : None |
---|
| 1251 | !! \n |
---|
| 1252 | !_ |
---|
| 1253 | !================================================================================================================================ |
---|
| 1254 | |
---|
| 1255 | SUBROUTINE thermosoil_cond (kjpindex, njsc, mcs, smc, qz, sh2o, cnd) |
---|
| 1256 | |
---|
| 1257 | !! 0. Variables and parameter declaration |
---|
| 1258 | |
---|
| 1259 | !! 0.1 Input variables |
---|
| 1260 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
| 1261 | INTEGER(i_std), DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) |
---|
| 1262 | REAL(r_std), DIMENSION (kjpindex), INTENT(IN) :: mcs !! Saturated moisture content (m3/m3) |
---|
| 1263 | REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: smc !! Volumetric Soil Moisture Content (m3/m3) |
---|
| 1264 | REAL(r_std), DIMENSION (nscm), INTENT(IN) :: qz !! Quartz Content (Soil Type Dependent) (0-1) |
---|
| 1265 | REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(IN) :: sh2o !! Unfrozen Soil Moisture Content; Frozen Soil Moisture = smc - sh2o |
---|
| 1266 | |
---|
| 1267 | !! 0.2 Output variables |
---|
| 1268 | REAL(r_std), DIMENSION (kjpindex,ngrnd), INTENT(OUT) :: cnd !! Soil Thermal Conductivity (W/m/k) |
---|
| 1269 | |
---|
| 1270 | !! 0.3 Modified variables |
---|
| 1271 | |
---|
| 1272 | !! 0.4 Local variables |
---|
| 1273 | REAL(r_std), DIMENSION (kjpindex,ngrnd) :: ake !! Kersten Number (unitless) |
---|
| 1274 | REAL(r_std), DIMENSION (kjpindex,ngrnd) :: thksat !! Saturated Thermal Conductivity (W/m/k) |
---|
| 1275 | REAL(r_std), DIMENSION (kjpindex,ngrnd) :: satratio !! Degree of Saturation (0-1) |
---|
| 1276 | REAL(r_std), DIMENSION (kjpindex,ngrnd) :: xu !! Unfrozen Volume For Saturation (0-1) |
---|
| 1277 | REAL(r_std), DIMENSION (kjpindex,ngrnd) :: xunfroz !! Unfrozon Volume Fraction (0-1) |
---|
| 1278 | REAL(r_std) :: thko !! Thermal Conductivity for Other Ssoil Components (W/m/k) |
---|
| 1279 | REAL(r_std) :: gammd !! Dry Dendity (kg/m3) |
---|
| 1280 | REAL(r_std) :: thkdry !! Dry Thermal Conductivity (W/m/k) |
---|
| 1281 | REAL(r_std) :: thks !! Thermal Conductivity for the Solids Combined (Quartz + Other) (W/m/k) |
---|
| 1282 | REAL(r_std), PARAMETER :: THKICE = 2.2 !! Ice Thermal Conductivity (W/m/k) |
---|
| 1283 | REAL(r_std), PARAMETER :: THKQTZ = 7.7 !! Thermal Conductivity for Quartz (W/m/k) |
---|
| 1284 | REAL(r_std), PARAMETER :: THKW = 0.57 !! Water Thermal Conductivity (W/m/k) |
---|
| 1285 | INTEGER(i_std) :: ji, jg, jst |
---|
| 1286 | |
---|
| 1287 | !_================================================================================================================================ |
---|
| 1288 | |
---|
| 1289 | !! 1. Dry and Saturated Thermal Conductivity. |
---|
| 1290 | |
---|
| 1291 | DO ji = 1,kjpindex |
---|
| 1292 | jst = njsc(ji) |
---|
| 1293 | |
---|
| 1294 | !! 1.1. Dry density (Kg/m3) and Dry thermal conductivity (W.M-1.K-1) |
---|
| 1295 | gammd = (1. - mcs(ji))*2700. |
---|
| 1296 | thkdry = (0.135* gammd+ 64.7)/ (2700. - 0.947* gammd) |
---|
| 1297 | |
---|
| 1298 | !! 1.2. thermal conductivity of "other" soil components |
---|
| 1299 | IF (qz(jst) > 0.2) THEN |
---|
| 1300 | thko = 2.0 |
---|
| 1301 | ELSEIF (qz(jst) <= 0.2) THEN |
---|
| 1302 | thko = 3.0 |
---|
| 1303 | ENDIF |
---|
| 1304 | |
---|
| 1305 | !! 1.3. Thermal conductivity of solids |
---|
| 1306 | thks = (THKQTZ ** qz(jst))* (thko ** (1. - qz(jst))) |
---|
| 1307 | |
---|
| 1308 | DO jg = 1,ngrnd |
---|
| 1309 | !! 1.4. saturation ratio |
---|
| 1310 | satratio(ji,jg) = smc(ji,jg) / mcs(ji) |
---|
| 1311 | |
---|
| 1312 | !! 1.5. Saturated Thermal Conductivity (thksat) |
---|
| 1313 | IF ( smc(ji,jg) > min_sechiba ) THEN |
---|
| 1314 | xunfroz(ji,jg) = sh2o(ji,jg) / smc(ji,jg) ! Unfrozen Fraction (From i.e., 100%Liquid, to 0. (100% Frozen)) |
---|
| 1315 | xu(ji,jg) = xunfroz(ji,jg) * mcs(ji) ! Unfrozen volume for saturation (porosity*xunfroz) |
---|
| 1316 | thksat(ji,jg) = thks ** (1. - mcs(ji))* THKICE ** (mcs(ji) - xu(ji,jg))* THKW ** (xu(ji,jg)) |
---|
| 1317 | ELSE |
---|
| 1318 | ! this value will not be used since ake=0 for this case |
---|
| 1319 | thksat(ji,jg)=0 |
---|
| 1320 | END IF |
---|
| 1321 | END DO ! DO jg = 1,ngrnd |
---|
| 1322 | |
---|
| 1323 | !! 2. Kersten Number (ake) |
---|
| 1324 | DO jg = 1,ngrnd |
---|
| 1325 | IF ( (sh2o(ji,jg) + 0.0005) < smc(ji,jg) ) THEN |
---|
| 1326 | ! Frozen |
---|
| 1327 | ake(ji,jg) = satratio(ji,jg) |
---|
| 1328 | ELSE |
---|
| 1329 | ! Unfrozen |
---|
| 1330 | ! Eq 11 in Peters-Lidard et al., 1998 |
---|
| 1331 | IF ( satratio(ji,jg) > 0.1 ) THEN |
---|
| 1332 | IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') ) THEN |
---|
| 1333 | ! Coarse |
---|
| 1334 | ake(ji,jg) = 0.7 * LOG10 (SATRATIO(ji,jg)) + 1.0 |
---|
| 1335 | ELSE |
---|
| 1336 | ! Fine |
---|
| 1337 | ake(ji,jg) = LOG10 (satratio(ji,jg)) + 1.0 |
---|
| 1338 | ENDIF |
---|
| 1339 | ELSEIF ( satratio(ji,jg) > 0.05 .AND. satratio(ji,jg) <= 0.1 ) THEN |
---|
| 1340 | IF ((jst < 4 .AND. soil_classif == 'usda') .OR. (jst == 1 .AND. soil_classif == 'zobler') ) THEN |
---|
| 1341 | ! Coarse |
---|
| 1342 | ake(ji,jg) = 0.7 * LOG10 (satratio(ji,jg)) + 1.0 |
---|
| 1343 | ELSE |
---|
| 1344 | ! Fine |
---|
| 1345 | ake(ji,jg) = 0.0 |
---|
| 1346 | ENDIF |
---|
| 1347 | ELSE |
---|
| 1348 | ake(ji,jg) = 0.0 ! use k = kdry |
---|
| 1349 | END IF |
---|
| 1350 | END IF |
---|
| 1351 | END DO ! DO jg = 1,ngrnd |
---|
| 1352 | |
---|
| 1353 | !! 3. Thermal conductivity (cnd) |
---|
| 1354 | DO jg = 1,ngrnd |
---|
| 1355 | cnd(ji,jg) = ake(ji,jg) * (thksat(ji,jg) - thkdry) + thkdry |
---|
| 1356 | END DO ! DO jg = 1,ngrnd |
---|
| 1357 | |
---|
| 1358 | END DO !DO ji = 1,kjpindex |
---|
| 1359 | |
---|
| 1360 | END SUBROUTINE thermosoil_cond |
---|
| 1361 | |
---|
| 1362 | |
---|
| 1363 | !! ================================================================================================================================ |
---|
| 1364 | !! SUBROUTINE : thermosoil_humlev |
---|
| 1365 | !! |
---|
| 1366 | !>\BRIEF Interpolate variables from the hydrology layers to the thermodynamic layers |
---|
| 1367 | !! |
---|
| 1368 | !! DESCRIPTION : Interpolate the volumetric soil moisture content from the node to the interface of the layer. |
---|
| 1369 | !! The values for the deep layers in thermosoil where hydrology is not existing are constant. |
---|
| 1370 | !! No interpolation is needed for the total soil moisture content and for the soil saturation degree. |
---|
| 1371 | !! |
---|
| 1372 | !! RECENT CHANGE(S) : None |
---|
| 1373 | !! |
---|
| 1374 | !! MAIN OUTPUT VARIABLE(S): mc_layt, mcl_layt, tmc_layt, shum_ngrnd_perma |
---|
| 1375 | !! |
---|
| 1376 | !! REFERENCE(S) : None |
---|
| 1377 | !! |
---|
| 1378 | !! FLOWCHART : None |
---|
| 1379 | !! \n |
---|
| 1380 | !_ ================================================================================================================================ |
---|
| 1381 | SUBROUTINE thermosoil_humlev(kjpindex, shumdiag_perma, mc_layh, mcl_layh, tmc_layh) |
---|
| 1382 | |
---|
| 1383 | !! 0. Variables and parameter declaration |
---|
| 1384 | |
---|
| 1385 | !! 0.1 Input variables |
---|
| 1386 | |
---|
| 1387 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
| 1388 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag_perma !! Soil saturation degree on the diagnostic axis (0-1, unitless) |
---|
| 1389 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mc_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid+ice) [m/s] |
---|
| 1390 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: mcl_layh !! Volumetric soil moisture content for each layer in hydrol at nodes(liquid) [m/s] |
---|
| 1391 | REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: tmc_layh !! Total soil moisture content for each layer in hydrol(liquid+ice) [mm] |
---|
| 1392 | |
---|
| 1393 | !! 0.2 Output variables |
---|
| 1394 | |
---|
| 1395 | !! 0.3 Modified variables |
---|
| 1396 | |
---|
| 1397 | !! 0.4 Local variables |
---|
| 1398 | INTEGER(i_std) :: ji, jd |
---|
| 1399 | |
---|
| 1400 | !_ ================================================================================================================================ |
---|
| 1401 | |
---|
| 1402 | IF (printlev >= 4) WRITE(numout,*) 'Start thermosoil_humlev' |
---|
| 1403 | |
---|
| 1404 | ! The values for the deep layers in thermosoil where hydrology is not existing are constant. |
---|
| 1405 | ! For exemple if thermosoil uses 8m, and hydrol uses 2m vertical discretization, |
---|
| 1406 | ! the values between 2m and 8m are constant. |
---|
| 1407 | ! The moisture computed in hydrol is at the nodes (except for the |
---|
| 1408 | ! top and bottom layer which are at interfaces) |
---|
| 1409 | ! A linear interpolation is applied to obtain the moisture values at |
---|
| 1410 | ! the interfaces (mc_layt), from the mc_layh at the nodes |
---|
| 1411 | |
---|
| 1412 | DO ji=1,kjpindex |
---|
| 1413 | DO jd = 1, nslm |
---|
| 1414 | IF(jd == 1) THEN ! the moisture at the 1st interface mc_layh(1) is at the surface, no interpolation |
---|
| 1415 | mc_layt(ji,jd) = mc_layh(ji,jd) |
---|
| 1416 | mcl_layt(ji,jd) = mcl_layh(ji,jd) |
---|
| 1417 | ELSEIF(jd == 2) THEN !! the mc_layt at the 2nd interface is interpolated using mc_layh(1) at surface and mc_layh(2) at the node |
---|
| 1418 | mc_layt(ji, jd) = mc_layh(ji,jd-1)*(znt(jd)-zlt(jd-1))/(znt(jd)-0.0) + & |
---|
| 1419 | mc_layh(ji, jd)*(zlt(jd-1)-0.0)/(znt(jd)-0.0) |
---|
| 1420 | mcl_layt(ji, jd) = mcl_layh(ji,jd-1)*(znt(jd)-zlt(jd-1))/(znt(jd)-0.0) + & |
---|
| 1421 | mcl_layh(ji, jd)*(zlt(jd-1)-0.0)/(znt(jd)-0.0) |
---|
| 1422 | ELSEIF(jd == nslm) THEN ! the mc_layt at the nslm interface is interpolated using mc_layh(nslm) and mc_layh(nslm-1) |
---|
| 1423 | mc_layt(ji, jd) = mc_layh(ji,jd-1)*(zlt(jd)-zlt(jd-1))/(zlt(jd)-znt(jd-1)) + & |
---|
| 1424 | mc_layh(ji,jd)*(zlt(jd-1)-znt(jd-1))/(zlt(jd)-znt(jd-1)) |
---|
| 1425 | mcl_layt(ji, jd) = mcl_layh(ji,jd-1)*(zlt(jd)-zlt(jd-1))/(zlt(jd)-znt(jd-1)) + & |
---|
| 1426 | mcl_layh(ji,jd)*(zlt(jd-1)-znt(jd-1))/(zlt(jd)-znt(jd-1)) |
---|
| 1427 | ELSE ! the mc_layt at the other interfaces are interpolated using mc_layh at adjacent nodes. |
---|
| 1428 | mc_layt(ji, jd) = mc_layh(ji, jd-1)*(1-dz5(jd-1)) + mc_layh(ji,jd)*dz5(jd-1) |
---|
| 1429 | mcl_layt(ji, jd) = mcl_layh(ji, jd-1)*(1-dz5(jd-1)) + mcl_layh(ji,jd)*dz5(jd-1) |
---|
| 1430 | ENDIF |
---|
| 1431 | |
---|
| 1432 | shum_ngrnd_perma(ji,jd) = shumdiag_perma(ji,jd) |
---|
| 1433 | tmc_layt(ji,jd) = tmc_layh(ji,jd) |
---|
| 1434 | ENDDO |
---|
| 1435 | |
---|
| 1436 | ! The deep layers in thermosoil where hydro is not existing |
---|
| 1437 | DO jd = nslm+1, ngrnd |
---|
| 1438 | shum_ngrnd_perma(ji,jd) = shumdiag_perma(ji,nslm) |
---|
| 1439 | mc_layt(ji,jd) = mc_layh(ji,nslm) |
---|
| 1440 | mcl_layt(ji,jd) = mcl_layh(ji,nslm) |
---|
| 1441 | tmc_layt(ji,jd) = tmc_layh(ji,nslm)/dlt(nslm) *dlt(jd) |
---|
| 1442 | ENDDO |
---|
| 1443 | ENDDO |
---|
| 1444 | |
---|
| 1445 | IF (printlev >= 4) WRITE(numout,*) 'thermosoil_humlev done' |
---|
| 1446 | |
---|
| 1447 | END SUBROUTINE thermosoil_humlev |
---|
| 1448 | |
---|
| 1449 | |
---|
| 1450 | !! ================================================================================================================================ |
---|
| 1451 | !! SUBROUTINE : thermosoil_energy_diag |
---|
| 1452 | !! |
---|
| 1453 | !>\BRIEF Calculate diagnostics |
---|
| 1454 | !! |
---|
| 1455 | !! DESCRIPTION : Calculate diagnostic variables coldcont_incr and coldcont_incr |
---|
| 1456 | !! |
---|
| 1457 | !! RECENT CHANGE(S) : None |
---|
| 1458 | !! |
---|
| 1459 | !! MAIN OUTPUT VARIABLE(S) : |
---|
| 1460 | !! |
---|
| 1461 | !! REFERENCE(S) : None |
---|
| 1462 | !! |
---|
| 1463 | !! FLOWCHART : None |
---|
| 1464 | !! \n |
---|
| 1465 | !_ ================================================================================================================================ |
---|
| 1466 | |
---|
| 1467 | SUBROUTINE thermosoil_energy_diag(kjpindex, temp_sol_new, soilcap) |
---|
| 1468 | |
---|
| 1469 | !! 0. Variables and parameter declaration |
---|
| 1470 | |
---|
| 1471 | !! 0.1 Input variables |
---|
| 1472 | |
---|
| 1473 | INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless) |
---|
| 1474 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! Surface temperature at the present time-step, Ts |
---|
| 1475 | !! @tex ($K$) @endtex |
---|
| 1476 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Apparent surface heat capacity |
---|
| 1477 | !! @tex ($J m^{-2} K^{-1}$) @endtex, |
---|
| 1478 | !! see eq. A29 of F. Hourdin\'s PhD thesis. |
---|
| 1479 | |
---|
| 1480 | !! 0.2 Output variables |
---|
| 1481 | |
---|
| 1482 | !! 0.3 Modified variables |
---|
| 1483 | |
---|
| 1484 | !! 0.4 Local variables |
---|
| 1485 | |
---|
| 1486 | INTEGER(i_std) :: ji, jg |
---|
| 1487 | !_ ================================================================================================================================ |
---|
| 1488 | |
---|
| 1489 | ! Sum up the energy content of all layers in the soil. |
---|
| 1490 | DO ji = 1, kjpindex |
---|
| 1491 | |
---|
| 1492 | IF (pcapa_en(ji,1) .LE. sn_capa) THEN |
---|
| 1493 | |
---|
| 1494 | ! Verify the energy conservation in the surface layer |
---|
| 1495 | coldcont_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji)) |
---|
| 1496 | surfheat_incr(ji) = zero |
---|
| 1497 | ELSE |
---|
| 1498 | |
---|
| 1499 | ! Verify the energy conservation in the surface layer |
---|
| 1500 | surfheat_incr(ji) = soilcap(ji) * (temp_sol_new(ji) - temp_sol_beg(ji)) |
---|
| 1501 | coldcont_incr(ji) = zero |
---|
| 1502 | ENDIF |
---|
| 1503 | ENDDO |
---|
| 1504 | |
---|
| 1505 | ! Save temp_sol_new to be used at next timestep |
---|
| 1506 | temp_sol_beg(:) = temp_sol_new(:) |
---|
| 1507 | |
---|
| 1508 | END SUBROUTINE thermosoil_energy_diag |
---|
| 1509 | |
---|
| 1510 | |
---|
| 1511 | |
---|
| 1512 | !! ================================================================================================================================ |
---|
| 1513 | !! SUBROUTINE : thermosoil_readjust |
---|
| 1514 | !! |
---|
| 1515 | !>\BRIEF |
---|
| 1516 | !! |
---|
| 1517 | !! DESCRIPTION : Energy conservation : Correction to make sure that the same latent heat is released and |
---|
| 1518 | !! consumed during freezing and thawing |
---|
| 1519 | !! |
---|
| 1520 | !! RECENT CHANGE(S) : None |
---|
| 1521 | !! |
---|
| 1522 | !! MAIN OUTPUT VARIABLE(S): ptn (soil temperature profile on the thermal axis), |
---|
| 1523 | !! |
---|
| 1524 | !! REFERENCE(S) : |
---|
| 1525 | !! |
---|
| 1526 | !! FLOWCHART : None |
---|
| 1527 | !! \n |
---|
| 1528 | !_ ================================================================================================================================ |
---|
| 1529 | |
---|
| 1530 | SUBROUTINE thermosoil_readjust(kjpindex, ptn) |
---|
| 1531 | |
---|
| 1532 | !! 0. Variables and parameter declaration |
---|
| 1533 | |
---|
| 1534 | !! 0.1 Input variables |
---|
| 1535 | INTEGER(i_std), INTENT(in) :: kjpindex |
---|
| 1536 | |
---|
| 1537 | !! 0.2 Modified variables |
---|
| 1538 | REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(inout) :: ptn |
---|
| 1539 | |
---|
| 1540 | !! 0.3 Local variables |
---|
| 1541 | INTEGER(i_std) :: ji, jg |
---|
| 1542 | INTEGER(i_std) :: lev3m !! Closest interface level to 3m |
---|
| 1543 | REAL(r_std) :: ptn_tmp |
---|
| 1544 | |
---|
| 1545 | ! The energy is spread over the layers down to approximatly 3m |
---|
| 1546 | ! Find the closest level to 3m. It can be below or above 3m. |
---|
| 1547 | lev3m=MINLOC(ABS(zlt(:)-3.0),dim=1) |
---|
| 1548 | IF (printlev >= 3) WRITE(numout,*) 'In thermosoil_adjust: lev3m=',lev3m, ' zlt(lev3m)=', zlt(lev3m) |
---|
| 1549 | |
---|
| 1550 | DO jg=1, ngrnd |
---|
| 1551 | DO ji=1, kjpindex |
---|
| 1552 | ! All soil latent energy is put into e_soil_lat(ji) |
---|
| 1553 | ! because the variable soil layers make it difficult to keep track of all |
---|
| 1554 | ! layers in this version |
---|
| 1555 | ! NOTE : pcapa has unit J/K/m3 and pcappa_supp has J/K |
---|
| 1556 | e_soil_lat(ji)=e_soil_lat(ji)+pcappa_supp(ji,jg)*(ptn(ji,jg)-ptn_beg(ji,jg)) |
---|
| 1557 | END DO |
---|
| 1558 | END DO |
---|
| 1559 | |
---|
| 1560 | DO ji=1, kjpindex |
---|
| 1561 | IF (e_soil_lat(ji).GT.min_sechiba.AND.MINVAL(ptn(ji,:)).GT.ZeroCelsius+fr_dT/2.) THEN |
---|
| 1562 | ! The soil is thawed: we spread the excess of energy over the uppermost lev3m levels |
---|
| 1563 | ! Here we increase the temperatures |
---|
| 1564 | DO jg=1, lev3m |
---|
| 1565 | ptn_tmp=ptn(ji,jg) |
---|
| 1566 | |
---|
| 1567 | ptn(ji,jg)=ptn(ji,jg)+MIN(e_soil_lat(ji)/pcapa(ji,jg)/zlt(lev3m), 0.5) |
---|
| 1568 | e_soil_lat(ji)=e_soil_lat(ji)-(ptn(ji,jg)-ptn_tmp)*pcapa(ji,jg)*dlt(jg) |
---|
| 1569 | ENDDO |
---|
| 1570 | ELSE IF (e_soil_lat(ji).LT.-min_sechiba.AND.MINVAL(ptn(ji,:)).GT.ZeroCelsius+fr_dT/2.) THEN |
---|
| 1571 | ! The soil is thawed |
---|
| 1572 | ! Here we decrease the temperatures |
---|
| 1573 | DO jg=1, lev3m |
---|
| 1574 | ptn_tmp=ptn(ji,jg) |
---|
| 1575 | ptn(ji,jg)=MAX(ZeroCelsius+fr_dT/2., ptn_tmp+e_soil_lat(ji)/pcapa(ji,jg)/zlt(lev3m)) |
---|
| 1576 | e_soil_lat(ji)=e_soil_lat(ji)+(ptn_tmp-ptn(ji,jg))*pcapa(ji,jg)*dlt(jg) |
---|
| 1577 | END DO |
---|
| 1578 | END IF |
---|
| 1579 | END DO |
---|
| 1580 | |
---|
| 1581 | END SUBROUTINE thermosoil_readjust |
---|
| 1582 | |
---|
| 1583 | !------------------------------------------------------------------- |
---|
| 1584 | |
---|
| 1585 | |
---|
| 1586 | |
---|
| 1587 | !! ================================================================================================================================ |
---|
| 1588 | !! SUBROUTINE : thermosoil_getdiff |
---|
| 1589 | !! |
---|
| 1590 | !>\BRIEF Computes soil and snow heat capacity and conductivity |
---|
| 1591 | !! |
---|
| 1592 | !! DESCRIPTION : Computation of the soil thermal properties; snow properties are also accounted for |
---|
| 1593 | !! |
---|
| 1594 | !! RECENT CHANGE(S) : None |
---|
| 1595 | !! |
---|
| 1596 | !! MAIN OUTPUT VARIABLE(S): |
---|
| 1597 | !! |
---|
| 1598 | !! REFERENCE(S) : |
---|
| 1599 | !! |
---|
| 1600 | !! FLOWCHART : None |
---|
| 1601 | !! \n |
---|
| 1602 | !_ ================================================================================================================================ |
---|
| 1603 | |
---|
| 1604 | SUBROUTINE thermosoil_getdiff( kjpindex, snow, ptn, mcs, njsc, snowrho, snowtemp, pb ) |
---|
| 1605 | |
---|
| 1606 | !! 0. Variables and parameter declaration |
---|
| 1607 | |
---|
| 1608 | !! 0.1 Input variables |
---|
| 1609 | INTEGER(i_std),INTENT(in) :: kjpindex |
---|
| 1610 | REAL(r_std),DIMENSION(kjpindex),INTENT (in) :: snow !! Snow mass |
---|
| 1611 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) |
---|
| 1612 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3) |
---|
| 1613 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density |
---|
| 1614 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K) |
---|
| 1615 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pressure (hPa) |
---|
| 1616 | REAL(r_std),DIMENSION(kjpindex,ngrnd),INTENT(in) :: ptn !! Soil temperature profile |
---|
| 1617 | |
---|
| 1618 | !! 0.3 Local variables |
---|
| 1619 | REAL :: xx !! Unfrozen fraction of the soil |
---|
| 1620 | REAL(r_std), DIMENSION(kjpindex) :: snow_h |
---|
| 1621 | REAL(r_std), DIMENSION(kjpindex,ngrnd) :: zx1, zx2 |
---|
| 1622 | INTEGER :: ji,jg |
---|
| 1623 | INTEGER :: jst |
---|
| 1624 | REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_tmp !! soil heat capacity (J/m3/K) |
---|
| 1625 | REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_spec !! SPECIFIC soil heat capacity (J/kg/K) |
---|
| 1626 | REAL(r_std) :: rho_tot !! Soil density (kg/m3) |
---|
| 1627 | |
---|
| 1628 | pcapa_tmp(:,:) = 0.0 |
---|
| 1629 | |
---|
| 1630 | !! Computes soil heat capacity and conductivity |
---|
| 1631 | DO ji = 1,kjpindex |
---|
| 1632 | |
---|
| 1633 | ! Since explicitsnow module is implemented set zx1=0 and zx2=1 |
---|
| 1634 | zx1(ji,:) = 0. |
---|
| 1635 | zx2(ji,:) = 1. |
---|
| 1636 | |
---|
| 1637 | DO jg = 1, ngrnd |
---|
| 1638 | jst = njsc(ji) |
---|
| 1639 | pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg) |
---|
| 1640 | ! |
---|
| 1641 | ! 2. Calculate volumetric heat capacity with allowance for permafrost |
---|
| 1642 | ! 2.1. soil heat capacity depending on temperature and humidity |
---|
| 1643 | ! For SP6MIP we also diagnose a specific heat capacity (pcapa_spec), |
---|
| 1644 | ! which requires to calculate the total density of the soil (rho_tot), always >> 0 |
---|
| 1645 | |
---|
| 1646 | IF (ptn(ji,jg) .LT. ZeroCelsius-fr_dT/2.) THEN |
---|
| 1647 | ! frozen soil |
---|
| 1648 | profil_froz(ji,jg) = 1. |
---|
| 1649 | pcappa_supp(ji,jg)= 0. |
---|
| 1650 | pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + so_capa_ice(ji) * tmc_layt(ji,jg) / mille / dlt(jg) |
---|
| 1651 | rho_tot = rho_soil * (1-mcs(ji)) + rho_ice * tmc_layt(ji,jg) / mille / dlt(jg) |
---|
| 1652 | pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot |
---|
| 1653 | |
---|
| 1654 | ELSEIF (ptn(ji,jg) .GT. ZeroCelsius+fr_dT/2.) THEN |
---|
| 1655 | ! unfrozen soil |
---|
| 1656 | pcapa(ji, jg) = pcapa_tmp(ji, jg) |
---|
| 1657 | profil_froz(ji,jg) = 0. |
---|
| 1658 | pcappa_supp(ji,jg)= 0. |
---|
| 1659 | rho_tot = rho_soil * (1-mcs(ji)) + rho_water * tmc_layt(ji,jg)/mille/dlt(jg) |
---|
| 1660 | pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot |
---|
| 1661 | ELSE |
---|
| 1662 | ! xx is the unfrozen fraction of soil water |
---|
| 1663 | xx = (ptn(ji,jg)-(ZeroCelsius-fr_dT/2.)) / fr_dT |
---|
| 1664 | profil_froz(ji,jg) = (1. - xx) |
---|
| 1665 | |
---|
| 1666 | IF (ok_freeze_thaw_latent_heat) THEN |
---|
| 1667 | pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + & |
---|
| 1668 | water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + & |
---|
| 1669 | so_capa_ice(ji) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) + & |
---|
| 1670 | shum_ngrnd_perma(ji,jg)*mcs(ji)*lhf*rho_water/fr_dT |
---|
| 1671 | ELSE |
---|
| 1672 | pcapa(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + & |
---|
| 1673 | water_capa * tmc_layt(ji,jg)/mille / dlt(jg) * xx + & |
---|
| 1674 | so_capa_ice(ji) * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) |
---|
| 1675 | ENDIF |
---|
| 1676 | |
---|
| 1677 | rho_tot = rho_soil* (1-mcs(ji)) + & |
---|
| 1678 | rho_water * tmc_layt(ji,jg)/mille / dlt(jg) * xx + & |
---|
| 1679 | rho_ice * tmc_layt(ji,jg) / mille/dlt(jg) * (1.-xx) |
---|
| 1680 | pcapa_spec(ji, jg) = pcapa(ji, jg) / rho_tot |
---|
| 1681 | |
---|
| 1682 | pcappa_supp(ji,jg)= shum_ngrnd_perma(ji,jg)*mcs(ji)*lhf*rho_water/fr_dT*zx2(ji,jg)*dlt(jg) |
---|
| 1683 | |
---|
| 1684 | ENDIF |
---|
| 1685 | |
---|
| 1686 | ! |
---|
| 1687 | ! 2.2. Take into account the snow and soil fractions in the layer |
---|
| 1688 | ! |
---|
| 1689 | pcapa(ji,jg) = zx1(ji,jg) * sn_capa + zx2(ji,jg) * pcapa(ji,jg) |
---|
| 1690 | |
---|
| 1691 | ! |
---|
| 1692 | ! 2.3. Calculate the heat capacity for energy conservation check |
---|
| 1693 | IF ( zx1(ji,jg).GT.0. ) THEN |
---|
| 1694 | pcapa_en(ji,jg) = sn_capa |
---|
| 1695 | ELSE |
---|
| 1696 | pcapa_en(ji,jg) = pcapa(ji,jg) |
---|
| 1697 | ENDIF |
---|
| 1698 | |
---|
| 1699 | END DO |
---|
| 1700 | ENDDO |
---|
| 1701 | |
---|
| 1702 | ! Output the specific heat capcaity for SP-MIP |
---|
| 1703 | CALL xios_orchidee_send_field("pcapa_spec",pcapa_spec) |
---|
| 1704 | |
---|
| 1705 | ! |
---|
| 1706 | ! 3. Calculate the heat conductivity with allowance for permafrost |
---|
| 1707 | ! |
---|
| 1708 | IF (ok_freeze_thaw_latent_heat) THEN |
---|
| 1709 | CALL thermosoil_cond (kjpindex, njsc, mcs, mc_layt, QZ, mcl_layt*(1-profil_froz), pkappa) |
---|
| 1710 | ELSE |
---|
| 1711 | CALL thermosoil_cond (kjpindex, njsc, mcs, mc_layt, QZ, mcl_layt, pkappa) |
---|
| 1712 | ENDIF |
---|
| 1713 | |
---|
| 1714 | !! Computes snow heat capacity and conductivity |
---|
| 1715 | DO ji = 1,kjpindex |
---|
| 1716 | pcapa_snow(ji,:) = snowrho(ji,:) * xci |
---|
| 1717 | pkappa_snow(ji,:) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2*snowrho(ji,:)*snowrho(ji,:)) + & |
---|
| 1718 | MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ & |
---|
| 1719 | ZSNOWTHRMCOND_CVAP)))*(XP00/(pb(ji)*100.))) |
---|
| 1720 | END DO |
---|
| 1721 | |
---|
| 1722 | END SUBROUTINE thermosoil_getdiff |
---|
| 1723 | |
---|
| 1724 | |
---|
| 1725 | !! ================================================================================================================================ |
---|
| 1726 | !! SUBROUTINE : thermosoil_getdiff_old_thermix_without_snow |
---|
| 1727 | !! |
---|
| 1728 | !>\BRIEF Computes soil and snow heat capacity and conductivity |
---|
| 1729 | !! |
---|
| 1730 | !! DESCRIPTION : Calculations of soil and snow thermal properties without effect of freezing. |
---|
| 1731 | !! |
---|
| 1732 | !! |
---|
| 1733 | !! RECENT CHANGE(S) : None |
---|
| 1734 | !! |
---|
| 1735 | !! MAIN OUTPUT VARIABLE(S): |
---|
| 1736 | !! |
---|
| 1737 | !! REFERENCE(S) : |
---|
| 1738 | !! |
---|
| 1739 | !! FLOWCHART : None |
---|
| 1740 | !! \n |
---|
| 1741 | !_ ================================================================================================================================ |
---|
| 1742 | |
---|
| 1743 | SUBROUTINE thermosoil_getdiff_old_thermix_without_snow( kjpindex, mcs, njsc, snowrho, snowtemp, pb ) |
---|
| 1744 | |
---|
| 1745 | !! 0. Variables and parameter declaration |
---|
| 1746 | |
---|
| 1747 | !! 0.1 Input variables |
---|
| 1748 | INTEGER(i_std), INTENT(in) :: kjpindex |
---|
| 1749 | INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) |
---|
| 1750 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: mcs !! Saturated moisture content (m3/m3) |
---|
| 1751 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density |
---|
| 1752 | REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature (K) |
---|
| 1753 | REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface pressure (hPa) |
---|
| 1754 | |
---|
| 1755 | |
---|
| 1756 | !! 0.1 Local variables |
---|
| 1757 | INTEGER(i_std) :: ji,jg, jst !! Index |
---|
| 1758 | REAL(r_std), DIMENSION (kjpindex,ngrnd) :: pcapa_tmp !! Soil heat capacity (J/m3/K) |
---|
| 1759 | |
---|
| 1760 | !! Computes soil heat capacity and conductivity |
---|
| 1761 | DO jg = 1,ngrnd |
---|
| 1762 | DO ji = 1,kjpindex |
---|
| 1763 | jst = njsc(ji) |
---|
| 1764 | pcapa_tmp(ji, jg) = so_capa_dry_ns(jst) * (1-mcs(ji)) + water_capa * tmc_layt(ji,jg)/mille/dlt(jg) |
---|
| 1765 | pcapa(ji,jg) = pcapa_tmp(ji, jg) |
---|
| 1766 | pcapa_en(ji,jg) = pcapa_tmp(ji, jg) |
---|
| 1767 | ENDDO |
---|
| 1768 | ENDDO |
---|
| 1769 | |
---|
| 1770 | CALL thermosoil_cond (kjpindex, njsc, mcs, mc_layt, QZ, mcl_layt, pkappa) |
---|
| 1771 | |
---|
| 1772 | IF (brk_flag == 1) THEN |
---|
| 1773 | ! Bedrock flag is activated |
---|
| 1774 | DO jg = ngrnd-1,ngrnd |
---|
| 1775 | DO ji = 1,kjpindex |
---|
| 1776 | pcapa(ji,jg) = brk_capa |
---|
| 1777 | pcapa_en(ji,jg) = brk_capa |
---|
| 1778 | pkappa(ji,jg) = brk_cond |
---|
| 1779 | ENDDO |
---|
| 1780 | ENDDO |
---|
| 1781 | ENDIF |
---|
| 1782 | |
---|
| 1783 | !! Computes snow heat capacity and conductivity |
---|
| 1784 | DO ji = 1,kjpindex |
---|
| 1785 | pcapa_snow(ji,:) = snowrho(ji,:) * xci |
---|
| 1786 | pkappa_snow(ji,:) = (ZSNOWTHRMCOND1 + ZSNOWTHRMCOND2*snowrho(ji,:)*snowrho(ji,:)) + & |
---|
| 1787 | MAX(0.0,(ZSNOWTHRMCOND_AVAP+(ZSNOWTHRMCOND_BVAP/(snowtemp(ji,:)+ & |
---|
| 1788 | ZSNOWTHRMCOND_CVAP)))*(XP00/(pb(ji)*100.))) |
---|
| 1789 | END DO |
---|
| 1790 | |
---|
| 1791 | END SUBROUTINE thermosoil_getdiff_old_thermix_without_snow |
---|
| 1792 | |
---|
| 1793 | |
---|
| 1794 | !! ================================================================================================================================ |
---|
| 1795 | !! SUBROUTINE : thermosoil_read_reftempfile |
---|
| 1796 | !! |
---|
| 1797 | !>\BRIEF |
---|
| 1798 | !! |
---|
| 1799 | !! DESCRIPTION : Read file with longterm soil temperature |
---|
| 1800 | !! |
---|
| 1801 | !! |
---|
| 1802 | !! RECENT CHANGE(S) : None |
---|
| 1803 | !! |
---|
| 1804 | !! MAIN OUTPUT VARIABLE(S): reftemp : Reference temerature |
---|
| 1805 | !! |
---|
| 1806 | !! REFERENCE(S) : |
---|
| 1807 | !! |
---|
| 1808 | !! FLOWCHART : None |
---|
| 1809 | !! \n |
---|
| 1810 | !_ ================================================================================================================================ |
---|
| 1811 | SUBROUTINE thermosoil_read_reftempfile(kjpindex,lalo,reftemp) |
---|
| 1812 | |
---|
| 1813 | USE interpweight |
---|
| 1814 | |
---|
| 1815 | IMPLICIT NONE |
---|
| 1816 | |
---|
| 1817 | !! 0. Variables and parameter declaration |
---|
| 1818 | |
---|
| 1819 | !! 0.1 Input variables |
---|
| 1820 | INTEGER(i_std), INTENT(in) :: kjpindex |
---|
| 1821 | REAL(r_std), DIMENSION(kjpindex,2), INTENT(in) :: lalo |
---|
| 1822 | |
---|
| 1823 | !! 0.2 Output variables |
---|
| 1824 | REAL(r_std), DIMENSION(kjpindex, ngrnd), INTENT(out) :: reftemp |
---|
| 1825 | |
---|
| 1826 | !! 0.3 Local variables |
---|
| 1827 | INTEGER(i_std) :: ib |
---|
| 1828 | CHARACTER(LEN=80) :: filename |
---|
| 1829 | REAL(r_std),DIMENSION(kjpindex) :: reftemp_file !! Horizontal temperature field interpolated from file [C] |
---|
| 1830 | INTEGER(i_std),DIMENSION(kjpindex,8) :: neighbours |
---|
| 1831 | REAL(r_std) :: vmin, vmax !! min/max values to use for the |
---|
| 1832 | !! renormalization |
---|
| 1833 | REAL(r_std), DIMENSION(kjpindex) :: areftemp !! Availability of data for the interpolation |
---|
| 1834 | CHARACTER(LEN=80) :: variablename !! Variable to interpolate |
---|
| 1835 | !! the file |
---|
| 1836 | CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file |
---|
| 1837 | REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable |
---|
| 1838 | !! (variabletypevals(1) = -un, not used) |
---|
| 1839 | CHARACTER(LEN=50) :: fractype !! method of calculation of fraction |
---|
| 1840 | !! 'XYKindTime': Input values are kinds |
---|
| 1841 | !! of something with a temporal |
---|
| 1842 | !! evolution on the dx*dy matrix' |
---|
| 1843 | LOGICAL :: nonegative !! whether negative values should be removed |
---|
| 1844 | CHARACTER(LEN=50) :: maskingtype !! Type of masking |
---|
| 1845 | !! 'nomask': no-mask is applied |
---|
| 1846 | !! 'mbelow': take values below maskvals(1) |
---|
| 1847 | !! 'mabove': take values above maskvals(1) |
---|
| 1848 | !! 'msumrange': take values within 2 ranges; |
---|
| 1849 | !! maskvals(2) <= SUM(vals(k)) <= maskvals(1) |
---|
| 1850 | !! maskvals(1) < SUM(vals(k)) <= maskvals(3) |
---|
| 1851 | !! (normalized by maskvals(3)) |
---|
| 1852 | !! 'var': mask values are taken from a |
---|
| 1853 | !! variable inside the file (>0) |
---|
| 1854 | REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to |
---|
| 1855 | !! `maskingtype') |
---|
| 1856 | CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask |
---|
| 1857 | REAL(r_std) :: reftemp_norefinf |
---|
| 1858 | REAL(r_std) :: reftemp_default !! Default value |
---|
| 1859 | |
---|
| 1860 | !Config Key = SOIL_REFTEMP_FILE |
---|
| 1861 | !Config Desc = File with climatological soil temperature |
---|
| 1862 | !Config If = READ_REFTEMP |
---|
| 1863 | !Config Def = reftemp.nc |
---|
| 1864 | !Config Help = |
---|
| 1865 | !Config Units = [FILE] |
---|
| 1866 | filename = 'reftemp.nc' |
---|
| 1867 | CALL getin_p('REFTEMP_FILE',filename) |
---|
| 1868 | |
---|
| 1869 | variablename = 'temperature' |
---|
| 1870 | |
---|
| 1871 | IF (printlev >= 1) WRITE(numout,*) "thermosoil_read_reftempfile: Read and interpolate file " & |
---|
| 1872 | // TRIM(filename) //" for variable " //TRIM(variablename) |
---|
| 1873 | |
---|
| 1874 | IF (xios_interpolation) THEN |
---|
| 1875 | |
---|
| 1876 | CALL xios_orchidee_recv_field('reftemp_interp',reftemp(:,1)) |
---|
| 1877 | |
---|
| 1878 | DO ib=1, kjpindex |
---|
| 1879 | reftemp(ib,:) = reftemp(ib,1) + ZeroCelsius |
---|
| 1880 | END DO |
---|
| 1881 | areftemp = 1.0 |
---|
| 1882 | ELSE |
---|
| 1883 | |
---|
| 1884 | |
---|
| 1885 | ! For this case there are not types/categories. We have 'only' a continuos field |
---|
| 1886 | ! Assigning values to vmin, vmax |
---|
| 1887 | |
---|
| 1888 | vmin = 0. |
---|
| 1889 | vmax = 9999. |
---|
| 1890 | |
---|
| 1891 | ! For this file we do not need neightbours! |
---|
| 1892 | neighbours = 0 |
---|
| 1893 | |
---|
| 1894 | !! Variables for interpweight |
---|
| 1895 | ! Type of calculation of cell fractions |
---|
| 1896 | fractype = 'default' |
---|
| 1897 | ! Name of the longitude and latitude in the input file |
---|
| 1898 | lonname = 'nav_lon' |
---|
| 1899 | latname = 'nav_lat' |
---|
| 1900 | ! Default value when no value is get from input file |
---|
| 1901 | reftemp_default = 1. |
---|
| 1902 | ! Reference value when no value is get from input file |
---|
| 1903 | reftemp_norefinf = 1. |
---|
| 1904 | ! Should negative values be set to zero from input file? |
---|
| 1905 | nonegative = .FALSE. |
---|
| 1906 | ! Type of mask to apply to the input data (see header for more details) |
---|
| 1907 | maskingtype = 'nomask' |
---|
| 1908 | ! Values to use for the masking (here not used) |
---|
| 1909 | maskvals = (/ undef_sechiba, undef_sechiba, undef_sechiba /) |
---|
| 1910 | ! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used) |
---|
| 1911 | namemaskvar = '' |
---|
| 1912 | |
---|
| 1913 | CALL interpweight_2Dcont(kjpindex, 0, 0, lalo, resolution, neighbours, & |
---|
| 1914 | contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, & |
---|
| 1915 | maskvals, namemaskvar, -1, fractype, reftemp_default, reftemp_norefinf, & |
---|
| 1916 | reftemp_file, areftemp) |
---|
| 1917 | IF (printlev >= 5) WRITE(numout,*)' thermosoil_read_reftempfile after interpweight_2Dcont' |
---|
| 1918 | |
---|
| 1919 | ! Copy reftemp_file temperature to all ground levels and transform into Kelvin |
---|
| 1920 | DO ib=1, kjpindex |
---|
| 1921 | reftemp(ib, :) = reftemp_file(ib)+ZeroCelsius |
---|
| 1922 | END DO |
---|
| 1923 | |
---|
| 1924 | END IF |
---|
| 1925 | |
---|
| 1926 | ! Write diagnostics |
---|
| 1927 | CALL xios_orchidee_send_field("areftemp",areftemp) |
---|
| 1928 | |
---|
| 1929 | END SUBROUTINE thermosoil_read_reftempfile |
---|
| 1930 | |
---|
| 1931 | END MODULE thermosoil |
---|