Ignore:
Timestamp:
Nov 19, 2024, 5:44:27 PM (2 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

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/changelog.txt

    r3518 r3525  
    489489== 14/11/2024 == JBC
    490490Committing the right correction in previous commit (r3516) would have been better...
     491
     492== 19/11/2024 == JBC
     493Computation 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.

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

    r3178 r3525  
    9090endif ! of if (firstcall)
    9191
    92 IF (.not. firstcall) THEN
     92if (.not. firstcall) THEN
    9393! 2. Compute soil temperatures
    9494! First layer:
    … …  
    120120END SUBROUTINE compute_tsoil_pem
    121121
    122 
     122!=======================================================================
    123123
    124124SUBROUTINE ini_tsoil_pem(ngrid,nsoil,therm_i,tsurf,tsoil)
    … …  
    221221END SUBROUTINE ini_tsoil_pem
    222222
    223 
    224 
    225223END MODULE compute_soiltemp_mod

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

    r3512 r3525  
    1818!-----------------------------------------------------------------------
    1919contains
     20
    2021!-----------------------------------------------------------------------
    2122SUBROUTINE ini_ice_table(ngrid,nslope,nsoil)
    … …  
    4546
    4647!------------------------------------------------------------------------------------------------------
    47 SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness)
     48SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness)
    4849!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    4950!!!
    … …  
    5152!!! density at the surface and at depth.
    5253!!! Computations are made following the methods in Schorgofer et al., 2005
    53 !!! This SUBROUTINE only gives the ice table at equilibrium and does not consider exchange with the atmosphere
     54!!! This subroutine only gives the ice table at equilibrium and does not consider exchange with the atmosphere
    5455!!!
    5556!!! Author: LL
    … …  
    6465! Inputs
    6566! ------
    66 integer,                                 intent(in) :: ngrid, nslope, nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM
    67 logical, dimension(ngrid),               intent(in) :: watercaptag              ! Boolean to check the presence of a perennial glacier
    68 real, dimension(ngrid,nslope),           intent(in) :: rhowatersurf_ave         ! Water density at the surface, yearly averaged [kg/m^3]
    69 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: rhowatersoil_ave         ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged [kg/m^3]
    70 real, dimension(ngrid,nslope),           intent(in) :: regolith_inertia         ! Thermal inertia of the regolith layer [SI]
     67integer,                             intent(in) :: ngrid, nslope, nsoil ! Size of the physical grid, number of subslope, number of soil layer in the PEM
     68logical, dimension(ngrid),           intent(in) :: watercaptag          ! Boolean to check the presence of a perennial glacier
     69real, dimension(ngrid,nslope),       intent(in) :: rhowatersurf_ave     ! Water density at the surface, yearly averaged [kg/m^3]
     70real, dimension(ngrid,nsoil,nslope), intent(in) :: rhowatersoil_ave     ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged [kg/m^3]
     71real, dimension(ngrid,nslope),       intent(in) :: regolith_inertia     ! Thermal inertia of the regolith layer [SI]
    7172! Ouputs
    7273! ------
    … …  
    7677! ------
    7778integer                       :: ig, islope, isoil, isoilend ! loop variables
    78 real, dimension(nsoil_PEM)    :: diff_rho                    ! difference of water vapor density between the surface and at depth [kg/m^3]
     79real, dimension(nsoil)        :: diff_rho                    ! difference of water vapor density between the surface and at depth [kg/m^3]
    7980real                          :: ice_table_end               ! depth of the end of the ice table  [m]
    8081real, dimension(ngrid,nslope) :: previous_icetable_depth     ! Ice table computed at previous ice depth [m]
    … …  
    8788    if (watercaptag(ig)) then
    8889        ice_table_beg(ig,:) = 0.
    89         ice_table_thickness(ig,:) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
     90        ice_table_thickness(ig,:) = layer_PEM(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
    9091    else
    9192        do islope = 1,nslope
    9293            ice_table_beg(ig,islope) = -1.
    9394            ice_table_thickness(ig,islope) = 0.
    94             do isoil = 1,nsoil_PEM
     95            do isoil = 1,nsoil
    9596                diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)
    9697            enddo
    9798            if (diff_rho(1) > 0) then ! ice is at the surface
    9899                ice_table_beg(ig,islope) = 0.
    99                 do isoilend = 2,nsoil_PEM - 1
     100                do isoilend = 2,nsoil - 1
    100101                    if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then
    101102                        call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end)
    … …  
    105106                enddo
    106107            else
    107                 do isoil = 1,nsoil_PEM - 1 ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
     108                do isoil = 1,nsoil - 1 ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
    108109                    if (diff_rho(isoil) < 0 .and. diff_rho(isoil + 1) > 0.) then
    109110                        call findroot(diff_rho(isoil),diff_rho(isoil + 1),mlayer_PEM(isoil),mlayer_PEM(isoil + 1),ice_table_beg(ig,islope))
    110111                        ! Now let's find the end of the ice table
    111                         ice_table_thickness(ig,islope) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
    112                         do isoilend = isoil + 1,nsoil_PEM - 1
     112                        ice_table_thickness(ig,islope) = layer_PEM(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
     113                        do isoilend = isoil + 1,nsoil - 1
    113114                            if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then
    114115                                call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end)
    … …  
    140141
    141142!-----------------------------------------------------------------------
    142 SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil_PEM,former_ice_table_thickness,new_ice_table_thickness,ice_table_depth,tsurf,tsoil,delta_m_h2o)
     143SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil,icetable_thickness_old,ice_porefilling_old,tsurf,tsoil,delta_m_h2o)
    143144!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    144145!!!
    … …  
    161162! Inputs
    162163! ------
    163 integer,                                 intent(in) :: ngrid, nslope, nsoil_PEM   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
    164 real, dimension(ngrid,nslope),           intent(in) :: former_ice_table_thickness ! ice table thickness at the former iteration [m]
    165 real, dimension(ngrid,nslope),           intent(in) :: new_ice_table_thickness    ! ice table thickness at the current iteration [m]
    166 real, dimension(ngrid,nslope),           intent(in) :: ice_table_depth            ! ice table depth [m]
    167 real, dimension(ngrid,nslope),           intent(in) :: tsurf                      ! Surface temperature [K]
    168 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil                      ! Soil temperature [K]
     164integer,                             intent(in) :: ngrid, nslope, nsoil   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
     165real, dimension(ngrid,nslope),       intent(in) :: icetable_thickness_old ! Ice table thickness at the previous iteration [m]
     166real, dimension(ngrid,nsoil,nslope), intent(in) :: ice_porefilling_old    ! Ice pore filling at the previous iteration [m]
     167real, dimension(ngrid,nslope),       intent(in) :: tsurf                  ! Surface temperature [K]
     168real, dimension(ngrid,nsoil,nslope), intent(in) :: tsoil                  ! Soil temperature [K]
    169169! Outputs
    170170! -------
    … …  
    172172! Locals
    173173!-------
    174 integer                       :: ig, islope, ilay, iref ! loop index
    175 real, dimension(ngrid,nslope) :: rho                    ! density of water ice [kg/m^3]
    176 real, dimension(ngrid,nslope) :: Tice                   ! ice temperature [k]
     174integer :: ig, islope, isoil ! loop index
     175real    :: rho               ! density of water ice [kg/m^3]
     176real    :: Tice              ! ice temperature [k]
    177177! Code
    178178! ----
    179 rho = 0.
    180 Tice = 0.
    181 !1. First let's compute Tice using a linear interpolation between the mlayer level
    182 do ig = 1,ngrid
    183     do islope = 1,nslope
    184         call compute_Tice_pem(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope),ice_table_depth(ig,islope),Tice(ig,islope))
    185         rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2 + 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012
    186     enddo
    187 enddo
    188 
    189 !2. Let's compute the amount of ice that has sublimated in each subslope
    190 do ig = 1,ngrid
    191     do islope = 1,nslope
    192         delta_m_h2o(ig) = delta_m_h2o(ig) + porosity*rho(ig,islope)*(new_ice_table_thickness(ig,islope) - former_ice_table_thickness(ig,islope)) & ! convention > 0. <=> it condenses
    193                           *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
    194     enddo
    195 enddo
     179
     180! Let's compute the amount of ice that has sublimated in each subslope
     181if (icetable_equilibrium) then
     182    delta_m_h2o = 0.
     183    do ig = 1,ngrid
     184        do islope = 1,nslope
     185            call compute_Tice_pem(nsoil,tsoil(ig,:,islope),tsurf(ig,islope),icetable_depth(ig,islope),Tice)
     186            rho = -3.5353e-4*Tice**2 + 0.0351*Tice + 933.5030 ! Rottgers, 2012
     187            delta_m_h2o(ig) = delta_m_h2o(ig) + porosity*rho*(icetable_thickness(ig,islope) - icetable_thickness_old(ig,islope)) & ! convention > 0. <=> it condenses
     188                              *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
     189        enddo
     190    enddo
     191else if (icetable_dynamic) then
     192    delta_m_h2o = 0.
     193    do ig = 1,ngrid
     194        do islope = 1,nslope
     195            do isoil = 1,nsoil
     196                call compute_Tice_pem(nsoil,tsoil(ig,:,islope),tsurf(ig,islope),mlayer_PEM(isoil - 1),Tice)
     197                rho = -3.5353e-4*Tice**2 + 0.0351*Tice + 933.5030 ! Rottgers, 2012
     198                delta_m_h2o(ig) = delta_m_h2o(ig) + rho*(ice_porefilling(ig,isoil,islope) - ice_porefilling_old(ig,isoil,islope)) & ! convention > 0. <=> it condenses
     199                                  *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
     200            enddo
     201        enddo
     202    enddo
     203endif
    196204
    197205END SUBROUTINE compute_massh2o_exchange_ssi
    … …  
    285293
    286294! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
    287 ! 2.0 preliminary: depth to shallowest  ice (discontinuity at interface)
     295! 2.0 preliminary: depth to shallowest ice (discontinuity at interface)
    288296index_firstice = -1
    289297do ilay = 1,nsoilmx_PEM
    … …  
    294302enddo
    295303
    296 ! 2.1: now we can computeCompute d_z (eta* d_z(rho))
     304! 2.1: now we can compute
    297305call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM - 1),dz_eta)
    298306
    … …  
    368376
    369377!-----------------------------------------------------------------------
    370 SUBROUTINE compute_Tice_pem(nsoil, ptsoil, ptsurf, ice_depth, Tice)
     378SUBROUTINE compute_Tice_pem(nsoil,ptsoil,ptsurf,ice_depth,Tice)
    371379!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    372380!!!
    … …  
    383391! Inputs
    384392! ------
    385 integer,                 intent(in) :: nsoil     ! Number of soil layers
    386 real, dimension(nsoil),  intent(in) :: ptsoil    ! Soil temperature (K)
    387 real,                    intent(in) :: ptsurf    ! Soil temperature (K)
    388 real,                    intent(in) :: ice_depth ! Ice depth (m)
     393integer,                intent(in) :: nsoil     ! Number of soil layers
     394real, dimension(nsoil), intent(in) :: ptsoil    ! Soil temperature (K)
     395real,                   intent(in) :: ptsurf    ! Surface temperature (K)
     396real,                   intent(in) :: ice_depth ! Ice depth (m)
    389397
    390398! Outputs
    391 ! ------
     399! -------
    392400real, intent(out) :: Tice ! Ice temperatures (K)
    393401

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

    r3516 r3525  
    223223logical                               :: bool_sublim                      ! logical to check if there is sublimation or not
    224224logical, dimension(:,:),  allocatable :: co2ice_disappeared               ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
    225 real, dimension(:,:),     allocatable :: icetable_thickness_prev_iter     ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
     225real, dimension(:,:),     allocatable :: icetable_thickness_old           ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
     226real, dimension(:,:,:),   allocatable :: ice_porefilling_old              ! ngrid x nslope: Ice pore filling at the previous iteration [m]
    226227real, dimension(:),       allocatable :: delta_h2o_icetablesublim         ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
    227228real, dimension(:),       allocatable :: porefill                         ! Pore filling (output) to compute the dynamic ice table
    … …  
    555556allocate(delta_co2_adsorbded(ngrid))
    556557allocate(co2ice_disappeared(ngrid,nslope))
    557 allocate(icetable_thickness_prev_iter(ngrid,nslope))
     558allocate(icetable_thickness_old(ngrid,nslope))
     559allocate(ice_porefilling_old(ngrid,nsoilmx_PEM,nslope))
    558560allocate(delta_h2o_icetablesublim(ngrid))
    559561allocate(delta_h2o_adsorbded(ngrid))
    … …  
    927929        if (icetable_equilibrium) then
    928930            write(*,*) "Compute ice table (equilibrium method)"
    929             icetable_thickness_prev_iter = icetable_thickness
     931            icetable_thickness_old = icetable_thickness
    930932            call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
    931             call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
     933            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
    932934        else if (icetable_dynamic) then
    933935            write(*,*) "Compute ice table (dynamic method)"
     936            ice_porefilling_old = ice_porefilling
    934937            allocate(porefill(nsoilmx_PEM))
    935938            do ig = 1,ngrid
    … …  
    941944            enddo
    942945            deallocate(porefill)
    943             !call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
     946            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
    944947        endif
    945948
    946949! II_d.4 Update the soil thermal properties
    947         call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,TI_PEM)
     950        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
    948951
    949952! II_d.5 Update the mass of the regolith adsorbed
    … …  
    12301233deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved)
    12311234deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
    1232 deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,icetable_thickness_prev_iter)
     1235deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim)
     1236deallocate(icetable_thickness_old,ice_porefilling_old)
    12331237deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
    12341238!----------------------------- END OUTPUT ------------------------------

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

    r3512 r3525  
    77!=======================================================================
    88
    9 SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table_depth,ice_table_thickness,      &
    10                     ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,d_h2oice,d_co2ice,co2_ice,h2o_ice, &
    11                     global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,                         &
     9SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table_depth,ice_table_thickness, &
     10                    ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,d_h2oice,d_co2ice,co2_ice,h2o_ice,    &
     11                    global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,                    &
    1212                    m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif)
    1313
    … …  
    1616use comsoil_h,                  only: volcapa, inertiedat
    1717use adsorption_mod,             only: regolith_adsorption, adsorption_pem
    18 use ice_table_mod,              only: computeice_table_equilibrium, icetable_equilibrium, icetable_dynamic
     18use ice_table_mod,              only: computeice_table_equilibrium, icetable_depth, icetable_thickness, icetable_equilibrium, icetable_dynamic
    1919use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
    2020use soil_thermalproperties_mod, only: update_soil_thermalproperties
    … …  
    315315                write(*,*)'will reconstruct the values of the ice table given the current state'
    316316                call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
    317                 call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
     317                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
    318318                do islope = 1,nslope
    319319                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
    … …  
    327327                write(*,*)'will reconstruct the values of the ice table given the current state'
    328328                ice_table_depth = -9999.
    329                 call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
     329                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
    330330                do islope = 1,nslope
    331331                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
    … …  
    504504        if (icetable_equilibrium) then
    505505            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
    506             call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
     506            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
    507507            do islope = 1,nslope
    508508                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
    … …  
    512512            ice_porefilling = 1.
    513513            ice_table_depth = -9999.
    514             call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
     514            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
    515515            do islope = 1,nslope
    516516                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))

  • 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.