Changeset 4071


Ignore:
Timestamp:
Feb 16, 2026, 10:28:56 AM (12 days ago)
Author:
jbclement
Message:

PEM:

  • Making the computation of ice tendencies more reliable by doing it after 'read_startpem' to know exactly where perennial ice is (no matter if there is a "startpem.nc" or not). Moreover, when there is no "startpem.nc", location of perennial ice depends now on 'watercaptag' and on huge amount of frost. This prevents negative ice tendency while there is no ice which can happen with weird PCM inputs (i.e. 'watercaptag' = F & 'watercap' /= 0 & no "stratpem.nc").
  • Few small safeguards throughout the code.

JBC

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

Legend:

Unmodified
Added
Removed

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

    r4068 r4071  
    868868== 12/02/2026 == JBC
    869869Following r4065, some safeguards forgotten in "io_netcdf.F90" + 'tsurf' is not needed anymore in the "start1D.txt" + few forgotten small updates.
     870
     871== 16/02/2026 == JBC
     872- Making the computation of ice tendencies more reliable by doing it after 'read_startpem' to know exactly where perennial ice is (no matter if there is a "startpem.nc" or not). Moreover, when there is no "startpem.nc", location of perennial ice depends now on 'watercaptag' and on huge amount of frost. This prevents negative ice tendency while there is no ice which can happen with weird PCM inputs (i.e. 'watercaptag' = F & 'watercap' /= 0 & no "stratpem.nc").
     873- Few small safeguards throughout the code.

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

    r4068 r4071  
    121121! DEPENDENCIES
    122122! ------------
    123 use geometry,   only: nlayer, nslope
     123use geometry,   only: nlayer
    124124use atmosphere, only: set_ap, set_bp, set_ps_PCM, set_teta_PCM, set_u_PCM, set_v_PCM, compute_alt_coord, set_ap, set_bp
    125125use tracers,    only: nq, set_q_PCM
    … …  
    134134! LOCAL VARIABLES
    135135! ---------------
    136 integer(di)                     :: i, j, k, ierr, funit
     136integer(di)                     :: i, k, ierr, funit
    137137character(30)                   :: header
    138138real(dp), dimension(1,1)        :: tmp
    … …  
    326326use planet,             only: TI_breccia, TI_bedrock, alpha_clap_h2o, beta_clap_h2o
    327327use frost,              only: h2ofrost_PCM, h2o_frost4PCM, co2frost_PCM, co2_frost4PCM
    328 use surf_ice,           only: is_h2o_perice_PCM, h2oice_huge_ini, co2_perice_PCM
     328use surf_ice,           only: is_h2o_perice_PCM, h2oice_huge_ini, co2_perice_PCM, threshold_h2oice_cap
    329329use soil,               only: do_soil, TI, depth_breccia, depth_bedrock, index_breccia, index_bedrock, inertiedat, layer, inertiedat_PCM
    330330use soil_temp,          only: ini_tsoil_pem, compute_tsoil
    … …  
    335335use layered_deposits,   only: layering, do_layering, array2map, ini_layering, add_stratum
    336336use surf_ice,           only: rho_co2ice, rho_h2oice
     337use slopes,             only: subslope_dist, def_slope_mean
     338use maths,              only: pi
    337339use display,            only: print_msg
    338340use utility,            only: int2str
    … …  
    344346! ARGUMENTS
    345347! ---------
    346 real(dp), dimension(:,:),       intent(in)    :: tsurf_avg_yr1, tsurf_avg_yr2 ! Surface temperature at the last but one PCM run and at the last PCM run
    347 real(dp),                       intent(in)    :: ps_avg_glob                  ! Mean average pressure on the planet
    348 real(dp), dimension(:,:),       intent(in)    :: ps_ts                        ! Surface pressure timeseries
    349 real(dp), dimension(:,:),       intent(in)    :: q_co2_ts, q_h2o_ts           ! MMR of CO2 and H2O
    350 real(dp), dimension(:,:),       intent(in)    :: h2o_surfdensity_avg          ! Average of surface water density
    351 real(dp), dimension(:,:),       intent(in)    :: d_h2oice, d_co2ice           ! Tendencies on ice
    352 real(dp), dimension(:,:),       intent(inout) :: h2o_ice, co2_ice             ! Surface ice
    353 real(dp), dimension(:,:,:),    intent(inout) :: tsoil_avg                    ! Soil temperature
    354 real(dp), dimension(:,:,:),    intent(inout) :: h2o_soildensity_avg          ! Average of soil water soil density
    355 real(dp), dimension(:,:),       intent(inout) :: icetable_depth               ! Ice table depth
    356 real(dp), dimension(:,:),       intent(inout) :: icetable_thickness           ! Ice table thickness
    357 real(dp), dimension(:,:,:),    intent(inout) :: ice_porefilling              ! Subsurface ice pore filling
    358 real(dp), dimension(:,:,:),    intent(inout) :: h2o_ads_reg, co2_ads_reg     ! Mass of H2O and CO2 adsorbed
    359 type(layering), dimension(:,:), intent(inout) :: layerings_map                ! Layerings
    360 real(dp), dimension(:),         intent(inout) :: delta_h2o_ads, delta_co2_ads ! Mass of H2O and CO2 exchanged due to adsorption/desorption
     348real(dp),       dimension(:,:),   intent(in)    :: tsurf_avg_yr1, tsurf_avg_yr2 ! Surface temperature at the last but one PCM run and at the last PCM run
     349real(dp),                         intent(in)    :: ps_avg_glob                  ! Mean average pressure on the planet
     350real(dp),       dimension(:,:),   intent(in)    :: ps_ts                        ! Surface pressure timeseries
     351real(dp),       dimension(:,:),   intent(in)    :: q_co2_ts, q_h2o_ts           ! MMR of CO2 and H2O
     352real(dp),       dimension(:,:),   intent(in)    :: h2o_surfdensity_avg          ! Average of surface water density
     353real(dp),       dimension(:,:),   intent(in)    :: d_h2oice, d_co2ice           ! Tendencies on ice
     354real(dp),       dimension(:,:),   intent(inout) :: h2o_ice, co2_ice             ! Surface ice
     355real(dp),       dimension(:,:,:), intent(inout) :: tsoil_avg                    ! Soil temperature
     356real(dp),       dimension(:,:,:), intent(inout) :: h2o_soildensity_avg          ! Average of soil water soil density
     357real(dp),       dimension(:,:),   intent(inout) :: icetable_depth               ! Ice table depth
     358real(dp),       dimension(:,:),   intent(inout) :: icetable_thickness           ! Ice table thickness
     359real(dp),       dimension(:,:,:), intent(inout) :: ice_porefilling              ! Subsurface ice pore filling
     360real(dp),       dimension(:,:,:), intent(inout) :: h2o_ads_reg, co2_ads_reg     ! Mass of H2O and CO2 adsorbed
     361type(layering), dimension(:,:),   intent(inout) :: layerings_map                ! Layerings
     362real(dp),       dimension(:),     intent(inout) :: delta_h2o_ads, delta_co2_ads ! Mass of H2O and CO2 exchanged due to adsorption/desorption
    361363
    362364! LOCAL VARIABLES
    … …  
    381383
    382384    ! H2O ice
    383     call print_msg("'h2o_ice' is initialized with default value 'h2oice_huge_ini' where 'watercaptag' is true.")
     385    call print_msg("'h2o_ice' is initialized with default value 'h2oice_huge_ini' where 'watercaptag' is true and where yearly minimum of frost can be considered as a huge reservoir ('threshold_h2oice_cap').")
    384386    h2o_ice(:,:) = 0._dp
    385387    do i = 1,ngrid
    386         if (is_h2o_perice_PCM(i)) h2o_ice(i,:) = h2oice_huge_ini + h2ofrost_PCM(i,:) - h2o_frost4PCM(i,:)
     388        if (is_h2o_perice_PCM(i)) then
     389            h2o_ice(i,:) = h2oice_huge_ini
     390        else if (sum((h2ofrost_PCM(i,:) - h2o_frost4PCM(i,:))*subslope_dist(i,:)/cos(pi*def_slope_mean(:)/180._dp)) > threshold_h2oice_cap) then
     391            h2o_ice(i,:) = h2oice_huge_ini
     392        end if
    387393    end do
    388394
    389395    ! CO2 ice
    390     call print_msg("'co2_ice' is initialized with 'perennial_co2ice' found in the PCM.")
     396    call print_msg("'co2_ice' is initialized with 'perennial_co2ice' and yearly minimum of frost found in the PCM.")
    391397    co2_ice(:,:) = co2_perice_PCM(:,:) + co2frost_PCM(:,:) - co2_frost4PCM(:,:)
    392398

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

    r4068 r4071  
    137137! DEPENDENCIES
    138138! ------------
    139 use geometry,   only: nlayer, nslope
     139use geometry,   only: nlayer
    140140use atmosphere, only: teta_PCM, u_PCM, v_PCM
    141141use tracers,    only: nq, qnames
    … …  
    155155! LOCAL VARIABLES
    156156! ---------------
    157 integer(di) :: funit, ierr, i, j, l
     157integer(di) :: funit, ierr, i, l
    158158
    159159! CODE
    … …  
    244244
    245245! Variables that have been modified
    246 call put_var_nc('watercaptag',merge(1.,0.,is_h2o_perice))
     246call put_var_nc('watercaptag',merge(1._dp,0._dp,is_h2o_perice))
    247247call put_var_nc('watercap',h2o_ice4PCM,1)
    248248call put_var_nc('h2o_ice',h2o_frost4PCM,1)

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

    r4065 r4071  
    164164
    165165!=======================================================================
    166 SUBROUTINE compute_frost4PCM(min_h2ofrost,min_co2frost)
     166SUBROUTINE compute_frost4PCM(minPCM_h2ofrost,minPCM_co2frost)
    167167!-----------------------------------------------------------------------
    168168! NAME
    … …  
    190190! ARGUMENTS
    191191! ---------
    192 real(dp), dimension(:,:), intent(in) :: min_h2ofrost, min_co2frost
     192real(dp), dimension(:,:), intent(in) :: minPCM_h2ofrost, minPCM_co2frost
    193193
    194194! CODE
    … …  
    198198h2o_frost4PCM(:,:) = 0._dp
    199199co2_frost4PCM(:,:) = 0._dp
    200 where (h2ofrost_PCM(:,:) > 0._dp) h2o_frost4PCM(:,:) = h2ofrost_PCM(:,:) - min_h2ofrost(:,:)
    201 where (co2frost_PCM(:,:) > 0._dp) co2_frost4PCM(:,:) = co2frost_PCM(:,:) - min_co2frost(:,:)
     200where (h2ofrost_PCM(:,:) > 0._dp) h2o_frost4PCM(:,:) = h2ofrost_PCM(:,:) - minPCM_h2ofrost(:,:)
     201where (co2frost_PCM(:,:) > 0._dp) co2_frost4PCM(:,:) = co2frost_PCM(:,:) - minPCM_co2frost(:,:)
    202202
    203203END SUBROUTINE compute_frost4PCM

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

    r4065 r4071  
    1818! DEPENDENCIES
    1919! ------------
    20 use numerics, only: dp, di, k4
     20use numerics, only: dp, di, k4, minieps
    2121
    2222! DECLARATION
    … …  
    108108! ARGUMENTS
    109109! ---------
    110 real(dp),    dimension(:,:),            intent(in)    :: vmr_co2_PEM     ! Grid points x Time field : VMR of CO2 in the first layer [mol/mol]
    111 real(dp),    dimension(:,:),            intent(in)    :: ps_PCM          ! Grid points x Time field: surface pressure given by the PCM [Pa]
    112 real(dp),                               intent(in)    :: ps_avg_glob_ini ! Global averaged surface pressure at the beginning [Pa]
    113 real(dp),                               intent(in)    :: ps_avg_global   ! Global averaged surface pressure during the PEM iteration [Pa]
    114 real(dp),    dimension(:,:),            intent(inout) :: co2ice          ! Grid points x Slope field: CO2 ice on the subgrid slopes [kg/m^2]
    115 logical(k4), dimension(:,:),            intent(out)   :: is_co2ice_flow  ! Flag to see if there is flow on the subgrid slopes
     110real(dp),    dimension(:,:), intent(in)    :: vmr_co2_PEM     ! Grid points x Time field : VMR of CO2 in the first layer [mol/mol]
     111real(dp),    dimension(:,:), intent(in)    :: ps_PCM          ! Grid points x Time field: surface pressure given by the PCM [Pa]
     112real(dp),                    intent(in)    :: ps_avg_glob_ini ! Global averaged surface pressure at the beginning [Pa]
     113real(dp),                    intent(in)    :: ps_avg_global   ! Global averaged surface pressure during the PEM iteration [Pa]
     114real(dp),    dimension(:,:), intent(inout) :: co2ice          ! Grid points x Slope field: CO2 ice on the subgrid slopes [kg/m^2]
     115logical(k4), dimension(:,:), intent(out)   :: is_co2ice_flow  ! Flag to see if there is flow on the subgrid slopes
    116116
    117117! LOCAL VARIABLES
    … …  
    301301                    iaval = islope + 1
    302302                end if
    303                 do while (iaval /= iflat .and. subslope_dist(ig,iaval) == 0)
     303                do while (iaval /= iflat .and. abs(subslope_dist(ig,iaval)) < minieps)
    304304                    if (iaval > iflat) then
    305305                        iaval = iaval - 1

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

    r4065 r4071  
    113113! ARGUMENTS
    114114! ---------
    115 integer(di), intent(in) :: stopPEM    ! Reason to stop
    116 real(dp),    intent(in) :: n_yr_run ! # of years
    117 real(dp),    intent(in) :: n_yr_sim ! Current simulated year
     115integer(di), intent(in) :: stopPEM     ! Reason to stop
     116real(dp),    intent(in) :: n_yr_run    ! # of years
     117real(dp),    intent(in) :: n_yr_sim    ! Current simulated year
    118118real(dp),    intent(in) :: nmax_yr_sim ! Maximum number of years to be simulated
    119119

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

    r4068 r4071  
    1616! DEPEDENCIES
    1717! -----------
    18 use numerics, only: dp, di, k4
     18use numerics, only: dp, di, k4, minieps
    1919use netcdf,   only: nf90_double, nf90_noerr, nf90_strerror, nf90_write, nf90_nowrite,                &
    2020                    nf90_open, nf90_close, nf90_redef, nf90_enddef, nf90_inquire, nf90_max_var_dims, &
    … …  
    959959if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
    960960if (has_fill) then
    961     if (var == fill_value) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
     961    if (abs(var - fill_value) < minieps) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
    962962end if
    963963
    … …  
    10441044if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
    10451045if (has_fill) then
    1046     if (any(var == fill_value)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
     1046    if (any(abs(var - fill_value) < minieps)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
    10471047end if
    10481048
    … …  
    11291129if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
    11301130if (has_fill) then
    1131     if (any(var == fill_value)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
     1131    if (any(abs(var - fill_value) < minieps)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
    11321132end if
    11331133
    … …  
    12141214if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
    12151215if (has_fill) then
    1216     if (any(var == fill_value)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
     1216    if (any(abs(var - fill_value) < minieps)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
    12171217end if
    12181218
    … …  
    12991299if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
    13001300if (has_fill) then
    1301     if (any(var == fill_value)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
     1301    if (any(abs(var - fill_value) < minieps)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
    13021302end if
    13031303

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

    r4065 r4071  
    406406! DEPENDENCIES
    407407! ------------
    408 use maths,     only: pi
    409408use evolution, only: pem_ini_date
    410409use display,   only: print_msg
    … …  
    572571
    573572use evolution,        only: pem_ini_date
    574 use maths,            only: pi
    575573use call_dayperi_mod, only: call_dayperi
    576574use stoppage,         only: stop_clean

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

    r4068 r4071  
    4848use surf_ice,           only: evolve_co2ice, evolve_h2oice, balance_h2oice_reservoirs, build4PCM_perice
    4949use surf_temp,          only: tsurf_PCM, adapt_tsurf2disappearedice, build4PCM_tsurf
    50 use tendencies,         only: evolve_tend_co2, evolve_tend_h2o
     50use tendencies,         only: compute_tendice, evolve_tend_co2ice, evolve_tend_h2oice
    5151use tracers,            only: adapt_tracers2pressure, build4PCM_tracers, nq
    5252use utility,            only: real2str
    … …  
    7777real(dp)                              :: preff4PCM       ! Reference pressure [Pa]
    7878! Ice-related:
    79 real(dp),    dimension(:,:), allocatable :: h2o_ice                    ! H2O ice [kg.m-2]
    80 real(dp),    dimension(:,:), allocatable :: co2_ice                    ! CO2 ice [kg.m-2]
    81 real(dp)                                 :: h2oice_sublim_coverage_ini ! Initial surface area of sublimating H2O ice [m2]
    82 real(dp)                                 :: co2ice_sublim_coverage_ini ! Initial surface area of sublimating CO2 ice [m2]
    83 logical(k4), dimension(:,:), allocatable :: is_h2oice_ini              ! Initial location of H2O ice
    84 logical(k4), dimension(:,:), allocatable :: is_co2ice_ini              ! Initial location of CO2 ice
    85 logical(k4), dimension(:,:), allocatable :: is_co2ice_flow             ! Flag for location of CO2 glacier flow
    86 logical(k4), dimension(:,:), allocatable :: is_h2oice_flow             ! Flag for location of H2O glacier flow
    87 real(dp),    dimension(:,:), allocatable :: h2o_ice4PCM                ! H2O ice reconstruction to feed back into PCM [kg.m-2]
    88 real(dp),    dimension(:,:), allocatable :: co2_ice4PCM                ! CO2 ice reconstruction to feed back into PCM [kg.m-2]
    89 logical(k4), dimension(:),   allocatable :: is_h2o_perice              ! Location of H2O infinite reservoir, called 'watercaptag' in PCM
    90 logical(k4), dimension(:,:), allocatable :: is_co2ice_disappeared      ! Flag to check if CO2 ice disappeared at the previous timestep
     79real(dp),    dimension(:,:),   allocatable :: h2o_ice                    ! H2O ice [kg.m-2]
     80real(dp),    dimension(:,:),   allocatable :: co2_ice                    ! CO2 ice [kg.m-2]
     81real(dp)                                   :: h2oice_sublim_coverage_ini ! Initial surface area of sublimating H2O ice [m2]
     82real(dp)                                   :: co2ice_sublim_coverage_ini ! Initial surface area of sublimating CO2 ice [m2]
     83logical(k4), dimension(:,:),   allocatable :: is_h2oice_ini              ! Initial location of H2O ice
     84logical(k4), dimension(:,:),   allocatable :: is_co2ice_ini              ! Initial location of CO2 ice
     85logical(k4), dimension(:,:),   allocatable :: is_co2ice_flow             ! Flag for location of CO2 glacier flow
     86logical(k4), dimension(:,:),   allocatable :: is_h2oice_flow             ! Flag for location of H2O glacier flow
     87real(dp),    dimension(:,:,:), allocatable :: minPCM_h2operice              ! Minimum of H2O perennial ice over the last PCM year [kg.m-2]
     88real(dp),    dimension(:,:,:), allocatable :: minPCM_co2perice              ! Minimum of CO2 perennial ice over the last PCM year [kg.m-2]
     89real(dp),    dimension(:,:,:), allocatable :: minPCM_h2ofrost               ! Minimum of H2O frost over the last PCM year [kg.m-2]
     90real(dp),    dimension(:,:,:), allocatable :: minPCM_co2frost               ! Minimum of CO2 frost over the last PCM year [kg.m-2]
     91real(dp),    dimension(:,:),   allocatable :: h2o_ice4PCM                ! H2O ice reconstruction to feed back into PCM [kg.m-2]
     92real(dp),    dimension(:,:),   allocatable :: co2_ice4PCM                ! CO2 ice reconstruction to feed back into PCM [kg.m-2]
     93logical(k4), dimension(:),     allocatable :: is_h2o_perice              ! Location of H2O infinite reservoir, called 'watercaptag' in PCM
     94logical(k4), dimension(:,:),   allocatable :: is_co2ice_disappeared      ! Flag to check if CO2 ice disappeared at the previous timestep
    9195! Surface-related:
    9296real(dp), dimension(:,:), allocatable :: tsurf_avg           ! Average surface temperature [K]
    … …  
    179183allocate(tsoil_avg(ngrid,nsoil,nslope),tsoil_ts(ngrid,nsoil,nslope,nday),h2o_soildensity_avg(ngrid,nsoil,nslope))
    180184allocate(q_h2o_ts(ngrid,nday),q_co2_ts(ngrid,nday))
    181 allocate(d_h2oice(ngrid,nslope),d_co2ice(ngrid,nslope))
     185allocate(minPCM_h2operice(ngrid,nslope,2),minPCM_co2perice(ngrid,nslope,2),minPCM_h2ofrost(ngrid,nslope,2),minPCM_co2frost(ngrid,nslope,2))
    182186
    183187call load_xios_data(ps_avg,ps_ts,tsurf_avg,tsurf_avg_yr1,tsoil_avg,tsoil_ts,h2o_surfdensity_avg,h2o_soildensity_avg, &
    184                     q_h2o_ts,q_co2_ts,d_h2oice,d_co2ice)
     188                    q_h2o_ts,q_co2_ts,minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost)
    185189
    186190! Initiate soil settings and TI
    … …  
    207211                   h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads)
    208212deallocate(tsurf_avg_yr1)
     213
     214! Compute ice tendencies from yearly minima
     215allocate(d_h2oice(ngrid,nslope),d_co2ice(ngrid,nslope))
     216call print_msg('> Computing surface ice tendencies')
     217call compute_tendice(minPCM_h2operice + minPCM_h2ofrost,h2o_ice(:,:) > 0._dp,d_h2oice)
     218call print_msg('H2O ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_h2oice))//' | '//real2str(maxval(d_h2oice)))
     219call compute_tendice(minPCM_co2perice + minPCM_co2frost,co2_ice(:,:) > 0._dp,d_co2ice)
     220call print_msg('CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)))
     221deallocate(minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost)
    209222
    210223! Save initial set-up useful for the next computations
    … …  
    290303        ! Conversion to surface ice
    291304        call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth)
     305        ! Balance H2O ice reservoirs
     306        allocate(d_h2oice_new(ngrid,nslope))
     307        call stopping_crit_h2o(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopcrit)
     308        call balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
     309        deallocate(d_h2oice_new)
    292310    else
    293311        zlag(:,:) = 0._dp
    294312        call evolve_h2oice(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,zshift_surf,stopcrit)
    295313        call evolve_co2ice(co2_ice,d_co2ice,zshift_surf)
    296     end if
    297 
    298     if (do_layering) then
    299         allocate(d_h2oice_new(ngrid,nslope))
    300         call stopping_crit_h2o(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopcrit)
    301         call balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
    302         deallocate(d_h2oice_new)
    303314    end if
    304315
    … …  
    406417
    407418    ! Evolve the tendencies
    408     call evolve_tend_co2(d_co2ice_ini,co2_ice,emissivity_PCM,q_co2_ts_ini,q_co2_ts,ps_ts,ps_avg_glob_ini,ps_avg_glob,d_co2ice)
    409 !~     call print_msg("> Updating the H2O sub-surface ice tendency due to lag layer")
     419    call evolve_tend_co2ice(d_co2ice_ini,co2_ice,emissivity_PCM,q_co2_ts_ini,q_co2_ts,ps_ts,ps_avg_glob_ini,ps_avg_glob,d_co2ice)
    410420!~     if (do_layering) then
    411421!~         do i = 1,ngrid
    412422!~             do islope = 1,nslope
    413 !~                 if (is_h2oice_sublim_ini(i,islope) .and. h2oice_depth(i,islope) > 0._dp) call evolve_tend_h2o(h2oice_depth_old(i,islope),h2oice_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope))
     423!~                 if (is_h2oice_sublim_ini(i,islope) .and. h2oice_depth(i,islope) > 0._dp) call evolve_tend_h2oice(h2oice_depth_old(i,islope),h2oice_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope))
    414424!~             end do
    415425!~         end do
    … …  
    417427!        do i = 1,ngrid
    418428!            do islope = 1,nslope
    419 !                call evolve_tend_h2o(icetable_depth_old(i,islope),icetable_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope))
     429!                call evolve_tend_h2oice(icetable_depth_old(i,islope),icetable_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope))
    420430!            end do
    421431!        end do

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

    r4065 r4071  
    200200
    201201!=======================================================================
    202 SUBROUTINE build4PCM_perice(h2o_ice,co2_ice,is_h2o_perice,h2oice_PCM,co2ice_PCM)
     202SUBROUTINE build4PCM_perice(h2o_ice,co2_ice,is_h2o_perice,h2o_ice4PCM,co2_ice4PCM)
    203203!-----------------------------------------------------------------------
    204204! NAME
    … …  
    230230! ---------
    231231real(dp),    dimension(:,:), intent(inout) :: h2o_ice, co2_ice
    232 logical(k4), dimension(:),   intent(out)   :: is_h2o_perice          ! H2O perennial ice flag
    233 real(dp),    dimension(:,:), intent(out)   :: h2oice_PCM, co2ice_PCM ! Ice for PCM
     232logical(k4), dimension(:),   intent(out)   :: is_h2o_perice            ! H2O perennial ice flag
     233real(dp),    dimension(:,:), intent(out)   :: h2o_ice4PCM, co2_ice4PCM ! Ice for PCM
    234234
    235235! LOCAL VARIABLES
    … …  
    240240! ----
    241241call print_msg('> Building surface ice/frost for the PCM')
    242 co2ice_PCM(:,:) = co2_ice(:,:)
    243 h2oice_PCM(:,:) = 0._dp ! Because in the Mars PCM, only the variation of perennial H2O ice is monitored, not the absolute quantity
     242co2_ice4PCM(:,:) = co2_ice(:,:)
     243h2o_ice4PCM(:,:) = 0._dp ! Because in the Mars PCM, only the variation of perennial H2O ice is monitored, not the absolute quantity
    244244do i = 1,ngrid
    245245    ! Is H2O ice still considered as an infinite reservoir for the PCM?

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

    r4065 r4071  
    1818! DEPENDENCIES
    1919! ------------
    20 use numerics, only: dp, di, minieps, tol
     20use numerics, only: dp, di, k4, minieps, tol
    2121
    2222! DECLARATION
    … …  
    2828
    2929!=======================================================================
    30 SUBROUTINE compute_tend(min_ice,d_ice)
    31 !-----------------------------------------------------------------------
    32 ! NAME
    33 !     compute_tend
     30SUBROUTINE compute_tendice(min_ice,is_perice,d_ice)
     31!-----------------------------------------------------------------------
     32! NAME
     33!     compute_tendice
    3434!
    3535! DESCRIPTION
    … …  
    5151! ARGUMENTS
    5252! ---------
    53 real(dp), dimension(:,:,:), intent(in)  :: min_ice
    54 real(dp), dimension(:,:),   intent(out) :: d_ice
     53real(dp),    dimension(:,:,:), intent(in)  :: min_ice
     54logical(k4), dimension(:,:),   intent(in)  :: is_perice
     55real(dp),    dimension(:,:),   intent(out) :: d_ice
    5556
    5657! CODE
    5758! ----
    58 ! We compute the difference
    59 d_ice(:,:) = min_ice(:,:,2) - min_ice(:,:,1)
     59! We compute the difference to get the tendency
     60d_ice = min_ice(:,:,2) - min_ice(:,:,1)
    6061
    6162! If the difference is too small, then there is no evolution
    6263where (abs(d_ice) < minieps) d_ice = 0._dp
    6364
    64 ! If the minimum over the last year is 0, then we have no perennial ice
    65 where (abs(min_ice(:,:,2)) < minieps) d_ice = 0._dp
    66 
    67 END SUBROUTINE compute_tend
    68 !=======================================================================
    69 
    70 !=======================================================================
    71 SUBROUTINE evolve_tend_co2(d_co2ice_ini,co2ice,emissivity,q_co2_ts_ini,q_co2_ts,ps_PCM,ps_avg_glob_ini,ps_avg_global,d_co2ice)
    72 !-----------------------------------------------------------------------
    73 ! NAME
    74 !     evolve_tend_co2
     65! If the tendency is negative but there is no ice reservoir for the PEM
     66where (abs(d_ice) < 0._dp .and. .not. is_perice) d_ice = 0._dp
     67
     68END SUBROUTINE compute_tendice
     69!=======================================================================
     70
     71!=======================================================================
     72SUBROUTINE 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)
     73!-----------------------------------------------------------------------
     74! NAME
     75!     evolve_tend_co2ice
    7576!
    7677! DESCRIPTION
    … …  
    135136call print_msg('New CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str( minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)))
    136137
    137 END SUBROUTINE evolve_tend_co2
    138 !=======================================================================
    139 
    140 !=======================================================================
    141 SUBROUTINE evolve_tend_h2o(h2oice_depth_old,h2oice_depth_new,tsurf,tsoil_ts_old,tsoil_ts_new,d_h2oice)
    142 !-----------------------------------------------------------------------
    143 ! NAME
    144 !     evolve_tend_h2o
     138END SUBROUTINE evolve_tend_co2ice
     139!=======================================================================
     140
     141!=======================================================================
     142SUBROUTINE evolve_tend_h2oice(h2oice_depth_old,h2oice_depth_new,tsurf,tsoil_ts_old,tsoil_ts_new,d_h2oice)
     143!-----------------------------------------------------------------------
     144! NAME
     145!     evolve_tend_h2oice
    145146!
    146147! DESCRIPTION
    … …  
    158159use soil_temp,      only: itp_tsoil
    159160use subsurface_ice, only: psv
     161use display,        only: print_msg
    160162
    161163! DECLARATION
    … …  
    181183! CODE
    182184! ----
     185call print_msg("> Updating the H2O sub-surface ice tendency due to lag layer")
     186
    183187! Higher resistance due to growing lag layer (higher depth)
    184188Rz_old = h2oice_depth_old*zcdv/coef_diff ! Old resistance from PCM
    … …  
    204208d_h2oice = d_h2oice*R_dec*hum_dec
    205209
    206 END SUBROUTINE evolve_tend_h2o
     210END SUBROUTINE evolve_tend_h2oice
    207211!=======================================================================
    208212

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

    r4065 r4071  
    2727!=======================================================================
    2828SUBROUTINE load_xios_data(ps_avg,ps_ts,tsurf_avg,tsurf_avg_yr1,tsoil_avg,tsoil_ts,h2o_surfdensity_avg,h2o_soildensity_avg, &
    29                           q_h2o_ts,q_co2_ts,d_h2oice,d_co2ice)
     29                          q_h2o_ts,q_co2_ts,minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost)
    3030!-----------------------------------------------------------------------
    3131! NAME
    … …  
    4747use io_netcdf,  only: xios_day_name2, xios_yr_name1, xios_yr_name2, open_nc, close_nc, get_var_nc, get_dim_nc
    4848use soil,       only: do_soil
    49 use tendencies, only: compute_tend
    5049use frost,      only: compute_frost4PCM
    5150use stoppage,   only: stop_clean
    … …  
    6059! ---------
    6160real(dp), dimension(:),       intent(out) :: ps_avg
    62 real(dp), dimension(:,:),     intent(out) :: tsurf_avg, tsurf_avg_yr1, h2o_surfdensity_avg, d_h2oice, d_co2ice, ps_ts, q_h2o_ts, q_co2_ts
     61real(dp), dimension(:,:),     intent(out) :: tsurf_avg, tsurf_avg_yr1, h2o_surfdensity_avg, ps_ts, q_h2o_ts, q_co2_ts
    6362real(dp), dimension(:,:,:),   intent(out) :: tsoil_avg, h2o_soildensity_avg
    6463real(dp), dimension(:,:,:,:), intent(out) :: tsoil_ts
     64real(dp), dimension(:,:,:),   intent(out) :: minPCM_h2operice, minPCM_co2perice, minPCM_h2ofrost, minPCM_co2frost
    6565
    6666! LOCAL VARIABLES
    … …  
    7272real(dp), dimension(:,:,:,:), allocatable    :: var_read_4d
    7373character(:),                 allocatable    :: num ! To read slope variables
    74 real(dp), dimension(ngrid,nslope,2)          :: min_h2operice, min_co2perice, min_h2ofrost, min_co2frost
    7574real(dp), dimension(ngrid,nsoil,nslope,nday) :: h2o_soildensity_ts
    7675! CODE
    7776! ----
    7877! Initialization
    79 min_h2operice(:,:,:) = 0._dp
    80 min_co2perice(:,:,:) = 0._dp
    81 min_h2ofrost(:,:,:) = 0._dp
    82 min_co2frost(:,:,:) = 0._dp
     78minPCM_h2operice(:,:,:) = 0._dp
     79minPCM_co2perice(:,:,:) = 0._dp
     80minPCM_h2ofrost(:,:,:) = 0._dp
     81minPCM_co2frost(:,:,:) = 0._dp
    8382if (nslope == 1) then ! No slope
    8483    allocate(character(0) :: num)
    … …  
    103102        num = '_slope'//num
    104103    end if
    105     call get_var_nc('co2ice'//num,var_read_2d)          ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2frost(:,islope,1))
    106     call get_var_nc('h2o_ice_s'//num,var_read_2d)       ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2ofrost(:,islope,1))
    107     call get_var_nc('watercap'//num,var_read_2d)        ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2operice(:,islope,1))
    108     call get_var_nc('perennial_co2ice'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2perice(:,islope,1))
     104    call get_var_nc('co2ice'//num,var_read_2d)          ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,minPCM_co2frost(:,islope,1))
     105    call get_var_nc('h2o_ice_s'//num,var_read_2d)       ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,minPCM_h2ofrost(:,islope,1))
     106    call get_var_nc('watercap'//num,var_read_2d)        ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,minPCM_h2operice(:,islope,1))
     107    call get_var_nc('perennial_co2ice'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,minPCM_co2perice(:,islope,1))
    109108    call get_var_nc('tsurf'//num,var_read_2d)           ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,tsurf_avg_yr1(:,islope))
    110109end do
    … …  
    129128    end if
    130129    call get_var_nc('tsurf'//num,var_read_2d)           ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,tsurf_avg(:,islope))
    131     call get_var_nc('co2ice'//num,var_read_2d)          ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2frost(:,islope,2))
    132     call get_var_nc('h2o_ice_s'//num,var_read_2d)       ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2ofrost(:,islope,2))
    133     call get_var_nc('watercap'//num,var_read_2d)        ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2operice(:,islope,2))
    134     call get_var_nc('perennial_co2ice'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2perice(:,islope,2))
     130    call get_var_nc('co2ice'//num,var_read_2d)          ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,minPCM_co2frost(:,islope,2))
     131    call get_var_nc('h2o_ice_s'//num,var_read_2d)       ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,minPCM_h2ofrost(:,islope,2))
     132    call get_var_nc('watercap'//num,var_read_2d)        ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,minPCM_h2operice(:,islope,2))
     133    call get_var_nc('perennial_co2ice'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,minPCM_co2perice(:,islope,2))
    135134    if (do_soil) then
    136135        call get_var_nc('soiltemp'//num,var_read_3d)
    … …  
    204203!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    205204! Compute frost from yearly minima
    206 call compute_frost4PCM(min_h2ofrost(:,:,2),min_co2frost(:,:,2))
    207 
    208 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    209 ! Compute ice tendencies from yearly minima
    210 call print_msg('> Computing surface ice tendencies')
    211 call compute_tend(min_h2operice + min_h2ofrost,d_h2oice)
    212 call print_msg('H2O ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_h2oice))//' | '//real2str(maxval(d_h2oice)))
    213 call compute_tend(min_co2perice + min_co2frost,d_co2ice)
    214 call print_msg('CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)))
     205call compute_frost4PCM(minPCM_h2ofrost(:,:,2),minPCM_co2frost(:,:,2))
    215206
    216207END SUBROUTINE load_xios_data
Note: See TracChangeset for help on using the changeset viewer.