Ignore:
Timestamp:
Feb 12, 2026, 9:09:12 AM (2 weeks ago)
Author:
jbclement
Message:

PEM:
Major refactor following the previous ones (r3989 and r3991) completing the large structural reorganization and cleanup of the PEM codebase. This revision introduces newly designed modules, standardizes interfaces with explicit ini/end APIs and adds native NetCDF I/O together with explicit PCM/PEM adapters. In detail:

  • Some PEM models were corrected or improved:
    • Frost/perennial ice semantics are clarified via renaming;
    • Soil temperature remapping clarified, notably by removing the rescaling of temperature deviation;
    • Geothermal flux for the PCM is computed based on the PEM state;
  • New explicit PEM/PCM adapters ("set_*"/"build4PCM_*") to decouple PEM internal representation from PCM file layouts and reconstruct consistent fields returned to the PCM;
  • New explicit build/teardown routines that centralize allocation and initialization ordering, reducing accidental use of uninitialized data and making the model lifecycle explicit;
  • Add native read/write helpers for NetCDF that centralize all low-level NetCDF interactions with major improvements (and more simplicity) compared to legacy PEM/PCM I/O (see the modules "io_netcdf" and "output"). They support reading, creation and writing of "diagevol.nc" (renamed from "diagpem.nc") and start/restart files;
  • Provide well-focused modules ("numerics"/"maths"/"utility"/"display") to host commonly-used primitives:
    • "numerics" defines numerical types and constants for reproducibility, portability across compilers and future transitions (e.g. quadruple precision experiments);
    • "display" provides a single controlled interface for runtime messages, status output and diagnostics, avoiding direct 'print'/'write' to enable silent mode, log redirection, and MPI-safe output in the future.
    • "utility" (new module) hosts generic helpers used throughout the code (e.g. "int2str" or "real2str");
  • Add modules "clim_state_init"/"clim_state_rec" which provide robust read/write logic for "start/startfi/startpem", including 1D fallbacks, mesh conversions and dimension checks. NetCDF file creation is centralized and explicit. Restart files are now self-consistent and future-proof, requiring changes only to affected variables;
  • Add module "atmosphere" which computes pressure fields, reconstructs potential temperature and air mass. It also holds the whole logic to define sigma or hybrid coordinates for altitudes;
  • Add module "geometry" to centrilize dimensions logic and grid conversions routines (including 2 new ones "dyngrd2vect"/"vect2dyngrd");
  • Add module "slopes" to isolate slopes handling;
  • Add module "surface" to isolate surface management. Notably, albedo and emissivity are now fully reconstructed following the PCM settings;
  • Add module "allocation" to check the dimension initialization and centrilize allocation/deallocation;
  • Finalize module decomposition and renaming to consolidate domain-based modules, purpose-based routines and physics/process-based variables;
  • The main program now drives a clearer sequence of conceptual steps (initialization / reading / evolution / update / build / writing) and fails explicitly instead of silently defaulting;
  • Ice table logic is made restart-safe;
  • 'Open'/'read' intrinsic logic is made safe and automatic;
  • Improve discoverability and standardize the data handling (private vs protected vs public);
  • Apply consistent documentation/header style already introduced;
  • Update deftank scripts to reflect new names and launch wrappers;

This revision is a structural milestone aiming to be behavior-preserving where possible. It has been tested via compilation and short integration runs. However, due to extensive renames, moves, and API changes, full validation is still ongoing.
Note: the revision includes one (possibly two) easter egg hidden in the code for future archaeologists of the PEM. No physics were harmed.
JBC

File:
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/soil_therm_inertia.F90

    r4064 r4065  
    1 MODULE soil_thermal_inertia
     1MODULE soil_therm_inertia
    22!-----------------------------------------------------------------------
    33! NAME
    4 !     soil_thermal_inertia
     4!     soil_therm_inertia
    55!
    66! DESCRIPTION
    … …  
    1616!-----------------------------------------------------------------------
    1717
     18! DEPENDENCIES
     19! ------------
     20use numerics, only: dp, di, k4, minieps
     21
    1822! DECLARATION
    1923! -----------
    2024implicit none
    2125
    22 ! MODULE PARAMETERS
    23 ! -----------------
    24 real, parameter :: reg_inertie_thresold = 370. ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2]
    25 real, parameter :: reg_inertie_minvalue = 50.  ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
    26 real, parameter :: reg_inertie_maxvalue = 370. ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
    27 real, parameter :: P610 = 610.0                ! current average pressure on Mars [Pa]
    28 real, parameter :: C = 0.0015                  ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [unitless]
    29 real, parameter :: K = 8.1*1e4                 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [Torr, or 133.3Pa]
    30 real, parameter :: Pa2Torr = 1./133.3          ! Conversion from Pa to tor [Pa/Torr]
     26! PARAMETERS
     27! ----------
     28real(dp), parameter :: reg_inertie_thresold = 370._dp ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2]
     29real(dp), parameter :: reg_inertie_minvalue = 50._dp  ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
     30real(dp), parameter :: reg_inertie_maxvalue = 370._dp ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
     31real(dp), parameter :: C = 0.0015_dp                  ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [unitless]
     32real(dp), parameter :: K = 8.1*1e4_dp                 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [Torr, or 133.3Pa]
     33real(dp), parameter :: Pa2Torr = 1./133.3_dp          ! Conversion from Pa to tor [Pa/Torr]
    3134
    3235contains
    … …  
    5356! DEPENDENCIES
    5457! ------------
    55 use constants_marspem_mod, only: porosity
     58use planet, only: porosity
    5659
    5760! DECLARATION
    … …  
    6164! ARGUMENTS
    6265! ---------
    63 logical, intent(in)  :: ispureice           ! Boolean to check if ice is massive or just pore filling
    64 real,    intent(in)  :: pore_filling        ! ice pore filling in each layer (1)
    65 real,    intent(in)  :: surf_thermalinertia ! surface thermal inertia (J/m^2/K/s^1/2)
    66 real,    intent(out) :: ice_thermalinertia  ! Thermal inertia of ice when present in the pore (J/m^2/K/s^1/2)
     66logical(k4), intent(in)  :: ispureice           ! Boolean to check if ice is massive or just pore filling
     67real(dp),    intent(in)  :: pore_filling        ! ice pore filling in each layer (1)
     68real(dp),    intent(in)  :: surf_thermalinertia ! surface thermal inertia (J/m^2/K/s^1/2)
     69real(dp),    intent(out) :: ice_thermalinertia  ! Thermal inertia of ice when present in the pore (J/m^2/K/s^1/2)
    6770
    6871! LOCAL VARIABLES
    6972! ---------------
    70 real :: inertie_purewaterice = 2100 ! 2050 is better according to my computations with the formula from Siegler et al., 2012, but let's follow Mellon et al. (2004)
     73real(dp) :: inertie_purewaterice = 2100._dp ! 2050 is better according to my computations with the formula from Siegler et al., 2012, but let's follow Mellon et al. (2004)
    7174
    7275! CODE
    … …  
    7679else
    7780    ice_thermalinertia = sqrt(surf_thermalinertia**2 + porosity*pore_filling*inertie_purewaterice**2) ! Siegler et al., 2012
    78 endif
     81end if
    7982
    8083END SUBROUTINE get_ice_TI
    … …  
    8285
    8386!=======================================================================
    84 SUBROUTINE update_soil_TI(ngrid,nslope,nsoil,tendencies_waterice,waterice,p_avg_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
     87SUBROUTINE update_soil_TI(tendencies_waterice,waterice,p_avg_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI)
    8588!-----------------------------------------------------------------------
    8689! NAME
    … …  
    101104! DEPENDENCIES
    102105! ------------
    103 use comsoil_h,             only: volcapa
    104 use soil,                  only: layer_PEM, inertiedat_PEM, depth_breccia, depth_bedrock, index_breccia, index_bedrock, reg_thprop_dependp
    105 use constants_marspem_mod, only: TI_breccia, TI_bedrock, TI_regolith_avg
     106use geometry, only: ngrid, nslope, nsoil
     107use soil,     only: volcapa, layer, inertiedat, depth_breccia, depth_bedrock, index_breccia, index_bedrock, reg_thprop_dependp
     108use planet,   only: TI_breccia, TI_bedrock, TI_regolith_avg, P610
     109use display,  only: print_msg
    106110
    107111! DECLARATION
    … …  
    111115! ARGUMENTS
    112116! ---------
    113 integer,                             intent(in)    :: ngrid, nslope, nsoil ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil
    114 real,                                intent(in)    :: p_avg_new            ! Global average surface pressure [Pa]
    115 real, dimension(ngrid,nslope),       intent(in)    :: tendencies_waterice  ! Tendencies on the water ice [kg/m^2/year]
    116 real, dimension(ngrid,nslope),       intent(in)    :: waterice             ! Surface Water ice [kg/m^2]
    117 real, dimension(ngrid,nslope),       intent(in)    :: icetable_depth       ! Ice table depth [m]
    118 real, dimension(ngrid,nslope),       intent(in)    :: icetable_thickness   ! Ice table thickness [m]
    119 real, dimension(ngrid,nsoil,nslope), intent(in)    :: ice_porefilling      ! Ice porefilling [m^3/m^3]
    120 logical,                             intent(in)    :: icetable_equilibrium, icetable_dynamic ! Computing method for ice table
    121 real, dimension(ngrid,nsoil,nslope), intent(inout) :: TI_PEM               ! Soil Thermal Inertia [J/m^2/K/s^1/2]
     117real(dp),                   intent(in)    :: p_avg_new            ! Global average surface pressure [Pa]
     118real(dp), dimension(:,:),   intent(in)    :: tendencies_waterice  ! Tendencies on the water ice [kg/m^2/year]
     119real(dp), dimension(:,:),   intent(in)    :: waterice             ! Surface Water ice [kg/m^2]
     120real(dp), dimension(:,:),   intent(in)    :: icetable_depth       ! Ice table depth [m]
     121real(dp), dimension(:,:),   intent(in)    :: icetable_thickness   ! Ice table thickness [m]
     122real(dp), dimension(:,:,:), intent(in)    :: ice_porefilling      ! Ice porefilling [m^3/m^3]
     123logical(k4),                intent(in)    :: icetable_equilibrium, icetable_dynamic ! Computing method for ice table
     124real(dp), dimension(:,:,:), intent(inout) :: TI                   ! Soil Thermal Inertia [J/m^2/K/s^1/2]
    122125
    123126! LOCAL VARIABLES
    124127! ---------------
    125 integer                       :: ig, islope, isoil, iref, iend ! Loop variables
    126 real, dimension(ngrid,nslope) :: regolith_inertia              ! Thermal inertia of the regolith (first layer in the GCM) [J/m^2/K/s^1/2]
    127 real                          :: delta                         ! Difference of depth to compute the  thermal inertia in a mixed layer [m]
    128 real                          :: ice_bottom_depth              ! Bottom depth of the subsurface ice [m]
    129 real                          :: d_part                        ! Regolith particle size [micrometer]
    130 real                          :: ice_inertia                   ! Inertia of water ice [SI]
     128integer(di)                       :: ig, islope, isoil, iref, iend ! Loop variables
     129real(dp), dimension(ngrid,nslope) :: regolith_inertia              ! Thermal inertia of the regolith (first layer in the GCM) [J/m^2/K/s^1/2]
     130real(dp)                          :: delta                         ! Difference of depth to compute the  thermal inertia in a mixed layer [m]
     131real(dp)                          :: ice_bottom_depth              ! Bottom depth of the subsurface ice [m]
     132real(dp)                          :: d_part                        ! Regolith particle size [micrometer]
     133real(dp)                          :: ice_inertia                   ! Inertia of water ice [SI]
    131134
    132135! CODE
    133136! ----
    134 write(*,*) "> Updating soil properties"
     137call print_msg("> Updating soil properties")
    135138
    136139! 1. Modification of the regolith thermal inertia.
    137140do islope = 1,nslope
    138     regolith_inertia(:,islope) = inertiedat_PEM(:,1)
     141    regolith_inertia(:,islope) = inertiedat(:,1)
    139142    do ig = 1,ngrid
    140         if (tendencies_waterice(ig,islope) < -1.e-5 .and. waterice(ig,islope) == 0) regolith_inertia(ig,islope) = TI_regolith_avg
     143        if (tendencies_waterice(ig,islope) < -1.e-5_dp .and. abs(waterice(ig,islope)) < minieps) regolith_inertia(ig,islope) = TI_regolith_avg
    141144        if (reg_thprop_dependp) then
    142             if (TI_PEM(ig,1,islope) < reg_inertie_thresold) then
    143                 d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Torr)**(0.6)))**(-1./(0.11*log10(P610*Pa2Torr/K))) ! compute particle size, in micrometer
    144                 regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Torr)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Torr/K)))
     145            if (TI(ig,1,islope) < reg_inertie_thresold) then
     146                d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Torr)**0.6))**(-1./(0.11_dp*log10(P610*Pa2Torr/K))) ! compute particle size, in micrometer
     147                regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Torr)**0.6*d_part**(-0.11_dp*log10(p_avg_new*Pa2Torr/K)))
    145148                if (regolith_inertia(ig,islope) > reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue
    146149                if (regolith_inertia(ig,islope) < reg_inertie_minvalue) regolith_inertia(ig,islope) = reg_inertie_minvalue
    147             endif
    148         endif
    149     enddo
    150 enddo
     150            end if
     151        end if
     152    end do
     153end do
    151154
    152155! 2. Build new Thermal Inertia
    … …  
    154157    do ig = 1,ngrid
    155158        do isoil = 1,index_breccia
    156             TI_PEM(ig,isoil,islope) = regolith_inertia(ig,islope)
    157         enddo
     159            TI(ig,isoil,islope) = regolith_inertia(ig,islope)
     160        end do
    158161        if (regolith_inertia(ig,islope) < TI_breccia) then
    159162!!! transition
    160163            delta = depth_breccia
    161             TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/              &
    162                                                   (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
    163                                                   ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
     164            TI(ig,index_breccia + 1,islope) = sqrt((layer(index_breccia + 1) - layer(index_breccia))/                  &
     165                                                  (((delta - layer(index_breccia))/(TI(ig,index_breccia,islope)**2)) + &
     166                                                  ((layer(index_breccia + 1) - delta)/(TI_breccia**2))))
    164167            do isoil = index_breccia + 2,index_bedrock
    165                 TI_PEM(ig,isoil,islope) = TI_breccia
    166             enddo
     168                TI(ig,isoil,islope) = TI_breccia
     169            end do
    167170        else ! we keep the high ti values
    168171            do isoil = index_breccia + 1,index_bedrock
    169                 TI_PEM(ig,isoil,islope) = TI_PEM(ig,index_breccia,islope)
    170             enddo
    171         endif ! TI PEM and breccia comparison
     172                TI(ig,isoil,islope) = TI(ig,index_breccia,islope)
     173            end do
     174        end if ! TI PEM and breccia comparison
    172175!!! transition
    173176        delta = depth_bedrock
    174         TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/              &
    175                                               (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
    176                                               ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
     177        TI(ig,index_bedrock + 1,islope) = sqrt((layer(index_bedrock + 1) - layer(index_bedrock))/                  &
     178                                              (((delta - layer(index_bedrock))/(TI(ig,index_bedrock,islope)**2)) + &
     179                                              ((layer(index_bedrock + 1) - delta)/(TI_bedrock**2))))
    177180        do isoil = index_bedrock + 2,nsoil
    178             TI_PEM(ig,isoil,islope) = TI_bedrock
    179         enddo
    180     enddo ! ig
    181 enddo ! islope
     181            TI(ig,isoil,islope) = TI_bedrock
     182        end do
     183    end do ! ig
     184end do ! islope
    182185
    183186! 3. Build new TI in case of ice table
    184187do ig = 1,ngrid
    185188    do islope = 1,nslope
    186         if (icetable_depth(ig,islope) > -1.e-6) then
     189        if (icetable_depth(ig,islope) > -1.e-6_dp) then
    187190        ! 3.0 Case where it is perennial ice
    188             if (icetable_depth(ig,islope) < 1.e-10) then
    189                 call get_ice_TI(.true.,1.,0.,ice_inertia)
     191            if (icetable_depth(ig,islope) < minieps) then
     192                call get_ice_TI(.true.,1._dp,0._dp,ice_inertia)
    190193                do isoil = 1,nsoil
    191                     TI_PEM(ig,isoil,islope) = ice_inertia
    192                 enddo
     194                    TI(ig,isoil,islope) = ice_inertia
     195                end do
    193196            else
    194197                if (icetable_equilibrium) then
    195                     call get_ice_TI(.false.,1.,regolith_inertia(ig,islope),ice_inertia)
     198                    call get_ice_TI(.false.,1._dp,regolith_inertia(ig,islope),ice_inertia)
    196199                    ! 3.1.1 find the index of the mixed layer
    197200                    iref = 0 ! initialize iref
    198201                    do isoil = 1,nsoil ! loop on layers to find the beginning of the ice table
    199                         if (icetable_depth(ig,islope) >= layer_PEM(isoil)) then
     202                        if (icetable_depth(ig,islope) >= layer(isoil)) then
    200203                            iref = isoil ! pure regolith layer up to here
    201204                        else
    202205                            exit ! correct iref was obtained in previous cycle
    203                         endif
    204                     enddo
     206                        end if
     207                    end do
    205208                    ! 3.1.2 find the index of the end of the ice table
    206209                    iend = 0 ! initialize iend
    207210                    ice_bottom_depth = icetable_depth(ig,islope) + icetable_thickness(ig,islope)
    208211                    do isoil = 1,nsoil ! loop on layers to find the end of the ice table
    209                         if (ice_bottom_depth >= layer_PEM(isoil)) then
     212                        if (ice_bottom_depth >= layer(isoil)) then
    210213                            iend = isoil ! pure regolith layer up to here
    211214                        else
    212215                            exit ! correct iref was obtained in previous cycle
    213                         endif
    214                     enddo
     216                        end if
     217                    end do
    215218                    ! 3.2 Build the new ti
    216219                    if (iref < nsoil) then
    … …  
    218221                            ! Ice table begins and end in the same mixture with three components
    219222                            if (iref /= 0) then ! mixed layer
    220                                 TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                                  &
    221                                                              (((icetable_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
    222                                                              ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) +            &
    223                                                              ((layer_PEM(iref + 1) - ice_bottom_depth)/(TI_PEM(ig,iref + 1,islope)**2))))
     223                                TI(ig,iref + 1,islope) = sqrt((layer(iref + 1) - layer(iref))/                                      &
     224                                                             (((icetable_depth(ig,islope) - layer(iref))/(TI(ig,iref,islope)**2)) + &
     225                                                             ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) +    &
     226                                                             ((layer(iref + 1) - ice_bottom_depth)/(TI(ig,iref + 1,islope)**2))))
    224227                            else ! first layer is already a mixed layer
    225228                                ! (ie: take layer(iref=0)=0)
    226                                 TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/                                                &
    227                                                       (((icetable_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) +           &
     229                                TI(ig,1,islope) = sqrt((layer(1))/                                                        &
     230                                                      (((icetable_depth(ig,islope))/(TI(ig,1,islope)**2)) +               &
    228231                                                      ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) + &
    229                                                       ((layer_PEM(2) - ice_bottom_depth)/(TI_PEM(ig,2,islope)**2))))
    230                             endif ! of if (iref /= 0)
     232                                                      ((layer(2) - ice_bottom_depth)/(TI(ig,2,islope)**2))))
     233                            end if ! of if (iref /= 0)
    231234                        else
    232235                            if (iref /= 0) then ! mixed layer
    233                                 TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                                  &
    234                                                              (((icetable_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
    235                                                              ((layer_PEM(iref + 1) - icetable_depth(ig,islope))/(ice_inertia**2))))
     236                                TI(ig,iref + 1,islope) = sqrt((layer(iref + 1) - layer(iref))/                                      &
     237                                                             (((icetable_depth(ig,islope) - layer(iref))/(TI(ig,iref,islope)**2)) + &
     238                                                             ((layer(iref + 1) - icetable_depth(ig,islope))/(ice_inertia**2))))
    236239                            else ! first layer is already a mixed layer
    237240                                ! (ie: take layer(iref=0)=0)
    238                                 TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/                                      &
    239                                                       (((icetable_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) + &
    240                                                       ((layer_PEM(1) - icetable_depth(ig,islope))/(ice_inertia**2))))
    241                             endif ! of if (iref /= 0)
    242                         endif ! iref == iend
    243 
    244                         TI_PEM(ig,iref + 2:iend,islope) = ice_inertia
     241                                TI(ig,1,islope) = sqrt((layer(1))/                                          &
     242                                                      (((icetable_depth(ig,islope))/(TI(ig,1,islope)**2)) + &
     243                                                      ((layer(1) - icetable_depth(ig,islope))/(ice_inertia**2))))
     244                            end if ! of if (iref /= 0)
     245                        end if ! iref == iend
     246
     247                        TI(ig,iref + 2:iend,islope) = ice_inertia
    245248                        if (iend < nsoil) then
    246                             TI_PEM(ig,iend + 1,islope) = sqrt((layer_PEM(iend + 1) - layer_PEM(iend))/                         &
    247                                                          (((ice_bottom_depth - layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + &
    248                                                          ((layer_PEM(iend + 1) - ice_bottom_depth)/(ice_inertia**2))))
    249                         endif
    250                     endif ! of if (iref < nsoilmx)
     249                            TI(ig,iend + 1,islope) = sqrt((layer(iend + 1) - layer(iend))/                             &
     250                                                         (((ice_bottom_depth - layer(iend))/(TI(ig,iend,islope)**2)) + &
     251                                                         ((layer(iend + 1) - ice_bottom_depth)/(ice_inertia**2))))
     252                        end if
     253                    end if ! of if (iref < nsoil)
    251254                else if (icetable_dynamic) then
    252255                    do  isoil = 1,nsoil
    253                         call get_ice_TI(.false.,ice_porefilling(ig,isoil,islope),regolith_inertia(ig,islope),TI_PEM(ig,isoil,islope))
    254                     enddo
    255                 endif ! of if icetable_equilibrium
    256             endif ! permanent glaciers
    257         endif ! icetable_depth(ig,islope) > 0.
    258     enddo !islope
    259 enddo !ig
     256                        call get_ice_TI(.false.,ice_porefilling(ig,isoil,islope),regolith_inertia(ig,islope),TI(ig,isoil,islope))
     257                    end do
     258                end if ! of if icetable_equilibrium
     259            end if ! permanent glaciers
     260        end if ! icetable_depth(ig,islope) > 0.
     261    end do !islope
     262end do !ig
    260263
    261264END SUBROUTINE update_soil_TI
    262265!=======================================================================
    263266
    264 END MODULE soil_thermal_inertia
     267END MODULE soil_therm_inertia
Note: See TracChangeset for help on using the changeset viewer.