Changeset 3143 for trunk/LMDZ.COMMON


Ignore:
Timestamp:
Nov 29, 2023, 9:42:08 AM (12 months ago)
Author:
jbclement
Message:

PEM:
'Watercap' has been removed from the water ice evolution since 'water_reservoir' does already the job + Some cleanings to simplify the code.

/!\ Commit for the PEM management of h2o ice before a rework of ice management in the PEM!
JBC

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

Legend:

Unmodified
Added
Removed

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

    r3136 r3143  
    145145== 22/11/2023 == JBC
    146146Update of files in the deftank: addition of variables in the xml definition files, inclusion of "callphys.def" into "run_PEM.def" and minor typo in "launch_pem.sh".
     147
     148== 29/11/2023 == JBC
     149'Watercap' has been removed from the water ice evolution since 'water_reservoir' does already the job + Some cleanings to simplify the code.
     150
     151/!\ Commit for the PEM management of h2o ice before a rework of ice management in the PEM!

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

    r3123 r3143  
    1 module comsoil_h_PEM
     1MODULE comsoil_h_PEM
    22
    33implicit none
    4   integer, parameter :: nsoilmx_PEM = 68                 ! number of layers in the PEM
    5   real,save,allocatable,dimension(:) :: layer_PEM        ! soil layer depths   [m]
    6   real,save,allocatable,dimension(:) :: mlayer_PEM       ! soil mid-layer depths [m]
    7   real,save,allocatable,dimension(:,:,:) :: TI_PEM       ! soil thermal inertia [SI]
    8   real,save,allocatable,dimension(:,:) :: inertiedat_PEM ! soil thermal inertia saved as reference for current climate [SI]
    9   ! variables (FC: built in firstcall in soil.F)
    10   REAL,SAVE,ALLOCATABLE :: tsoil_PEM(:,:,:)              ! sub-surface temperatures [K]
    11   real,save,allocatable :: mthermdiff_PEM(:,:)           ! (FC) mid-layer thermal diffusivity [SI]
    12   real,save,allocatable :: thermdiff_PEM(:,:)            ! (FC) inter-layer thermal diffusivity [SI]
    13   real,save,allocatable :: coefq_PEM(:)                  ! (FC) q_{k+1/2} coefficients [SI]
    14   real,save,allocatable :: coefd_PEM(:,:)                ! (FC) d_k coefficients [SI]
    15   real,save,allocatable :: alph_PEM(:,:)                 ! (FC) alpha_k coefficients [SI]
    16   real,save,allocatable :: beta_PEM(:,:)                 ! beta_k coefficients [SI]
    17   real,save :: mu_PEM                                    ! mu coefficient [SI]
    18   real,save :: fluxgeo                                   ! Geothermal flux [W/m^2]
    19   real,save :: depth_breccia                             ! Depth at which we have breccia [m]
    20   real,save :: depth_bedrock                             ! Depth at which we have bedrock [m]
    21   integer,save :: index_breccia                          ! last index of the depth grid before having breccia
    22   integer,save :: index_bedrock                          ! last index of the depth grid before having bedrock
    23   LOGICAL soil_pem  ! True by default, to run with the subsurface physic. Read in pem.def
    24   real,save,allocatable,dimension(:) :: water_reservoir  ! Large reserve of water   [kg/m^2]
    25   real,save :: water_reservoir_nom                       ! Nominal value for the large reserve of water   [kg/m^2]
    26   logical, save :: reg_thprop_dependp                    ! thermal properites of the regolith vary with the surface pressure
    274
     5integer, parameter                        :: nsoilmx_PEM = 68    ! number of layers in the PEM
     6real, allocatable, dimension(:),     save :: layer_PEM           ! soil layer depths [m]
     7real, allocatable, dimension(:),     save :: mlayer_PEM          ! soil mid-layer depths [m]
     8real, allocatable, dimension(:,:,:), save :: TI_PEM              ! soil thermal inertia [SI]
     9real, allocatable, dimension(:,:),   save :: inertiedat_PEM      ! soil thermal inertia saved as reference for current climate [SI]
     10! Variables (FC: built in firstcall in soil.F)
     11real, allocatable, dimension(:,:,:), save :: tsoil_PEM           ! sub-surface temperatures [K]
     12real, allocatable, dimension(:,:),   save :: mthermdiff_PEM      ! (FC) mid-layer thermal diffusivity [SI]
     13real, allocatable, dimension(:,:),   save :: thermdiff_PEM       ! (FC) inter-layer thermal diffusivity [SI]
     14real, allocatable, dimension(:),     save :: coefq_PEM           ! (FC) q_{k+1/2} coefficients [SI]
     15real, allocatable, dimension(:,:),   save :: coefd_PEM           ! (FC) d_k coefficients [SI]
     16real, allocatable, dimension(:,:),   save :: alph_PEM            ! (FC) alpha_k coefficients [SI]
     17real, allocatable, dimension(:,:),   save :: beta_PEM            ! beta_k coefficients [SI]
     18real,                                save :: mu_PEM              ! mu coefficient [SI]
     19real,                                save :: fluxgeo             ! Geothermal flux [W/m^2]
     20real,                                save :: depth_breccia       ! Depth at which we have breccia [m]
     21real,                                save :: depth_bedrock       ! Depth at which we have bedrock [m]
     22integer,                             save :: index_breccia       ! last index of the depth grid before having breccia
     23integer,                             save :: index_bedrock       ! last index of the depth grid before having bedrock
     24logical                                   :: soil_pem            ! True by default, to run with the subsurface physic. Read in pem.def
     25real, allocatable, dimension(:),     save :: water_reservoir     ! Large reserve of water [kg/m^2]
     26real,                                save :: water_reservoir_nom ! Nominal value for the large reserve of water [kg/m^2]
     27logical,                             save :: reg_thprop_dependp  ! thermal properites of the regolith vary with the surface pressure
     28
     29!=======================================================================
    2830contains
     31!=======================================================================
    2932
    30   subroutine ini_comsoil_h_PEM(ngrid,nslope)
    31  
    32   implicit none
    33   integer,intent(in) :: ngrid ! number of atmospheric columns
    34   integer,intent(in) :: nslope ! number of slope within a mesh
     33SUBROUTINE ini_comsoil_h_PEM(ngrid,nslope)
    3534
    36     allocate(layer_PEM(nsoilmx_PEM))
    37     allocate(mlayer_PEM(0:nsoilmx_PEM-1))
    38     allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope))
    39     allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope))
    40     allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM-1))
    41     allocate(thermdiff_PEM(ngrid,nsoilmx_PEM-1))
    42     allocate(coefq_PEM(0:nsoilmx_PEM-1))
    43     allocate(coefd_PEM(ngrid,nsoilmx_PEM-1))
    44     allocate(alph_PEM(ngrid,nsoilmx_PEM-1))
    45     allocate(beta_PEM(ngrid,nsoilmx_PEM-1))
    46     allocate(inertiedat_PEM(ngrid,nsoilmx_PEM))
    47     allocate(water_reservoir(ngrid))
    48   end subroutine ini_comsoil_h_PEM
     35implicit none
    4936
     37integer, intent(in) :: ngrid  ! number of atmospheric columns
     38integer, intent(in) :: nslope ! number of slope within a mesh
    5039
    51   subroutine end_comsoil_h_PEM
     40allocate(layer_PEM(nsoilmx_PEM))
     41allocate(mlayer_PEM(0:nsoilmx_PEM-1))
     42allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope))
     43allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope))
     44allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM-1))
     45allocate(thermdiff_PEM(ngrid,nsoilmx_PEM-1))
     46allocate(coefq_PEM(0:nsoilmx_PEM-1))
     47allocate(coefd_PEM(ngrid,nsoilmx_PEM-1))
     48allocate(alph_PEM(ngrid,nsoilmx_PEM-1))
     49allocate(beta_PEM(ngrid,nsoilmx_PEM-1))
     50allocate(inertiedat_PEM(ngrid,nsoilmx_PEM))
     51allocate(water_reservoir(ngrid))
    5252
    53   implicit none
     53END SUBROUTINE ini_comsoil_h_PEM
    5454
    55     if (allocated(layer_PEM)) deallocate(layer_PEM)
    56     if (allocated(mlayer_PEM)) deallocate(mlayer_PEM)
    57     if (allocated(TI_PEM)) deallocate(TI_PEM)
    58     if (allocated(tsoil_PEM)) deallocate(tsoil_PEM)
    59     if (allocated(mthermdiff_PEM)) deallocate(mthermdiff_PEM)
    60     if (allocated(thermdiff_PEM)) deallocate(thermdiff_PEM)
    61     if (allocated(coefq_PEM)) deallocate(coefq_PEM)
    62     if (allocated(coefd_PEM)) deallocate(coefd_PEM)
    63     if (allocated(alph_PEM)) deallocate(alph_PEM)
    64     if (allocated(beta_PEM)) deallocate(beta_PEM)
    65     if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)
    66     if (allocated(water_reservoir)) deallocate(water_reservoir)
    67   end subroutine end_comsoil_h_PEM
     55!=======================================================================
    6856
    69 end module comsoil_h_PEM
     57SUBROUTINE end_comsoil_h_PEM
     58
     59implicit none
     60
     61if (allocated(layer_PEM)) deallocate(layer_PEM)
     62if (allocated(mlayer_PEM)) deallocate(mlayer_PEM)
     63if (allocated(TI_PEM)) deallocate(TI_PEM)
     64if (allocated(tsoil_PEM)) deallocate(tsoil_PEM)
     65if (allocated(mthermdiff_PEM)) deallocate(mthermdiff_PEM)
     66if (allocated(thermdiff_PEM)) deallocate(thermdiff_PEM)
     67if (allocated(coefq_PEM)) deallocate(coefq_PEM)
     68if (allocated(coefd_PEM)) deallocate(coefd_PEM)
     69if (allocated(alph_PEM)) deallocate(alph_PEM)
     70if (allocated(beta_PEM)) deallocate(beta_PEM)
     71if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)
     72if (allocated(water_reservoir)) deallocate(water_reservoir)
     73
     74END SUBROUTINE end_comsoil_h_PEM
     75
     76END MODULE comsoil_h_PEM

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

    r3130 r3143  
    3838
    3939! Conversion H2O/CO2 frost to perennial frost and vice versa
    40 real, parameter :: threshold_water_frost2perennial = 1000. !~ 1 m
    41 real, parameter :: threshold_co2_frost2perennial = 16000.   !~ 10 m
     40real, parameter :: threshold_h2o_frost2perennial = 1000. !~ 1 m
     41real, parameter :: threshold_co2_frost2perennial = 16000. !~ 10 m
    4242
    4343END MODULE constants_marspem_mod
    44 

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

    r3130 r3143  
    33implicit none
    44
     5! ----------------------------------------------------------------------
    56contains
     7! ----------------------------------------------------------------------
    68
    79!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    4446STOPPING = .false.
    4547
    46 ! Computation of the present surface of water ice sublimating
     48! Computation of the present surface of h2o ice
    4749present_surf = 0.
    4850do i = 1,ngrid
    4951    do islope = 1,nslope
    50         !if (initial_h2o_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
    51         if (initial_h2o_ice(i,islope) > 0.5) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
     52        if (initial_h2o_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
    5253    enddo
    5354enddo
    5455
    5556! Check of the criterion
    56 if (present_surf < ini_surf*(1. - water_ice_criterion) .or. present_surf > ini_surf*(1. + water_ice_criterion)) then
     57if (present_surf < ini_surf*(1. - water_ice_criterion)) then
    5758    STOPPING = .true.
    5859    write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold"
     60    write(*,*) "present_surf < ini_surf*(1. - water_ice_criterion)", present_surf < ini_surf*(1. - water_ice_criterion)
    5961    write(*,*) "Current surface of water ice sublimating =", present_surf
    6062    write(*,*) "Initial surface of water ice sublimating =", ini_surf
    6163    write(*,*) "Percentage of change accepted =", water_ice_criterion*100
    62     write(*,*) "present_surf < ini_surf*(1. - water_ice_criterion)", (present_surf < ini_surf*(1. - water_ice_criterion))
     64endif
     65if (present_surf > ini_surf*(1. + water_ice_criterion)) then
     66    STOPPING = .true.
     67    write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold"
     68    write(*,*) "present_surf > ini_surf*(1. + water_ice_criterion)", present_surf > ini_surf*(1. + water_ice_criterion)
     69    write(*,*) "Current surface of water ice sublimating =", present_surf
     70    write(*,*) "Initial surface of water ice sublimating =", ini_surf
     71    write(*,*) "Percentage of change accepted =", water_ice_criterion*100
    6372endif
    6473
    … …  
    105114STOPPING_ps = .false.
    106115
    107 ! Computation of the actual surface
     116! Computation of the present surface of co2 ice
    108117present_surf = 0.
    109118do i = 1,ngrid
    … …  
    114123
    115124! Check of the criterion
    116 if (present_surf < ini_surf*(1. - co2_ice_criterion) .or. present_surf > ini_surf*(1. + co2_ice_criterion)) then
     125if (present_surf < ini_surf*(1. - co2_ice_criterion)) then
    117126    STOPPING_ice = .true.
    118     write(*,*) "Reason of stopping: the surface of co2 ice sublimating reach the threshold"
     127    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reached the threshold"
     128    write(*,*) "present_surf < ini_surf*(1. - co2_ice_criterion)", present_surf < ini_surf*(1. - co2_ice_criterion)
    119129    write(*,*) "Current surface of co2 ice sublimating =", present_surf
    120130    write(*,*) "Initial surface of co2 ice sublimating =", ini_surf
    121131    write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
    122     write(*,*) "present_surf < ini_surf*(1. - co2_ice_criterion)", (present_surf < ini_surf*(1. - co2_ice_criterion))
     132endif
     133if (present_surf > ini_surf*(1. + co2_ice_criterion)) then
     134    STOPPING_ice = .true.
     135    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reach the threshold"
     136    write(*,*) "present_surf > ini_surf*(1. + co2_ice_criterion)", present_surf > ini_surf*(1. + co2_ice_criterion)
     137    write(*,*) "Current surface of co2 ice sublimating =", present_surf
     138    write(*,*) "Initial surface of co2 ice sublimating =", ini_surf
     139    write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
    123140endif
    124141
    125142if (abs(ini_surf) < 1.e-5) STOPPING_ice = .false.
    126143
    127 if (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion) .or. global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)) then
     144if (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)) then
    128145    STOPPING_ps = .true.
    129146    write(*,*) "Reason of stopping: the global pressure reach the threshold"
     147    write(*,*) "global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)", global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)
    130148    write(*,*) "Current global pressure =", global_ave_press_new
    131149    write(*,*) "PCM global pressure =", global_ave_press_GCM
    132150    write(*,*) "Percentage of change accepted =", ps_criterion*100.
    133     write(*,*) "global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)", (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion))
     151endif
     152if (global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)) then
     153    STOPPING_ps = .true.
     154    write(*,*) "Reason of stopping: the global pressure reach the threshold"
     155    write(*,*) "global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)", global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)
     156    write(*,*) "Current global pressure =", global_ave_press_new
     157    write(*,*) "PCM global pressure =", global_ave_press_GCM
     158    write(*,*) "Percentage of change accepted =", ps_criterion*100.
    134159endif
    135160

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

    r3130 r3143  
    11!------------------------
    22! I   Initialization
    3 !    I_a READ run.def
    4 !    I_b READ of start_evol.nc and starfi_evol.nc
     3!    I_a Read run.def
     4!    I_b Read of start_evol.nc and starfi_evol.nc
    55!    I_c Subslope parametrisation
    6 !    I_d READ PCM data and convert to the physical grid
     6!    I_d Read PCM data and convert to the physical grid
    77!    I_e Initialization of the PEM variable and soil
    8 !    I_f Compute tendencies & Save initial situation
     8!    I_f Compute tendencies
    99!    I_g Save initial PCM situation
    1010!    I_h Read the startpem.nc
    … …  
    4242use criterion_pem_stop_mod,       only: criterion_waterice_stop, criterion_co2_stop
    4343use constants_marspem_mod,        only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
    44                                         m_noco2, threshold_water_frost2perennial, threshold_co2_frost2perennial
     44                                        m_noco2, threshold_h2o_frost2perennial, threshold_co2_frost2perennial
    4545use evol_co2_ice_s_mod,           only: evol_co2_ice_s
    4646use evol_h2o_ice_s_mod,           only: evol_h2o_ice_s
    … …  
    8686    use paleoclimate_mod,   only: albedo_perennialco2
    8787    use comcstfi_h,         only: pi, rad, g, cpp, mugaz, r
     88    use surfini_mod,        only: surfini
    8889#else
    8990    use tracer_h,           only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names
    … …  
    131132
    132133! Variables to read start.nc
    133 character(len = *), parameter :: FILE_NAME_start = "start_evol.nc" ! Name of the file used for initialsing the PEM
     134character(*), parameter :: start_name = "start_evol.nc" ! Name of the file used to initialize the PEM
    134135
    135136! Dynamic variables
    … …  
    138139real, dimension(ip1jmp1,llm)        :: teta          ! temperature potentielle
    139140real, dimension(:,:,:), allocatable :: q             ! champs advectes
    140 real, dimension(ip1jmp1)            :: ps            ! pression  au sol
     141real, dimension(ip1jmp1)            :: ps            ! pression au sol
    141142real, dimension(:),     allocatable :: ps_start_GCM  ! (ngrid) pression  au sol
    142 real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) ! pression  au sol instantannées
     143real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) ! pression au sol instantannées
    143144real, dimension(ip1jmp1,llm)        :: masse         ! masse d'air
    144145real, dimension(ip1jmp1)            :: phis          ! geopotentiel au sol
    … …  
    146147
    147148! Variables to read starfi.nc
    148 character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" ! Name of the file used for initialsing the PEM
    149 character(2)                   :: str2
    150 integer                        :: ncid, varid, status                      ! Variable for handling opening of files
    151 integer                        :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf files
    152 integer                        :: lonvarid, latvarid, areavarid, sdvarid   ! Variable ID for Netcdf files
    153 integer                        :: apvarid, bpvarid                         ! Variable ID for Netcdf files
     149character(*), parameter :: startfi_name = "startfi_evol.nc" ! Name of the file used to initialize the PEM
     150character(2)            :: str2
     151integer                 :: ncid, varid, status                      ! Variable for handling opening of files
     152integer                 :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf files
     153integer                 :: lonvarid, latvarid, areavarid, sdvarid   ! Variable ID for Netcdf files
     154integer                 :: apvarid, bpvarid                         ! Variable ID for Netcdf files
    154155
    155156! Variables to read starfi.nc and write restartfi.nc
    156 real, dimension(:), allocatable :: longitude     ! Longitude read in FILE_NAME and written in restartfi
    157 real, dimension(:), allocatable :: latitude      ! Latitude read in FILE_NAME and written in restartfi
    158 real, dimension(:), allocatable :: cell_area     ! Cell_area read in FILE_NAME and written in restartfi
     157real, dimension(:), allocatable :: longitude     ! Longitude read in startfi_name and written in restartfi
     158real, dimension(:), allocatable :: latitude      ! Latitude read in startfi_name and written in restartfi
     159real, dimension(:), allocatable :: cell_area     ! Cell_area read in startfi_name and written in restartfi
    159160real                            :: Total_surface ! Total surface of the planet
    160161
    … …  
    187188real, dimension(:,:), allocatable :: co2frost_ini         ! Initial amount of co2 frost (at the start of the PEM run)
    188189real, dimension(:,:), allocatable :: perennial_co2ice_ini ! Initial amoun of perennial co2 ice (at the start of the PEM run)
    189 
    190 
    191190
    192191! Variables for slopes
    … …  
    243242    real                                :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice
    244243    real                                :: albedo_h2o_frost              ! albedo of h2o frost
    245     real, dimension(:),     allocatable :: tsurf_read_generic            ! Temporary variable to do the subslope transfert dimensiion when reading form generic
    246     real, dimension(:,:),   allocatable :: qsurf_read_generic            ! Temporary variable to do the subslope transfert dimensiion when reading form generic
    247     real, dimension(:,:),   allocatable :: tsoil_read_generic            ! Temporary variable to do the subslope transfert dimensiion when reading form generic
    248     real, dimension(:),     allocatable :: emis_read_generic             ! Temporary variable to do the subslope transfert dimensiion when reading form generic
    249     real, dimension(:,:),   allocatable :: albedo_read_generic           ! Temporary variable to do the subslope transfert dimensiion when reading form generic
     244    real, dimension(:),     allocatable :: tsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
     245    real, dimension(:,:),   allocatable :: qsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
     246    real, dimension(:,:),   allocatable :: tsoil_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
     247    real, dimension(:),     allocatable :: emis_read_generic             ! Temporary variable to do the subslope transfert dimension when reading form generic
     248    real, dimension(:,:),   allocatable :: albedo_read_generic           ! Temporary variable to do the subslope transfert dimension when reading form generic
    250249    real, dimension(:,:),   allocatable :: tsurf                         ! Subslope variable, only needed in the GENERIC case
    251250    real, dimension(:,:,:), allocatable :: qsurf                         ! Subslope variable, only needed in the GENERIC case
    … …  
    259258
    260259#ifdef CPP_1D
    261     integer                           :: nsplit_phys
    262     integer, parameter                :: jjm_value = jjm - 1
    263     integer                           :: ierr
     260    integer            :: nsplit_phys
     261    integer, parameter :: jjm_value = jjm - 1
     262    integer            :: ierr
    264263
    265264    ! Dummy variables to use the subroutine 'init_testphys1d'
    … …  
    272271    real, dimension(nlayer + 1)         :: plev
    273272#else
    274     integer, parameter                :: jjm_value = jjm
    275     real, dimension(:), allocatable   :: ap ! Coefficient ap read in FILE_NAME_start and written in restart
    276     real, dimension(:), allocatable   :: bp ! Coefficient bp read in FILE_NAME_start and written in restart
     273    integer, parameter              :: jjm_value = jjm
     274    real, dimension(:), allocatable :: ap ! Coefficient ap read in start_name and written in restart
     275    real, dimension(:), allocatable :: bp ! Coefficient bp read in start_name and written in restart
    277276#endif
    278277
    … …  
    302301!------------------------
    303302! I   Initialization
    304 !    I_a READ run.def
     303!    I_a Read run.def
    305304!------------------------
    306305#ifndef CPP_1D
    … …  
    315314    allocate(longitude(1),latitude(1),cell_area(1))
    316315
    317     therestart1D = .false.
     316    therestart1D = .false. ! Default value
    318317    inquire(file = 'start1D_evol.txt',exist = therestart1D)
    319318    if (.not. therestart1D) then
    … …  
    321320        error stop 'Initialization cannot be done for the 1D PEM.'
    322321    endif
    323     therestartfi = .false.
     322    therestartfi = .false. ! Default value
    324323    inquire(file = 'startfi_evol.nc',exist = therestartfi)
    325324    if (.not. therestartfi) then
    … …  
    339338!------------------------
    340339! I   Initialization
    341 !    I_b READ of start_evol.nc and starfi_evol.nc
    342 !------------------------
    343 ! I_b.1 READ start_evol.nc
     340!    I_b Read of start_evol.nc and starfi_evol.nc
     341!------------------------
     342! I_b.1 Read start_evol.nc
    344343allocate(ps_start_GCM(ngrid))
    345344#ifndef CPP_1D
    346     call dynetat0(FILE_NAME_start,vcov,ucov,teta,q,masse,ps,phis,time_0)
     345    call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0)
    347346
    348347    call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM)
    … …  
    353352    allocate(ap(nlayer + 1))
    354353    allocate(bp(nlayer + 1))
    355     status = nf90_open(FILE_NAME_start,NF90_NOWRITE,ncid)
     354    status = nf90_open(start_name,NF90_NOWRITE,ncid)
    356355    status = nf90_inq_varid(ncid,"ap",apvarid)
    357356    status = nf90_get_var(ncid,apvarid,ap)
    … …  
    368367! In the PCM, these values are given to the physic by the dynamic.
    369368! Here we simply read them in the startfi_evol.nc file
    370 status = nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
     369status = nf90_open(startfi_name, NF90_NOWRITE, ncid)
    371370
    372371status = nf90_inq_varid(ncid,"longitude",lonvarid)
    … …  
    384383status = nf90_close(ncid)
    385384
    386 ! I_b.2 READ startfi_evol.nc
     385! I_b.2 Read startfi_evol.nc
    387386! First we read the initial state (starfi.nc)
    388387#ifndef CPP_STD
    389     call phyetat0(FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
     388    call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
    390389                  tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,       &
    391390                  watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
    … …  
    394393    where (qsurf < 0.) qsurf = 0.
    395394
    396     call surfini(ngrid,qsurf)
     395    call surfini(ngrid,nslope,qsurf)
    397396#else
    398397    call phys_state_var_init(nq)
    … …  
    414413    allocate(albedo(ngrid,2,1))
    415414    allocate(inertiesoil(ngrid,nsoilmx,1))
    416     call phyetat0(.true.,ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
     415    call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
    417416                  tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,            &
    418417                  qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice,     &
    … …  
    449448    if (noms(nnq) == "h2o_vap") then
    450449        igcm_h2o_vap = nnq
    451         mmol(igcm_h2o_vap)=18.
     450        mmol(igcm_h2o_vap) = 18.
    452451    endif
    453452    if (noms(nnq) == "co2") igcm_co2 = nnq
    … …  
    480479!------------------------
    481480! I   Initialization
    482 !    I_d READ PCM data and convert to the physical grid
     481!    I_d Read PCM data and convert to the physical grid
    483482!------------------------
    484483! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the PCM run, saving only the minimum value
    … …  
    552551!------------------------
    553552! I   Initialization
    554 !    I_f Compute tendencies & Save initial situation
     553!    I_f Compute tendencies
    555554!------------------------
    556555allocate(tendencies_co2_ice(ngrid,nslope))
    … …  
    603602write(*,*) "Total initial surface of co2 ice sublimating (slope) =", ini_surf_co2
    604603write(*,*) "Total initial surface of h2o ice sublimating (slope) =", ini_surf_h2o
    605 write(*,*) "Total surface of the planet=", Total_surface
     604write(*,*) "Total surface of the planet =", Total_surface
    606605allocate(zplev_gcm(ngrid,nlayer + 1))
    607606
    … …  
    617616global_ave_press_GCM = global_ave_press_old
    618617global_ave_press_new = global_ave_press_old
    619 write(*,*) "Initial global average pressure=", global_ave_press_GCM
     618write(*,*) "Initial global average pressure =", global_ave_press_GCM
    620619
    621620!------------------------
    … …  
    638637do ig = 1,ngrid
    639638    do islope = 1,nslope
    640         qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
     639        qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
    641640        qsurf(ig,igcm_co2,islope) = qsurf(ig,igcm_co2,islope) + perennial_co2ice(ig,islope)
    642641    enddo
    … …  
    657656    enddo
    658657
    659     write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
    660     write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
     658    write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded
     659    write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded
    661660endif ! adsorption
    662661deallocate(tsurf_ave_yr1)
    … …  
    735734            q_h2o_PEM_phys(ig,t) = q_h2o_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ &
    736735                                   (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
    737             if (q_h2o_PEM_phys(ig,t) < 0) q_h2o_PEM_phys(ig,t) = 1.e-30
    738             if (q_h2o_PEM_phys(ig,t) > 1) q_h2o_PEM_phys(ig,t) = 1.
     736            if (q_h2o_PEM_phys(ig,t) < 0) then
     737                q_h2o_PEM_phys(ig,t) = 1.e-30
     738            else if (q_h2o_PEM_phys(ig,t) > 1) then
     739                q_h2o_PEM_phys(ig,t) = 1.
     740            endif
    739741        enddo
    740742    enddo
    … …  
    749751            if (q_co2_PEM_phys(ig,t) < 0) then
    750752                q_co2_PEM_phys(ig,t) = 1.e-30
    751             elseif (q_co2_PEM_phys(ig,t) > 1) then
     753            else if (q_co2_PEM_phys(ig,t) > 1) then
    752754                q_co2_PEM_phys(ig,t) = 1.
    753755            endif
    … …  
    898900        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
    899901
    900 ! II_d.5 Update the mass of the regolith adsorbded
     902! II_d.5 Update the mass of the regolith adsorbed
    901903        if (adsorption_pem) then
    902904            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, &
    … …  
    985987        exit
    986988    else
    987         write(*,*) "We continue!"
     989        write(*,*) "The PEM can continue!"
    988990        write(*,*) "Number of iterations of the PEM: year_iter =", year_iter
    989991        write(*,*) "Number of simulated Martian years: i_myear =", i_myear
    … …  
    10071009        watercap_sum = 0.
    10081010        do islope = 1,nslope
    1009             if (qsurf(ig,igcm_h2o_ice,islope) > (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then ! water_reservoir and water cap have not changed since PCM call: here we check if we have accumulate frost or not. 1st case we have more ice than initialy
    1010                 qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) ! put back ancien frost
     1011            ! water_reservoir and watercap have not changed since PCM call: here we check if we have accumulate frost or not.
     1012            if (qsurf(ig,igcm_h2o_ice,islope) > water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) then
     1013                ! 1st case: more ice than initialy. Put back ancien frost
     1014                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
    10111015            else
    1012 !               2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap
    1013                 watercap(ig,islope) = watercap(ig,islope) + qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
     1016                ! 2nd case: ice sublimated. Then let's put qsurf = 0 and add the difference in watercap
     1017                watercap(ig,islope) = qsurf(ig,igcm_h2o_ice,islope)
    10141018                qsurf(ig,igcm_h2o_ice,islope) = 0.
    10151019            endif
    … …  
    10271031    enddo
    10281032    if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming.
    1029         if (water_sum > threshold_water_frost2perennial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1
     1033        if (water_sum > threshold_h2o_frost2perennial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1
    10301034            watercaptag(ig) = .true.
    1031             water_reservoir(ig) = water_reservoir(ig) + threshold_water_frost2perennial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost
     1035            water_reservoir(ig) = water_reservoir(ig) + threshold_h2o_frost2perennial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost
    10321036            do islope = 1,nslope
    1033                 qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perennial/2.*cos(pi*def_slope_mean(islope)/180.)
     1037                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_h2o_frost2perennial/2.*cos(pi*def_slope_mean(islope)/180.)
    10341038            enddo
    10351039        endif
    10361040    else ! let's check that the infinite source of water disapear
    1037         if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perennial) then
     1041        if ((water_sum + water_reservoir(ig)) < threshold_h2o_frost2perennial) then
    10381042            watercaptag(ig) = .false.
    10391043            do islope = 1,nslope

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

    r3123 r3143  
    104104!0. Check if the start_PEM exist.
    105105
    106 inquire(file=filename,exist = startpem_file)
     106inquire(file = filename,exist = startpem_file)
    107107
    108108write(*,*)'Is start PEM?',startpem_file
    … …  
    312312!. 5 water reservoir
    313313#ifndef CPP_STD
     314   water_reservoir = 0.
    314315   call get_field("water_reservoir",water_reservoir,found)
    315    if(.not.found) then
    316       write(*,*)'Pemetat0: failed loading <water_reservoir>'
    317       write(*,*)'will reconstruct the values from watercaptag'
    318       do ig=1,ngrid
    319         if(watercaptag(ig)) then
    320            water_reservoir(ig)=water_reservoir_nom
    321         else
    322            water_reservoir(ig)=0.
    323         endif
    324       enddo
     316   if (.not. found) then
     317       write(*,*)'Pemetat0: failed loading <water_reservoir>'
     318       write(*,*)'will reconstruct the values from watercaptag'
     319       where (watercaptag) water_reservoir = water_reservoir_nom
    325320   endif
    326321#endif
    … …  
    481476!. e) water reservoir
    482477#ifndef CPP_STD
    483     do ig = 1,ngrid
    484         if (watercaptag(ig)) then
    485             water_reservoir(ig) = water_reservoir_nom
    486         else
    487             water_reservoir(ig) = 0.
    488         endif
    489     enddo
     478    water_reservoir = 0.
     479    where (watercaptag) water_reservoir = water_reservoir_nom
    490480#endif
    491481

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

    r3130 r3143  
    104104    write(*,*) "Data for h2o_ice_s downloaded!"
    105105
     106    call get_var3("tsurf",tsurf_gcm_dyn(:,:,1,:))
     107    write(*,*) "Data for tsurf downloaded!"
     108
    106109#ifndef CPP_STD
    107110    call get_var3("watercap",watercap(:,:,1,:))
    … …  
    110113    call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:))
    111114    write(*,*) "Data for perennial_co2ice downloaded!"
    112 #endif
    113 
    114     call get_var3("tsurf",tsurf_gcm_dyn(:,:,1,:))
    115     write(*,*) "Data for tsurf downloaded!"
    116 
    117 #ifndef CPP_STD
     115
    118116    if (soil_pem) then
    119117        call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
    … …  
    140138    write(*,*) "Data for h2o_ice_s downloaded!"
    141139
     140    do islope = 1,nslope
     141        write(num,'(i2.2)') islope
     142        call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
     143    enddo
     144    write(*,*) "Data for tsurf downloaded!"
     145
    142146#ifndef CPP_STD
    143147    do islope = 1,nslope
    … …  
    152156    enddo
    153157    write(*,*) "Data for perennial_co2ice downloaded!"
    154 #endif
    155 
    156     do islope = 1,nslope
    157         write(num,'(i2.2)') islope
    158         call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
    159     enddo
    160     write(*,*) "Data for tsurf downloaded!"
    161 
    162 #ifndef CPP_STD
     158
    163159    if (soil_pem) then
    164160        do islope = 1,nslope
Note: See TracChangeset for help on using the changeset viewer.