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

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/y] (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/y] (min|max): '//real2str( minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO)

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

!=======================================================================
SUBROUTINE evolve_flux_ssice(h2oice_depth_old,h2oice_depth_new,tsurf,tsoil_ts_old,tsoil_ts_new,flux_ssice_avg)
!-----------------------------------------------------------------------
! NAME
!     evolve_flux_ssice
!
! 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 geometry,       only: ngrid, nslope
use soil_temp,      only: itp_tsoil
use subsurface_ice, only: psv
use ice_table,      only: coef_ssdif_PCM
use display,        only: print_msg, LVL_NFO

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:),     intent(in) :: h2oice_depth_old ! Old H2O ice depth
real(dp), dimension(:,:),     intent(in) :: h2oice_depth_new ! New H2O ice depth
real(dp), dimension(:,:),     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), dimension(:,:),     intent(inout) :: flux_ssice_avg       ! Tendency of sub-surface ice

! LOCAL VARIABLES
! ---------------
real(dp), parameter :: zcdv = 0.0325_dp ! Drag coefficient
real(dp)            :: Rz_old, Rz_new, R_dec, hum_dec, psv_max_old, psv_max_new
integer(di)         :: i, islope, iday

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

do i = 1,ngrid
    do islope = 1,nslope
        ! Higher resistance due to growing lag layer (higher depth)
        Rz_old = h2oice_depth_old(i,islope)*zcdv/coef_ssdif_PCM(i,islope) ! Old resistance from PCM
        Rz_new = h2oice_depth_new(i,islope)*zcdv/coef_ssdif_PCM(i,islope) ! 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(i,:,islope,iday),tsurf(i,islope),h2oice_depth_old(i,islope))))
            psv_max_new = max(psv_max_new,psv(itp_tsoil(tsoil_ts_new(i,:,islope,iday),tsurf(i,islope),h2oice_depth_new(i,islope))))
        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)
        flux_ssice_avg(i,islope) = flux_ssice_avg(i,islope)*R_dec*hum_dec
    enddo
enddo

END SUBROUTINE evolve_flux_ssice
!=======================================================================

END MODULE tendencies