Ignore:
Timestamp:
Feb 12, 2026, 9:09:12 AM (2 weeks ago)
Author:
jbclement
Message:

PEM:
Major refactor following the previous ones (r3989 and r3991) completing the large structural reorganization and cleanup of the PEM codebase. This revision introduces newly designed modules, standardizes interfaces with explicit ini/end APIs and adds native NetCDF I/O together with explicit PCM/PEM adapters. In detail:

  • Some PEM models were corrected or improved:
    • Frost/perennial ice semantics are clarified via renaming;
    • Soil temperature remapping clarified, notably by removing the rescaling of temperature deviation;
    • Geothermal flux for the PCM is computed based on the PEM state;
  • New explicit PEM/PCM adapters ("set_*"/"build4PCM_*") to decouple PEM internal representation from PCM file layouts and reconstruct consistent fields returned to the PCM;
  • New explicit build/teardown routines that centralize allocation and initialization ordering, reducing accidental use of uninitialized data and making the model lifecycle explicit;
  • Add native read/write helpers for NetCDF that centralize all low-level NetCDF interactions with major improvements (and more simplicity) compared to legacy PEM/PCM I/O (see the modules "io_netcdf" and "output"). They support reading, creation and writing of "diagevol.nc" (renamed from "diagpem.nc") and start/restart files;
  • Provide well-focused modules ("numerics"/"maths"/"utility"/"display") to host commonly-used primitives:
    • "numerics" defines numerical types and constants for reproducibility, portability across compilers and future transitions (e.g. quadruple precision experiments);
    • "display" provides a single controlled interface for runtime messages, status output and diagnostics, avoiding direct 'print'/'write' to enable silent mode, log redirection, and MPI-safe output in the future.
    • "utility" (new module) hosts generic helpers used throughout the code (e.g. "int2str" or "real2str");
  • Add modules "clim_state_init"/"clim_state_rec" which provide robust read/write logic for "start/startfi/startpem", including 1D fallbacks, mesh conversions and dimension checks. NetCDF file creation is centralized and explicit. Restart files are now self-consistent and future-proof, requiring changes only to affected variables;
  • Add module "atmosphere" which computes pressure fields, reconstructs potential temperature and air mass. It also holds the whole logic to define sigma or hybrid coordinates for altitudes;
  • Add module "geometry" to centrilize dimensions logic and grid conversions routines (including 2 new ones "dyngrd2vect"/"vect2dyngrd");
  • Add module "slopes" to isolate slopes handling;
  • Add module "surface" to isolate surface management. Notably, albedo and emissivity are now fully reconstructed following the PCM settings;
  • Add module "allocation" to check the dimension initialization and centrilize allocation/deallocation;
  • Finalize module decomposition and renaming to consolidate domain-based modules, purpose-based routines and physics/process-based variables;
  • The main program now drives a clearer sequence of conceptual steps (initialization / reading / evolution / update / build / writing) and fails explicitly instead of silently defaulting;
  • Ice table logic is made restart-safe;
  • 'Open'/'read' intrinsic logic is made safe and automatic;
  • Improve discoverability and standardize the data handling (private vs protected vs public);
  • Apply consistent documentation/header style already introduced;
  • Update deftank scripts to reflect new names and launch wrappers;

This revision is a structural milestone aiming to be behavior-preserving where possible. It has been tested via compilation and short integration runs. However, due to extensive renames, moves, and API changes, full validation is still ongoing.
Note: the revision includes one (possibly two) easter egg hidden in the code for future archaeologists of the PEM. No physics were harmed.
JBC

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r3991 r4065  
    1616!-----------------------------------------------------------------------
    1717
     18! DEPENDENCIES
     19! ------------
     20use numerics, only: dp, di, k4
     21
    1822! DECLARATION
    1923! -----------
    … …  
    2226! PARAMETERS
    2327! ----------
    24 logical :: h2oice_flow ! True by default, flag to compute H2O ice flow
    25 logical :: co2ice_flow ! True by default, flag to compute CO2 ice flow
     28logical(k4), protected :: h2oice_flow ! Flag to compute H2O ice flow
     29logical(k4), protected :: co2ice_flow ! Flag to compute CO2 ice flow
    2630
    2731contains
    28 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    29 
    30 !=======================================================================
    31 SUBROUTINE flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,co2ice,flag_co2flow)
     32!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     33
     34!=======================================================================
     35SUBROUTINE set_glaciers_config(h2oice_flow_in,co2ice_flow_in)
     36!-----------------------------------------------------------------------
     37! NAME
     38!     set_glaciers_config
     39!
     40! DESCRIPTION
     41!     Setter for 'glaciers' configuration parameters.
     42!
     43! AUTHORS & DATE
     44!     JB Clement, 02/2026
     45!
     46! NOTES
     47!
     48!-----------------------------------------------------------------------
     49
     50! DEPENDENCIES
     51! ------------
     52use utility,  only: bool2str
     53use geometry, only: nslope
     54use display,  only: print_msg
     55
     56! DECLARATION
     57! -----------
     58implicit none
     59
     60! ARGUMENTS
     61! ---------
     62logical(k4), intent(in) :: h2oice_flow_in, co2ice_flow_in
     63
     64! CODE
     65! ----
     66h2oice_flow = h2oice_flow_in
     67co2ice_flow = co2ice_flow_in
     68if (h2oice_flow .and. nslope == 1) then
     69    call print_msg('Warning: h2oice_flow = .true. but nslope = 1. So there will be no flow!')
     70    h2oice_flow = .false.
     71end if
     72if (co2ice_flow .and. nslope == 1) then
     73    call print_msg('Warning: co2ice_flow = .true. but nslope = 1. So there will be no flow!')
     74    co2ice_flow = .false.
     75end if
     76call print_msg('h2oice_flow = '//bool2str(h2oice_flow))
     77call print_msg('co2ice_flow = '//bool2str(co2ice_flow))
     78
     79END SUBROUTINE set_glaciers_config
     80!=======================================================================
     81
     82!=======================================================================
     83SUBROUTINE flow_co2glaciers(vmr_co2_PEM,ps_PCM,ps_avg_glob_ini,ps_avg_global,co2ice,is_co2ice_flow)
    3284!-----------------------------------------------------------------------
    3385! NAME
    … …  
    4597!-----------------------------------------------------------------------
    4698
    47 ! DECLARATION
    48 ! -----------
    49 implicit none
    50 
    51 ! ARGUMENTS
    52 ! ---------
    53 integer,                                  intent(in)    :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
    54 real, dimension(ngrid,nslope),            intent(in)    :: subslope_dist                 ! Grid points x Slopes: Distribution of the subgrid slopes
    55 real, dimension(ngrid),                   intent(in)    :: def_slope_mean                ! Grid points: values of the sub grid slope angles
    56 real, dimension(ngrid,timelen),           intent(in)    :: vmr_co2_PEM                   ! Grid points x Time field : VMR of CO2 in the first layer [mol/mol]
    57 real, dimension(ngrid,timelen),           intent(in)    :: ps_PCM                        ! Grid points x Time field: surface pressure given by the PCM [Pa]
    58 real,                                     intent(in)    :: ps_avg_global_ini             ! Global averaged surface pressure at the beginning [Pa]
    59 real,                                     intent(in)    :: ps_avg_global                 ! Global averaged surface pressure during the PEM iteration [Pa]
    60 real, dimension(ngrid,nslope),            intent(inout) :: co2ice                        ! Grid points x Slope field: CO2 ice on the subgrid slopes [kg/m^2]
    61 integer(kind=1), dimension(ngrid,nslope), intent(out)   :: flag_co2flow                  ! Flag to see if there is flow on the subgrid slopes
    62 
    63 ! LOCAL VARIABLES
    64 ! ---------------
    65 real, dimension(ngrid,nslope) :: Tcond ! CO2 condensation temperature [K]
    66 real, dimension(ngrid,nslope) :: hmax  ! Maximum thickness before flow
    67 
    68 ! CODE
    69 ! ----
    70 write(*,*) "> Flow of CO2 glaciers"
    71 call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,Tcond)
    72 call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
    73 call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tcond,co2ice,flag_co2flow)
     99! DEPENDENCIES
     100! ------------
     101use geometry, only: ngrid, nslope
     102use display,  only: print_msg
     103
     104! DECLARATION
     105! -----------
     106implicit none
     107
     108! ARGUMENTS
     109! ---------
     110real(dp),    dimension(:,:),            intent(in)    :: vmr_co2_PEM     ! Grid points x Time field : VMR of CO2 in the first layer [mol/mol]
     111real(dp),    dimension(:,:),            intent(in)    :: ps_PCM          ! Grid points x Time field: surface pressure given by the PCM [Pa]
     112real(dp),                               intent(in)    :: ps_avg_glob_ini ! Global averaged surface pressure at the beginning [Pa]
     113real(dp),                               intent(in)    :: ps_avg_global   ! Global averaged surface pressure during the PEM iteration [Pa]
     114real(dp),    dimension(:,:),            intent(inout) :: co2ice          ! Grid points x Slope field: CO2 ice on the subgrid slopes [kg/m^2]
     115logical(k4), dimension(:,:),            intent(out)   :: is_co2ice_flow  ! Flag to see if there is flow on the subgrid slopes
     116
     117! LOCAL VARIABLES
     118! ---------------
     119real(dp), dimension(ngrid,nslope) :: Tcond ! CO2 condensation temperature [K]
     120real(dp), dimension(ngrid,nslope) :: hmax  ! Maximum thickness before flow
     121
     122! CODE
     123! ----
     124call print_msg("> Flow of CO2 glaciers")
     125call computeTcondCO2(vmr_co2_PEM,ps_PCM,ps_avg_glob_ini,ps_avg_global,Tcond)
     126call compute_hmaxglaciers(Tcond,"co2",hmax)
     127call transfer_ice_duringflow(hmax,Tcond,co2ice,is_co2ice_flow)
    74128
    75129END SUBROUTINE flow_co2glaciers
    … …  
    77131
    78132!=======================================================================
    79 SUBROUTINE flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow)
     133SUBROUTINE flow_h2oglaciers(Tice,h2oice,is_h2oice_flow)
    80134!-----------------------------------------------------------------------
    81135! NAME
    … …  
    93147!-----------------------------------------------------------------------
    94148
    95 ! DECLARATION
    96 ! -----------
    97 implicit none
    98 
    99 ! ARGUMENTS
    100 ! ---------
    101 integer,                                  intent(in)    ::  ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
    102 real, dimension(ngrid,nslope),            intent(in)    ::  subslope_dist        ! Distribution of the subgrid slopes
    103 real, dimension(ngrid),                   intent(in)    ::  def_slope_mean       ! Values of the sub grid slope angles
    104 real, dimension(ngrid,nslope),            intent(in)    ::  Tice                 ! Ice Temperature [K]
    105 real, dimension(ngrid,nslope),            intent(inout) ::  h2oice               ! H2O ice on the subgrid slopes [kg/m^2]
    106 integer(kind=1), dimension(ngrid,nslope), intent(out)   ::  flag_h2oflow         ! Flow flag on subgrid slopes
    107 
    108 ! LOCAL VARIABLES
    109 ! ---------------
    110 real, dimension(ngrid,nslope) :: hmax ! Maximum thickness before flow
    111 
    112 ! CODE
    113 ! ----
    114 write(*,*) "> Flow of H2O glaciers"
    115 call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
    116 call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,h2oice,flag_h2oflow)
     149! DEPENDENCIES
     150! ------------
     151use geometry, only: ngrid, nslope
     152use display,  only: print_msg
     153
     154! DECLARATION
     155! -----------
     156implicit none
     157
     158! ARGUMENTS
     159! ---------
     160real(dp),    dimension(:,:), intent(in)    ::  Tice           ! Ice Temperature [K]
     161real(dp),    dimension(:,:), intent(inout) ::  h2oice         ! H2O ice on the subgrid slopes [kg/m^2]
     162logical(k4), dimension(:,:), intent(out)   ::  is_h2oice_flow ! Flow flag on subgrid slopes
     163
     164! LOCAL VARIABLES
     165! ---------------
     166real(dp), dimension(ngrid,nslope) :: hmax ! Maximum thickness before flow
     167
     168! CODE
     169! ----
     170call print_msg("> Flow of H2O glaciers")
     171call compute_hmaxglaciers(Tice,"h2o",hmax)
     172call transfer_ice_duringflow(hmax,Tice,h2oice,is_h2oice_flow)
    117173
    118174END SUBROUTINE flow_h2oglaciers
    … …  
    120176
    121177!=======================================================================
    122 SUBROUTINE compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
     178SUBROUTINE compute_hmaxglaciers(Tice,name_ice,hmax)
    123179!-----------------------------------------------------------------------
    124180! NAME
    … …  
    140196! DEPENDENCIES
    141197! ------------
    142 use ice_table,      only: rho_ice
    143 use stop_pem,       only: stop_clean
    144 use phys_constants, only: pi, g
    145 
    146 ! DECLARATION
    147 ! -----------
    148 implicit none
    149 
    150 ! ARGUMENTS
    151 ! ---------
    152 integer,                       intent(in)  :: ngrid, nslope  ! # of grid points and subslopes
    153 integer,                       intent(in)  :: iflat          ! index of the flat subslope
    154 real, dimension(nslope),       intent(in)  :: def_slope_mean ! Values of the subgrid slope angles [deg]
    155 real, dimension(ngrid,nslope), intent(in)  :: Tice           ! Ice temperature [K]
    156 character(3),                  intent(in)  :: name_ice      ! Nature of ice
    157 real, dimension(ngrid,nslope), intent(out) :: hmax           ! Maximum thickness before flow [m]
    158 
    159 ! LOCAL VARIABLES
    160 ! ---------------
    161 real    :: tau_d ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
    162 integer :: ig, islope
    163 real    :: slo_angle
     198use geometry,  only: ngrid, nslope
     199use slopes,    only: iflat, def_slope_mean
     200use ice_table, only: rho_ice
     201use stoppage,  only: stop_clean
     202use maths,     only: pi
     203use physics,   only: g
     204
     205! DECLARATION
     206! -----------
     207implicit none
     208
     209! ARGUMENTS
     210! ---------
     211real(dp), dimension(:,:), intent(in)  :: Tice     ! Ice temperature [K]
     212character(3),             intent(in)  :: name_ice ! Nature of ice
     213real(dp), dimension(:,:), intent(out) :: hmax     ! Maximum thickness before flow [m]
     214
     215! LOCAL VARIABLES
     216! ---------------
     217real(dp)    :: tau_d ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
     218integer(di) :: ig, islope
     219real(dp)    :: slo_angle
    164220
    165221! CODE
    … …  
    167223select case (trim(adjustl(name_ice)))
    168224    case('h2o')
    169         tau_d = 1.e5
     225        tau_d = 1.e5_dp
    170226    case('co2')
    171         tau_d = 5.e3
     227        tau_d = 5.e3_dp
    172228    case default
    173         call stop_clean("compute_hmaxglaciers","Type of ice not known!",1)
     229        call stop_clean(__FILE__,__LINE__,"type of ice unknown!",1)
    174230end select
    175231
    … …  
    177233    do islope = 1,nslope
    178234        if (islope == iflat) then
    179             hmax(ig,islope) = 1.e8
     235            hmax(ig,islope) = 1.e8_dp
    180236        else
    181             slo_angle = abs(def_slope_mean(islope)*pi/180.)
     237            slo_angle = abs(def_slope_mean(islope)*pi/180._dp)
    182238            hmax(ig,islope) = tau_d/(rho_ice(Tice(ig,islope),name_ice)*g*slo_angle)
    183         endif
    184     enddo
    185 enddo
     239        end if
     240    end do
     241end do
    186242
    187243END SUBROUTINE compute_hmaxglaciers
    … …  
    189245
    190246!=======================================================================
    191 SUBROUTINE transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,qice,flag_flow)
     247SUBROUTINE transfer_ice_duringflow(hmax,Tice,qice,flag_flow)
    192248!-----------------------------------------------------------------------
    193249! NAME
    … …  
    208264! DEPENDENCIES
    209265! ------------
    210 use ice_table,      only: rho_ice
    211 use stop_pem,       only: stop_clean
    212 use phys_constants, only: pi
    213 
    214 ! DECLARATION
    215 ! -----------
    216 implicit none
    217 
    218 ! ARGUMENTS
    219 ! ---------
    220 integer,                                  intent(in)    :: ngrid, nslope  ! # of physical points and subslope
    221 integer,                                  intent(in)    :: iflat          ! Index of the flat subslope
    222 real, dimension(ngrid,nslope),            intent(in)    :: subslope_dist  ! Distribution of the subgrid slopes
    223 real, dimension(nslope),                  intent(in)    :: def_slope_mean ! Values of the subgrid slopes
    224 real, dimension(ngrid,nslope),            intent(in)    :: hmax           ! Maximum height before initiating flow [m]
    225 real, dimension(ngrid,nslope),            intent(in)    :: Tice           ! Ice temperature [K]
    226 real, dimension(ngrid,nslope),            intent(inout) :: qice           ! Ice in the subslope [kg/m^2]
    227 integer(kind=1), dimension(ngrid,nslope), intent(out)   :: flag_flow      ! Flow flag on subgrid slopes
    228 
    229 ! LOCAL VARIABLES
    230 ! ---------------
    231 integer :: ig, islope
    232 integer :: iaval ! index where ice is transferred
    233 
    234 ! CODE
    235 ! ----
    236 flag_flow = 0
     266use geometry,  only: ngrid, nslope
     267use slopes,    only: iflat, subslope_dist, def_slope_mean
     268use ice_table, only: rho_ice
     269use stoppage,  only: stop_clean
     270use maths,     only: pi
     271
     272! DECLARATION
     273! -----------
     274implicit none
     275
     276! ARGUMENTS
     277! ---------
     278real(dp),    dimension(:,:), intent(in)    :: hmax      ! Maximum height before initiating flow [m]
     279real(dp),    dimension(:,:), intent(in)    :: Tice      ! Ice temperature [K]
     280real(dp),    dimension(:,:), intent(inout) :: qice      ! Ice in the subslope [kg/m^2]
     281logical(k4), dimension(:,:), intent(out)   :: flag_flow ! Flow flag on subgrid slopes
     282
     283! LOCAL VARIABLES
     284! ---------------
     285integer(di) :: ig, islope
     286integer(di) :: iaval ! Index where ice is transferred
     287
     288! CODE
     289! ----
     290flag_flow = .false.
    237291
    238292do ig = 1,ngrid
    … …  
    240294        if (islope /= iflat) then ! ice can be infinite on flat ground
    241295! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely  on flat ground
    242             if (qice(ig,islope) >= rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)) then
     296            if (qice(ig,islope) >= rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180._dp)) then
    243297! Second: determine the flatest slopes possible:
    244298                if (islope > iflat) then
    … …  
    246300                else
    247301                    iaval = islope + 1
    248                 endif
     302                end if
    249303                do while (iaval /= iflat .and. subslope_dist(ig,iaval) == 0)
    250304                    if (iaval > iflat) then
    … …  
    252306                    else
    253307                        iaval = iaval + 1
    254                     endif
    255                 enddo
    256                 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.)) &
    257                                *subslope_dist(ig,islope)/subslope_dist(ig,iaval)*cos(pi*def_slope_mean(iaval)/180.)/cos(pi*def_slope_mean(islope)/180.)
    258 
    259                 qice(ig,islope) = rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)
    260                 flag_flow(ig,islope) = 1
    261             endif ! co2ice > hmax
    262         endif ! iflat
    263     enddo !islope
    264 enddo !ig
     308                    end if
     309                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)) &
     311                               *subslope_dist(ig,islope)/subslope_dist(ig,iaval)*cos(pi*def_slope_mean(iaval)/180._dp)/cos(pi*def_slope_mean(islope)/180._dp)
     312
     313                qice(ig,islope) = rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180._dp)
     314                flag_flow(ig,islope) = .true.
     315            end if ! co2ice > hmax
     316        end if ! iflat
     317    end do !islope
     318end do !ig
    265319
    266320END SUBROUTINE
    267321
    268322!=======================================================================
    269 SUBROUTINE computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,Tcond)
     323SUBROUTINE computeTcondCO2(vmr_co2_PEM,ps_PCM,ps_avg_glob_ini,ps_avg_global,Tcond)
    270324!-----------------------------------------------------------------------
    271325! NAME
    … …  
    285339! DEPENDENCIES
    286340! ------------
    287 use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2
    288 
    289 ! DECLARATION
    290 ! -----------
    291 implicit none
    292 
    293 ! ARGUMENTS
    294 ! ---------
    295 integer,                        intent(in)  :: timelen, ngrid, nslope
    296 real, dimension(ngrid,timelen), intent(in)  :: vmr_co2_PEM       ! VMR of CO2 in the first layer [mol/mol]
    297 real, dimension(ngrid,timelen), intent(in)  :: ps_PCM            ! Surface pressure in the PCM [Pa]
    298 real,                           intent(in)  :: ps_avg_global_ini ! Global averaged surfacepressure in the PCM [Pa]
    299 real,                           intent(in)  :: ps_avg_global     ! Global averaged surface pressure computed during the PEM iteration
    300 real, dimension(ngrid,nslope), intent(out) :: Tcond             ! Condensation temperature of CO2, yearly averaged
    301 
    302 ! LOCAL VARIABLES
    303 ! ---------------
    304 integer :: ig, it
     341use geometry, only: ngrid, nday
     342use planet,   only: alpha_clap_co2, beta_clap_co2
     343
     344! DECLARATION
     345! -----------
     346implicit none
     347
     348! ARGUMENTS
     349! ---------
     350real(dp), dimension(:,:), intent(in)  :: vmr_co2_PEM       ! VMR of CO2 in the first layer [mol/mol]
     351real(dp), dimension(:,:), intent(in)  :: ps_PCM            ! Surface pressure in the PCM [Pa]
     352real(dp),                 intent(in)  :: ps_avg_glob_ini ! Global averaged surfacepressure in the PCM [Pa]
     353real(dp),                 intent(in)  :: ps_avg_global     ! Global averaged surface pressure computed during the PEM iteration
     354real(dp), dimension(:,:), intent(out) :: Tcond             ! Condensation temperature of CO2, yearly averaged
     355
     356! LOCAL VARIABLES
     357! ---------------
     358integer(di) :: ig, it
    305359
    306360! CODE
    307361! ----
    308362do ig = 1,ngrid
    309     Tcond(ig,:) = 0
    310     do it = 1,timelen
    311         Tcond(ig,:) = Tcond(ig,:) + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM(ig,it)*ps_PCM(ig,it)*ps_avg_global_ini/ps_avg_global/100))
    312     enddo
    313     Tcond(ig,:) = Tcond(ig,:)/timelen
    314 enddo
     363    Tcond(ig,:) = 0._dp
     364    do it = 1,nday
     365        Tcond(ig,:) = Tcond(ig,:) + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM(ig,it)*ps_PCM(ig,it)*ps_avg_glob_ini/ps_avg_global/100._dp))
     366    end do
     367    Tcond(ig,:) = Tcond(ig,:)/nday
     368end do
    315369
    316370END SUBROUTINE computeTcondCO2
Note: See TracChangeset for help on using the changeset viewer.