Changeset 4074


Ignore:
Timestamp:
Feb 17, 2026, 2:45:53 PM (8 weeks ago)
Author:
jbclement
Message:

PEM:

  • Correct management of H2O ice tendency in 1D when there is not enough ice anymore.
  • Clean initialization of allocatable module arrays (especially needed when no slope)
  • One more renaming for consistency + few small updates thoughout the code.

JBC

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
26 edited
1 moved

Legend:

Unmodified
Added
Removed

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

    r4065 r4074  
    6969if (.not. allocated(u_PCM)) allocate(u_PCM(ngrid,nlayer))
    7070if (.not. allocated(v_PCM)) allocate(v_PCM(ngrid,nlayer))
     71ap(:) = 0._dp
     72bp(:) = 0._dp
     73ps_PCM(:) = 0._dp
     74teta_PCM(:,:) = 0._dp
     75u_PCM(:,:) = 0._dp
     76v_PCM(:,:) = 0._dp
    7177
    7278END SUBROUTINE ini_atmosphere
    … …  
    662668! The pressure deviation is rescaled to avoid disproportionate oscillations in case of huge change of average pressure during the PEM run
    663669ps4PCM(:) = ps_avg(:) + ps_dev(:)*ps_avg_glob/ps_avg_glob_ini
    664 pa4PCM = ps_avg_glob_ini/30.  ! For now the altitude grid is not changed
    665 preff4PCM = ps_avg_glob_ini   ! For now the altitude grid is not changed
    666 
     670pa4PCM = ps_avg_glob_ini/30. ! For now the altitude grid is not changed
     671preff4PCM = ps_avg_glob_ini  ! For now the altitude grid is not changed
     672
     673! Correction on teta due to surface pressure changes
    667674call print_msg('> Building potential temperature for the PCM')
    668675do l = 1,nlayer
    669     ! Correction on teta due to surface pressure changes
    670676    teta4PCM(:,l) = teta_PCM(:,l)*ps4PCM(:)**rcp
    671677end do
    672678
     679! Compute atmospheric pressure
    673680call print_msg('> Building air mass for the PCM')
    674 ! Compute atmospheric pressure
    675681do l = 1,nlayer + 1
    676682    p(:,l) = ap(l) + bp(l)*ps4PCM(:)

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

    r4072 r4074  
    875875== 16/02/2026 == JBC
    876876Modify the PEM semantics to standardize file/functions/variable names and avoid ambiguity: workflow = repetion of cycles; cycle = PCM phase + PEM phase; PCM phase = 2+ PCM runs of 1 year each, PEM phase = 1 PEM run of 'x' years.
     877
     878== 17/02/2026 == JBC
     879- Correct management of H2O ice tendency in 1D when there is not enough ice anymore.
     880- Clean initialization of allocatable module arrays (especially needed when no slope)
     881- One more renaming for consistency + few small updates thoughout the code.

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

    r4072 r4074  
    2424
    2525character(7),  parameter :: rundef_name = 'run.def'
    26 character(11), parameter :: runPCMdef_name = 'run_PCM.def'
     26character(11), parameter :: runPCMdef_name = 'run_pcm.def'
    2727character(12), parameter :: callphys_name = 'callphys.def'
    2828
    … …  
    203203call set_layered_deposits_config(do_layering_l,impose_dust_ratio_l,d_dust_l,dust2ice_ratio_l)
    204204
    205 ! Read "run_PCM.def" parameters
     205! Read "run_pcm.def" parameters
    206206hybrid = .true. ! Default
    207207call get_hybrid(hybrid)
    … …  
    274274!
    275275! NOTES
    276 !     To work, it needs that "run_PEM.def" hols a line with
     276!     To work, it needs that "run_pem.def" hols a line with
    277277!     "INCLUDEDEF=callphys.def".
    278278!-----------------------------------------------------------------------
    … …  
    318318!
    319319! DESCRIPTION
    320 !     Get the key definition in "run_PCM.def".
     320!     Get the key definition in "run_pcm.def".
    321321!
    322322! AUTHORS & DATE

  • trunk/LMDZ.COMMON/libf/evolution/deftank/README

    r4072 r4074  
    4444    > the xml files for XIOS which can be found in the PCM deftank folder: "iodef.xml", "context_lmdz_physics.xml", "file_def_physics_mars.xml" and "field_def_physics_mars.xml";
    4545    > the def files you want to run the PCM: "run.def", "callphys.def", "traceur.def", etc.
    46       /!\ Do not forget to rename the PCM "run.def" into "run_PCM.def";
     46      /!\ Do not forget to rename the PCM "run.def" into "run_pcm.def";
    4747    > the starting files you want to run the PCM: "startfi.nc", "start.nc"/"start1D.txt"/profiles;
    48     > the necessary PEM files: "pem_workflow.sh", "pem_workflow_lib.sh", "pcm_run.job", "pem_run.job", "run_PEM.def" and "obl_ecc_lsp.asc";
     48    > the necessary PEM files: "pem_workflow.sh", "pem_workflow_lib.sh", "pcm_run.job", "pem_run.job", "run_pem.def" and "obl_ecc_lsp.asc";
    4949    > the optional PEM files: "diagevol.def" to define the PEM variables to be ouputted and "startpem.nc" to set the initial state of the PEM.
    5050
    … …  
    5555    > "pcm_run.job";
    5656    > "pem_run.job";
    57     > def files, especially for "run_PEM.def", "run_PCM.def", "callphys.def".
    58 In addition, the user has to provide a "startfi.nc" whose orbital parameters are consistent with the initial date set in "run_PEM.def". The script "ini_pem_orbit.sh" can do it automatically with "obl_ecc_lsp.asc".
     57    > def files, especially for "run_pem.def", "run_pcm.def", "callphys.def".
     58In addition, the user has to provide a "startfi.nc" whose orbital parameters are consistent with the initial date set in "run_pem.def". The script "ini_pem_orbit.sh" can do it automatically with "obl_ecc_lsp.asc".
    5959
    6060Outputs:
    … …  
    100100      > The PEM executable can have an optional argument which should be specified according to the set-up ("--auto-exit" for SLURM and PBS/TORQUE | "" when the script is not run as a job).
    101101
    102 # run_PEM.def
     102# run_pem.def
    103103  All the possible parameters to define a PEM run (read in "config.F90").
    104   It needs to include in "run_PCM.def" with "INCLUDEDEF=run_PEM.def".
     104  It needs to include in "run_pcm.def" with "INCLUDEDEF=run_pem.def".
    105105
    106106# obl_ecc_lsp.asc [default], obl_ecc_lsp_pos.asc [future years]
    … …  
    114114
    115115# ini_pem_orbit.sh:
    116   Bash script file to set the orbital parameters of a file "startfi.nc" from Laskar's data contained in "obl_ecc_lsp.asc" according to the initial date 'pem_ini_earth_date' defined in "run_PEM.def". See also "modify_startfi_orbit.sh".
     116  Bash script file to set the orbital parameters of a file "startfi.nc" from Laskar's data contained in "obl_ecc_lsp.asc" according to the initial date 'pem_ini_earth_date' defined in "run_pem.def". See also "modify_startfi_orbit.sh".
    117117
    118118# concat_pem.py:

  • trunk/LMDZ.COMMON/libf/evolution/deftank/ini_pem_orbit.sh

    r4072 r4074  
    22######################################################################
    33### Script to modify the orbital parameters of a file "startfi.nc" ###
    4 ### according to the date set in the file "run_PEM.def"            ###
     4### according to the date set in the file "run_pem.def"            ###
    55######################################################################
    66set -e
    … …  
    1717
    1818# Name of the file containing the starting date (Earth years)
    19 date_file="run_PEM.def"
     19date_file="run_pem.def"
    2020######################################################################
    2121

  • trunk/LMDZ.COMMON/libf/evolution/deftank/pcm_run.job

    r4072 r4074  
    2424
    2525# Name of executable for the PCM:
    26 exePCM="gcm_64x48x32_phymars_para.e"
     26pcm_exe="gcm_64x48x32_phymars_para.e"
    2727
    2828# Execution command:
    29 exe_cmd="srun --ntasks-per-node=${SLURM_NTASKS_PER_NODE} --cpu-bind=none --mem-bind=none --label -- ./adastra_cpu_binding.sh ./$exePCM"
     29exec_cmd="srun --ntasks-per-node=${SLURM_NTASKS_PER_NODE} --cpu-bind=none --mem-bind=none --label -- ./adastra_cpu_binding.sh ./$pcm_exe"
    3030########################################################################
    3131
    … …  
    3636read n_yr_sim ntot_yr_sim r_plnt2earth_yr i_pcm_run i_pem_run n_pcm_runs n_pcm_runs_ini < pem_workflow.sts
    3737echo "Run \"PCM $i_pcm_run\" is starting."
    38 cp run_PCM.def run.def
    39 eval "$exe_cmd > run.log 2>&1"
     38cp run_pcm.def run.def
     39eval "$exec_cmd > run.log 2>&1"
    4040if [ ! -f "restartfi.nc" ] || ! (tail -n 100 run.log | grep -iq "everything is cool!"); then # Check if it ended abnormally
    4141    echo "Error: the run \"PCM $i_pcm_run\" crashed!"

  • trunk/LMDZ.COMMON/libf/evolution/deftank/pem_run.job

    r4072 r4074  
    2121
    2222# Name of executable for the PEM:
    23 exePEM="pem_64x48x32_phymars_seq.e"
     23pem_exe="pem_64x48x32_phymars_seq.e"
    2424
    2525# Argument for the PEM execution ("--auto-exit" for SLURM and PBS/TORQUE | "" when the script is not run as a job):
    26 arg_pem="--auto-exit"
     26pem_arg="--auto-exit"
    2727########################################################################
    2828
    … …  
    3333read n_yr_sim ntot_yr_sim r_plnt2earth_yr i_pcm_run i_pem_run n_pcm_runs n_pcm_runs_ini < pem_workflow.sts
    3434echo "Run \"PEM $i_pem_run\" is starting."
    35 cp run_PEM.def run.def
    36 eval "./$exePEM $arg_pem > run.log 2>&1"
     35cp run_pem.def run.def
     36eval "./$pem_exe $pem_arg > run.log 2>&1"
    3737if [ ! -f "restartfi.nc" ] || ! (tail -n 100 run.log | grep -iq "so far, so good!"); then # Check if it ended abnormally
    3838    echo "Error: the run \"PEM $i_pem_run\" crashed!"

  • trunk/LMDZ.COMMON/libf/evolution/deftank/pem_workflow.sh

    r4072 r4074  
    6868        submit_cycle $exec_mode $n_pcm_runs
    6969
    70     # Starting a resume
     70    # Starting a resumption of the workflow
    7171    elif [ $1 = "re" ]; then
    7272        if [ ! -f "pem_workflow.sts" ]; then
    … …  
    132132        fi
    133133
    134     # Continuing the PEM run
    135     elif [ $1 = "cont" ]; then
    136         exec >> pem_workflow.log 2>&1
    137         echo
    138         echo "This is a continuation of the previous PEM run."
    139         date
    140         submit_pem_phase $exec_mode
     134    # Continuing a PEM run
     135    # CANNOT BE DONE FOR NOW BECAUSE THE PEM DOES NOT SAVE ITS STATE TO BE ABLE TO RECOVER
     136    #elif [ $1 = "cont" ]; then
     137    #    exec >> pem_workflow.log 2>&1
     138    #    echo
     139    #    echo "This is a continuation of the previous PEM run."
     140    #    date
     141    #    submit_pem_phase $exec_mode
    141142
    142143    # Default case: error

  • trunk/LMDZ.COMMON/libf/evolution/deftank/pem_workflow_lib.sh

    r4072 r4074  
    4848        ns=$(ncdump -h startfi.nc | sed -n 's/.*nslope = \([0-9]*\) ;.*/\1/p')
    4949    else
    50         for f in run_PCM.def callphys.def; do
     50        for f in run_pcm.def callphys.def; do
    5151            if [[ -f "$f" ]]; then
    5252                while IFS= read -r line; do
    … …  
    162162                  print val
    163163                  exit
    164                   }' run_PCM.def)
     164                  }' run_pcm.def)
    165165
    166166    if [ -z "$sol_in_file" ]; then
    167         echo "Error: no length of year found in \"run_PCM.def\"!"
     167        echo "Error: no length of year found in \"run_pcm.def\"!"
    168168        abort_workflow
    169169    elif [ "$sol_in_file" -eq "$year_sol" ]; then
    … …  
    171171        :
    172172    else
    173         echo "Error: length of year mismatch between \"run_PCM.def\" ($sol_in_file) and \"startfi.nc\" ($year_sol)!"
     173        echo "Error: length of year mismatch between \"run_pcm.def\" ($sol_in_file) and \"startfi.nc\" ($year_sol)!"
    174174        abort_workflow
    175175    fi
    … …  
    212212        abort_workflow
    213213    fi
    214     if [ ! -f "run_PCM.def" ]; then
    215         echo "Error: file \"run_PCM.def\" does not exist in $dir!"
    216         abort_workflow
    217     fi
    218     if [ ! -f "run_PEM.def" ]; then
    219         echo "Error: file \"run_PEM.def\" does not exist in $dir!"
     214    if [ ! -f "run_pcm.def" ]; then
     215        echo "Error: file \"run_pcm.def\" does not exist in $dir!"
     216        abort_workflow
     217    fi
     218    if [ ! -f "run_pem.def" ]; then
     219        echo "Error: file \"run_pem.def\" does not exist in $dir!"
    220220        abort_workflow
    221221    fi

  • trunk/LMDZ.COMMON/libf/evolution/deftank/run_pem.def

    r4073 r4074  
    1717# pem_ini_earth_date=0
    1818
    19 # Do you want the obliquity to eveolve when following "obl_ecc_lsp.asc"? Default = .true.
     19# Do you want the obliquity to evolve when following "obl_ecc_lsp.asc"? Default = .true.
    2020# evol_obl=.true.
    2121
    … …  
    5656
    5757# Do you want to run with adsorption/desorption in the PEM? Default = .false.
    58 # sorption=.false.
     58# do_sorption=.false.
    5959
    6060# Do you want to modify the soil thermal properties with the pressure? Default = .false.
    … …  
    102102# dust2ice_ratio=0.1
    103103
    104 # Some definitions for the physics in file "run_PCM.def"
    105 INCLUDEDEF=run_PCM.def
     104# Some definitions for the physics in file "run_pcm.def"
     105INCLUDEDEF=run_pcm.def

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

    r4071 r4074  
    6767if (.not. allocated(h2o_frost4PCM)) allocate(h2o_frost4PCM(ngrid,nslope))
    6868if (.not. allocated(co2_frost4PCM)) allocate(co2_frost4PCM(ngrid,nslope))
     69h2ofrost_PCM(:,:) = 0._dp
     70co2frost_PCM(:,:) = 0._dp
     71h2o_frost4PCM(:,:) = 0._dp
     72co2_frost4PCM(:,:) = 0._dp
    6973
    7074END SUBROUTINE ini_frost

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

    r4071 r4074  
    124124call print_msg("> Flow of CO2 glaciers")
    125125call computeTcondCO2(vmr_co2_PEM,ps_PCM,ps_avg_glob_ini,ps_avg_global,Tcond)
    126 call compute_hmaxglaciers(Tcond,"co2",hmax)
     126call compute_hmaxglaciers(Tcond,"CO2",hmax)
    127127call transfer_ice_duringflow(hmax,Tcond,co2ice,is_co2ice_flow)
    128128
    … …  
    169169! ----
    170170call print_msg("> Flow of H2O glaciers")
    171 call compute_hmaxglaciers(Tice,"h2o",hmax)
     171call compute_hmaxglaciers(Tice,"H2O",hmax)
    172172call transfer_ice_duringflow(hmax,Tice,h2oice,is_h2oice_flow)
    173173
    … …  
    222222! ----
    223223select case (trim(adjustl(name_ice)))
    224     case('h2o')
     224    case('H2O')
    225225        tau_d = 1.e5_dp
    226     case('co2')
     226    case('CO2')
    227227        tau_d = 5.e3_dp
    228228    case default
    … …  
    294294        if (islope /= iflat) then ! ice can be infinite on flat ground
    295295! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely  on flat ground
    296             if (qice(ig,islope) >= rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180._dp)) then
     296            if (qice(ig,islope) >= rho_ice(Tice(ig,islope),'H2O')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180._dp)) then
    297297! Second: determine the flatest slopes possible:
    298298                if (islope > iflat) then
    … …  
    308308                    end if
    309309                end do
    310                 qice(ig,iaval) = qice(ig,iaval) + (qice(ig,islope) - rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180._dp)) &
     310                qice(ig,iaval) = qice(ig,iaval) + (qice(ig,islope) - rho_ice(Tice(ig,islope),'H2O')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180._dp)) &
    311311                               *subslope_dist(ig,islope)/subslope_dist(ig,iaval)*cos(pi*def_slope_mean(iaval)/180._dp)/cos(pi*def_slope_mean(islope)/180._dp)
    312312
    313                 qice(ig,islope) = rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180._dp)
     313                qice(ig,islope) = rho_ice(Tice(ig,islope),'H2O')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180._dp)
    314314                flag_flow(ig,islope) = .true.
    315315            end if ! co2ice > hmax

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

    r4065 r4074  
    6767allocate(icetable_thickness(ngrid,nslope))
    6868allocate(ice_porefilling(ngrid,nsoil,nslope))
     69icetable_depth(:,:) = 0._dp
     70icetable_thickness(:,:) = 0._dp
     71ice_porefilling(:,:,:) = 0._dp
    6972
    7073END SUBROUTINE ini_ice_table
    … …  
    303306        do islope = 1,nslope
    304307            call compute_Tice(tsoil(ig,:,islope),tsurf(ig,islope),icetable_depth(ig,islope),Tice)
    305             delta_icetable(ig) = delta_icetable(ig) + porosity*rho_ice(Tice,'h2o')*(icetable_thickness(ig,islope) - icetable_thickness_old(ig,islope)) & ! convention > 0. <=> it condenses
     308            delta_icetable(ig) = delta_icetable(ig) + porosity*rho_ice(Tice,'H2O')*(icetable_thickness(ig,islope) - icetable_thickness_old(ig,islope)) & ! convention > 0. <=> it condenses
    306309                              *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180._dp)
    307310        end do
    … …  
    313316            do isoil = 1,nsoil
    314317                call compute_Tice(tsoil(ig,:,islope),tsurf(ig,islope),mlayer(isoil - 1),Tice)
    315                 delta_icetable(ig) = delta_icetable(ig) + rho_ice(Tice,'h2o')*(ice_porefilling(ig,isoil,islope) - ice_porefilling_old(ig,isoil,islope)) & ! convention > 0. <=> it condenses
     318                delta_icetable(ig) = delta_icetable(ig) + rho_ice(Tice,'H2O')*(ice_porefilling(ig,isoil,islope) - ice_porefilling_old(ig,isoil,islope)) & ! convention > 0. <=> it condenses
    316319                                  *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180._dp)
    317320            end do
    … …  
    586589! ----
    587590select case (trim(adjustl(ice)))
    588     case('h2o')
     591    case('H2O')
    589592        rho = -3.5353e-4_dp*T**2 + 0.0351_dp*T + 933.5030_dp ! Rottgers 2012
    590     case('co2')
     593    case('CO2')
    591594        rho = (1.72391_dp - 2.53e-4_dp*T - 2.87_dp*1.e-7_dp*T**2)*1.e3_dp ! Mangan et al. 2017
    592595    case default

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

    r4071 r4074  
    283283if (.not. allocated(ecc_lask)) allocate(ecc_lask(n))
    284284if (.not. allocated(lsp_lask)) allocate(lsp_lask(n))
     285year_lask(:) = 0._dp
     286obl_lask(:) = 0._dp
     287ecc_lask(:) = 0._dp
     288lsp_lask(:) = 0._dp
    285289
    286290END SUBROUTINE ini_orbit

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

    r4072 r4074  
    246246        if (h2o_ice(i,islope) > 0._dp) then
    247247            is_h2oice_ini(i,islope) = .true.
    248             if (d_co2ice(i,islope) < 0._dp) h2oice_sublim_coverage_ini = h2oice_sublim_coverage_ini + cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp)
     248            if (d_h2oice(i,islope) < 0._dp) h2oice_sublim_coverage_ini = h2oice_sublim_coverage_ini + cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp)
    249249        end if
    250250    end do

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

    r4065 r4074  
    2828real(dp), parameter :: m_co2 = 44.01e-3_dp               ! CO2 molecular mass [kg/mol]
    2929real(dp), parameter :: m_noco2 = 33.37e-3_dp             ! Non condensible mol mass [kg/mol]
    30 real(dp), parameter :: m_h2o = 18.01528e-3_dp            ! Molecular weight of h2o [kg/mol]
     30real(dp), parameter :: m_h2o = 18.01528e-3_dp            ! Molecular weight of H2O [kg/mol]
    3131real(dp), parameter :: A = (1._dp/m_co2 - 1._dp/m_noco2) ! Intermediate parameter for computation
    3232real(dp), parameter :: B = 1._dp/m_noco2                 ! Intermediate parameter for computation

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

    r4065 r4074  
    5959if (.not. allocated(def_slope_mean)) allocate(def_slope_mean(nslope))
    6060if (.not. allocated(subslope_dist)) allocate(subslope_dist(ngrid,nslope))
     61def_slope_mean(:) = 0._dp
     62subslope_dist(:,:) = 1._dp
    6163
    6264END SUBROUTINE ini_slopes

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

    r4065 r4074  
    9292if (.not. allocated(inertiedat_PCM)) allocate(inertiedat_PCM(ngrid,nsoil_PCM))
    9393if (.not. allocated(TI_PCM)) allocate(TI_PCM(ngrid,nsoil_PCM,nslope))
     94inertiedat_PCM(:,:) = 0._dp
     95TI_PCM(:,:,:) = 0._dp
    9496
    9597END SUBROUTINE ini_soil

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

    r4065 r4074  
    5959if (.not. allocated(tsoil_PCM)) allocate(tsoil_PCM(ngrid,nsoil_PCM,nslope))
    6060if (.not. allocated(flux_geo_PCM)) allocate(flux_geo_PCM(ngrid,nslope))
     61tsoil_PCM(:,:,:) = 0._dp
     62flux_geo_PCM(:,:) = 0._dp
    6163
    6264END SUBROUTINE ini_soil_temp

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

    r4065 r4074  
    102102real(dp), dimension(:,:,:), intent(in)    :: tsoil              ! Soil temperature [K]
    103103real(dp), dimension(:,:),   intent(in)    :: ps                 ! Average surface pressure [Pa]
    104 real(dp), dimension(:,:),   intent(in)    :: q_co2_ts           ! Mass mixing ratio of co2 in the first layer [kg/kg]
    105 real(dp), dimension(:,:),   intent(in)    :: q_h2o_ts           ! Mass mixing ratio of H2o in the first layer [kg/kg]
    106 real(dp), dimension(:),     intent(out)   :: delta_h2o_ads      ! Difference density of h2o adsorbed [kg/m^3]
    107 real(dp), dimension(:),     intent(out)   :: delta_co2_ads      ! Difference density of co2 adsorbed [kg/m^3]
     104real(dp), dimension(:,:),   intent(in)    :: q_co2_ts           ! Mass mixing ratio of CO2 in the first layer [kg/kg]
     105real(dp), dimension(:,:),   intent(in)    :: q_h2o_ts           ! Mass mixing ratio of H2O in the first layer [kg/kg]
     106real(dp), dimension(:),     intent(out)   :: delta_h2o_ads      ! Difference density of H2O adsorbed [kg/m^3]
     107real(dp), dimension(:),     intent(out)   :: delta_co2_ads      ! Difference density of CO2 adsorbed [kg/m^3]
    108108real(dp), dimension(:,:,:), intent(inout) :: h2o_ads_reg, co2_ads_reg
    109109
    … …  
    156156real(dp), dimension(:,:),   intent(in)    :: co2ice             ! CO2 ice at the surface [kg/m^2]
    157157real(dp), dimension(:,:),   intent(in)    :: ps                 ! Surface pressure [Pa]
    158 real(dp), dimension(:,:),   intent(in)    :: q_co2_ts           ! Mass mixing ratio of co2 in the first layer [kg/kg]
    159 real(dp), dimension(:,:),   intent(in)    :: q_h2o_ts           ! Mass mixing ratio of H2o in the first layer [kg/kg]
     158real(dp), dimension(:,:),   intent(in)    :: q_co2_ts           ! Mass mixing ratio of CO2 in the first layer [kg/kg]
     159real(dp), dimension(:,:),   intent(in)    :: q_h2o_ts           ! Mass mixing ratio of H2O in the first layer [kg/kg]
    160160real(dp), dimension(:,:,:), intent(in)    :: TI                 ! Soil Thermal inertia [J/K/^2/s^1/2]
    161161real(dp), dimension(:,:,:), intent(in)    :: tsoil              ! Soil temperature [K]
    162 real(dp), dimension(:,:,:), intent(inout) :: h2o_ads_reg        ! Density of h2o adsorbed [kg/m^3]
     162real(dp), dimension(:,:,:), intent(inout) :: h2o_ads_reg        ! Density of H2O adsorbed [kg/m^3]
    163163real(dp), dimension(:,:,:), intent(out)   :: theta_h2o_ads      ! Fraction of the pores occupied by H2O molecules
    164 real(dp), dimension(:),     intent(out)   :: delta_h2o_ads      ! Difference density of h2o adsorbed [kg/m^3]
     164real(dp), dimension(:),     intent(out)   :: delta_h2o_ads      ! Difference density of H2O adsorbed [kg/m^3]
    165165
    166166! LOCAL VARIABLES
    … …  
    170170real(dp) :: e = 2573.9_dp               ! Jackosky et al. 1997
    171171real(dp) :: mu = 0.48_dp                ! Jackosky et al. 1997
    172 real(dp) :: m_theta = 2.84e-7_dp        ! Mass of h2o per m^2 absorbed Jackosky et al. 1997
     172real(dp) :: m_theta = 2.84e-7_dp        ! Mass of H2O per m^2 absorbed Jackosky et al. 1997
    173173! real(dp) :: as = 18.9e3_dp              ! Specific area, Buhler & Piqueux 2021
    174174real(dp) :: as = 9.48e4_dp              ! Specific area, Zent
    … …  
    177177real(dp)                                        :: K                       ! Used to compute theta
    178178integer(di)                                     :: ig, iloop, islope, it   ! For loops
    179 logical(k4),  dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent
    180 logical(k4),  dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent
    181 real(dp), dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference density of h2o adsorbed per slope [kg/m^3]
    182 real(dp), dimension(ngrid,nsoil,nslope)         :: dm_h2o_regolith_slope   ! Elementary h2o mass adsorded per mesh per slope
     179logical(k4),  dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the CO2 glacier is permanent
     180logical(k4),  dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the H2O glacier is permanent
     181real(dp), dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference density of H2O adsorbed per slope [kg/m^3]
     182real(dp), dimension(ngrid,nsoil,nslope)         :: dm_h2o_regolith_slope   ! Elementary H2O mass adsorded per mesh per slope
    183183real(dp)                                        :: A, B                    ! Used to compute the mean mass above the surface
    184184real(dp)                                        :: p_sat                   ! Saturated vapor pressure of ice
    … …  
    299299real(dp), dimension(:,:,:), intent(in)    :: TI                 ! Mean Thermal Inertia [USI]
    300300real(dp), dimension(:,:),   intent(in)    :: d_h2oice, d_co2ice ! Tendencies on the glaciers ()
    301 real(dp), dimension(:,:),   intent(in)    :: q_co2_ts, q_h2o_ts ! Mass mixing ratio of co2 and h2o in the first layer [kg/kg]
     301real(dp), dimension(:,:),   intent(in)    :: q_co2_ts, q_h2o_ts ! Mass mixing ratio of CO2 and H2O in the first layer [kg/kg]
    302302real(dp), dimension(:,:),   intent(in)    :: h2oice             ! Water ice at the surface [kg/m^2]
    303303real(dp), dimension(:,:),   intent(in)    :: co2ice             ! CO2 ice at the surface [kg/m^2]
    304 real(dp), dimension(:,:,:), intent(inout) :: co2_ads_reg        ! Density of co2 adsorbed [kg/m^3]
    305 real(dp), dimension(:),     intent(out)   :: delta_co2_ads      ! Difference density of co2 adsorbed [kg/m^3]
     304real(dp), dimension(:,:,:), intent(inout) :: co2_ads_reg        ! Density of CO2 adsorbed [kg/m^3]
     305real(dp), dimension(:),     intent(out)   :: delta_co2_ads      ! Difference density of CO2 adsorbed [kg/m^3]
    306306
    307307! LOCAL VARIABLES
    … …  
    311311real(dp) :: beta = -1541.5_dp          ! Zent & Quinn 1995
    312312real(dp) :: inertie_thresold = 800._dp ! TI > 800 means cementation
    313 real(dp) :: m_theta = 4.27e-7_dp       ! Mass of co2 per m^2 absorbed
     313real(dp) :: m_theta = 4.27e-7_dp       ! Mass of CO2 per m^2 absorbed
    314314! real(dp) :: as = 18.9e3_dp           ! Specific area, Buhler & Piqueux 2021
    315315real(dp) :: as = 9.48e4_dp             ! Same as previous but from zent
    … …  
    317317integer(di)                                        :: ig, islope, iloop, it   ! For loops
    318318real(dp),    dimension(ngrid,nsoil,nslope)         :: dm_co2_regolith_slope   ! Elementary mass adsorded per mesh per slope
    319 logical(k4), dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent
    320 logical(k4), dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent
     319logical(k4), dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the CO2 glacier is permanent
     320logical(k4), dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the H2O glacier is permanent
    321321real(dp),    dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer [kg/m^3]
    322322real(dp),    dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference in the mass per slope  [kg/m^3]
    323323real(dp),    dimension(ngrid,nsoil,nslope)         :: m_h2o_adsorbed          ! Density of CO2 adsorbed [kg/m^3]
    324324real(dp),    dimension(ngrid,nsoil,nslope)         :: theta_h2o_ads           ! Fraction of the pores occupied by H2O molecules
    325 real(dp),    dimension(ngrid)                      :: delta_mh2o              ! Difference density of h2o adsorbed [kg/m^3]
     325real(dp),    dimension(ngrid)                      :: delta_mh2o              ! Difference density of H2O adsorbed [kg/m^3]
    326326! nday array are allocated because heavy...
    327327real(dp),    dimension(:,:), allocatable           :: mass_mean               ! Mean mass above the surface
    … …  
    368368call evolve_h2o_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_co2_ts,q_h2o_ts,tsoil,TI,theta_h2o_ads,m_h2o_adsorbed,delta_mh2o)
    369369
    370 ! 2. we compute the mass of co2 adsorded in each layer of the meshes
     370! 2. we compute the mass of CO2 adsorded in each layer of the meshes
    371371do ig = 1,ngrid
    372372    do islope = 1,nslope

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

    r4065 r4074  
    2424! PARAMETERS
    2525! ----------
    26 real(dp), protected :: h2oice_crit ! Change of the surface of h2o ice sublimating before stopping the PEM
    27 real(dp), protected :: co2ice_crit ! Change of the surface of co2 ice sublimating before stopping the PEM
     26real(dp), protected :: h2oice_crit ! Change of the surface of H2O ice sublimating before stopping the PEM
     27real(dp), protected :: co2ice_crit ! Change of the surface of CO2 ice sublimating before stopping the PEM
    2828real(dp), protected :: ps_crit     ! Change of averaged surface pressure before stopping the PEM
    2929
    … …  
    196196! CODE
    197197! ----
    198 if     (this%surface_h2oice_change)  then; msg = "**** STOPPING because surface of h2o ice has changed too much (see message above)."
    199 elseif (this%zero_dh2oice)           then; msg = "**** STOPPING because tendencies on h2o ice = 0 (see message above)."
    200 elseif (this%surface_co2ice_change)  then; msg = "**** STOPPING because surface of co2 ice has changed too much (see message above)."
     198if     (this%surface_h2oice_change)  then; msg = "**** STOPPING because surface of H2O ice has changed too much (see message above)."
     199elseif (this%zero_dh2oice)           then; msg = "**** STOPPING because tendencies on H2O ice = 0 (see message above)."
     200elseif (this%surface_co2ice_change)  then; msg = "**** STOPPING because surface of CO2 ice has changed too much (see message above)."
    201201elseif (this%pressure_change)        then; msg = "**** STOPPING because surface global pressure has changed too much (see message above)."
    202202elseif (this%nmax_yr_run_reached)    then; msg = "**** STOPPING because maximum number of years is reached."
    … …  
    253253! CODE
    254254! ----
    255 ! Check to escape the subroutine (not relevant in 1D or if the PEM should stop already)
    256 if (ngrid == 1 .or. stopcrit%stop_code() > 0) return
     255! Check to escape the subroutine (if the PEM should stop already)
     256if (stopcrit%stop_code() > 0) return
    257257
    258258! Computation of the present surface of sublimating H2O ice
    … …  
    313313! ---------------
    314314integer(di) :: i, islope
    315 real(dp)    :: co2ice_now_surf ! Current surface of co2 ice
     315real(dp)    :: co2ice_now_surf ! Current surface of CO2 ice
    316316
    317317! CODE
    … …  
    320320if (ngrid == 1 .or. stopcrit%stop_code() > 0) return
    321321
    322 ! Computation of the present surface of co2 ice still sublimating
     322! Computation of the present surface of CO2 ice still sublimating
    323323co2ice_now_surf = 0.
    324324do i = 1,ngrid
    … …  
    331331if (co2ice_now_surf < co2ice_sublim_coverage_ini*(1. - co2ice_crit) .or. co2ice_now_surf > co2ice_sublim_coverage_ini*(1._dp + co2ice_crit)) then
    332332    stopcrit%surface_co2ice_change = .true.
    333     call print_msg("Initial surface of co2 ice sublimating (ini|now) [m2] = "//real2str(co2ice_sublim_coverage_ini)//' | '//real2str(co2ice_now_surf))
     333    call print_msg("Initial surface of CO2 ice sublimating (ini|now) [m2] = "//real2str(co2ice_sublim_coverage_ini)//' | '//real2str(co2ice_now_surf))
    334334    call print_msg("Change (accepted|current) [%] = +/- "//real2str(co2ice_crit*100._dp)//' | '//real2str(100._dp*(co2ice_now_surf - co2ice_sublim_coverage_ini)/co2ice_sublim_coverage_ini))
    335335end if

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

    r4065 r4074  
    379379     fracIR = 0.04*p0(k)/600.; fracDust = 0.02*p0(k)/600.
    380380     B = Diff*bigstep*86400.*365.24/(porosity*927.)
    381      !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(Tb(),'h2o'))
     381     !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(Tb(),'H2O'))
    382382
    383383     typeT = -9
    … …  
    869869
    870870  B = Diff*bigstep*86400.*365.24/(porosity*927.)
    871   !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(T,'h2o'))
     871  !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(T,'H2O'))
    872872
    873873  ! advance ice table, avdrho>0 is retreat
    … …  
    889889        beta = (1-icefrac)/(1-porosity)/icefrac
    890890        beta = Diff*bigstep*beta*86400*365.24/927.
    891         !beta = Diff*bigstep*beta*86400*365.24/rho_ice(T,'h2o')
     891        !beta = Diff*bigstep*beta*86400*365.24/rho_ice(T,'H2O')
    892892        zdepthT = sqrt(2*beta*avdrho*18./8314. + zdepthT**2)
    893893     endif

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

    r4071 r4074  
    6262if (.not. allocated(is_h2o_perice_PCM)) allocate(is_h2o_perice_PCM(ngrid))
    6363if (.not. allocated(co2_perice_PCM)) allocate(co2_perice_PCM(ngrid,nslope))
     64is_h2o_perice_PCM(:) = .false.
     65co2_perice_PCM(:,:) = 0._dp
    6466
    6567END SUBROUTINE ini_surf_ice
    … …  
    359361
    360362if (ngrid == 1) then ! In 1D
    361     h2o_ice = h2o_ice + d_h2oice*dt
    362     zshift_surf = d_h2oice*dt/rho_h2oice
     363    where (d_h2oice(:,:) < 0._dp .and. abs(d_h2oice(:,:)*dt) > h2o_ice(:,:)) ! Not enough ice to sublimate everything
     364        ! We sublimate what we can
     365        d_h2oice_new(:,:) = h2o_ice(:,:)/dt
     366        ! It means the tendency is zero next time
     367        d_h2oice(:,:) = 0._dp
     368    else where
     369        d_h2oice_new(:,:) = d_h2oice(:,:)
     370    end where
    363371else ! In 3D
    364372    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)
    365373    if (stopcrit%stop_code() > 0) return
    366 
    367374    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)
    368     h2o_ice = h2o_ice + d_h2oice_new*dt
    369     zshift_surf = d_h2oice_new*dt/rho_h2oice
    370375end if
     376
     377h2o_ice = h2o_ice + d_h2oice_new*dt
     378zshift_surf = d_h2oice_new*dt/rho_h2oice
    371379
    372380END SUBROUTINE evolve_h2oice
    … …  
    408416! ---------------
    409417integer(di) :: i, islope
    410 real(qp)    :: S_target, S_target_subl_h2oice, S_target_cond_h2oice, S_ghostice, d_target ! Balance variables
     418real(qp)    :: S_target, S_target_subl_h2oice, S_target_cond_h2oice, S_ghostice ! Balance variables
     419real(dp)    :: d_target
    411420
    412421! CODE
    … …  
    421430    do islope = 1,nslope
    422431        if (d_h2oice(i,islope) > 0._dp) then ! Condensing
    423             d_h2oice_new(i,islope) = d_h2oice_new(i,islope)*real(S_target_cond_h2oice/S_atm_2_h2oice,dp)
     432            d_h2oice_new(i,islope) = d_h2oice(i,islope)*real(S_target_cond_h2oice/S_atm_2_h2oice,dp)
    424433        else if (d_h2oice(i,islope) < 0._dp) then ! Sublimating
    425434            d_target = d_h2oice(i,islope)*real(S_target_subl_h2oice/S_h2oice_2_atm,dp)
    426             if (abs(d_target*dt) < h2o_ice(i,islope)) then ! Enough ice to sublimate everything
    427                 d_h2oice_new(i,islope) = real(d_target,dp)
     435            if (abs(d_target*dt) <= h2o_ice(i,islope)) then ! Enough ice to sublimate everything
     436                d_h2oice_new(i,islope) = d_target
    428437            else ! Not enough ice to sublimate everything
    429438                ! We sublimate what we can
    … …  
    438447end do
    439448
    440 ! We need to remove this ice unable to sublimate from places where ice condenseds in order to keep balance
     449! We need to remove this ice unable to sublimate from places where ice condensed in order to keep balance
    441450where (d_h2oice_new > 0._dp) d_h2oice_new = d_h2oice_new*real((S_target_cond_h2oice - S_ghostice)/S_target_cond_h2oice,dp)
    442451

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

    r4065 r4074  
    5656! ----
    5757if (.not. allocated(tsurf_PCM)) allocate(tsurf_PCM(ngrid,nslope))
     58tsurf_PCM(:,:) = 0._dp
    5859
    5960END SUBROUTINE ini_surf_temp

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

    r4065 r4074  
    6060if (.not. allocated(albedodat_PCM)) allocate(albedodat_PCM(ngrid))
    6161if (.not. allocated(emissivity_PCM)) allocate(emissivity_PCM(ngrid,nslope))
     62albedo_PCM(:,:) = 0._dp
     63albedodat_PCM(:) = 0._dp
     64emissivity_PCM(:,:) = 1._dp
    6265
    6366END SUBROUTINE ini_surface

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

    r4071 r4074  
    178178real(dp)            :: Rz_old, Rz_new, R_dec, hum_dec, psv_max_old, psv_max_new
    179179integer(di)         :: iday
    180 real(dp), parameter :: coef_diff = 4.e-4 ! Diffusion coefficient
    181 real(dp), parameter :: zcdv = 0.0325     ! Drag coefficient
     180real(dp), parameter :: coef_diff = 4.e-4_dp ! Diffusion coefficient
     181real(dp), parameter :: zcdv = 0.0325_dp     ! Drag coefficient
    182182
    183183! CODE

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

    r4065 r4074  
    8080if (.not. allocated(qnames)) allocate(qnames(nq))
    8181if (.not. allocated(q_PCM)) allocate(q_PCM(ngrid,nlayer,nq))
     82mmol(:) = 0._dp
     83q_PCM(:,:,:) = 0._dp
    8284
    8385! Initialize the names of tracers
Note: See TracChangeset for help on using the changeset viewer.