Changeset 4170


Ignore:
Timestamp:
Apr 3, 2026, 4:34:51 PM (7 days ago)
Author:
jbclement
Message:

PEM:

  • Deletion of 'flux_ssice' ('zqdsdif_tot') from the "startfi.nc" as it is not needed.
  • Using the yearly average flux for the sublimating subsurface ice instead of the value at last timestep.
  • Keeping a clear separation between the subsurface ice flux and the surface ice tendency.
  • Making sure that subsurface ice depth is well given to the PCM at the end of the PEM (ice table VS layering).

JBC

Location:
trunk
Files:
26 edited

Legend:

Unmodified
Added
Removed

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

    r4157 r4170  
    2828use orbit,      only: ini_orbit, end_orbit
    2929use output,     only: ini_output, end_output
     30use ice_table,  only: ini_ice_table, end_ice_table
    3031use planet,     only: ini_planet, end_planet
    3132
    … …  
    8081call ini_orbit()
    8182call ini_output()
     83call ini_ice_table()
    8284
    8385END SUBROUTINE ini_allocation
    … …  
    106108! CODE
    107109! ----
     110call end_ice_table()
    108111call end_output()
    109112call end_orbit()

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

    r4157 r4170  
    118118!
    119119! DESCRIPTION
    120 !     Setter for 'atmosphere' configuration parameters.
     120!     Setter for "atmosphere" configuration parameters.
    121121!
    122122! AUTHORS & DATE
    … …  
    155155!
    156156! DESCRIPTION
    157 !     Setter for 'set_CO2cond_ps_PCM' parameter.
     157!     Setter for 'CO2cond_ps_PCM' parameter.
    158158!
    159159! AUTHORS & DATE

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

    r4157 r4170  
    3636!
    3737! DESCRIPTION
    38 !     Setter for backup configuration parameters.
     38!     Setter for "backup" configuration parameters.
    3939!
    4040! AUTHORS & DATE
    … …  
    7070!=======================================================================
    7171SUBROUTINE save_clim_state(h2o_ice,co2_ice,tsurf_avg,tsurf_dev,tsoil_avg,tsoil_dev,ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini, &
    72                            icetable_depth,icetable_thickness,ice_porefilling,h2o_ads_reg,co2_ads_reg,layerings_map,backup_idt)
     72                           icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map,backup_idt)
    7373!-----------------------------------------------------------------------
    7474! NAME
    … …  
    9595use atmosphere,       only: build4PCM_atmosphere
    9696use tracers,          only: build4PCM_tracers, nq
     97use ice_table,        only: build4PCM_ssice
    9798use clim_state_rec,   only: write_restart, write_restartfi, write_restartevo
    9899use layered_deposits, only: layering
    … …  
    104105! ARGUMENTS
    105106! ---------
    106 real(dp),       dimension(:),     intent(in)    :: ps_avg, ps_dev
    107 real(dp),                         intent(in)    :: ps_avg_glob, ps_avg_glob_ini
    108 real(dp),       dimension(:,:),   intent(in)    :: tsurf_avg, tsurf_dev, icetable_depth, icetable_thickness
    109 real(dp),       dimension(:,:,:), intent(in)    :: tsoil_avg, tsoil_dev, ice_porefilling, h2o_ads_reg, co2_ads_reg
    110 type(layering), dimension(:,:),   intent(in)    :: layerings_map
    111 integer(di), optional,            intent(in)    :: backup_idt
     107real(dp),       dimension(:),     intent(in)           :: ps_avg, ps_dev
     108real(dp),                         intent(in)           :: ps_avg_glob, ps_avg_glob_ini
     109real(dp),       dimension(:,:),   intent(in)           :: tsurf_avg, tsurf_dev, icetable_depth, icetable_thickness, h2oice_depth, flux_ssice_avg
     110real(dp),       dimension(:,:,:), intent(in)           :: tsoil_avg, tsoil_dev, ice_porefilling, h2o_ads_reg, co2_ads_reg
     111type(layering), dimension(:,:),   intent(in)           :: layerings_map
     112integer(di),                      intent(in), optional :: backup_idt
    112113real(dp),       dimension(:,:),   intent(inout) :: h2o_ice, co2_ice
    113114
    114115! LOCAL VARIABLES
    115116! ---------------
    116 real(dp),    dimension(ngrid,nslope)           :: h2o_ice4PCM, co2_ice4PCM, tsurf4PCM, flux_geo4PCM, albedo4PCM, emissivity4PCM
     117real(dp),    dimension(ngrid,nslope)           :: h2o_ice4PCM, co2_ice4PCM, tsurf4PCM, flux_geo4PCM, albedo4PCM, emissivity4PCM, flux_ssice4PCM, h2oice_depth4PCM
    117118real(dp),    dimension(ngrid,nlayer)           :: teta4PCM, air_mass4PCM
    118119real(dp),    dimension(ngrid,nsoil_PCM,nslope) :: tsoil4PCM, inertiesoil4PCM
    … …  
    130131call build4PCM_tsurf(tsurf_avg,tsurf_dev,tsurf4PCM)
    131132
    132 ! Build soil for the PCM
    133 if (do_soil) call build4PCM_soil(tsoil_avg,tsoil_dev,inertiesoil4PCM,tsoil4PCM,flux_geo4PCM)
     133if (do_soil) then
     134    ! Build soil for the PCM
     135    call build4PCM_soil(tsoil_avg,tsoil_dev,inertiesoil4PCM,tsoil4PCM,flux_geo4PCM)
     136
     137    ! Build subsurface ice for the PCM
     138    call build4PCM_ssice(icetable_depth,h2oice_depth,h2oice_depth4PCM)
     139end if
    134140
    135141! Build atmosphere for the PCM
    … …  
    144150! Write restart files
    145151call write_restartevo(h2o_ice,co2_ice,tsoil_avg,TI,icetable_depth,icetable_thickness,ice_porefilling,h2o_ads_reg,co2_ads_reg,layerings_map)
    146 call write_restartfi(is_h2o_perice,h2o_ice4PCM,co2_ice4PCM,tsurf4PCM,tsoil4PCM,inertiesoil4PCM,albedo4PCM,emissivity4PCM,flux_geo4PCM)
     152call write_restartfi(is_h2o_perice,h2o_ice4PCM,co2_ice4PCM,tsurf4PCM,tsoil4PCM,inertiesoil4PCM,albedo4PCM,emissivity4PCM,flux_geo4PCM,h2oice_depth4PCM)
    147153call write_restart(ps4PCM,pa4PCM,preff4PCM,q4PCM,teta4PCM,air_mass4PCM)
    148154
    … …  
    191197! CODE
    192198! ----
    193 suffix = ''
    194 suffix = suffix//'_ts'//int2str(backup_idt)
     199suffix = '_ts'//int2str(backup_idt)
    195200
    196201call print_msg('> Backup of "restart" files at dt = '//int2str(backup_idt),LVL_NFO)

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

    r4157 r4170  
    950950- Add lifecycle helpers for clear allocation/deallocation logic.
    951951- Simplify string suffix for slopes variables.
     952
     953== 03/04/2026 == JBC
     954- Deletion of 'flux_ssice' ('zqdsdif_tot') from the "startfi.nc" as it is not needed.
     955- Using the yearly average flux for the sublimating subsurface ice instead of the value at last timestep.
     956- Keeping a clear separation between the subsurface ice flux and the surface ice tendency.
     957- Making sure that subsurface ice depth is well given to the PCM at the end of the PEM (ice table VS layering).

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

    r4157 r4170  
    205205use surf_ice,  only: set_is_h2o_perice_PCM, set_co2_perice_PCM
    206206use soil,      only: set_TI_PCM, set_inertiedat_PCM
     207use ice_table, only: set_coef_ssdif_PCM
    207208use display,   only: print_msg, LVL_NFO
    208209
    … …  
    276277call set_tsoil_PCM(tmp3d)
    277278
     279call get_var_nc('coef_ssdif',tmp2d)
     280call set_coef_ssdif_PCM(tmp2d)
     281
    278282deallocate(tmp2d); allocate(tmp2d(ngrid,nsoil_PCM))
    279283call get_var_nc('inertiedat',tmp2d)
    … …  
    282286call get_var_nc('inertiesoil',tmp3d)
    283287call set_TI_PCM(tmp3d)
    284 
    285 ! To do?
    286 !   h2oice_depth
    287 !   d_coef
    288 !   lag_co2_ice
    289288
    290289! Close/Deallocate
    … …  
    398397                    !!! transition
    399398                    delta = depth_breccia
    400                     TI(i,index_breccia + 1,islope) = sqrt((layer(index_breccia + 1) - layer(index_breccia))/                       &
    401                                                                (((delta - layer(index_breccia))/(TI(i,index_breccia,islope)**2)) + &
    402                                                                ((layer(index_breccia + 1) - delta)/(TI_breccia**2))))
     399                    TI(i,index_breccia + 1,islope) = sqrt((layer(index_breccia + 1) - layer(index_breccia))/                  &
     400                                                          (((delta - layer(index_breccia))/(TI(i,index_breccia,islope)**2)) + &
     401                                                          ((layer(index_breccia + 1) - delta)/(TI_breccia**2))))
    403402                    do k = index_breccia + 2,index_bedrock
    404403                        TI(i,k,islope) = TI_breccia
    … …  
    411410                !! transition
    412411                delta = depth_bedrock
    413                 TI(i,index_bedrock + 1,islope) = sqrt((layer(index_bedrock + 1) - layer(index_bedrock))/                   &
    414                                                        (((delta - layer(index_bedrock))/(TI(i,index_bedrock,islope)**2)) + &
    415                                                        ((layer(index_bedrock + 1) - delta)/(TI_breccia**2))))
     412                TI(i,index_bedrock + 1,islope) = sqrt((layer(index_bedrock + 1) - layer(index_bedrock))/                  &
     413                                                      (((delta - layer(index_bedrock))/(TI(i,index_bedrock,islope)**2)) + &
     414                                                      ((layer(index_bedrock + 1) - delta)/(TI_breccia**2))))
    416415                do k = index_bedrock + 2,nsoil
    417416                    TI(i,k,islope) = TI_bedrock
    … …  
    427426        do i = 1,ngrid
    428427            if (inertiedat(i,index_breccia) < TI_breccia) then
    429                 inertiedat(i,index_breccia + 1) = sqrt((layer(index_breccia + 1) - layer(index_breccia))/                    &
    430                                                         (((delta - layer(index_breccia))/(inertiedat(i,index_breccia)**2)) + &
    431                                                         ((layer(index_breccia + 1) - delta)/(TI_breccia**2))))
     428                inertiedat(i,index_breccia + 1) = sqrt((layer(index_breccia + 1) - layer(index_breccia))/                   &
     429                                                       (((delta - layer(index_breccia))/(inertiedat(i,index_breccia)**2)) + &
     430                                                       ((layer(index_breccia + 1) - delta)/(TI_breccia**2))))
    432431                do k = index_breccia + 2,index_bedrock
    433432                    inertiedat(i,k) = TI_breccia
    … …  
    442441        delta = depth_bedrock
    443442        do i = 1,ngrid
    444             inertiedat(i,index_bedrock + 1) = sqrt((layer(index_bedrock + 1) - layer(index_bedrock))/                    &
    445                                                     (((delta - layer(index_bedrock))/(inertiedat(i,index_bedrock)**2)) + &
    446                                                     ((layer(index_bedrock + 1) - delta)/(TI_bedrock**2))))
     443            inertiedat(i,index_bedrock + 1) = sqrt((layer(index_bedrock + 1) - layer(index_bedrock))/                   &
     444                                                   (((delta - layer(index_bedrock))/(inertiedat(i,index_bedrock)**2)) + &
     445                                                   ((layer(index_bedrock + 1) - delta)/(TI_bedrock**2))))
    447446        end do
    448447

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

    r4157 r4170  
    180180
    181181!=======================================================================
    182 SUBROUTINE write_restartfi(is_h2o_perice,h2o_ice4PCM,co2_ice4PCM,tsurf4PCM,tsoil4PCM,inertiesoil4PCM,albedo4PCM,emissivity4PCM,flux_geo4PCM)
     182SUBROUTINE write_restartfi(is_h2o_perice,h2o_ice4PCM,co2_ice4PCM,tsurf4PCM,tsoil4PCM,inertiesoil4PCM,albedo4PCM,emissivity4PCM,flux_geo4PCM,h2oice_depth4PCM)
    183183!-----------------------------------------------------------------------
    184184! NAME
    … …  
    208208! ARGUMENTS
    209209! ---------
    210 real(dp),    dimension(:,:),   intent(in) :: h2o_ice4PCM, co2_ice4PCM, albedo4PCM, emissivity4PCM, tsurf4PCM, flux_geo4PCM
     210real(dp),    dimension(:,:),   intent(in) :: h2o_ice4PCM, co2_ice4PCM, albedo4PCM, emissivity4PCM, tsurf4PCM, flux_geo4PCM, h2oice_depth4PCM
    211211real(dp),    dimension(:,:,:), intent(in) :: tsoil4PCM, inertiesoil4PCM
    212212logical(k4), dimension(:),     intent(in) :: is_h2o_perice
    … …  
    240240
    241241! Orbital parameters (controle)
    242 controle(18) = obliquity ! Obliquity
     242controle(18) = obliquity  ! Obliquity
    243243controle(15) = perihelion ! Perihelion
    244 controle(16) = aphelion ! Aphelion
    245 controle(17) = date_peri ! Date of perihelion
     244controle(16) = aphelion   ! Aphelion
     245controle(17) = date_peri  ! Date of perihelion
    246246call put_var_nc('controle',controle)
    247247deallocate(controle)
    … …  
    259259call put_var_nc('emis',emissivity4PCM,1)
    260260call put_var_nc('flux_geo',flux_geo4PCM,1)
    261 
    262 ! h2oice_depth
    263 ! lag_co2_ice
    264 ! zdqsdif_ssi_tot
    265 ! d_coef
    266 
     261call put_var_nc('h2oice_depth',h2oice_depth4PCM)
     262
     263! Close
    267264call close_nc('re'//startfi_name)
    268265

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

    r4152 r4170  
    172172call getin('h2oice_huge_ini',h2oice_huge_ini_l)
    173173
    174 threshold_h2oice_cap_l = 460._dp ! kg.m-2 (= 0.5 m)
     174threshold_h2oice_cap_l = 460._dp ! kg/m2 (= 0.5 m)
    175175call getin('threshold_h2oice_cap',threshold_h2oice_cap_l)
    176176
    … …  
    186186call getin('do_layering',do_layering_l)
    187187
    188 d_dust_l = 5.78e-2_dp ! kg.m-2.y-1 (= 1.e-9 kg.m-2.s-1)
     188d_dust_l = 5.78e-2_dp ! kg/m2/y (= 1.e-9 kg/m2/s)
    189189call getin('d_dust',d_dust_l)
    190190

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

    r4117 r4170  
    5050!
    5151! DESCRIPTION
    52 !     Setter for 'display' configuration parameters.
     52!     Setter for "display" configuration parameters.
    5353!
    5454! AUTHORS & DATE

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

    r4147 r4170  
    7777!
    7878! DESCRIPTION
    79 !     Setter for 'evolution' configuration parameters.
     79!     Setter for "evolution" configuration parameters.
    8080!
    8181! AUTHORS & DATE

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

    r4135 r4170  
    3939!
    4040! DESCRIPTION
    41 !     Setter for 'glaciers' configuration parameters.
     41!     Setter for "glaciers" configuration parameters.
    4242!
    4343! AUTHORS & DATE

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

    r4138 r4170  
    2525! PARAMETERS
    2626! ----------
    27 logical(k4), protected :: icetable_equilibrium ! Flag to compute the icetable depth according to equilibrium
    28 logical(k4), protected :: icetable_dynamic     ! Flag to compute the icetable depth according to the dynamic method
     27logical(k4),                           protected :: icetable_equilibrium ! Flag to compute the icetable depth according to equilibrium
     28logical(k4),                           protected :: icetable_dynamic     ! Flag to compute the icetable depth according to the dynamic method
     29real(dp), dimension(:,:), allocatable, protected :: coef_ssdif_PCM       ! Diffusion coefficient for subsurface water
    2930
    3031contains
    … …  
    3233
    3334!=======================================================================
     35SUBROUTINE ini_ice_table()
     36!-----------------------------------------------------------------------
     37! NAME
     38!     ini_ice_table
     39!
     40! DESCRIPTION
     41!     Initialize the parameters of module 'ice_table'.
     42!
     43! AUTHORS & DATE
     44!     JB Clement, 03/2026
     45!
     46! NOTES
     47!
     48!-----------------------------------------------------------------------
     49
     50! DEPENDENCIES
     51! ------------
     52use geometry, only: ngrid, nslope
     53
     54! DECLARATION
     55! -----------
     56implicit none
     57
     58! CODE
     59! ----
     60allocate(coef_ssdif_PCM(ngrid,nslope))
     61coef_ssdif_PCM(:,:) = 0._dp
     62
     63END SUBROUTINE ini_ice_table
     64!=======================================================================
     65
     66!=======================================================================
     67SUBROUTINE end_ice_table()
     68!-----------------------------------------------------------------------
     69! NAME
     70!     end_ice_table
     71!
     72! DESCRIPTION
     73!     Deallocate ice_table arrays.
     74!
     75! AUTHORS & DATE
     76!     JB Clement, 03/2026
     77!
     78! NOTES
     79!
     80!-----------------------------------------------------------------------
     81
     82! DECLARATION
     83! -----------
     84implicit none
     85
     86! CODE
     87! ----
     88if (allocated(coef_ssdif_PCM)) deallocate(coef_ssdif_PCM)
     89
     90END SUBROUTINE end_ice_table
     91!=======================================================================
     92
     93!=======================================================================
    3494SUBROUTINE set_ice_table_config(icetable_equilibrium_in,icetable_dynamic_in)
    3595!-----------------------------------------------------------------------
    … …  
    3898!
    3999! DESCRIPTION
    40 !     Setter for 'ice_table' configuration parameters.
     100!     Setter for "ice_table" configuration parameters.
    41101!
    42102! AUTHORS & DATE
    … …  
    72132
    73133END SUBROUTINE set_ice_table_config
     134!=======================================================================
     135
     136!=======================================================================
     137SUBROUTINE set_coef_ssdif_PCM(coef_ssdif_PCM_in)
     138!-----------------------------------------------------------------------
     139! NAME
     140!     set_coef_ssdif_PCM
     141!
     142! DESCRIPTION
     143!     Setter for 'coef_ssdif_PCM' parameter.
     144!
     145! AUTHORS & DATE
     146!     JB Clement, 03/2026
     147!
     148! NOTES
     149!
     150!-----------------------------------------------------------------------
     151
     152! DECLARATION
     153! -----------
     154implicit none
     155
     156! ARGUMENTS
     157! ---------
     158real(dp), dimension(:,:), intent(in) :: coef_ssdif_PCM_in
     159
     160! CODE
     161! ----
     162coef_ssdif_PCM(:,:) = coef_ssdif_PCM_in(:,:)
     163
     164END SUBROUTINE set_coef_ssdif_PCM
    74165!=======================================================================
    75166
    … …  
    609700!=======================================================================
    610701
     702!=======================================================================
     703SUBROUTINE build4PCM_ssice(icetable_depth,h2oice_depth,h2oice_depth4PCM)
     704!-----------------------------------------------------------------------
     705! NAME
     706!     build4PCM_ssice
     707!
     708! DESCRIPTION
     709!     Reconstructs subsurface ice for the PCM.
     710!
     711! AUTHORS & DATE
     712!     JB Clement, 03/2026
     713!
     714! NOTES
     715!
     716!-----------------------------------------------------------------------
     717
     718! DEPENDENCIES
     719! ------------
     720use layered_deposits, only: do_layering
     721use display,          only: print_msg, LVL_NFO
     722
     723! DECLARATION
     724! -----------
     725implicit none
     726
     727! ARGUMENTS
     728! ---------
     729real(dp), dimension(:,:), intent(in)  :: icetable_depth, h2oice_depth
     730real(dp), dimension(:,:), intent(out) :: h2oice_depth4PCM
     731
     732! CODE
     733! ----
     734call print_msg('> Building subsurface ice for the PCM',LVL_NFO)
     735h2oice_depth4PCM(:,:) = -999. ! Default initialization (no subsurface ice)
     736if (do_layering) then
     737    h2oice_depth4PCM(:,:) = h2oice_depth(:,:)
     738else if (icetable_equilibrium .or. icetable_dynamic) then
     739    h2oice_depth4PCM(:,:) = icetable_depth(:,:)
     740end if
     741
     742END SUBROUTINE build4PCM_ssice
     743!=======================================================================
     744
    611745END MODULE ice_table

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

    r4135 r4170  
    3333logical(k4), protected, private :: impose_dust_ratio            ! Impose dust-to-ice ratio instead of dust tendency (see 'dust2ice_ratio')
    3434real(dp),    protected, private :: dust2ice_ratio               ! Dust-to-ice ratio when ice condenses (10% by default)
    35 real(dp),    protected, private :: d_dust                       ! Tendency of dust [kg.m-2.y-1]
     35real(dp),    protected, private :: d_dust                       ! Tendency of dust [kg/m2/y]
    3636real(dp),    parameter, private :: dry_lag_porosity  = 0.2_dp   ! Porosity of dust lag
    3737real(dp),    parameter, private :: wet_lag_porosity  = 0.15_dp  ! Porosity of water ice lag
    … …  
    8585!
    8686! DESCRIPTION
    87 !     Setter for 'layered_deposits' configuration parameters.
     87!     Setter for "layered_deposits" configuration parameters.
    8888!
    8989! AUTHORS & DATE
    … …  
    749749            h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice*rho_h2oice
    750750        else
    751             call subsurface_ice_layering(layerings_map(ig,islope),h2oice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
     751            call get_ssice_layering(layerings_map(ig,islope),h2oice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
    752752        end if
    753753    end do
    … …  
    10341034
    10351035!=======================================================================
    1036 SUBROUTINE subsurface_ice_layering(this,h_toplag,h2o_ice,co2_ice)
    1037 !-----------------------------------------------------------------------
    1038 ! NAME
    1039 !     subsurface_ice_layering
     1036SUBROUTINE get_ssice_layering(this,h_toplag,h2o_ice,co2_ice)
     1037!-----------------------------------------------------------------------
     1038! NAME
     1039!     get_ssice_layering
    10401040!
    10411041! DESCRIPTION
    … …  
    10951095end if
    10961096
    1097 END SUBROUTINE subsurface_ice_layering
     1097END SUBROUTINE get_ssice_layering
    10981098!=======================================================================
    10991099

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

    r4138 r4170  
    5959!
    6060! DESCRIPTION
    61 !     Setter for 'orbit' configuration parameters.
     61!     Setter for "orbit" configuration parameters.
    6262!
    6363! AUTHORS & DATE

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

    r4152 r4170  
    100100!
    101101! DESCRIPTION
    102 !     Setter for 'output' configuration parameters.
     102!     Setter for "output" configuration parameters.
    103103!
    104104! AUTHORS & DATE

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

    r4157 r4170  
    5252use surf_ice,           only: evolve_co2ice, evolve_h2oice, balance_h2oice_reservoirs
    5353use surf_temp,          only: tsurf_PCM, adapt_tsurf2disappearedice
    54 use tendencies,         only: compute_tendice, evolve_tend_co2ice, evolve_tend_h2oice
     54use tendencies,         only: compute_tendice, evolve_tend_co2ice, evolve_flux_ssice
    5555use tracers,            only: adapt_tracers2pressure
    5656use utility,            only: real2str
    … …  
    7272logical(k4), dimension(:,:),   allocatable :: is_co2ice_flow   ! Flag for location of CO2 glacier flow
    7373logical(k4), dimension(:,:),   allocatable :: is_h2oice_flow   ! Flag for location of H2O glacier flow
    74 real(dp),    dimension(:,:,:), allocatable :: minPCM_h2operice ! Minimum of H2O perennial ice over the last PCM year [kg.m-2]
    75 real(dp),    dimension(:,:,:), allocatable :: minPCM_co2perice ! Minimum of CO2 perennial ice over the last PCM year [kg.m-2]
    76 real(dp),    dimension(:,:,:), allocatable :: minPCM_h2ofrost  ! Minimum of H2O frost over the last PCM year [kg.m-2]
    77 real(dp),    dimension(:,:,:), allocatable :: minPCM_co2frost  ! Minimum of CO2 frost over the last PCM year [kg.m-2]
     74real(dp),    dimension(:,:,:), allocatable :: minPCM_h2operice ! Minimum of H2O perennial ice over the last PCM year [kg/m2]
     75real(dp),    dimension(:,:,:), allocatable :: minPCM_co2perice ! Minimum of CO2 perennial ice over the last PCM year [kg/m2]
     76real(dp),    dimension(:,:,:), allocatable :: minPCM_h2ofrost  ! Minimum of H2O frost over the last PCM year [kg/m2]
     77real(dp),    dimension(:,:,:), allocatable :: minPCM_co2frost  ! Minimum of CO2 frost over the last PCM year [kg/m2]
    7878! Surface-related:
    7979real(dp), dimension(:,:), allocatable :: tsurf_avg_yr1 ! Average surface temperature of the second to last PCM run [K]
    … …  
    141141
    142142call load_xios_data(ps_avg,ps_ts,tsurf_avg,tsurf_avg_yr1,tsoil_avg,tsoil_ts,h2o_surfdensity_avg,h2o_soildensity_avg, &
    143                     q_h2o_ts,q_co2_ts,minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost)
     143                    q_h2o_ts,q_co2_ts,minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost,flux_ssice_avg)
    144144
    145145! Initiate soil settings and TI
    … …  
    169169call print_msg('> Computing surface ice tendencies',LVL_NFO)
    170170call compute_tendice(minPCM_h2operice,minPCM_h2ofrost,h2o_ice,d_h2oice)
    171 call print_msg('H2O ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_h2oice))//' | '//real2str(maxval(d_h2oice)),LVL_NFO)
     171call print_msg('H2O ice tendencies [kg/m2/y] (min|max): '//real2str(minval(d_h2oice))//' | '//real2str(maxval(d_h2oice)),LVL_NFO)
    172172call compute_tendice(minPCM_co2perice,minPCM_co2frost,co2_ice,d_co2ice)
    173 call print_msg('CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO)
     173call print_msg('CO2 ice tendencies [kg/m2/y] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO)
    174174deallocate(minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost)
    175175
    … …  
    220220    ! Conversion to surface ice
    221221    call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth)
    222     ! We put the sublimating tendency coming from subsurface ice into the overall tendency
    223     !where (h2oice_depth > 0. .and. zdqsdif_ssi_tot < -minieps) d_h2oice = zdqsdif_ssi_tot
    224222end if
    225223call allocate_loop_state()
    … …  
    227225idt = 0
    228226do while (n_yr_run < min(nmax_yr_runorb,nmax_yr_run) .and. n_yr_sim < ntot_yr_sim)
    229     call print_msg('**** Iteration of the PEM run [plnt yr]: '//real2str(n_yr_run + dt),LVL_NFO)
     227    call print_msg('**** Iteration of the PEM run [plnt y]: '//real2str(n_yr_run + dt),LVL_NFO)
    230228    ! Evolve global surface pressure
    231229    call evolve_pressure(d_co2ice,delta_co2_ads,do_sorption,ps_avg_glob_old,ps_avg_glob,ps_avg)
    … …  
    237235    allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
    238236    if (do_layering) then
     237        allocate(d_h2oice_new(ngrid,nslope))
    239238        h2oice_depth_old(:,:) = h2oice_depth(:,:)
     239        ! The sublimating tendency coming from subsurface ice is given through the surface ice tendency
     240        where (h2oice_depth(:,:) > 0. .and. flux_ssice_avg(:,:) < 0._dp) d_h2oice(:,:) = flux_ssice_avg(:,:)
    240241        do islope = 1,nslope
    241242            do i = 1,ngrid
    … …  
    247248        call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth)
    248249        ! Balance H2O ice reservoirs
    249         allocate(d_h2oice_new(ngrid,nslope))
    250250        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)
    251251        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)
    … …  
    253253    else
    254254        zlag(:,:) = 0._dp
     255        zshift_surf(:,:) = 0._dp
    255256        call evolve_h2oice(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,zshift_surf,stopcrit)
    256257        call evolve_co2ice(co2_ice,d_co2ice,zshift_surf)
    … …  
    306307            write(num,'("_slope",i2.2)') islope
    307308        end if
    308         call write_diagevo('h2oice'//trim(num),'H2O ice','kg.m-2',h2o_ice(:,islope),(/dim_ngrid/))
    309         call write_diagevo('co2ice'//trim(num),'CO2 ice','kg.m-2',co2_ice(:,islope),(/dim_ngrid/))
    310         call write_diagevo('d_h2oice'//trim(num),'H2O ice tendency','kg.m-2.yr-1',d_h2oice(:,islope),(/dim_ngrid/))
    311         call write_diagevo('d_co2ice'//trim(num),'CO2 ice tendency','kg.m-2.yr-1',d_co2ice(:,islope),(/dim_ngrid/))
     309        call write_diagevo('h2oice'//trim(num),'H2O ice','kg/m2',h2o_ice(:,islope),(/dim_ngrid/))
     310        call write_diagevo('co2ice'//trim(num),'CO2 ice','kg/m2',co2_ice(:,islope),(/dim_ngrid/))
     311        call write_diagevo('d_h2oice'//trim(num),'H2O ice tendency','kg/m2/y',d_h2oice(:,islope),(/dim_ngrid/))
     312        call write_diagevo('d_co2ice'//trim(num),'CO2 ice tendency','kg/m2/y',d_co2ice(:,islope),(/dim_ngrid/))
    312313        if (co2ice_flow) then
    313314            call write_diagevo('Flow_co2ice'//trim(num),'CO2 ice flow location','T/F',merge(1._dp,0._dp,is_co2ice_flow(:,islope)),(/dim_ngrid/))
    … …  
    328329            call write_diagevo('inertiesoil'//trim(num),'Thermal inertia','SI',TI(:,:,islope),(/dim_ngrid,dim_nsoil/))
    329330            if (do_sorption) then
    330                 call write_diagevo('co2_ads_reg'//trim(num),'CO2 adsorbed in regolith','kg.m-2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
    331                 call write_diagevo('h2o_ads_reg'//trim(num),'H2O adsorbed in regolith','kg.m-2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
     331                call write_diagevo('co2_ads_reg'//trim(num),'CO2 adsorbed in regolith','kg/m2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
     332                call write_diagevo('h2o_ads_reg'//trim(num),'H2O adsorbed in regolith','kg/m2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
    332333            end if
    333334        end if
    … …  
    361362    ! Evolve the tendencies
    362363    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)
    363 !~     if (do_layering) then
    364 !~         do i = 1,ngrid
    365 !~             do islope = 1,nslope
    366 !~                 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))
    367 !~             end do
    368 !~         end do
    369 !    else
    370 !        do i = 1,ngrid
    371 !            do islope = 1,nslope
    372 !                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))
    373 !            end do
    374 !        end do
    375 !~     end if
     364    if (do_layering) then
     365        call evolve_flux_ssice(h2oice_depth_old,h2oice_depth,tsurf_avg,tsoil_ts_old,tsoil_ts,flux_ssice_avg)
     366    else if (icetable_equilibrium .or. icetable_dynamic) then
     367        call evolve_flux_ssice(icetable_depth_old,icetable_depth,tsurf_avg,tsoil_ts_old,tsoil_ts,flux_ssice_avg)
     368    end if
    376369
    377370    ! Increment time
    … …  
    383376    if (backup_rate > 0) then
    384377        if (mod(idt,backup_rate) == 0) call save_clim_state(h2o_ice,co2_ice,tsurf_avg,tsurf_dev,tsoil_avg,tsoil_dev,ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini, &
    385                                                             icetable_depth,icetable_thickness,ice_porefilling,h2o_ads_reg,co2_ads_reg,layerings_map,idt)
     378                                                            icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map,idt)
    386379    end if
    387380
    … …  
    415408
    416409call save_clim_state(h2o_ice,co2_ice,tsurf_avg,tsurf_dev,tsoil_avg,tsoil_dev,ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini, &
    417                      icetable_depth,icetable_thickness,ice_porefilling,h2o_ads_reg,co2_ads_reg,layerings_map)
     410                     icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map)
    418411
    419412! Update the duration information of the workflow

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

    r4157 r4170  
    3838
    3939! Ice-related:
    40 real(dp),    dimension(:,:), allocatable :: h2o_ice                    ! H2O ice [kg.m-2]
    41 real(dp),    dimension(:,:), allocatable :: co2_ice                    ! CO2 ice [kg.m-2]
     40real(dp),    dimension(:,:), allocatable :: h2o_ice                    ! H2O ice [kg/m2]
     41real(dp),    dimension(:,:), allocatable :: co2_ice                    ! CO2 ice [kg/m2]
    4242real(dp)                                 :: h2oice_sublim_coverage_ini ! Initial surface area of sublimating H2O ice [m2]
    4343real(dp)                                 :: co2ice_sublim_coverage_ini ! Initial surface area of sublimating CO2 ice [m2]
    … …  
    7373real(dp), dimension(:,:),   allocatable :: icetable_depth_old ! Old depth of the ice table [m]
    7474real(dp), dimension(:),     allocatable :: delta_icetable     ! Total mass of the H2O exchanged with the ice table [kg]
     75real(dp), dimension(:,:),   allocatable :: flux_ssice_avg     ! Average of total flux exchanged with subsurface ice [kg/m2/y]
    7576
    7677! Tracer-related:
    … …  
    153154allocate(q_h2o_ts(ngrid,nday))
    154155allocate(q_co2_ts(ngrid,nday))
     156allocate(flux_ssice_avg(ngrid,nslope))
    155157
    156158ps_avg(:) = 0._dp
    … …  
    163165q_h2o_ts(:,:) = 0._dp
    164166q_co2_ts(:,:) = 0._dp
     167flux_ssice_avg(:,:) = 0._dp
    165168
    166169END SUBROUTINE allocate_xios_state
    … …  
    422425if (allocated(d_h2oice)) deallocate(d_h2oice)
    423426if (allocated(layerings_map)) deallocate(layerings_map)
     427if (allocated(flux_ssice_avg)) deallocate(flux_ssice_avg)
    424428
    425429END SUBROUTINE end_planet

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

    r4136 r4170  
    148148!
    149149! DESCRIPTION
    150 !     Setter for 'soil' configuration parameters.
     150!     Setter for "soil" configuration parameters.
    151151!
    152152! AUTHORS & DATE

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

    r4135 r4170  
    3838!
    3939! DESCRIPTION
    40 !     Setter for 'sorption' configuration parameters.
     40!     Setter for "sorption" configuration parameters.
    4141!
    4242! AUTHORS & DATE

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

    r4152 r4170  
    5555!
    5656! DESCRIPTION
    57 !     Setter for 'stopping_crit' configuration parameters.
     57!     Setter for "stopping_crit" configuration parameters.
    5858!
    5959! AUTHORS & DATE

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

    r4152 r4170  
    2626real(dp),                                 parameter :: rho_co2ice = 1650._dp ! Density of CO2 ice [kg.m-3]
    2727real(dp),                                 parameter :: rho_h2oice = 920._dp  ! Density of H2O ice [kg.m-3]
    28 real(dp),                                 protected :: threshold_h2oice_cap  ! Threshold to consider H2O ice as infinite reservoir [kg.m-2]
     28real(dp),                                 protected :: threshold_h2oice_cap  ! Threshold to consider H2O ice as infinite reservoir [kg/m2]
    2929real(dp),                                 protected :: h2oice_huge_ini       ! Initial value for huge reservoirs of H2O ice [kg/m^2]
    30 logical(k4), dimension(:),   allocatable, protected :: is_h2o_perice_PCM     ! Flag for the presence of perennial H2O ice in the PCM at the beginning [kg.m-2]
    31 real(dp),    dimension(:,:), allocatable, protected :: co2_perice_PCM        ! Perennial CO2 ice in the PCM at the beginning [kg.m-2]
     30logical(k4), dimension(:),   allocatable, protected :: is_h2o_perice_PCM     ! Flag for the presence of perennial H2O ice in the PCM at the beginning [kg/m2]
     31real(dp),    dimension(:,:), allocatable, protected :: co2_perice_PCM        ! Perennial CO2 ice in the PCM at the beginning [kg/m2]
    3232
    3333contains
    … …  
    103103!
    104104! DESCRIPTION
    105 !     Setter for 'surf_ice' configuration parameters.
     105!     Setter for "surf_ice" configuration parameters.
    106106!
    107107! AUTHORS & DATE
    … …  
    306306call print_msg('> Evolution of CO2 ice',LVL_NFO)
    307307
    308 co2_ice_old = co2_ice
    309 co2_ice = co2_ice + d_co2ice*dt
     308co2_ice_old(:,:) = co2_ice(:,:)
     309co2_ice(:,:) = co2_ice(:,:) + d_co2ice(:,:)*dt
    310310where (co2_ice < 0._dp)
    311     co2_ice = 0._dp
    312     d_co2ice = -co2_ice_old/dt
     311    co2_ice(:,:) = 0._dp
     312    d_co2ice(:,:) = -co2_ice_old(:,:)/dt
    313313end where
    314314
    315 zshift_surf = d_co2ice*dt/rho_co2ice
     315zshift_surf(:,:) = zshift_surf(:,:) + d_co2ice(:,:)*dt/rho_co2ice
    316316
    317317END SUBROUTINE evolve_co2ice
    … …  
    379379
    380380h2o_ice(:,:) = h2o_ice(:,:) + d_h2oice_new(:,:)*dt
    381 zshift_surf(:,:) = d_h2oice_new(:,:)*dt/rho_h2oice
     381zshift_surf(:,:) = zshift_surf(:,:) + d_h2oice_new(:,:)*dt/rho_h2oice
    382382
    383383END SUBROUTINE evolve_h2oice

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

    r4145 r4170  
    9191if (allocated(albedo_PCM)) deallocate(albedo_PCM)
    9292if (allocated(albedodat_PCM)) deallocate(albedodat_PCM)
    93 if (allocated(emissivity_PCM)) deallocate(emissivity_PCM )
     93if (allocated(emissivity_PCM)) deallocate(emissivity_PCM)
    9494
    9595END SUBROUTINE end_surface

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

    r4152 r4170  
    114114! ----
    115115call print_msg("> Evolving the CO2 ice tendency according to the new pressure",LVL_NFO)
    116 call print_msg('Old CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO)
     116call print_msg('Old CO2 ice tendencies [kg/m2/y] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO)
    117117
    118118! Compute volume mixing ratio
    … …  
    135135    end do
    136136end do
    137 call print_msg('New CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str( minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO)
     137call print_msg('New CO2 ice tendencies [kg/m2/y] (min|max): '//real2str( minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO)
    138138
    139139END SUBROUTINE evolve_tend_co2ice
    … …  
    141141
    142142!=======================================================================
    143 SUBROUTINE evolve_tend_h2oice(h2oice_depth_old,h2oice_depth_new,tsurf,tsoil_ts_old,tsoil_ts_new,d_h2oice)
    144 !-----------------------------------------------------------------------
    145 ! NAME
    146 !     evolve_tend_h2oice
     143SUBROUTINE evolve_flux_ssice(h2oice_depth_old,h2oice_depth_new,tsurf,tsoil_ts_old,tsoil_ts_new,flux_ssice_avg)
     144!-----------------------------------------------------------------------
     145! NAME
     146!     evolve_flux_ssice
    147147!
    148148! DESCRIPTION
    … …  
    158158! DEPENDENCIES
    159159! ------------
     160use geometry,       only: ngrid, nslope
    160161use soil_temp,      only: itp_tsoil
    161162use subsurface_ice, only: psv
     163use ice_table,      only: coef_ssdif_PCM
    162164use display,        only: print_msg, LVL_NFO
    163165
    … …  
    168170! ARGUMENTS
    169171! ---------
    170 real(dp),                 intent(in)    :: h2oice_depth_old ! Old H2O ice depth
    171 real(dp),                 intent(in)    :: h2oice_depth_new ! New H2O ice depth
    172 real(dp),                 intent(in)    :: tsurf            ! Surface temperature
    173 real(dp), dimension(:,:), intent(in)    :: tsoil_ts_old     ! Old soil temperature time series
    174 real(dp), dimension(:,:), intent(in)    :: tsoil_ts_new     ! New soil temperature time series
    175 real(dp),                 intent(inout) :: d_h2oice         ! Evolution of perennial ice
     172real(dp), dimension(:,:),     intent(in) :: h2oice_depth_old ! Old H2O ice depth
     173real(dp), dimension(:,:),     intent(in) :: h2oice_depth_new ! New H2O ice depth
     174real(dp), dimension(:,:),     intent(in) :: tsurf            ! Surface temperature
     175real(dp), dimension(:,:,:,:), intent(in) :: tsoil_ts_old     ! Old soil temperature time series
     176real(dp), dimension(:,:,:,:), intent(in) :: tsoil_ts_new     ! New soil temperature time series
     177real(dp), dimension(:,:),     intent(inout) :: flux_ssice_avg       ! Tendency of sub-surface ice
    176178
    177179! LOCAL VARIABLES
    178180! ---------------
     181real(dp), parameter :: zcdv = 0.0325_dp ! Drag coefficient
    179182real(dp)            :: Rz_old, Rz_new, R_dec, hum_dec, psv_max_old, psv_max_new
    180 integer(di)         :: iday
    181 real(dp), parameter :: coef_diff = 4.e-4_dp ! Diffusion coefficient
    182 real(dp), parameter :: zcdv = 0.0325_dp     ! Drag coefficient
     183integer(di)         :: i, islope, iday
    183184
    184185! CODE
    185186! ----
    186 call print_msg("> Updating the H2O sub-surface ice tendency due to lag layer",LVL_NFO)
    187 
    188 ! Higher resistance due to growing lag layer (higher depth)
    189 Rz_old = h2oice_depth_old*zcdv/coef_diff ! Old resistance from PCM
    190 Rz_new = h2oice_depth_new*zcdv/coef_diff ! New resistance based on new depth
    191 R_dec = Rz_old/Rz_new ! Decrease because of resistance
    192 
    193 ! The maxmimum of the daily averages over one year for the saturation vapor pressure at the ice table location
    194 psv_max_old = 0._dp
    195 psv_max_new = 0._dp
    196 do iday = 1,size(tsoil_ts_old,2)
    197     psv_max_old = max(psv_max_old,psv(itp_tsoil(tsoil_ts_old(:,iday),tsurf,h2oice_depth_old)))
    198     psv_max_new = max(psv_max_new,psv(itp_tsoil(tsoil_ts_new(:,iday),tsurf,h2oice_depth_new)))
    199 end do
    200 
    201 ! Lower humidity due to growing lag layer (higher depth)
    202 if (abs(psv_max_old) < tol) then
    203     hum_dec = 1._dp
    204 else
    205     hum_dec = psv_max_new/psv_max_old ! Decrease because of lower water vapor pressure at the new depth
    206 end if
    207 
    208 ! Flux correction (decrease)
    209 d_h2oice = d_h2oice*R_dec*hum_dec
    210 
    211 END SUBROUTINE evolve_tend_h2oice
     187call print_msg("> Updating the H2O sub-surface ice tendency due to dust lag layer thickness",LVL_NFO)
     188
     189do i = 1,ngrid
     190    do islope = 1,nslope
     191        ! Higher resistance due to growing lag layer (higher depth)
     192        Rz_old = h2oice_depth_old(i,islope)*zcdv/coef_ssdif_PCM(i,islope) ! Old resistance from PCM
     193        Rz_new = h2oice_depth_new(i,islope)*zcdv/coef_ssdif_PCM(i,islope) ! New resistance based on new depth
     194        R_dec = Rz_old/Rz_new ! Decrease because of resistance
     195
     196        ! The maxmimum of the daily averages over one year for the saturation vapor pressure at the ice table location
     197        psv_max_old = 0._dp
     198        psv_max_new = 0._dp
     199        do iday = 1,size(tsoil_ts_old,2)
     200            psv_max_old = max(psv_max_old,psv(itp_tsoil(tsoil_ts_old(i,:,islope,iday),tsurf(i,islope),h2oice_depth_old(i,islope))))
     201            psv_max_new = max(psv_max_new,psv(itp_tsoil(tsoil_ts_new(i,:,islope,iday),tsurf(i,islope),h2oice_depth_new(i,islope))))
     202        end do
     203
     204        ! Lower humidity due to growing lag layer (higher depth)
     205        if (abs(psv_max_old) < tol) then
     206            hum_dec = 1._dp
     207        else
     208            hum_dec = psv_max_new/psv_max_old ! Decrease because of lower water vapor pressure at the new depth
     209        end if
     210
     211        ! Flux correction (decrease)
     212        flux_ssice_avg(i,islope) = flux_ssice_avg(i,islope)*R_dec*hum_dec
     213    enddo
     214enddo
     215
     216END SUBROUTINE evolve_flux_ssice
    212217!=======================================================================
    213218

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

    r4157 r4170  
    152152    n_id_PCM = n_pcm_runs
    153153end if
    154 headline = ''
    155 headline = headline//'cycle=PCM('//int2str(id_PCM_1)
     154headline = 'cycle=PCM('//int2str(id_PCM_1)
    156155do i = 1,n_id_PCM - 1
    157156    headline = headline//'+'//int2str(id_PCM_1 + i)

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

    r4157 r4170  
    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,minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost)
     29                          q_h2o_ts,q_co2_ts,minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost,flux_ssice_avg)
    3030!-----------------------------------------------------------------------
    3131! NAME
    … …  
    5959! ---------
    6060real(dp), dimension(:),       intent(out) :: ps_avg
    61 real(dp), dimension(:,:),     intent(out) :: tsurf_avg, tsurf_avg_yr1, h2o_surfdensity_avg, 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, flux_ssice_avg
    6262real(dp), dimension(:,:,:),   intent(out) :: tsoil_avg, h2o_soildensity_avg
    6363real(dp), dimension(:,:,:,:), intent(out) :: tsoil_ts
    … …  
    130130    call get_var_nc('perennial_co2ice'//trim(num),var_read_2d); call lonlat2vect(var_read_2d,minPCM_co2perice(:,islope,2))
    131131    if (do_soil) then
     132        call get_var_nc('flux_ssice'//trim(num),var_read_2d); call lonlat2vect(var_read_2d,flux_ssice_avg(:,islope))
    132133        call get_var_nc('soiltemp'//trim(num),var_read_3d)
    133134        do isoil = 1,nsoil_PCM

  • trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90

    r4166 r4170  
    3737use comsoil_h,           only: flux_geo
    3838use comslope_mod,        only: nslope, major_slope
    39 use paleoclimate_mod,    only: paleoclimate, h2oice_depth, co2ice_depth, coef_ssdif, zdqsdif_ssi_tot
     39use paleoclimate_mod,    only: paleoclimate, h2oice_depth, co2ice_depth, coef_ssdif
    4040use comcstfi_h,          only: pi
    4141use geometry_mod,        only: latitude
    … …  
    914914     h2oice_depth(:,:) = -999.
    915915   endif
    916    
    917 ! Total flux with SSI
    918    call get_field("flux_ssice",zdqsdif_ssi_tot,found,indextime)
    919    if (.not.found) then
    920      write(*,*) "phyetat0: Failed loading <flux_ssice> : ", &
    921                           "<flux_ssice> is set as 0 (no subsurface ice)"
    922      zdqsdif_ssi_tot(:,:) = 0.
    923    endif
    924916
    925917! Diffusion coeficent
    … …  
    950942    h2oice_depth(:,:) = -999.
    951943    co2ice_depth(:,:) = -999.
    952     zdqsdif_ssi_tot(:,:) = 0.
    953944    coef_ssdif(:,:) = 4.e-4
    954945    perennial_co2ice(:,:) = 0.
    … …  
    958949   h2oice_depth(:,:) = -999.
    959950   co2ice_depth(:,:) = -999.
    960    zdqsdif_ssi_tot(:,:) = 0.
    961951   coef_ssdif(:,:) = 4.e-4
    962952   perennial_co2ice(:,:) = 0.

  • trunk/LMDZ.MARS/libf/phymars/phyredem.F90

    r4160 r4170  
    210210  use callkeys_mod,        only: calltherm, dustinjection, calllott_nonoro
    211211  use callkeys_mod, only: CLFvarying
    212   use paleoclimate_mod,    only: paleoclimate, h2oice_depth, co2ice_depth, coef_ssdif, zdqsdif_ssi_tot
     212  use paleoclimate_mod,    only: paleoclimate, h2oice_depth, co2ice_depth, coef_ssdif
    213213
    214214  implicit none
    … …  
    390390    call put_field("h2oice_depth","Depth of the shallowest H2O ice layer",h2oice_depth)
    391391    call put_field("co2ice_depth","Depth of the shallowest CO2 ice layer",co2ice_depth)
    392     call put_field("flux_ssice","Total flux exchanged with subsurface water ice",zdqsdif_ssi_tot)
    393392    call put_field("coef_ssdif","Diffusion coefficent for subsurface water",coef_ssdif)
    394393  endif
Note: See TracChangeset for help on using the changeset viewer.