Ignore:
Timestamp:
Nov 19, 2024, 5:44:27 PM (5 days ago)
Author:
jbclement
Message:

PEM:
Computation of <soil thermal intertia> and <H2O mass subsurface/surface exchange> according to the presence of subsurface ice provided by the (Norbert's) dynamic method + few cleanings.
JBC

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r3347 r3525  
    99!!!
    1010!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     11
     12! Constants:
     13real, 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]
     14real, parameter :: reg_inertie_minvalue = 50.  ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
     15real, parameter :: reg_inertie_maxvalue = 370. ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
     16real, parameter :: P610 = 610.0                ! current average pressure on Mars [Pa]
     17real, parameter :: C = 0.0015                  ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [unitless]
     18real, parameter :: K = 8.1*1e4                 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [Torr, or 133.3Pa]
     19real, parameter :: Pa2Torr = 1./133.3          ! Conversion from Pa to tor [Pa/Torr]
    1120
    1221!=======================================================================
    … …  
    5059    ice_thermalinertia = inertie_purewaterice
    5160else
    52     ice_thermalinertia = sqrt(surf_thermalinertia**2 + porosity*pore_filling*inertie_purewaterice**2)  ! Siegler et al., 2012
     61    ice_thermalinertia = sqrt(surf_thermalinertia**2 + porosity*pore_filling*inertie_purewaterice**2) ! Siegler et al., 2012
    5362endif
    5463
    … …  
    5665!=======================================================================
    5766
    58 SUBROUTINE update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tendencies_waterice,waterice,p_avg_new,ice_depth,ice_thickness,TI_PEM)
     67SUBROUTINE update_soil_thermalproperties(ngrid,nslope,nsoil,tendencies_waterice,waterice,p_avg_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
    5968
    6069use comsoil_h,             only: volcapa
    … …  
    6574
    6675! Input:
    67 integer,                       intent(in) :: ngrid, nslope, nsoil_PEM ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil
    68 real,                          intent(in) :: p_avg_new                ! Global average surface pressure [Pa]
    69 real, dimension(ngrid,nslope), intent(in) :: tendencies_waterice      ! Tendencies on the water ice [kg/m^2/year]
    70 real, dimension(ngrid,nslope), intent(in) :: waterice                 ! Surface Water ice [kg/m^2]
    71 real, dimension(ngrid,nslope), intent(in) :: ice_depth                ! Ice table depth [m]
    72 real, dimension(ngrid,nslope), intent(in) :: ice_thickness            ! Ice table thickness [m]
     76integer,                             intent(in) :: ngrid, nslope, nsoil ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil
     77real,                                intent(in) :: p_avg_new            ! Global average surface pressure [Pa]
     78real, dimension(ngrid,nslope),       intent(in) :: tendencies_waterice  ! Tendencies on the water ice [kg/m^2/year]
     79real, dimension(ngrid,nslope),       intent(in) :: waterice             ! Surface Water ice [kg/m^2]
     80real, dimension(ngrid,nslope),       intent(in) :: icetable_depth       ! Ice table depth [m]
     81real, dimension(ngrid,nslope),       intent(in) :: icetable_thickness   ! Ice table thickness [m]
     82real, dimension(ngrid,nsoil,nslope), intent(in) :: ice_porefilling      ! Ice porefilling [m^3/m^3]
     83logical,                             intent(in) :: icetable_equilibrium, icetable_dynamic ! Computing method for ice table
    7384
    7485! Outputs:
    75 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! Soil Thermal Inertia [J/m^2/K/s^1/2]
    76 
    77 ! Constants:
    78 real :: 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]
    79 real :: reg_inertie_minvalue = 50.  ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
    80 real :: reg_inertie_maxvalue = 370. ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
    81 real :: ice_inertia                 ! Inertia of water ice [SI]
    82 real :: P610 = 610.0                ! current average pressure on Mars [Pa]
    83 real :: C = 0.0015                  ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [uniteless]
    84 real :: K = 8.1*1e4                 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [tor, or 133.3Pa]
    85 real :: Pa2Tor = 1./133.3           ! Conversion from Pa to tor [Pa/tor]
     86real, dimension(ngrid,nsoil,nslope), intent(inout) :: TI_PEM ! Soil Thermal Inertia [J/m^2/K/s^1/2]
    8687
    8788! Local variables:
    … …  
    9192real    :: ice_bottom_depth                       ! Bottom depth of the subsurface ice [m]
    9293real    :: d_part                                 ! Regolith particle size [micrometer]
     94real    :: ice_inertia                            ! Inertia of water ice [SI]
    9395
    9496write(*,*) "Update soil properties"
    … …  
    101103        if (reg_thprop_dependp) then
    102104            if (TI_PEM(ig,1,islope) < reg_inertie_thresold) then
    103                 d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Tor)**(0.6)))**(-1./(0.11*log10(P610*Pa2Tor/K))) ! compute particle size, in micrometer
    104                 regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Tor)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Tor/K)))
     105                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
     106                regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Torr)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Torr/K)))
    105107                if (regolith_inertia(ig,islope) > reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue
    106108                if (regolith_inertia(ig,islope) < reg_inertie_minvalue) regolith_inertia(ig,islope) = reg_inertie_minvalue
    … …  
    135137                                              (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
    136138                                              ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
    137         do isoil = index_bedrock + 2,nsoil_PEM
     139        do isoil = index_bedrock + 2,nsoil
    138140            TI_PEM(ig,isoil,islope) = TI_bedrock
    139141        enddo
    … …  
    141143enddo ! islope
    142144
    143 !  3. Build new TI in case of ice table
     145! 3. Build new TI in case of ice table
    144146do ig = 1,ngrid
    145147    do islope = 1,nslope
    146         if (ice_depth(ig,islope) > -1.e-6) then
     148        if (icetable_depth(ig,islope) > -1.e-6) then
    147149        ! 3.0 Case where it is perennial ice
    148             if (ice_depth(ig,islope) < 1.e-10) then
     150            if (icetable_depth(ig,islope) < 1.e-10) then
    149151                call ice_thermal_properties(.true.,1.,0.,ice_inertia)
    150                 do isoil = 1,nsoil_PEM
     152                do isoil = 1,nsoil
    151153                    TI_PEM(ig,isoil,islope) = ice_inertia
    152154                enddo
    153155            else
    154                 call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),ice_inertia)
    155         ! 3.1.1 find the index of the mixed layer
    156                 iref = 0 ! initialize iref
    157                 do isoil = 1,nsoil_PEM ! loop on layers to find the beginning of the ice table
    158                     if (ice_depth(ig,islope) >= layer_PEM(isoil)) then
    159                         iref = isoil ! pure regolith layer up to here
    160                     else
    161                         exit ! correct iref was obtained in previous cycle
    162                     endif
    163                 enddo
    164         ! 3.1.2 find the index of the end of the ice table
    165                 iend = 0 ! initialize iend
    166                 ice_bottom_depth = ice_depth(ig,islope) + ice_thickness(ig,islope)
    167                 do isoil = 1,nsoil_PEM ! loop on layers to find the end of the ice table
    168                     if (ice_bottom_depth >= layer_PEM(isoil)) then
    169                         iend = isoil ! pure regolith layer up to here
    170                     else
    171                         exit ! correct iref was obtained in previous cycle
    172                     endif
    173                 enddo
    174         ! 3.2 Build the new ti
    175                 if (iref < nsoil_PEM) then
    176                     if (iref == iend) then
    177                         ! Ice table begins and end in the same mixture Mixtures with three components
    178                         if (iref /= 0) then ! mixed layer
    179                             TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                             &
    180                                                          (((ice_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
    181                                                          ((ice_bottom_depth - ice_depth(ig,islope))/(ice_inertia**2)) +            &
    182                                                          ((layer_PEM(iref + 1) - ice_bottom_depth)/(TI_PEM(ig,iref + 1,islope)**2))))
    183                         else ! first layer is already a mixed layer
    184                             ! (ie: take layer(iref=0)=0)
    185                             TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/ &
    186                                                   (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) +           &
    187                                                   ((ice_bottom_depth - ice_depth(ig,islope))/(ice_inertia**2)) + &
    188                                                   ((layer_PEM(2) - ice_bottom_depth)/(TI_PEM(ig,2,islope)**2))))
    189                         endif ! of if (iref /= 0)
    190                     else
    191                         if (iref /= 0) then ! mixed layer
    192                             TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                             &
    193                                                          (((ice_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
    194                                                          ((layer_PEM(iref + 1) - ice_depth(ig,islope))/(ice_inertia**2))))
    195                         else ! first layer is already a mixed layer
    196                             ! (ie: take layer(iref=0)=0)
    197                             TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/                                 &
    198                                                   (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) + &
    199                                                   ((layer_PEM(1) - ice_depth(ig,islope))/(ice_inertia**2))))
    200                         endif ! of if (iref /= 0)
    201                     endif ! iref == iend
    202         ! 3.3 Build the new ti
    203                     do isoil = iref + 2,iend
    204                         TI_PEM(ig,isoil,islope) = ice_inertia
     156                if (icetable_equilibrium) then
     157                    call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),ice_inertia)
     158                    ! 3.1.1 find the index of the mixed layer
     159                    iref = 0 ! initialize iref
     160                    do isoil = 1,nsoil ! loop on layers to find the beginning of the ice table
     161                        if (icetable_depth(ig,islope) >= layer_PEM(isoil)) then
     162                            iref = isoil ! pure regolith layer up to here
     163                        else
     164                            exit ! correct iref was obtained in previous cycle
     165                        endif
    205166                    enddo
    206                     if (iend < nsoil_PEM) then
    207                         TI_PEM(ig,iend + 1,islope) = sqrt((layer_PEM(iend + 1) - layer_PEM(iend))/                         &
    208                                                      (((ice_bottom_depth - layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + &
    209                                                      ((layer_PEM(iend + 1) - ice_bottom_depth)/(ice_inertia**2))))
    210                     endif
    211                 endif ! of if (iref < nsoilmx)
     167                    ! 3.1.2 find the index of the end of the ice table
     168                    iend = 0 ! initialize iend
     169                    ice_bottom_depth = icetable_depth(ig,islope) + icetable_thickness(ig,islope)
     170                    do isoil = 1,nsoil ! loop on layers to find the end of the ice table
     171                        if (ice_bottom_depth >= layer_PEM(isoil)) then
     172                            iend = isoil ! pure regolith layer up to here
     173                        else
     174                            exit ! correct iref was obtained in previous cycle
     175                        endif
     176                    enddo
     177                    ! 3.2 Build the new ti
     178                    if (iref < nsoil) then
     179                        if (iref == iend) then
     180                            ! Ice table begins and end in the same mixture with three components
     181                            if (iref /= 0) then ! mixed layer
     182                                TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                                  &
     183                                                             (((icetable_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
     184                                                             ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) +            &
     185                                                             ((layer_PEM(iref + 1) - ice_bottom_depth)/(TI_PEM(ig,iref + 1,islope)**2))))
     186                            else ! first layer is already a mixed layer
     187                                ! (ie: take layer(iref=0)=0)
     188                                TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/                                                &
     189                                                      (((icetable_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) +           &
     190                                                      ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) + &
     191                                                      ((layer_PEM(2) - ice_bottom_depth)/(TI_PEM(ig,2,islope)**2))))
     192                            endif ! of if (iref /= 0)
     193                        else
     194                            if (iref /= 0) then ! mixed layer
     195                                TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/                                  &
     196                                                             (((icetable_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + &
     197                                                             ((layer_PEM(iref + 1) - icetable_depth(ig,islope))/(ice_inertia**2))))
     198                            else ! first layer is already a mixed layer
     199                                ! (ie: take layer(iref=0)=0)
     200                                TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/                                      &
     201                                                      (((icetable_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) + &
     202                                                      ((layer_PEM(1) - icetable_depth(ig,islope))/(ice_inertia**2))))
     203                            endif ! of if (iref /= 0)
     204                        endif ! iref == iend
     205
     206                        TI_PEM(ig,iref + 2:iend,islope) = ice_inertia
     207                        if (iend < nsoil) then
     208                            TI_PEM(ig,iend + 1,islope) = sqrt((layer_PEM(iend + 1) - layer_PEM(iend))/                         &
     209                                                         (((ice_bottom_depth - layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + &
     210                                                         ((layer_PEM(iend + 1) - ice_bottom_depth)/(ice_inertia**2))))
     211                        endif
     212                    endif ! of if (iref < nsoilmx)
     213                else if (icetable_dynamic) then
     214                    do  isoil = 1,nsoil
     215                        call ice_thermal_properties(.false.,ice_porefilling(ig,isoil,islope),regolith_inertia(ig,islope),TI_PEM(ig,isoil,islope))
     216                    enddo
     217                endif ! of if icetable_equilibrium
    212218            endif ! permanent glaciers
    213         endif ! ice_depth(ig,islope) > 0.
    214 !        write(*,*) 'TI=', TI_PEM(ig,:,islope)
     219        endif ! icetable_depth(ig,islope) > 0.
    215220    enddo !islope
    216221enddo !ig
Note: See TracChangeset for help on using the changeset viewer.