Ignore:
Timestamp:
Oct 29, 2025, 10:02:05 AM (5 weeks ago)
Author:
jbclement
Message:

PEM:

  • Correction of the process to balance the H2O flux from and into the atmosphere accross reservoirs: (i) computation of H2O amount going from/in the atmosphere is corrected (missing dt and wrong condition); (ii) computation of the scaling factor is corrected because we need to balance all the icetable/adsorbed/surface ice flux but it affects only sublimating/condensing surface ice flux; (iii) process of balancing is modified to be applied to every flux by scaling instead of applying it only to positive or negative flux by trimming; (iv) balance is now fully processed by the tendency before evolving the ice. This is done through dedicated subroutines.
  • Correction of ice conversion between the layering algorithm and PEM ice variables (density factor was missing).
  • Addition of H2O flux balance and related stopping criterion for the layering algorithm.
  • Few updates for files in the deftank to be in line with the wiki Planeto pages.

JBC

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

Legend:

Unmodified
Added
Removed

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

    r3933 r3938  
    781781
    782782== 23/10/2025 == JBC
    783 Correction on dimensions for the new algoirthm introduced in r3907 to update of surface temperatue when ice disappeared.
     783Correction on dimensions for the new algorithm introduced in r3907 to update the surface temperature when ice disappeared.
     784
     785== 29/10/2025 == JBC
     786- Correction of the process to balance the H2O flux from and into the atmosphere accross reservoirs: (i) computation of H2O amount going from/in the atmosphere is corrected (missing dt and wrong condition); (ii) computation of the scaling factor is corrected because we need to balance all the icetable/adsorbed/surface ice flux but it affects only sublimating/condensing surface ice flux; (iii) process of balancing is modified to be applied to every flux by scaling instead of applying it only to positive or negative flux by trimming; (iv) balance is now fully processed by the tendency before evolving the ice. This is done through dedicated subroutines.
     787- Correction of ice conversion between the layering algorithm and PEM ice variables (density factor was missing).
     788- Addition of H2O flux balance and related stopping criterion for the layering algorithm.
     789- Few updates for files in the deftank to be in line with the wiki Planeto pages.

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

    r3933 r3938  
    2222------------
    2323To compile the PEM, in "LMDZ.COMMON", do: ./makelmdz_fcm -arch [local] -p [planet] -d [dimensions] -j 8 pem
    24 Options:
     24Options with example:
    2525    1) [local]     : root name of arch files, assuming that they have been set up for your configuration;
    26     2) [planet]    : mars to use the planet physics package;
    27     3) [dimensions]: 64x48x54 to define the grid you want to use.
     26    2) [planet]    : mars to use the Mars planet physics package;
     27    3) [dimensions]: 64x48x54 to define the grid you want to use  (longitude x latitude x atmospheric layers).
    2828To run the PEM, you need a dedicated reshaping tool with consistent options. To compile it, in "LMDZ_COMMON", do: ./makelmdz_fcm -arch [local] -p [planet] -d [dimensions] -j 8 reshape_XIOS_output
    29 To run the PEM, you also need a PCM working with XIOS and consistent options. To compile it, in "LMDZ.COMMON", do for example: ./makelmdz_fcm -arch [local] -p [planet] -parallel mpi_omp -io XIOS -d [dimensions] -j 8 gcm
     29To run the PEM, you also need a PCM working with XIOS and consistent options. To compile it, in "LMDZ.COMMON", do: ./makelmdz_fcm -arch [local] -p [planet] -parallel mpi_omp -io XIOS -d [dimensions] -j 8 gcm
    3030After compilation, the executable file can be found in the "bin" sub-directory.
    3131
    3232Usage:
    3333------
    34   To run a PEM simulation, do: ./launchPEM.sh [options]
    35   Options:
    36       1) None: to start a simulation from scratch;
    37       2) 're': to relaunch a simulation from a starting point (interactive prompt).
     34To run a PEM simulation, do: ./launchPEM.sh [options]
     35Options:
     36    1) None: to start a simulation from scratch;
     37    2) 're': to relaunch a simulation from a starting point (interactive prompt).
     38
     39The Bash file ''launchPEM.sh'' is the master script to launch the PEM chained simulation. It checks if necessary files and required options for your simulation are ok.
    3840
    3941Requirements:
    4042-------------
    41 To run the PEM, you need the following files:
    42     > your executable files for the PCM, the PEM and the reshaping tool with consistent dimensions;
    43     > the xml files for XIOS ("iodef.xml", "context_lmdz_physics.xml", "file_def_physics_mars.xml" and "field_def_physics_mars.xml") which can be found in the PCM deftank folder;
    44     > the def files you want to run the PCM ("run.def", "callphys.def", "traceur.def", etc).
     43To run the PEM, you can create a folder in which you need the following files:
     44    > your executable files for the PCM, the PEM and the reshaping tool with consistent options;
     45    > 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";
     46    > the def files you want to run the PCM: "run.def", "callphys.def", "traceur.def", etc.
    4547      /!\ Do not forget to rename the PCM "run.def" into "run_PCM.def";
    46     > the starting files you want to run the PCM ("startfi.nc", "start.nc"/"start1D.txt"/profiles);
    47     > the necessary PEM files ("launchPEM.sh", "lib_launchPEM.sh", "PCMrun.job", "PEMrun.job", "run_PEM.def" and "obl_ecc_lsp.asc";
    48     > the optional PEM files "diagpem.def" to define the PEM variables to be ouputted and "startpem.nc" to set the initial state of the PEM.
     48    > the starting files you want to run the PCM: "startfi.nc", "start.nc"/"start1D.txt"/profiles;
     49    > the necessary PEM files: "launchPEM.sh", "lib_launchPEM.sh", "PCMrun.job", "PEMrun.job", "run_PEM.def" and "obl_ecc_lsp.asc";
     50    > the optional PEM files: "diagpem.def" to define the PEM variables to be ouputted and "startpem.nc" to set the initial state of the PEM.
     51
     52The PEM files can be found in the deftank folder.
     53
     54Before a simulation, you have to set up some parameters/options in:
     55    > "launchPEM.sh";
     56    > "PCMrun.job";
     57    > "PEMrun.job";
     58    > def files, especially for "run_PEM.def", "run_PCM.def", "callphys.def".
     59In 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 "inipem_orbit.sh" can do it automatically with "obl_ecc_lsp.asc".
    4960
    5061Outputs:
    5162--------
    5263The PEM simulation generates the following files:
    53     > the usual outputs of the PCM ("restartfi.nc", "restart.nc", "diagfi.nc", etc);
    54     > the XIOS outputs of the PCM, then reshaped ("Xdiurnalave.nc"/"data2reshape*.nc"/"data_PCM_Y*.nc");
    55     > the outputs of the chained simulation ("launchPEM.log", "info_PEM.txt" and possibly "kill_launchPEM.sh");
    56     > the usual outputs of the PEM ("restartfi.nc", "restart.nc"/"restart1D.txt" and "diagpem.nc").
     64    > the usual outputs of the PCM: "restartfi.nc", "restart.nc", "diagfi.nc", etc;
     65    > the XIOS outputs of the PCM, then reshaped: "Xdiurnalave.nc"/"data2reshape*.nc"/"data_PCM_Y*.nc";
     66    > the outputs of the chained simulation: "launchPEM.log", "info_PEM.txt" and possibly "kill_launchPEM.sh";
     67    > the usual outputs of the PEM: "restartfi.nc", "restart.nc"/"restart1D.txt" and "diagpem.nc".
    5768During the simulation, the PCM/PEM run files are renamed conveniently and stored in the sub-directories "logs" (log files), "starts" (starting files) and "diags" (diagnostic files).
    5869If you run a simulation by submitting jobs, the script "kill_launchPEM.sh" is automatically generated. It can be used to kill in the queue of the job scheduler the jobs related to your chained simulation.
    … …  
    6980      > mode -> the launching mode (0 = "processing scripts"; any other values = "submitting jobs"). The former option is usually used to process the script on a local machine while the latter is used to submit jobs on a supercomputer with SLURM or PBS/TORQUE.
    7081  The script can take an argument:
    71       1) If there is no argument, then the script initiates a PEM simulation from scratch.
    72       2) If the argument is 're', then the script relaunches an existing PEM simulation. It will ask for parameters to know the starting point that you want to.
     82      1) None: to start a simulation from scratch;
     83      2) 're': to relaunch a simulation from a starting point (interactive prompt).
    7384
    7485# liblaunchPEM.sh:
    … …  
    7687
    7788# PCMrun.job:
    78   Bash script file to submit a PCM job (with SLURM or PBS/TORQUE). The headers correspond to the ADASTRA supercomputer and should be changed for other machines and job schedulers. In case of "processing scripts" launching mode, the headers are naturally omitted.
    79   The path to source the arch file should be adapted to the machine.
    80   The name of the PCM executable file should be adapted.
    81   The execution command should also be adapted according to the set-up.
     89  Bash script file to submit a PCM job (with SLURM or PBS/TORQUE).
     90  The user has to specify:
     91      > The headers correspond to the ADASTRA supercomputer and should be changed for other machines and job schedulers. In case of "processing scripts" launching mode, the headers are naturally omitted.
     92      > The path to source the arch file should be adapted to the machine.
     93      > The name of the PCM executable file should be adapted.
     94      > The execution command should also be adapted according to the set-up.
    8295
    8396# PEMrun.job:
    84   Bash script file to submit PEM job (with SLURM or PBS/TORQUE).The headers correspond to the ADASTRA supercomputer and should be changed for other machines and job schedulers. In case of "processing scripts" launching mode, the headers are naturally omitted.
    85   The path to source the arch file should be adapted to the machine.
    86   The name of the PEM executable file and Reshaping executable file should be adapted.
    87   The PEM executable can have an optional argument which should be specified according to the set-up. Especially, the value of '--jobid' which is the job ID to make the PEM detect the job time limit.
     97  Bash script file to submit PEM job (with SLURM or PBS/TORQUE).
     98  The user has to specify:
     99      > The headers correspond to the ADASTRA supercomputer and should be changed for other machines and job schedulers. In case of "processing scripts" launching mode, the headers are naturally omitted.
     100      > The path to source the arch file should be adapted to the machine.
     101      > The name of the PEM and Reshaping executable files should be adapted.
     102      > 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).
    88103
    89104# run_PEM.def

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

    r3926 r3938  
    44########################################################################
    55# This script can take an argument:
    6 #   - If there is no argument, then the script initiates a PEM simulation from scratch.
    7 #   - If the argument is 're', then the script relaunches an existing PEM simulation.
    8 #     It will ask for parameters to know the starting point that you want to.
     6#     1) None: to start a simulation from scratch;
     7#     2) 're': to relaunch a simulation from a starting point (interactive prompt).
    98########################################################################
    109set -e
    … …  
    2524nPCM=2
    2625
    27 # Set the counting method for the number of years to be simulated (0 = "only PEM runs count"; any other values = "PCM runs are taken into account"):
     26# Set the counting method for the number of years to be simulated (0 = "only PEM runs count"; any other values = "PCM runs are taken into account"). The former option is the usual one:
    2827counting=0
    2928

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

    r3603 r3938  
    5151SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
    5252
    53 use time_evol_mod, only: dt
    54 use comslope_mod,  only: subslope_dist, def_slope_mean
    55 use glaciers_mod,  only: rho_h2oice
    56 
    57 #ifndef CPP_STD
    58     use comcstfi_h,   only: pi
    59 #else
    60     use comcstfi_mod, only: pi
    61 #endif
     53use time_evol_mod,     only: dt
     54use glaciers_mod,      only: rho_h2oice
     55use stopping_crit_mod, only: stopping_crit_h2o
    6256
    6357implicit none
    … …  
    7771
    7872! OUTPUT
    79 real, dimension(ngrid,nslope), intent(inout) :: h2o_ice  ! Previous and actual density of h2o ice (kg/m^2)
    80 real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! Evolution of perennial ice over one year (kg/m^2/year)
     73real, dimension(ngrid,nslope), intent(inout) :: h2o_ice  ! H2O ice (kg/m^2)
     74real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! Tendency of H2O ice (kg/m^2/year)
    8175integer,                       intent(inout) :: stopPEM  ! Stopping criterion code
    8276real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m]
    … …  
    8478! local:
    8579! ------
    86 integer                       :: i, islope                                           ! Loop variables
    87 real                          :: pos_tend, neg_tend, real_coefficient, negative_part ! Variable to conserve h2o
    88 real, dimension(ngrid,nslope) :: new_tend                                            ! Tendencies computed in order to conserve h2o ice on the surface, only exchange between surface are done
     80integer                       :: i, islope                                                ! Loop variables
     81real                          :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to balance H2O
     82real, dimension(ngrid,nslope) :: d_h2oice_new                                             ! Tendencies computed to keep balance
    8983!=======================================================================
    9084write(*,*) '> Evolution of H2O ice'
    91 if (ngrid /= 1) then ! Not in 1D
    92     ! We compute the amount of condensing and sublimating h2o ice
    93     pos_tend = 0.
    94     neg_tend = 0.
    95     do i = 1,ngrid
    96         if (delta_h2o_adsorbed(i) > 0.) then
    97             pos_tend = pos_tend + delta_h2o_adsorbed(i)*cell_area(i)
    98         else
    99             neg_tend = neg_tend + delta_h2o_adsorbed(i)*cell_area(i)
    100         endif
    101         if (delta_h2o_icetablesublim(i) > 0.) then
    102             pos_tend = pos_tend + delta_h2o_icetablesublim(i)*cell_area(i)
    103         else
    104             neg_tend = neg_tend + delta_h2o_icetablesublim(i)*cell_area(i)
    105         endif
    106         do islope = 1,nslope
    107             if (h2o_ice(i,islope) > 0.) then
    108                 if (d_h2oice(i,islope) > 0.) then
    109                     pos_tend = pos_tend + d_h2oice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
    110                 else
    111                     neg_tend = neg_tend - d_h2oice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
    112                 endif
    113             endif
    114         enddo
    115     enddo
    11685
    117     new_tend = 0.
    118     if (abs(pos_tend) < 1.e-10 .or. abs(neg_tend) < 1.e-10) then
    119         write(*,*) "Reason of stopping: there is no sublimating or condensing h2o ice!"
    120         write(*,*) "Tendencies on ice sublimating =", neg_tend
    121         write(*,*) "Tendencies on ice increasing =", pos_tend
    122         write(*,*) "This can be due to the absence of h2o ice in the PCM run!"
    123         stopPEM = 2
    124     else
    125         ! We adapt the tendencies to conserve h2o and do only exchange between grid points
    126         if (neg_tend > pos_tend .and. pos_tend > 0.) then ! More sublimation on the planet than condensation
    127             where (d_h2oice < 0.) ! We lower the sublimating rate by a coefficient
    128                 new_tend = d_h2oice*pos_tend/neg_tend
    129             elsewhere ! We don't change the condensing rate
    130                 new_tend = d_h2oice
    131             end where
    132         else if (neg_tend < pos_tend .and. neg_tend > 0.) then ! More condensation on the planet than sublimation
    133             where (d_h2oice < 0.) ! We don't change the sublimating rate
    134                 new_tend = d_h2oice
    135             elsewhere ! We lower the condensing rate by a coefficient
    136                 new_tend = d_h2oice*neg_tend/pos_tend
    137             end where
    138         endif
    139     endif
     86if (ngrid == 1) then ! In 1D
     87    h2o_ice = h2o_ice + d_h2oice*dt
     88    zshift_surf = d_h2oice*dt/rho_h2oice
     89else ! In 3D
     90    call stopping_crit_h2o(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopPEM)
     91    if (stopPEM > 0) return
    14092
    141     ! Evolution of the h2o ice for each physical point
    142     h2o_ice = h2o_ice + new_tend*dt
    143 
    144     ! We compute the amount of h2o that is sublimated in excess
    145     negative_part = 0.
    146     do i = 1,ngrid
    147         do islope = 1,nslope
    148             if (h2o_ice(i,islope) < 0.) then
    149                 negative_part = negative_part - h2o_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
    150                 h2o_ice(i,islope) = 0.
    151                 d_h2oice(i,islope) = 0.
    152             endif
    153         enddo
    154     enddo
    155 
    156     ! We compute a coefficient by which we should remove the ice that has been added to places
    157     ! even if this ice was contributing to an unphysical negative amount of ice at other places
    158     if (abs(pos_tend) < 1.e-10) then
    159         real_coefficient = 0.
    160     else
    161         real_coefficient = negative_part/pos_tend
    162     endif
    163     ! In the place of accumulation of ice, we remove a bit of ice in order to conserve h2o
    164     do islope = 1,nslope
    165         do i = 1,ngrid
    166              if (new_tend(i,islope) > 0.) h2o_ice(i,islope) = h2o_ice(i,islope) - new_tend(i,islope)*real_coefficient*dt*cos(def_slope_mean(islope)*pi/180.)
    167         enddo
    168      enddo
    169 else ! ngrid == 1, i.e. in 1D
    170     h2o_ice = h2o_ice + d_h2oice*dt
     93    call balance_h2oice_reservoirs(ngrid,nslope,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
     94    h2o_ice = h2o_ice + d_h2oice_new*dt
     95    zshift_surf = d_h2oice_new*dt/rho_h2oice
    17196endif
    172 
    173 zshift_surf = d_h2oice*dt/rho_h2oice
    17497
    17598END SUBROUTINE evol_h2o_ice
    17699
     100!=======================================================================
     101
     102SUBROUTINE balance_h2oice_reservoirs(ngrid,nslope,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
     103
     104use time_evol_mod, only: dt
     105
     106implicit none
     107
     108!=======================================================================
     109!
     110! Routine to balance the H2O flux from and into the atmosphere accross reservoirs
     111!
     112!=======================================================================
     113! arguments:
     114! ----------
     115! INPUT
     116integer,                       intent(in) :: ngrid, nslope ! # of grid points, # of subslopes
     117real,                          intent(in) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! ! Variables to conserve H2O
     118real, dimension(ngrid,nslope), intent(in) :: h2o_ice ! H2O ice (kg/m^2)
     119
     120! OUTPUT
     121real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! Tendency of H2O ice (kg/m^2/year)
     122real, dimension(ngrid,nslope), intent(out) :: d_h2oice_new ! Adjusted tendencies to keep balance between donor and recipient reservoirs
     123
     124! local:
     125! ------
     126integer :: i, islope
     127real    :: S_target, S_target_subl_h2oice, S_target_cond_h2oice, S_ghostice, d_target
     128!=======================================================================
     129S_target = (S_atm_2_h2o + S_h2o_2_atm)/2.
     130S_target_cond_h2oice = S_atm_2_h2oice + S_target - S_atm_2_h2o
     131S_target_subl_h2oice = S_h2oice_2_atm + S_target - S_h2o_2_atm
     132
     133d_h2oice_new = 0.
     134S_ghostice = 0.
     135do i = 1,ngrid
     136    do islope = 1,nslope
     137        if (d_h2oice(i,islope) > 0.) then ! Condensing
     138            d_h2oice_new(i,islope) = d_h2oice_new(i,islope)*S_target_cond_h2oice/S_atm_2_h2oice
     139        else if (d_h2oice(i,islope) < 0.) then ! Sublimating
     140            d_target = d_h2oice(i,islope)*S_target_subl_h2oice/S_h2oice_2_atm
     141            if (abs(d_target*dt) < h2o_ice(i,islope)) then ! Enough ice to sublimate everything
     142                d_h2oice_new(i,islope) = d_target
     143            else ! Not enough ice to sublimate everything
     144                ! We sublimate what we can
     145                d_h2oice_new(i,islope) = h2o_ice(i,islope)/dt
     146                ! It means the tendency is zero next time
     147                d_h2oice(i,islope) = 0.
     148                ! We compute the amount of H2O ice that we could not make sublimate
     149                S_ghostice = S_ghostice + abs(d_target*dt) - h2o_ice(i,islope)
     150            endif
     151        endif
     152    enddo
     153enddo
     154
     155! We need to remove this ice unable to sublimate from places where ice condenseds in order to keep balance
     156where (d_h2oice_new > 0.) d_h2oice_new = d_h2oice_new*(S_target_cond_h2oice - S_ghostice)/S_target_cond_h2oice
     157
     158END SUBROUTINE balance_h2oice_reservoirs
     159
    177160END MODULE evol_ice_mod

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

    r3926 r3938  
    594594    if (h_toplag < h_patchy_dust) then ! But the dust lag is too thin to considered ice as subsurface ice
    595595        h_toplag = 0.
    596         h2o_ice = current%h_h2oice
     596        h2o_ice = current%h_h2oice*rho_h2oice
    597597    endif
    598598else if (current%h_co2ice > 0. .and. current%h_co2ice > current%h_h2oice) then ! No, there is more CO2 ice than H2O ice
    599599    h_toplag = 0. ! Because it matters only for H2O ice
    600     if (h_toplag < h_patchy_ice) co2_ice = current%h_co2ice ! But the dust lag is too thin to considered ice as subsurface ice
     600    if (h_toplag < h_patchy_ice) co2_ice = current%h_co2ice*rho_co2ice ! But the dust lag is too thin to considered ice as subsurface ice
    601601endif
    602602

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

    r3933 r3938  
    4141use glaciers_mod,               only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, &
    4242                                      metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice, computeTcondCO2
    43 use stopping_crit_mod,          only: stopping_crit_h2o_ice, stopping_crit_co2
     43use stopping_crit_mod,          only: stopping_crit_h2o_ice, stopping_crit_co2, stopping_crit_h2o
    4444use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
    45 use evol_ice_mod,               only: evol_co2_ice, evol_h2o_ice
     45use evol_ice_mod,               only: evol_co2_ice, evol_h2o_ice, balance_h2oice_reservoirs
    4646use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
    4747                                      TI_PEM,               & ! Soil thermal inertia
    … …  
    7070use co2condens_mod,             only: CO2cond_ps
    7171use layering_mod,               only: layering, del_layering, make_layering, layering_algo, subsurface_ice_layering, &
    72                                       ptrarray, stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str
     72                                      rho_co2ice, rho_h2oice, ptrarray,                                              &
     73                                      stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str
    7374use dyn_ss_ice_m_mod,           only: dyn_ss_ice_m
    7475use parse_args_mod,             only: parse_args
    … …  
    182183integer                               :: timelen              ! # time samples
    183184real                                  :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
     185real, dimension(:,:),    allocatable  :: d_h2oice_new         ! Adjusted tendencies to keep balance between donor and recipient
     186real                                  :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to balance H2O
    184187
    185188! Variables for co2 ice evolution
    … …  
    615618        do islope = 1,nslope
    616619            if (is_co2ice_str(layerings_map(ig,islope)%top)) then
    617                 co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
     620                co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice*rho_co2ice
    618621            else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
    619                 h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
     622                h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice*rho_h2oice
    620623            else
    621624                call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
    … …  
    782785    if (layering_algo) then
    783786        h2o_ice_depth_old = h2o_ice_depth
     787
     788        allocate(d_h2oice_new(ngrid,nslope))
     789        call stopping_crit_h2o(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopPEM)
     790        call balance_h2oice_reservoirs(ngrid,nslope,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
     791   
    784792        do islope = 1,nslope
    785793            do ig = 1,ngrid
    786                 call make_layering(layerings_map(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),new_str(ig,islope),zshift_surf(ig,islope),new_lag(ig,islope),zlag(ig,islope),current(ig,islope)%p)
     794                call make_layering(layerings_map(ig,islope),d_co2ice(ig,islope),d_h2oice_new(ig,islope),new_str(ig,islope),zshift_surf(ig,islope),new_lag(ig,islope),zlag(ig,islope),current(ig,islope)%p)
    787795                !call print_layering(layerings_map(ig,islope))
    788796                co2_ice(ig,islope) = 0.
    … …  
    790798                h2o_ice_depth(ig,islope) = 0.
    791799                if (is_co2ice_str(layerings_map(ig,islope)%top)) then
    792                     co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
     800                    co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice*rho_co2ice
    793801                else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
    794                     h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
     802                    h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice*rho_h2oice
    795803                else
    796804                    call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
    … …  
    798806            enddo
    799807        enddo
     808        deallocate(d_h2oice_new)
    800809    else
    801810        zlag = 0.
    … …  
    815824            do islope = 1,nslope
    816825                do ig = 1,ngrid
    817                     layerings_map(ig,islope)%top%h_co2ice = co2_ice(ig,islope)
     826                    layerings_map(ig,islope)%top%h_co2ice = co2_ice(ig,islope)/rho_co2ice
    818827                enddo
    819828            enddo
    … …  
    825834            do islope = 1,nslope
    826835                do ig = 1,ngrid
    827 !~                     layerings_map(ig,islope)%top%h_h2oice = h2o_ice(ig,islope)
     836                    layerings_map(ig,islope)%top%h_h2oice = h2o_ice(ig,islope)/rho_h2oice
    828837                enddo
    829838            enddo
    … …  
    968977        totmass_ini = max(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini,1.e-10)
    969978        write(*,'(a,f8.3,a)') " > Relative total CO2 mass balance = ", 100.*(totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini, ' %'
    970         if ((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini > 0.01) then
    971             write(*,*) '  /!\ Warning: mass balance is not conseved!'
     979        if (abs((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini) > 0.01) then
     980            write(*,*) '  /!\ Warning: mass balance is not conserved!'
    972981            totmass_ini = max(totmass_atmco2_ini,1.e-10)
    973             write(*,'(a,f8.3,a)') '       Atmospheric CO2 mass balance = ', 100.*(totmass_atmco2 - totmass_atmco2_ini)/totmass_ini, ' %'
     982            write(*,*) '       Atmospheric CO2 mass balance = ', 100.*(totmass_atmco2 - totmass_atmco2_ini)/totmass_ini, ' %'
    974983            totmass_ini = max(totmass_co2ice_ini,1.e-10)
    975             write(*,'(a,f8.3,a)') '       CO2 ice mass balance         = ', 100.*(totmass_co2ice - totmass_co2ice_ini)/totmass_ini, ' %'
     984            write(*,*) '       CO2 ice mass balance         = ', 100.*(totmass_co2ice - totmass_co2ice_ini)/totmass_ini, ' %'
    976985            totmass_ini = max(totmass_adsco2_ini,1.e-10)
    977             write(*,'(a,f8.3,a)') '       Adsorbed CO2 mass balance    = ', 100.*(totmass_adsco2 - totmass_adsco2_ini)/totmass_ini, ' %'
     986            write(*,*) '       Adsorbed CO2 mass balance    = ', 100.*(totmass_adsco2 - totmass_adsco2_ini)/totmass_ini, ' %'
    978987        endif
    979988    endif

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

    r3770 r3938  
    149149END SUBROUTINE stopping_crit_co2
    150150
     151!=======================================================================
     152
     153SUBROUTINE stopping_crit_h2o(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopPEM)
     154
     155use time_evol_mod, only: dt
     156use comslope_mod,  only: subslope_dist, def_slope_mean
     157#ifndef CPP_STD
     158    use comcstfi_h,   only: pi
     159#else
     160    use comcstfi_mod, only: pi
     161#endif
     162
     163implicit none
     164
     165!=======================================================================
     166!
     167! Routine to check if the h2o is only exchanged between grid points
     168!
     169!=======================================================================
     170! Inputs
     171!-------
     172integer,                       intent(in) :: ngrid, nslope            ! # of grid points, # of subslopes
     173real, dimension(ngrid),        intent(in) :: cell_area                ! Area of each mesh grid (m^2)
     174real, dimension(ngrid),        intent(in) :: delta_h2o_adsorbed       ! Mass of H2O adsorbed/desorbded in the soil (kg/m^2)
     175real, dimension(ngrid),        intent(in) :: delta_h2o_icetablesublim ! Mass of H2O that condensed/sublimated at the ice table
     176real, dimension(ngrid,nslope), intent(in) :: h2o_ice                  ! H2O ice (kg/m^2)
     177real, dimension(ngrid,nslope), intent(in) :: d_h2oice                 ! Tendency of H2O ice (kg/m^2/year)
     178! Outputs
     179!--------
     180integer, intent(inout) :: stopPEM ! Stopping criterion code
     181real,    intent(out) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to conserve H2O
     182
     183! Locals
     184! ------
     185integer :: i, islope ! Loop
     186
     187!=======================================================================
     188if (stopPEM > 0) return
     189
     190! We compute the amount of water going out from the atmosphere (S_atm_2_h2o) and going into the atmophere (S_h2o_2_atm)
     191S_atm_2_h2o = 0.
     192S_h2o_2_atm = 0.
     193S_atm_2_h2oice = 0.
     194S_h2oice_2_atm = 0.
     195do i = 1,ngrid
     196    if (delta_h2o_adsorbed(i) > 0.) then
     197        S_atm_2_h2o = S_atm_2_h2o + delta_h2o_adsorbed(i)*cell_area(i)
     198    else
     199        S_h2o_2_atm = S_h2o_2_atm + delta_h2o_adsorbed(i)*cell_area(i)
     200    endif
     201    if (delta_h2o_icetablesublim(i) > 0.) then
     202        S_atm_2_h2o = S_atm_2_h2o + delta_h2o_icetablesublim(i)*cell_area(i)
     203    else
     204        S_h2o_2_atm = S_h2o_2_atm + delta_h2o_icetablesublim(i)*cell_area(i)
     205    endif
     206    do islope = 1,nslope
     207        if (d_h2oice(i,islope) > 0.) then
     208            S_atm_2_h2o = S_atm_2_h2o + d_h2oice(i,islope)*dt*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
     209            S_atm_2_h2oice = S_atm_2_h2oice + d_h2oice(i,islope)*dt*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
     210        else if (d_h2oice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then
     211            S_h2o_2_atm = S_h2o_2_atm - d_h2oice(i,islope)*dt*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
     212            S_h2oice_2_atm = S_h2oice_2_atm - d_h2oice(i,islope)*dt*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
     213        endif
     214    enddo
     215enddo
     216
     217! Since relative atmospheric water is kept constant, we need to equate condensing reservoirs to the sublimating ones.
     218! It is not possible if one of them is missing.
     219if (abs(S_atm_2_h2oice) < 1.e-10 .or. abs(S_h2oice_2_atm) < 1.e-10) then
     220    write(*,*) "Reason of stopping: there is no sublimating or condensing h2o ice!"
     221    write(*,*) "This can be due to the absence of h2o ice in the PCM run."
     222    write(*,*) "Amount of condensing ice  =", S_atm_2_h2oice
     223    write(*,*) "Amount of sublimating ice =", S_h2oice_2_atm
     224    stopPEM = 2
     225endif
     226
     227END SUBROUTINE stopping_crit_h2o
     228
    151229END MODULE
Note: See TracChangeset for help on using the changeset viewer.