source: trunk/LMDZ.COMMON/libf/evolution/tendencies.F90 @ 3991

MODULE tendencies
!-----------------------------------------------------------------------
! NAME
!     tendencies
!
! DESCRIPTION
!     Computation and update of PEM ice evolution tendencies.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!=======================================================================
SUBROUTINE compute_tend(ngrid,nslope,min_ice,d_ice)
!-----------------------------------------------------------------------
! NAME
!     compute_tend
!
! DESCRIPTION
!     Compute initial ice evolution tendencies from PCM data.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!     Based on minima of ice at each point for the PCM years.
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
integer,                         intent(in)  :: ngrid
integer,                         intent(in)  :: nslope
real, dimension(ngrid,nslope,2), intent(in)  :: min_ice ! Minima of ice at each point for the PCM years
real, dimension(ngrid,nslope),   intent(out) :: d_ice   ! Evolution of perennial ice

! CODE
! ----
! We compute the difference
d_ice = min_ice(:,:,2) - min_ice(:,:,1)

! If the difference is too small, then there is no evolution
where (abs(d_ice) < 1.e-10) d_ice = 0.

! If the minimum over the last year is 0, then we have no perennial ice
where (abs(min_ice(:,:,2)) < 1.e-10) d_ice = 0.

END SUBROUTINE compute_tend
!=======================================================================

!=======================================================================
SUBROUTINE recomp_tend_co2(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_ini,co2ice,emissivity, &
                           vmr_co2_PCM,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global)
!-----------------------------------------------------------------------
! NAME
!     recomp_tend_co2
!
! DESCRIPTION
!     Recompute CO2 ice tendency based on pressure and atmospheric changes.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!     Adjusts CO2 ice evolution based on Clausius-Clapeyron changes.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB, Lco2, sols_per_my, sec_per_sol

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
integer,                        intent(in)    :: timelen, ngrid, nslope ! Time length, # of grid points and slopes
real, dimension(ngrid,timelen), intent(in)    :: vmr_co2_PCM            ! CO2 VMR in PCM first layer
real, dimension(ngrid,timelen), intent(in)    :: vmr_co2_PEM            ! CO2 VMR in PEM first layer
real, dimension(ngrid,timelen), intent(in)    :: ps_PCM                 ! Surface pressure in PCM
real,                           intent(in)    :: ps_avg_global_ini      ! Global average pressure (initial)
real,                           intent(in)    :: ps_avg_global          ! Global average pressure (current)
real, dimension(ngrid,nslope),  intent(in)    :: d_co2ice_ini           ! Initial CO2 ice evolution
real, dimension(ngrid,nslope),  intent(in)    :: co2ice                 ! CO2 ice surface [kg/m^2]
real, dimension(ngrid,nslope),  intent(in)    :: emissivity             ! Emissivity
real, dimension(ngrid,nslope),  intent(inout) :: d_co2ice_phys          ! Updated CO2 ice evolution

! LOCAL VARIABLES
! ---------------
integer :: i, t, islope
real    :: coef, avg

! CODE
! ----
write(*,*) "> Updating the CO2 ice tendency for the new pressure"

! Evolution of the water ice for each physical point
do i = 1,ngrid
    do islope = 1,nslope
        coef = sols_per_my*sec_per_sol*emissivity(i,islope)*sigmaB/Lco2
        avg = 0.
        if (co2ice(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then
            do t = 1,timelen
                avg = avg + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM(i,t)/100.)))**4 &
                      - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(i,t)*ps_PCM(i,t)*(ps_avg_global/ps_avg_global_ini)/100.)))**4
            enddo
            if (avg < 1.e-4) avg = 0.
            d_co2ice_phys(i,islope) = d_co2ice_ini(i,islope) - coef*avg/timelen
        endif
    enddo
enddo

END SUBROUTINE recomp_tend_co2
!=======================================================================

!=======================================================================
SUBROUTINE recomp_tend_h2o(h2oice_depth_old,h2oice_depth_new,tsurf,tsoil_PEM_timeseries_old,tsoil_PEM_timeseries_new,d_h2oice)
!-----------------------------------------------------------------------
! NAME
!     recomp_tend_h2o
!
! DESCRIPTION
!     Recompute H2O ice tendency based on soil depth and temperature changes.
!
! AUTHORS & DATE
!     JB Clement, 2025 (following E. Vos's work)
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use soil_temp,      only: itp_tsoil
use subsurface_ice, only: psv

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real,                 intent(in)    :: h2oice_depth_old ! Old H2O ice depth
real,                 intent(in)    :: h2oice_depth_new ! New H2O ice depth
real,                 intent(in)    :: tsurf            ! Surface temperature
real, dimension(:,:), intent(in)    :: tsoil_PEM_timeseries_old ! Old soil temperature time series
real, dimension(:,:), intent(in)    :: tsoil_PEM_timeseries_new ! New soil temperature time series
real,                 intent(inout) :: d_h2oice         ! Evolution of perennial ice

! LOCAL VARIABLES
! ---------------
real            :: Rz_old, Rz_new, R_dec, hum_dec, psv_max_old, psv_max_new
integer         :: t
real, parameter :: coef_diff = 4.e-4 ! Diffusion coefficient
real, parameter :: zcdv = 0.0325     ! Drag coefficient

! CODE
! ----
! Higher resistance due to growing lag layer (higher depth)
Rz_old = h2oice_depth_old*zcdv/coef_diff ! Old resistance from PCM
Rz_new = h2oice_depth_new*zcdv/coef_diff ! New resistance based on new depth
R_dec = Rz_old/Rz_new ! Decrease because of resistance

! The maxmimum of the daily averages over one year for the saturation vapor pressure at the ice table location
psv_max_old = 0.
psv_max_new = 0.
do t = 1,size(tsoil_PEM_timeseries_old,2)
    psv_max_old = max(psv_max_old,psv(itp_tsoil(tsoil_PEM_timeseries_old(:,t),tsurf,h2oice_depth_old)))
    psv_max_new = max(psv_max_new,psv(itp_tsoil(tsoil_PEM_timeseries_new(:,t),tsurf,h2oice_depth_new)))
enddo

! Lower humidity due to growing lag layer (higher depth)
if (abs(psv_max_old) < 1.e2*epsilon(1.)) then
    hum_dec = 1.
else
    hum_dec = psv_max_new/psv_max_old ! Decrease because of lower water vapor pressure at the new depth
endif

! Flux correction (decrease)
d_h2oice = d_h2oice*R_dec*hum_dec

END SUBROUTINE recomp_tend_h2o
!=======================================================================

END MODULE tendencies