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