[7541] | 1 | ! ================================================================================================================================ |
---|
| 2 | ! MODULE : vertical_soil |
---|
| 3 | ! |
---|
| 4 | ! CONTACT : orchidee-help _at_ listes.ipsl.fr |
---|
| 5 | ! |
---|
| 6 | ! LICENCE : IPSL (2006) |
---|
| 7 | ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
---|
| 8 | ! |
---|
| 9 | !>\BRIEF Initialization of vertical discretization in the soil for hydrology and thermodynamics |
---|
| 10 | !! |
---|
| 11 | !!\n DESCRIPTION: Initialization of vertical discretization in the soil for hydrology and thermodynamics |
---|
| 12 | !! |
---|
| 13 | !! RECENT CHANGES: |
---|
| 14 | !! |
---|
| 15 | !! REFERENCE(S) : None |
---|
| 16 | !! |
---|
| 17 | !! SVN : |
---|
| 18 | !! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/vertical_soil.f90 $ |
---|
| 19 | !! $Date: 2018-10-03 10:22:28 +0200 (Wed, 03 Oct 2018) $ |
---|
| 20 | !! $Revision: 5454 $ |
---|
| 21 | !! \n |
---|
| 22 | !_ ================================================================================================================================ |
---|
| 23 | MODULE vertical_soil |
---|
| 24 | |
---|
| 25 | USE defprec |
---|
| 26 | USE ioipsl_para |
---|
| 27 | USE vertical_soil_var |
---|
| 28 | USE constantes_var |
---|
| 29 | |
---|
| 30 | IMPLICIT NONE |
---|
| 31 | |
---|
| 32 | |
---|
| 33 | PUBLIC :: vertical_soil_init |
---|
| 34 | |
---|
| 35 | CONTAINS |
---|
| 36 | !! ================================================================================================================================ |
---|
| 37 | !! SUBROUTINE : vertical_soil_init |
---|
| 38 | !! |
---|
| 39 | !>\BRIEF Initialization of vertical discretization in the soil for hydrology and thermodynamics |
---|
| 40 | !! |
---|
| 41 | !! DESCRIPTION : To define the total number of layers |
---|
| 42 | !! for hydrology model (nslm), |
---|
| 43 | !! and for thermodynamics model ngrnd. |
---|
| 44 | !! |
---|
| 45 | !! We have 2 type of levels in the hydrology and the thermodynamics : |
---|
| 46 | !! (1) the nodes where the scalar variables are computed. |
---|
| 47 | !! (2) the intermediate levels: these are the bounds of the layers and thus the |
---|
| 48 | !! contact point between layer i and i+1. |
---|
| 49 | !! |
---|
| 50 | !! Situation in hydrol.f90 |
---|
| 51 | !! The CWRR model only uses two variables to describe the vertical discretization : |
---|
| 52 | !! (1) znh : the depth of the nodes where soil moisture is computed. The nodes are not |
---|
| 53 | !! supposed to be in the middle of the layer. For the first and last layers the |
---|
| 54 | !! nodes are put onto the intermediate level. At 0 for the first layer and |
---|
| 55 | !! zmaxh for the last layer. |
---|
| 56 | !! (2) dnh : the distance between two nodes |
---|
| 57 | !! |
---|
| 58 | !! Situation in thermosoil.f90 : |
---|
| 59 | !! For the thermodynamics other variables are needed to describe the vertical structure : |
---|
| 60 | !! (1) znt(znh) : the depth of the nodes are the same as in hydrol.f90 except for |
---|
| 61 | !! the first and last levels. In older versions than revision XXXX of thermosoil it is |
---|
| 62 | !! computed with fz(i+1/2) |
---|
| 63 | !! (2) zlt (Zint) (zz_coef in thermosoil) : is the depth of the intermediate levels. |
---|
| 64 | !! In revision before XXXX of thermosoil.f90 it's computed by fz(i). |
---|
| 65 | !! Since revision XXXX, in the new formulation (vertical.f90) it is computed as |
---|
| 66 | !! zint(i) = (znh(i)+znh(i+1))/2 |
---|
| 67 | !! (3) dlt(dz_tmp) (dz2 in thermosoil): This is the thickness of the layer, i.e the |
---|
| 68 | !! distance between the top and bottom of the layer. Thus it is given by |
---|
| 69 | !! zint(i+1)-zint(i). It should not be confused with dnh in hydrology which is the |
---|
| 70 | !! distance between nodes. |
---|
| 71 | !! |
---|
| 72 | !! in standard hydrology model (de Rosnay Patricia, 2000; depth_topthickness = 9.77517107e-04, zmaxh=2.0, refinebottom = .FALSE.): |
---|
| 73 | !! |
---|
| 74 | !! Example for standard depth of the hydrology |
---|
| 75 | !! Number of layers in hydrol = 11 |
---|
| 76 | !! znh = depth of nodes |
---|
| 77 | !! [ 0.00000000e+00 1.95503421e-03 5.86510264e-03 1.36852395e-02 |
---|
| 78 | !! 2.93255132e-02 6.06060606e-02 1.23167155e-01 2.48289345e-01 |
---|
| 79 | !! 4.98533725e-01 9.99022483e-01 2.00000000e+00 ] |
---|
| 80 | !! dnh = internode distance |
---|
| 81 | !! [ 0. 0.00195503 0.00391007 0.00782014 0.01564027 0.03128055 |
---|
| 82 | !! 0.06256109 0.12512219 0.25024438 0.50048876 1.00097752 ] |
---|
| 83 | !! dlh = soil layer thickness |
---|
| 84 | !! [ 0.00097752 0.00293255 0.0058651 0.01173021 0.02346041 0.04692082 |
---|
| 85 | !! 0.09384164 0.18768328 0.37536657 0.75073314 0.50048876 ] |
---|
| 86 | !! hcum = depth of soil layer bottom |
---|
| 87 | !! [ 9.77517107e-04 3.91006843e-03 9.77517107e-03 2.15053764e-02 |
---|
| 88 | !! 4.49657869e-02 9.18866081e-02 1.85728250e-01 3.73411535e-01 |
---|
| 89 | !! 7.48778104e-01 1.49951124e+00 2.00000000e+00 ] |
---|
| 90 | !! |
---|
| 91 | !! |
---|
| 92 | !! In thermal model (an example of: depth_topthickness = 9.77517107e-04, zmaxt=10, depth_geom=10, refinebottom = .FALSE.): |
---|
| 93 | !! Number of layers in thermosoil = 19 |
---|
| 94 | !! The approximate maximal depth for thermosoil = 10.0078201332 |
---|
| 95 | !! The actual maximal depth for thermosoil = 10.0 |
---|
| 96 | !! znt= |
---|
| 97 | !! [ 4.88758554e-04 1.95503421e-03 5.86510264e-03 1.36852395e-02 |
---|
| 98 | !! 2.93255132e-02 6.06060606e-02 1.23167155e-01 2.48289345e-01 |
---|
| 99 | !! 4.98533725e-01 9.99022483e-01 1.74975562e+00 2.50048876e+00 |
---|
| 100 | !! 3.50146627e+00 4.50244379e+00 5.50342131e+00 6.50439882e+00 |
---|
| 101 | !! 7.50537634e+00 8.50635386e+00 9.50733137e+00 ] |
---|
| 102 | !! dlt= |
---|
| 103 | !! [ 9.77517107e-04 2.93255132e-03 5.86510264e-03 1.17302053e-02 |
---|
| 104 | !! 2.34604106e-02 4.69208211e-02 9.38416423e-02 1.87683285e-01 |
---|
| 105 | !! 3.75366569e-01 7.50733138e-01 5.00488758e-01 1.00097752e+00 |
---|
| 106 | !! 1.00097752e+00 1.00097752e+00 1.00097752e+00 1.00097752e+00 |
---|
| 107 | !! 1.00097752e+00 1.00097752e+00 9.93157383e-01 ] |
---|
| 108 | !! zlt= |
---|
| 109 | !! [ 9.77517107e-04 3.91006843e-03 9.77517107e-03 2.15053764e-02 |
---|
| 110 | !! 4.49657869e-02 9.18866081e-02 1.85728250e-01 3.73411535e-01 |
---|
| 111 | !! 7.48778104e-01 1.49951124e+00 2.00000000e+00 3.00097752e+00 |
---|
| 112 | !! 4.00195503e+00 5.00293255e+00 6.00391007e+00 7.00488758e+00 |
---|
| 113 | !! 8.00586510e+00 9.00684262e+00 1.00000000e+01 ] |
---|
| 114 | !! |
---|
| 115 | !! |
---|
| 116 | !! A simple figure below shows the discretization. ^ |
---|
| 117 | !! '------' means interface; 'X' means node; '...' means many layers; | means distance. |
---|
| 118 | !! v |
---|
| 119 | !! Keep in mind that the nodes are not necessarly in the middle of the layers. |
---|
| 120 | !! |
---|
| 121 | !! Hydrology Thermodynamics |
---|
| 122 | !! |
---|
| 123 | !! --------X------- znh(1) ^ ------------------- 1st ^ |
---|
| 124 | !! | | |
---|
| 125 | !! | X znt(1) | depth_topthickness, dlt(1) |
---|
| 126 | !! |dnh(2) | |
---|
| 127 | !! ----------------- | ------------------- 2nd zlt(1) v ^ |
---|
| 128 | !! | | |
---|
| 129 | !! X znh(2) v ^ X znt(2) |dlt(2) |
---|
| 130 | !! | | |
---|
| 131 | !! ----------------- |dnh(3) ------------------- 3rd zlt(2) v |
---|
| 132 | !! | |
---|
| 133 | !! X znh(3) v X znt(3) |
---|
| 134 | !! |
---|
| 135 | !! ----------------- ... ------------------- ... |
---|
| 136 | !! |
---|
| 137 | !! ... X ... |
---|
| 138 | !! |
---|
| 139 | !! --------X------- znh(nslm) ------------------- nslm zlt(nslm) |
---|
| 140 | !! |
---|
| 141 | !! X znt(nslm) |
---|
| 142 | !! |
---|
| 143 | !! ------------------- nslm+1 zlt(nslm+1) |
---|
| 144 | !! |
---|
| 145 | !! X znt(nslm+1) |
---|
| 146 | !! |
---|
| 147 | !! ------------------- ... |
---|
| 148 | !! |
---|
| 149 | !! X znt(ngrnd) |
---|
| 150 | !! |
---|
| 151 | !! ------------------- ngrnd zlt(ngrnd) |
---|
| 152 | !! |
---|
| 153 | !! |
---|
| 154 | !! RECENT CHANGE(S): None |
---|
| 155 | !! |
---|
| 156 | !! MAIN OUTPUT VARIABLE(S): |
---|
| 157 | !! |
---|
| 158 | !! REFERENCE(S) : |
---|
| 159 | !! |
---|
| 160 | !! FLOWCHART : |
---|
| 161 | !! \n |
---|
| 162 | !_ ================================================================================================================================ |
---|
| 163 | |
---|
| 164 | SUBROUTINE vertical_soil_init |
---|
| 165 | |
---|
| 166 | !! 0. Variable and parameter declaration |
---|
| 167 | !! 0.1 Local variables |
---|
| 168 | INTEGER(i_std), PARAMETER :: nblayermax=500 !! Preset 500 for max. number of soil layers. |
---|
| 169 | INTEGER(i_std) :: i, irev, iref, ntmp, ier |
---|
| 170 | REAL(r_std) :: zgeoend !! The node position where the geometrical increases of nodes stopped. (m) |
---|
| 171 | REAL(r_std) :: dgeoend !! The distance of the node at zgeoend and the node above. (m) |
---|
| 172 | REAL(r_std) :: dcst !! Related to refine at bottom. Work in progress... |
---|
| 173 | REAL(r_std) :: cstint !! Related to refine at bottom. Work in progress... |
---|
| 174 | REAL(r_std) :: hh !! Temporay variable, to calculate the layer thickness for temperature at zmaxh. |
---|
| 175 | REAL(r_std) :: ratio !! The ratio of geometric increas. |
---|
| 176 | INTEGER(i_std) :: nbrefine !! Related to refine at bottom. Work in progress... |
---|
| 177 | INTEGER(i_std) :: nbcst !! Related to refine at bottom. Work in progress... |
---|
| 178 | INTEGER(i_std) :: igeo !! Related to refine at bottom. Work in progress... |
---|
| 179 | REAL(r_std), DIMENSION(nblayermax+1) :: ztmp !! Depth at the node of each layer (m) |
---|
| 180 | REAL(r_std), DIMENSION(nblayermax) :: zint !! Depth at the interface of each layer (m) |
---|
| 181 | REAL(r_std), DIMENSION(nblayermax+1) :: dtmp !! Distance between the current node and the one above (m) |
---|
| 182 | REAL(r_std), DIMENSION(nblayermax+1) :: drefine !! Related to refine at bottom. Work in progress... |
---|
| 183 | INTEGER(i_std), DIMENSION(1) :: imin !! Related to refine at bottom. Work in progress... |
---|
| 184 | |
---|
| 185 | !! Variables controling the vertical discretiztaion |
---|
| 186 | REAL(r_std) :: depth_topthickness !! Thickness of the top Layer (m) |
---|
| 187 | REAL(r_std) :: depth_cstthickness !! Depth at which constant layer thickness start (m) |
---|
| 188 | REAL(r_std) :: depth_geom !! Depth at which we resume geometrical increases for temperature (m) |
---|
| 189 | REAL(r_std) :: ratio_geom_below !! Ratio of the geometrical series defining the thickness below DEPTH_GEOM. |
---|
| 190 | LOGICAL :: refinebottom !! Whether or not the hydrology layers will be refined towards the bottom. |
---|
| 191 | |
---|
| 192 | |
---|
| 193 | !! 1. Read parameters from run.def file |
---|
| 194 | |
---|
| 195 | !Config Key = DEPTH_MAX_T |
---|
| 196 | !Config Desc = Maximum depth of the soil thermodynamics |
---|
| 197 | !Config If = |
---|
| 198 | !Config Def = 90.0 |
---|
| 199 | !Config Help = Maximum depth of soil for temperature. |
---|
| 200 | !Config Units = m |
---|
| 201 | zmaxt = 90.0 |
---|
| 202 | CALL getin_p("DEPTH_MAX_T",zmaxt) |
---|
| 203 | |
---|
| 204 | !Config Key = DEPTH_MAX_H |
---|
| 205 | !Config Desc = Maximum depth of soil moisture |
---|
| 206 | !Config If = |
---|
| 207 | !Config Def = 2.0 |
---|
| 208 | !Config Help = Maximum depth of soil for soil moisture (CWRR). |
---|
| 209 | !Config Units = m |
---|
| 210 | zmaxh=2.0 |
---|
| 211 | CALL getin_p("DEPTH_MAX_H",zmaxh) |
---|
| 212 | |
---|
| 213 | ! Verification |
---|
| 214 | IF ( zmaxh > zmaxt) THEN |
---|
| 215 | CALL ipslerr_p(3,'vertical_soil_init',"ERROR : zmaxh needs to be smaller than zmaxt.", & |
---|
| 216 | "Correction : zmaxh set to zmaxt", " ") |
---|
| 217 | ENDIF |
---|
| 218 | |
---|
| 219 | !Config Key = DEPTH_TOPTHICK |
---|
| 220 | !Config Desc = Thickness of upper most Layer |
---|
| 221 | !Config If = |
---|
| 222 | !Config Def = 9.77517107e-04 |
---|
| 223 | !Config Help = Thickness of top hydrology layer for soil moisture (CWRR). |
---|
| 224 | !Config Units = m |
---|
| 225 | depth_topthickness = 9.77517107e-04 |
---|
| 226 | CALL getin_p("DEPTH_TOPTHICK",depth_topthickness) |
---|
| 227 | |
---|
| 228 | !Config Key = DEPTH_CSTTHICK |
---|
| 229 | !Config Desc = Depth at which constant layer thickness start |
---|
| 230 | !Config If = |
---|
| 231 | !Config Def = DEPTH_MAX_H |
---|
| 232 | !Config Help = Depth at which constant layer thickness start (smaller than zmaxh/2) |
---|
| 233 | !Config Units = m |
---|
| 234 | depth_cstthickness=zmaxh |
---|
| 235 | CALL getin_p("DEPTH_CSTTHICK",depth_cstthickness) |
---|
| 236 | |
---|
| 237 | IF ( (depth_cstthickness /= zmaxh) .AND. (depth_cstthickness > zmaxh/2.0) ) THEN |
---|
| 238 | CALL ipslerr_p(2,'vertical_soil_init',"ERROR : depth_cstthickness needs to be smaller than zmaxh/2.0", & |
---|
| 239 | "Correction : Constant thickness disabled", " ") |
---|
| 240 | depth_cstthickness = zmaxh |
---|
| 241 | ENDIF |
---|
| 242 | |
---|
| 243 | !Config Key = REFINEBOTTOM |
---|
| 244 | !Config Desc = Depth at which the hydrology layers will be refined towards the bottom. |
---|
| 245 | !Config If = |
---|
| 246 | !Config Def = .FALSE. |
---|
| 247 | !Config Help = Refinebottom is important when lower boundary conditions is different from a free drainage. |
---|
| 248 | !Config Units = - |
---|
| 249 | refinebottom = .FALSE. |
---|
| 250 | CALL getin_p("REFINEBOTTOM",refinebottom) |
---|
| 251 | |
---|
| 252 | !Config Key = DEPTH_GEOM |
---|
| 253 | !Config Desc = Depth at which we resume geometrical increases for temperature |
---|
| 254 | !Config If = |
---|
| 255 | !Config Def = DEPTH_MAX_H |
---|
| 256 | !Config Help = Depth at which the thickness increases again for temperature. |
---|
| 257 | ! This depth has to be deeper than the bottom of the hydrology. |
---|
| 258 | !Config Units = m |
---|
| 259 | depth_geom=zmaxh |
---|
| 260 | CALL getin_p("DEPTH_GEOM",depth_geom) |
---|
| 261 | |
---|
| 262 | IF ( depth_geom < zmaxh ) THEN |
---|
| 263 | CALL ipslerr_p(3,'vertical_soil_init',"ERROR : depth_geom needs to be larger than zmaxh", & |
---|
| 264 | "Correction : setting depth_geom to zmaxh", " ") |
---|
| 265 | ENDIF |
---|
| 266 | |
---|
| 267 | !Config Key = RATIO_GEOM_BELOW |
---|
| 268 | !Config Desc = Ratio of the geometrical series defining the thickness below DEPTH_GEOM |
---|
| 269 | !Config If = |
---|
| 270 | !Config Def = 2 |
---|
| 271 | !Config Help = Ratio of the geometrical series defining the thickness below DEPTH_GEOM. |
---|
| 272 | ! This parameter allows to cover the depth needed for temperature with fewer layers. |
---|
| 273 | !Config Units = - |
---|
| 274 | ratio_geom_below=2 |
---|
| 275 | CALL getin_p("RATIO_GEOM_BELOW",ratio_geom_below) |
---|
| 276 | |
---|
| 277 | ! |
---|
| 278 | ! Computing the layer depth for soil moisture. This defines the number of layers needed. |
---|
| 279 | ! |
---|
| 280 | ztmp(:) = 0.0 |
---|
| 281 | dtmp(:) = 0.0 |
---|
| 282 | DO i=1,nblayermax |
---|
| 283 | IF ( ztmp(i) < depth_cstthickness ) THEN |
---|
| 284 | ztmp(i+1) = depth_topthickness*2.0*((2**i)-1) |
---|
| 285 | dtmp(i+1) = ztmp(i+1)-ztmp(i) |
---|
| 286 | igeo=i+1 |
---|
| 287 | zgeoend=ztmp(i+1) |
---|
| 288 | dgeoend=dtmp(i+1) |
---|
| 289 | ELSE |
---|
| 290 | ztmp(i+1) = ztmp(i)+dtmp(i) |
---|
| 291 | dtmp(i+1) = dtmp(i) |
---|
| 292 | ENDIF |
---|
| 293 | ENDDO |
---|
| 294 | ! |
---|
| 295 | ! refine at bottom. Work in progress... |
---|
| 296 | ! |
---|
| 297 | nbrefine = 1 |
---|
| 298 | drefine(:) = 0.0 |
---|
| 299 | IF (refinebottom) THEN |
---|
| 300 | ! |
---|
| 301 | ! Compute parameters for the constant increment interval before refining, |
---|
| 302 | ! If needed ! |
---|
| 303 | cstint=zmaxh-(2.0*zgeoend) |
---|
| 304 | nbcst=MAX(INT(cstint/dgeoend), 0) |
---|
| 305 | IF ( nbcst > 0 ) THEN |
---|
| 306 | dcst=cstint/nbcst |
---|
| 307 | ELSE |
---|
| 308 | dcst=dgeoend |
---|
| 309 | ENDIF |
---|
| 310 | ! |
---|
| 311 | ! If we have to add constant increments |
---|
| 312 | ! |
---|
| 313 | IF ( nbcst > 0 ) THEN |
---|
| 314 | ! Add linear levels |
---|
| 315 | DO i=igeo,igeo+nbcst-1 |
---|
| 316 | ztmp(i+1) = ztmp(i)+dcst |
---|
| 317 | dtmp(i+1) = dcst |
---|
| 318 | ENDDO |
---|
| 319 | ! Refine the levels toward the bottom |
---|
| 320 | DO i=igeo+nbcst,2*igeo+nbcst-2 |
---|
| 321 | irev=(2*igeo+nbcst-1)-i |
---|
| 322 | ztmp(i+1) = ztmp(i)+dtmp(irev+1) |
---|
| 323 | dtmp(i+1) = ztmp(i+1)-ztmp(i) |
---|
| 324 | drefine(nbrefine)=dtmp(i+1) |
---|
| 325 | nbrefine=nbrefine+1 |
---|
| 326 | ENDDO |
---|
| 327 | ! Without constant increments |
---|
| 328 | ELSE |
---|
| 329 | imin=MINLOC(ABS(ztmp(:)-(zmaxh/2.0))) |
---|
| 330 | igeo=imin(1) |
---|
| 331 | ztmp(igeo)=zmaxh/2.0 |
---|
| 332 | dtmp(igeo) = ztmp(igeo)-ztmp(igeo-1) |
---|
| 333 | DO i=igeo,2*igeo |
---|
| 334 | irev=(2*igeo)-i |
---|
| 335 | ztmp(i+1) = ztmp(i)+dtmp(irev) |
---|
| 336 | dtmp(i+1) = ztmp(i-1)-ztmp(i) |
---|
| 337 | drefine(nbrefine)=dtmp(i+1) |
---|
| 338 | nbrefine=nbrefine+1 |
---|
| 339 | ENDDO |
---|
| 340 | ENDIF |
---|
| 341 | nbrefine = nbrefine-1 |
---|
| 342 | ENDIF |
---|
| 343 | |
---|
| 344 | |
---|
| 345 | |
---|
| 346 | |
---|
| 347 | ! Find the index (nslm) of the node closest to the zmaxh |
---|
| 348 | imin=MINLOC(ABS(ztmp(:)-zmaxh)) |
---|
| 349 | nslm=imin(1) |
---|
| 350 | ! |
---|
| 351 | ! ALLOCATE the arrays we need to keep |
---|
| 352 | ! |
---|
| 353 | ALLOCATE(znh(nslm), stat=ier) |
---|
| 354 | IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable znh','','') |
---|
| 355 | ALLOCATE(dnh(nslm), stat=ier) |
---|
| 356 | IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dnh','','') |
---|
| 357 | ALLOCATE(dlh(nslm), stat=ier) |
---|
| 358 | IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dlh','','') |
---|
| 359 | ALLOCATE(zlh(nslm), stat=ier) |
---|
| 360 | IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable zlh','','') |
---|
| 361 | ! |
---|
| 362 | ! Extract now the values we need to reach zmaxh |
---|
| 363 | ! |
---|
| 364 | znh(:)=ztmp(1:nslm) |
---|
| 365 | |
---|
| 366 | ! Force the last hydrological layer and node to be zmaxh |
---|
| 367 | znh(nslm)=zmaxh |
---|
| 368 | dnh(:)=dtmp(1:nslm) |
---|
| 369 | ! Recalculate the distance between the last 2 nodes |
---|
| 370 | dnh(nslm) = ztmp(nslm)-ztmp(nslm-1) |
---|
| 371 | |
---|
| 372 | ! Calculate the thickness of the layers |
---|
| 373 | DO i = 1, nslm-1 |
---|
| 374 | dlh(i) = (dnh(i) + dnh(i+1)) / 2.0 |
---|
| 375 | ENDDO |
---|
| 376 | dlh(nslm) = dnh(nslm)/2 |
---|
| 377 | |
---|
| 378 | ! |
---|
| 379 | ! Extension for the thermodynamics below zmaxh |
---|
| 380 | ! |
---|
| 381 | ntmp = nslm |
---|
| 382 | ztmp(:)=0.0 |
---|
| 383 | ! Exception for the first layer. It is at the top in the hydrology and in |
---|
| 384 | ! the middle for temperature. |
---|
| 385 | ztmp(1)=depth_topthickness/2 |
---|
| 386 | DO i=2,ntmp |
---|
| 387 | ztmp(i)=znh(i) |
---|
| 388 | ENDDO |
---|
| 389 | |
---|
| 390 | ! Exception for the last layer where again temperature needs to be in |
---|
| 391 | ! the middle of the layer. Also add a layer with the same thickness |
---|
| 392 | hh=dnh(ntmp)/2.0 |
---|
| 393 | ztmp(ntmp)=ztmp(ntmp)-hh/2.0 |
---|
| 394 | ztmp(ntmp+1)=ztmp(ntmp)+hh*1.5 |
---|
| 395 | ztmp(ntmp+2)=ztmp(ntmp+1)+hh*2.0 |
---|
| 396 | ntmp=ntmp+2 |
---|
| 397 | |
---|
| 398 | ! If we have created a refined region at the bottom of the soil moisture. We need to |
---|
| 399 | ! unwinde it for temperature below zmaxh. Only done for option refinebottom. |
---|
| 400 | IF ( nbrefine > 1 ) THEN |
---|
| 401 | DO i=ntmp,ntmp+nbrefine-1 |
---|
| 402 | iref=nbrefine-(i-ntmp) |
---|
| 403 | ztmp(i+1)=ztmp(i)+drefine(iref) |
---|
| 404 | ENDDO |
---|
| 405 | ntmp=ntmp+nbrefine |
---|
| 406 | ENDIF |
---|
| 407 | ! |
---|
| 408 | ! Resume the geometric increas of thickness |
---|
| 409 | DO i=ntmp,nblayermax |
---|
| 410 | IF ( ztmp(i) < depth_geom ) THEN |
---|
| 411 | ratio = 1.0 |
---|
| 412 | ELSE |
---|
| 413 | ratio = ratio_geom_below |
---|
| 414 | ENDIF |
---|
| 415 | ztmp(i+1)=ztmp(i)+ratio*(ztmp(i)-ztmp(i-1)) |
---|
| 416 | ENDDO |
---|
| 417 | |
---|
| 418 | ! Compute the depth of the lower interface of each layer. |
---|
| 419 | ! |
---|
| 420 | zint(1) = depth_topthickness |
---|
| 421 | DO i=2,nblayermax-1 |
---|
| 422 | zint(i) = (ztmp(i)+ztmp(i+1))/2.0 |
---|
| 423 | ENDDO |
---|
| 424 | zint(nslm-1) = (znh(nslm-1) + znh(nslm))/2.0 |
---|
| 425 | zint(nslm) = (znh(nslm)) |
---|
| 426 | zint(nblayermax) = ztmp(nblayermax)+(ztmp(nblayermax)-ztmp(nblayermax-1))/2.0 |
---|
| 427 | |
---|
| 428 | ! Determine the total number of layers for thermal (ngrnd). |
---|
| 429 | ! |
---|
| 430 | imin=MINLOC(ABS(zint(:)-zmaxt)) |
---|
| 431 | ngrnd=imin(1) |
---|
| 432 | ALLOCATE(znt(ngrnd), stat=ier) |
---|
| 433 | IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable znt','','') |
---|
| 434 | ALLOCATE(dlt(ngrnd), stat=ier) |
---|
| 435 | IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dlt','','') |
---|
| 436 | ALLOCATE(zlt(ngrnd), stat=ier) |
---|
| 437 | IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable zlt','','') |
---|
| 438 | |
---|
| 439 | ! Assign values for znt, zlt, dlt |
---|
| 440 | ! |
---|
| 441 | znt(:)=ztmp(1:ngrnd) |
---|
| 442 | zlt(:) = zint(1:ngrnd) |
---|
| 443 | dlt(1) = zint(1) |
---|
| 444 | DO i=2,ngrnd |
---|
| 445 | dlt(i) = zint(i)-zint(i-1) |
---|
| 446 | ENDDO |
---|
| 447 | ! Depth of layers for the hydrology are the same as for thermosoil but only the upper nslm layers are used |
---|
| 448 | zlh(:) = zlt(1:nslm) |
---|
| 449 | |
---|
| 450 | ! Force the last thermal layer and node to be zmaxt |
---|
| 451 | zlt(ngrnd) = zmaxt |
---|
| 452 | dlt(ngrnd) = zmaxt-zint(ngrnd-1) |
---|
| 453 | |
---|
| 454 | |
---|
| 455 | |
---|
| 456 | |
---|
| 457 | IF (printlev >=1) THEN |
---|
| 458 | WRITE(numout,*) "=== Vertical discretization for hydrology ===================================" |
---|
| 459 | WRITE(numout,*) "Number of layers in hydrology: nslm=", nslm |
---|
| 460 | WRITE(numout,*) "Total depth in hydrology (DEPTH_MAX_H) (m): znh(nslm)=", znh(nslm) |
---|
| 461 | WRITE(numout,*) "Depth of nodes in hydrology (m): znh=", znh |
---|
| 462 | WRITE(numout,*) "Depth of lower layer-interface for hydrology (m): zlh=", zlh |
---|
| 463 | WRITE(numout,*) "Layer thickness for hydrology (m): dlh=",dlh |
---|
| 464 | WRITE(numout,*) "Internode distance for hydrology (m): dnh=",dnh |
---|
| 465 | WRITE(numout,*) "=== Vertical discretization for thermosoil ==================================" |
---|
| 466 | WRITE(numout,*) "Number of layers in thermosoil: ngrnd=", ngrnd |
---|
| 467 | WRITE(numout,*) "Total depth for thermosoil (DEPTH_MAX_T) (m): zlt=", zlt(ngrnd) |
---|
| 468 | WRITE(numout,*) "to be compared with calculated total depth for thermosoil (not used) zint=", zint(ngrnd) |
---|
| 469 | WRITE(numout,*) "Depth of the nodes in thermosoil (m): znt=",znt |
---|
| 470 | WRITE(numout,*) "Depth of lower layer-interface for thermosoil (m): zlt=", zlt(1:ngrnd) |
---|
| 471 | WRITE(numout,*) "Layer thickness for thermosoil (m): dlt=", dlt |
---|
| 472 | WRITE(numout,*) "=============================================================================" |
---|
| 473 | END IF |
---|
| 474 | END SUBROUTINE vertical_soil_init |
---|
| 475 | |
---|
| 476 | END MODULE vertical_soil |
---|