! ================================================================================================================================ ! MODULE : vertical_soil ! ! CONTACT : orchidee-help _at_ listes.ipsl.fr ! ! LICENCE : IPSL (2006) ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC ! !>\BRIEF Initialization of vertical discretization in the soil for hydrology and thermodynamics !! !!\n DESCRIPTION: Initialization of vertical discretization in the soil for hydrology and thermodynamics !! !! RECENT CHANGES: !! !! REFERENCE(S) : None !! !! SVN : !! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/vertical_soil.f90 $ !! $Date: 2018-10-03 10:22:28 +0200 (Wed, 03 Oct 2018) $ !! $Revision: 5454 $ !! \n !_ ================================================================================================================================ MODULE vertical_soil USE defprec USE ioipsl_para USE vertical_soil_var USE constantes_var IMPLICIT NONE PUBLIC :: vertical_soil_init CONTAINS !! ================================================================================================================================ !! SUBROUTINE : vertical_soil_init !! !>\BRIEF Initialization of vertical discretization in the soil for hydrology and thermodynamics !! !! DESCRIPTION : To define the total number of layers !! for hydrology model (nslm), !! and for thermodynamics model ngrnd. !! !! We have 2 type of levels in the hydrology and the thermodynamics : !! (1) the nodes where the scalar variables are computed. !! (2) the intermediate levels: these are the bounds of the layers and thus the !! contact point between layer i and i+1. !! !! Situation in hydrol.f90 !! The CWRR model only uses two variables to describe the vertical discretization : !! (1) znh : the depth of the nodes where soil moisture is computed. The nodes are not !! supposed to be in the middle of the layer. For the first and last layers the !! nodes are put onto the intermediate level. At 0 for the first layer and !! zmaxh for the last layer. !! (2) dnh : the distance between two nodes !! !! Situation in thermosoil.f90 : !! For the thermodynamics other variables are needed to describe the vertical structure : !! (1) znt(znh) : the depth of the nodes are the same as in hydrol.f90 except for !! the first and last levels. In older versions than revision XXXX of thermosoil it is !! computed with fz(i+1/2) !! (2) zlt (Zint) (zz_coef in thermosoil) : is the depth of the intermediate levels. !! In revision before XXXX of thermosoil.f90 it's computed by fz(i). !! Since revision XXXX, in the new formulation (vertical.f90) it is computed as !! zint(i) = (znh(i)+znh(i+1))/2 !! (3) dlt(dz_tmp) (dz2 in thermosoil): This is the thickness of the layer, i.e the !! distance between the top and bottom of the layer. Thus it is given by !! zint(i+1)-zint(i). It should not be confused with dnh in hydrology which is the !! distance between nodes. !! !! in standard hydrology model (de Rosnay Patricia, 2000; depth_topthickness = 9.77517107e-04, zmaxh=2.0, refinebottom = .FALSE.): !! !! Example for standard depth of the hydrology !! Number of layers in hydrol = 11 !! znh = depth of nodes !! [ 0.00000000e+00 1.95503421e-03 5.86510264e-03 1.36852395e-02 !! 2.93255132e-02 6.06060606e-02 1.23167155e-01 2.48289345e-01 !! 4.98533725e-01 9.99022483e-01 2.00000000e+00 ] !! dnh = internode distance !! [ 0. 0.00195503 0.00391007 0.00782014 0.01564027 0.03128055 !! 0.06256109 0.12512219 0.25024438 0.50048876 1.00097752 ] !! dlh = soil layer thickness !! [ 0.00097752 0.00293255 0.0058651 0.01173021 0.02346041 0.04692082 !! 0.09384164 0.18768328 0.37536657 0.75073314 0.50048876 ] !! hcum = depth of soil layer bottom !! [ 9.77517107e-04 3.91006843e-03 9.77517107e-03 2.15053764e-02 !! 4.49657869e-02 9.18866081e-02 1.85728250e-01 3.73411535e-01 !! 7.48778104e-01 1.49951124e+00 2.00000000e+00 ] !! !! !! In thermal model (an example of: depth_topthickness = 9.77517107e-04, zmaxt=10, depth_geom=10, refinebottom = .FALSE.): !! Number of layers in thermosoil = 19 !! The approximate maximal depth for thermosoil = 10.0078201332 !! The actual maximal depth for thermosoil = 10.0 !! znt= !! [ 4.88758554e-04 1.95503421e-03 5.86510264e-03 1.36852395e-02 !! 2.93255132e-02 6.06060606e-02 1.23167155e-01 2.48289345e-01 !! 4.98533725e-01 9.99022483e-01 1.74975562e+00 2.50048876e+00 !! 3.50146627e+00 4.50244379e+00 5.50342131e+00 6.50439882e+00 !! 7.50537634e+00 8.50635386e+00 9.50733137e+00 ] !! dlt= !! [ 9.77517107e-04 2.93255132e-03 5.86510264e-03 1.17302053e-02 !! 2.34604106e-02 4.69208211e-02 9.38416423e-02 1.87683285e-01 !! 3.75366569e-01 7.50733138e-01 5.00488758e-01 1.00097752e+00 !! 1.00097752e+00 1.00097752e+00 1.00097752e+00 1.00097752e+00 !! 1.00097752e+00 1.00097752e+00 9.93157383e-01 ] !! zlt= !! [ 9.77517107e-04 3.91006843e-03 9.77517107e-03 2.15053764e-02 !! 4.49657869e-02 9.18866081e-02 1.85728250e-01 3.73411535e-01 !! 7.48778104e-01 1.49951124e+00 2.00000000e+00 3.00097752e+00 !! 4.00195503e+00 5.00293255e+00 6.00391007e+00 7.00488758e+00 !! 8.00586510e+00 9.00684262e+00 1.00000000e+01 ] !! !! !! A simple figure below shows the discretization. ^ !! '------' means interface; 'X' means node; '...' means many layers; | means distance. !! v !! Keep in mind that the nodes are not necessarly in the middle of the layers. !! !! Hydrology Thermodynamics !! !! --------X------- znh(1) ^ ------------------- 1st ^ !! | | !! | X znt(1) | depth_topthickness, dlt(1) !! |dnh(2) | !! ----------------- | ------------------- 2nd zlt(1) v ^ !! | | !! X znh(2) v ^ X znt(2) |dlt(2) !! | | !! ----------------- |dnh(3) ------------------- 3rd zlt(2) v !! | !! X znh(3) v X znt(3) !! !! ----------------- ... ------------------- ... !! !! ... X ... !! !! --------X------- znh(nslm) ------------------- nslm zlt(nslm) !! !! X znt(nslm) !! !! ------------------- nslm+1 zlt(nslm+1) !! !! X znt(nslm+1) !! !! ------------------- ... !! !! X znt(ngrnd) !! !! ------------------- ngrnd zlt(ngrnd) !! !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCE(S) : !! !! FLOWCHART : !! \n !_ ================================================================================================================================ SUBROUTINE vertical_soil_init !! 0. Variable and parameter declaration !! 0.1 Local variables INTEGER(i_std), PARAMETER :: nblayermax=500 !! Preset 500 for max. number of soil layers. INTEGER(i_std) :: i, irev, iref, ntmp, ier REAL(r_std) :: zgeoend !! The node position where the geometrical increases of nodes stopped. (m) REAL(r_std) :: dgeoend !! The distance of the node at zgeoend and the node above. (m) REAL(r_std) :: dcst !! Related to refine at bottom. Work in progress... REAL(r_std) :: cstint !! Related to refine at bottom. Work in progress... REAL(r_std) :: hh !! Temporay variable, to calculate the layer thickness for temperature at zmaxh. REAL(r_std) :: ratio !! The ratio of geometric increas. INTEGER(i_std) :: nbrefine !! Related to refine at bottom. Work in progress... INTEGER(i_std) :: nbcst !! Related to refine at bottom. Work in progress... INTEGER(i_std) :: igeo !! Related to refine at bottom. Work in progress... REAL(r_std), DIMENSION(nblayermax+1) :: ztmp !! Depth at the node of each layer (m) REAL(r_std), DIMENSION(nblayermax) :: zint !! Depth at the interface of each layer (m) REAL(r_std), DIMENSION(nblayermax+1) :: dtmp !! Distance between the current node and the one above (m) REAL(r_std), DIMENSION(nblayermax+1) :: drefine !! Related to refine at bottom. Work in progress... INTEGER(i_std), DIMENSION(1) :: imin !! Related to refine at bottom. Work in progress... !! Variables controling the vertical discretiztaion REAL(r_std) :: depth_topthickness !! Thickness of the top Layer (m) REAL(r_std) :: depth_cstthickness !! Depth at which constant layer thickness start (m) REAL(r_std) :: depth_geom !! Depth at which we resume geometrical increases for temperature (m) REAL(r_std) :: ratio_geom_below !! Ratio of the geometrical series defining the thickness below DEPTH_GEOM. LOGICAL :: refinebottom !! Whether or not the hydrology layers will be refined towards the bottom. !! 1. Read parameters from run.def file !Config Key = DEPTH_MAX_T !Config Desc = Maximum depth of the soil thermodynamics !Config If = !Config Def = 90.0 !Config Help = Maximum depth of soil for temperature. !Config Units = m zmaxt = 90.0 CALL getin_p("DEPTH_MAX_T",zmaxt) !Config Key = DEPTH_MAX_H !Config Desc = Maximum depth of soil moisture !Config If = !Config Def = 2.0 !Config Help = Maximum depth of soil for soil moisture (CWRR). !Config Units = m zmaxh=2.0 CALL getin_p("DEPTH_MAX_H",zmaxh) ! Verification IF ( zmaxh > zmaxt) THEN CALL ipslerr_p(3,'vertical_soil_init',"ERROR : zmaxh needs to be smaller than zmaxt.", & "Correction : zmaxh set to zmaxt", " ") ENDIF !Config Key = DEPTH_TOPTHICK !Config Desc = Thickness of upper most Layer !Config If = !Config Def = 9.77517107e-04 !Config Help = Thickness of top hydrology layer for soil moisture (CWRR). !Config Units = m depth_topthickness = 9.77517107e-04 CALL getin_p("DEPTH_TOPTHICK",depth_topthickness) !Config Key = DEPTH_CSTTHICK !Config Desc = Depth at which constant layer thickness start !Config If = !Config Def = DEPTH_MAX_H !Config Help = Depth at which constant layer thickness start (smaller than zmaxh/2) !Config Units = m depth_cstthickness=zmaxh CALL getin_p("DEPTH_CSTTHICK",depth_cstthickness) IF ( (depth_cstthickness /= zmaxh) .AND. (depth_cstthickness > zmaxh/2.0) ) THEN CALL ipslerr_p(2,'vertical_soil_init',"ERROR : depth_cstthickness needs to be smaller than zmaxh/2.0", & "Correction : Constant thickness disabled", " ") depth_cstthickness = zmaxh ENDIF !Config Key = REFINEBOTTOM !Config Desc = Depth at which the hydrology layers will be refined towards the bottom. !Config If = !Config Def = .FALSE. !Config Help = Refinebottom is important when lower boundary conditions is different from a free drainage. !Config Units = - refinebottom = .FALSE. CALL getin_p("REFINEBOTTOM",refinebottom) !Config Key = DEPTH_GEOM !Config Desc = Depth at which we resume geometrical increases for temperature !Config If = !Config Def = DEPTH_MAX_H !Config Help = Depth at which the thickness increases again for temperature. ! This depth has to be deeper than the bottom of the hydrology. !Config Units = m depth_geom=zmaxh CALL getin_p("DEPTH_GEOM",depth_geom) IF ( depth_geom < zmaxh ) THEN CALL ipslerr_p(3,'vertical_soil_init',"ERROR : depth_geom needs to be larger than zmaxh", & "Correction : setting depth_geom to zmaxh", " ") ENDIF !Config Key = RATIO_GEOM_BELOW !Config Desc = Ratio of the geometrical series defining the thickness below DEPTH_GEOM !Config If = !Config Def = 2 !Config Help = Ratio of the geometrical series defining the thickness below DEPTH_GEOM. ! This parameter allows to cover the depth needed for temperature with fewer layers. !Config Units = - ratio_geom_below=2 CALL getin_p("RATIO_GEOM_BELOW",ratio_geom_below) ! ! Computing the layer depth for soil moisture. This defines the number of layers needed. ! ztmp(:) = 0.0 dtmp(:) = 0.0 DO i=1,nblayermax IF ( ztmp(i) < depth_cstthickness ) THEN ztmp(i+1) = depth_topthickness*2.0*((2**i)-1) dtmp(i+1) = ztmp(i+1)-ztmp(i) igeo=i+1 zgeoend=ztmp(i+1) dgeoend=dtmp(i+1) ELSE ztmp(i+1) = ztmp(i)+dtmp(i) dtmp(i+1) = dtmp(i) ENDIF ENDDO ! ! refine at bottom. Work in progress... ! nbrefine = 1 drefine(:) = 0.0 IF (refinebottom) THEN ! ! Compute parameters for the constant increment interval before refining, ! If needed ! cstint=zmaxh-(2.0*zgeoend) nbcst=MAX(INT(cstint/dgeoend), 0) IF ( nbcst > 0 ) THEN dcst=cstint/nbcst ELSE dcst=dgeoend ENDIF ! ! If we have to add constant increments ! IF ( nbcst > 0 ) THEN ! Add linear levels DO i=igeo,igeo+nbcst-1 ztmp(i+1) = ztmp(i)+dcst dtmp(i+1) = dcst ENDDO ! Refine the levels toward the bottom DO i=igeo+nbcst,2*igeo+nbcst-2 irev=(2*igeo+nbcst-1)-i ztmp(i+1) = ztmp(i)+dtmp(irev+1) dtmp(i+1) = ztmp(i+1)-ztmp(i) drefine(nbrefine)=dtmp(i+1) nbrefine=nbrefine+1 ENDDO ! Without constant increments ELSE imin=MINLOC(ABS(ztmp(:)-(zmaxh/2.0))) igeo=imin(1) ztmp(igeo)=zmaxh/2.0 dtmp(igeo) = ztmp(igeo)-ztmp(igeo-1) DO i=igeo,2*igeo irev=(2*igeo)-i ztmp(i+1) = ztmp(i)+dtmp(irev) dtmp(i+1) = ztmp(i-1)-ztmp(i) drefine(nbrefine)=dtmp(i+1) nbrefine=nbrefine+1 ENDDO ENDIF nbrefine = nbrefine-1 ENDIF ! Find the index (nslm) of the node closest to the zmaxh imin=MINLOC(ABS(ztmp(:)-zmaxh)) nslm=imin(1) ! ! ALLOCATE the arrays we need to keep ! ALLOCATE(znh(nslm), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable znh','','') ALLOCATE(dnh(nslm), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dnh','','') ALLOCATE(dlh(nslm), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dlh','','') ALLOCATE(zlh(nslm), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable zlh','','') ! ! Extract now the values we need to reach zmaxh ! znh(:)=ztmp(1:nslm) ! Force the last hydrological layer and node to be zmaxh znh(nslm)=zmaxh dnh(:)=dtmp(1:nslm) ! Recalculate the distance between the last 2 nodes dnh(nslm) = ztmp(nslm)-ztmp(nslm-1) ! Calculate the thickness of the layers DO i = 1, nslm-1 dlh(i) = (dnh(i) + dnh(i+1)) / 2.0 ENDDO dlh(nslm) = dnh(nslm)/2 ! ! Extension for the thermodynamics below zmaxh ! ntmp = nslm ztmp(:)=0.0 ! Exception for the first layer. It is at the top in the hydrology and in ! the middle for temperature. ztmp(1)=depth_topthickness/2 DO i=2,ntmp ztmp(i)=znh(i) ENDDO ! Exception for the last layer where again temperature needs to be in ! the middle of the layer. Also add a layer with the same thickness hh=dnh(ntmp)/2.0 ztmp(ntmp)=ztmp(ntmp)-hh/2.0 ztmp(ntmp+1)=ztmp(ntmp)+hh*1.5 ztmp(ntmp+2)=ztmp(ntmp+1)+hh*2.0 ntmp=ntmp+2 ! If we have created a refined region at the bottom of the soil moisture. We need to ! unwinde it for temperature below zmaxh. Only done for option refinebottom. IF ( nbrefine > 1 ) THEN DO i=ntmp,ntmp+nbrefine-1 iref=nbrefine-(i-ntmp) ztmp(i+1)=ztmp(i)+drefine(iref) ENDDO ntmp=ntmp+nbrefine ENDIF ! ! Resume the geometric increas of thickness DO i=ntmp,nblayermax IF ( ztmp(i) < depth_geom ) THEN ratio = 1.0 ELSE ratio = ratio_geom_below ENDIF ztmp(i+1)=ztmp(i)+ratio*(ztmp(i)-ztmp(i-1)) ENDDO ! Compute the depth of the lower interface of each layer. ! zint(1) = depth_topthickness DO i=2,nblayermax-1 zint(i) = (ztmp(i)+ztmp(i+1))/2.0 ENDDO zint(nslm-1) = (znh(nslm-1) + znh(nslm))/2.0 zint(nslm) = (znh(nslm)) zint(nblayermax) = ztmp(nblayermax)+(ztmp(nblayermax)-ztmp(nblayermax-1))/2.0 ! Determine the total number of layers for thermal (ngrnd). ! imin=MINLOC(ABS(zint(:)-zmaxt)) ngrnd=imin(1) ALLOCATE(znt(ngrnd), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable znt','','') ALLOCATE(dlt(ngrnd), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable dlt','','') ALLOCATE(zlt(ngrnd), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'vertical_soil_init','Problem in allocate of variable zlt','','') ! Assign values for znt, zlt, dlt ! znt(:)=ztmp(1:ngrnd) zlt(:) = zint(1:ngrnd) dlt(1) = zint(1) DO i=2,ngrnd dlt(i) = zint(i)-zint(i-1) ENDDO ! Depth of layers for the hydrology are the same as for thermosoil but only the upper nslm layers are used zlh(:) = zlt(1:nslm) ! Force the last thermal layer and node to be zmaxt zlt(ngrnd) = zmaxt dlt(ngrnd) = zmaxt-zint(ngrnd-1) IF (printlev >=1) THEN WRITE(numout,*) "=== Vertical discretization for hydrology ===================================" WRITE(numout,*) "Number of layers in hydrology: nslm=", nslm WRITE(numout,*) "Total depth in hydrology (DEPTH_MAX_H) (m): znh(nslm)=", znh(nslm) WRITE(numout,*) "Depth of nodes in hydrology (m): znh=", znh WRITE(numout,*) "Depth of lower layer-interface for hydrology (m): zlh=", zlh WRITE(numout,*) "Layer thickness for hydrology (m): dlh=",dlh WRITE(numout,*) "Internode distance for hydrology (m): dnh=",dnh WRITE(numout,*) "=== Vertical discretization for thermosoil ==================================" WRITE(numout,*) "Number of layers in thermosoil: ngrnd=", ngrnd WRITE(numout,*) "Total depth for thermosoil (DEPTH_MAX_T) (m): zlt=", zlt(ngrnd) WRITE(numout,*) "to be compared with calculated total depth for thermosoil (not used) zint=", zint(ngrnd) WRITE(numout,*) "Depth of the nodes in thermosoil (m): znt=",znt WRITE(numout,*) "Depth of lower layer-interface for thermosoil (m): zlt=", zlt(1:ngrnd) WRITE(numout,*) "Layer thickness for thermosoil (m): dlt=", dlt WRITE(numout,*) "=============================================================================" END IF END SUBROUTINE vertical_soil_init END MODULE vertical_soil