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

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

! DEPENDENCIES
! ------------
use numerics, only: dp, di, k4, minieps, tol

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

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

!=======================================================================
SUBROUTINE compute_tendice(min_perice,min_frost,perice,d_ice)
!-----------------------------------------------------------------------
! NAME
!     compute_tendice
!
! 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
! ---------
real(dp), dimension(:,:,:), intent(in)  :: min_perice
real(dp), dimension(:,:,:), intent(in)  :: min_frost
real(dp), dimension(:,:),   intent(in)  :: perice
real(dp), dimension(:,:),   intent(out) :: d_ice

! CODE
! ----
! We compute the difference to get the tendency
d_ice(:,:) = min_perice(:,:,2) + min_frost(:,:,2) - (min_perice(:,:,1) + min_frost(:,:,1))

! If the difference is too small, then there is no evolution
where (abs(d_ice) < minieps) d_ice = 0._dp

! If the tendency is negative but there is no ice reservoir for the PEM
where (d_ice < 0._dp .and. abs(perice) < minieps) d_ice = 0._dp

END SUBROUTINE compute_tendice
!=======================================================================

!=======================================================================
SUBROUTINE evolve_tend_co2ice(d_co2ice_ini,co2ice,emissivity,q_co2_ts_ini,q_co2_ts,ps_PCM,ps_avg_glob_ini,ps_avg_global,d_co2ice)
!-----------------------------------------------------------------------
! NAME
!     evolve_tend_co2ice
!
! 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 geometry,   only: ngrid, nslope, nday
use physics,    only: sigmaB, alpha_clap_co2, beta_clap_co2, Lco2, m_co2, A, B
use orbit,      only: yr_len_sol, sol_len_s
use display,    only: print_msg, LVL_NFO
use utility,    only: real2str

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(in)    :: q_co2_ts, q_co2_ts_ini, ps_PCM, d_co2ice_ini, co2ice, emissivity
real(dp),                 intent(in)    :: ps_avg_global, ps_avg_glob_ini
real(dp), dimension(:,:), intent(inout) :: d_co2ice

! LOCAL VARIABLES
! ---------------
integer(di)                     :: i, iday, islope
real(dp)                        :: coef, avg
real(dp), dimension(ngrid,nday) :: vmr_co2_ts, vmr_co2_ts_ini

! CODE
! ----
call print_msg("> Evolving the CO2 ice tendency according to the new pressure",LVL_NFO)
call print_msg('Old CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO)

! Compute volume mixing ratio
vmr_co2_ts(:,:) = q_co2_ts(:,:)/(A*q_co2_ts(:,:) + B)/m_co2
vmr_co2_ts_ini(:,:) = q_co2_ts_ini(:,:)/(A*q_co2_ts_ini(:,:) + B)/m_co2

! Recompute the tendency
do i = 1,ngrid
    do islope = 1,nslope
        coef = yr_len_sol*sol_len_s*emissivity(i,islope)*sigmaB/Lco2
        avg = 0._dp
        if (co2ice(i,islope) > 1.e-4_dp .and. abs(d_co2ice(i,islope)) > 1.e-5_dp) then
            do iday = 1,nday
                avg = avg + (beta_clap_co2/(alpha_clap_co2 - log(q_co2_ts_ini(i,iday)*ps_PCM(i,iday)/100._dp)))**4 &
                      - (beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_ts(i,iday)*ps_PCM(i,iday)*(ps_avg_global/ps_avg_glob_ini)/100._dp)))**4
            end do
            if (avg < 1.e-4_dp) avg = 0._dp
            d_co2ice(i,islope) = d_co2ice_ini(i,islope) - coef*avg/nday
        end if
    end do
end do
call print_msg('New CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str( minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO)

END SUBROUTINE evolve_tend_co2ice
!=======================================================================

!=======================================================================
SUBROUTINE evolve_tend_h2oice(h2oice_depth_old,h2oice_depth_new,tsurf,tsoil_ts_old,tsoil_ts_new,d_h2oice)
!-----------------------------------------------------------------------
! NAME
!     evolve_tend_h2oice
!
! 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
use display,        only: print_msg, LVL_NFO

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

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

! LOCAL VARIABLES
! ---------------
real(dp)            :: Rz_old, Rz_new, R_dec, hum_dec, psv_max_old, psv_max_new
integer(di)         :: iday
real(dp), parameter :: coef_diff = 4.e-4_dp ! Diffusion coefficient
real(dp), parameter :: zcdv = 0.0325_dp     ! Drag coefficient

! CODE
! ----
call print_msg("> Updating the H2O sub-surface ice tendency due to lag layer",LVL_NFO)

! 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._dp
psv_max_new = 0._dp
do iday = 1,size(tsoil_ts_old,2)
    psv_max_old = max(psv_max_old,psv(itp_tsoil(tsoil_ts_old(:,iday),tsurf,h2oice_depth_old)))
    psv_max_new = max(psv_max_new,psv(itp_tsoil(tsoil_ts_new(:,iday),tsurf,h2oice_depth_new)))
end do

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

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

END SUBROUTINE evolve_tend_h2oice
!=======================================================================

END MODULE tendencies