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