Changeset 3039


Ignore:
Timestamp:
Sep 11, 2023, 12:41:50 PM (17 months ago)
Author:
jbclement
Message:

Mars PEM:
Big cleaning of main program pem.F90 (indentation, declarations, comments, simplification of conditions/loops, etc).
JBC

Location:
trunk
Files:
12 edited
2 moved

Legend:

Unmodified
Added
Removed

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

    r3002 r3039  
    55! Purpose: Read the parameter for a PEM run in run_pem.def
    66!
    7 ! Author: RV 
     7! Author: RV, JBC
    88!=======================================================================
    99
    … …  
    1212CONTAINS
    1313
    14   SUBROUTINE conf_pem
     14  SUBROUTINE conf_pem(i_myear,n_myears)
    1515
    1616#ifdef CPP_IOIPSL
    … …  
    2121#endif
    2222 
    23   USE temps_mod_evol, ONLY: year_bp_ini, dt_pem, water_ice_criterion, co2_ice_criterion, ps_criterion, &
    24                 Max_iter_pem, evol_orbit_pem, var_obl, var_ex, var_lsp
    25   USE comsoil_h_pem, only: soil_pem,fluxgeo,water_reservoir_nom,depth_breccia,depth_bedrock,reg_thprop_dependp
    26   USE adsorption_mod,only: adsorption_pem
    27   USE glaciers_mod, only: co2glaciersflow,h2oglaciersflow
    28   use ice_table_mod, only: icetable_equilibrium, icetable_dynamic
     23  use time_evol_mod,  only: year_bp_ini, dt_pem, water_ice_criterion, co2_ice_criterion, ps_criterion, Max_iter_pem, &
     24                            evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years
     25  use comsoil_h_pem,  only: soil_pem, fluxgeo, water_reservoir_nom, depth_breccia, depth_bedrock, reg_thprop_dependp
     26  use adsorption_mod, only: adsorption_pem
     27  use glaciers_mod,   only: co2glaciersflow, h2oglaciersflow
     28  use ice_table_mod,  only: icetable_equilibrium, icetable_dynamic
    2929
    30   CHARACTER(len=20),parameter :: modname ='conf_pem'
     30  integer, intent(out) :: i_myear, n_myears ! Current simulated Martian year and maximum number of Martian years to be simulated
     31
     32  character(len=20), parameter :: modname ='conf_pem'
     33  integer                      :: ierr
     34  integer                      :: year_earth_bp_ini ! Initial year (in Earth years) of the simulation of the PEM defined in run.def
    3135
    3236!PEM parameters
    3337
    34 ! #---------- ORBITAL parameters --------------#
     38! #---------- Martian years parameters from launching script
     39open(100,file = 'tmp_PEMyears.txt',status = 'old',form = 'formatted',iostat = ierr)
     40if (ierr /= 0) then
     41    write(*,*) 'Cannot find required file "tmp_PEMyears.txt"!'
     42    write(*,*) 'It should be created by the launching script...'
     43    stop
     44else
     45    read(100,*) i_myear, n_myears, convert_years
     46endif
     47close(100)
    3548
     49! #---------- Orbital parameters
    3650  evol_orbit_pem=.false.
    37   CALL getin('evol_orbit_pem', evol_orbit_pem)
     51  call getin('evol_orbit_pem',evol_orbit_pem)
    3852
    39   year_bp_ini=0.
    40   CALL getin('year_bp_ini', year_bp_ini)
     53  year_bp_ini = 0.
     54  call getin('year_earth_bp_ini',year_earth_bp_ini)
     55  year_bp_ini = year_earth_bp_ini/convert_years
    4156
    4257  var_obl = .true.
    43   CALL getin('var_obl',var_obl)
    44   print*,'Does obliquity vary ?',var_obl
     58  call getin('var_obl',var_obl)
     59  write(*,*) 'Does obliquity vary ?',var_obl
    4560
    46   var_ex = .true.
    47   CALL getin('var_ex',var_ex)
    48   print*,'Does excentricity vary ?',var_ex
     61  var_ecc = .true.
     62  call getin('var_ecc',var_ecc)
     63  write(*,*) 'Does excentricity vary ?',var_ecc
    4964
    5065  var_lsp = .true.
    51   CALL getin('var_lsp',var_lsp)
    52   print*,'Does Ls peri vary ?',var_lsp
     66  call getin('var_lsp',var_lsp)
     67  write(*,*) 'Does Ls peri vary ?',var_lsp
    5368
    54 ! #---------- Stopping criterion parameters --------------#
    55 
    56   Max_iter_pem=99999999
    57   CALL getin('Max_iter_pem', Max_iter_pem)
     69! #---------- Stopping criteria parameters
     70  Max_iter_pem=100000000
     71  call getin('Max_iter_pem', Max_iter_pem)
    5872
    5973  water_ice_criterion=0.2
    60   CALL getin('water_ice_criterion', water_ice_criterion)
     74  call getin('water_ice_criterion', water_ice_criterion)
    6175
    6276  co2_ice_criterion=0.2
    63   CALL getin('co2_ice_criterion', co2_ice_criterion)
     77  call getin('co2_ice_criterion', co2_ice_criterion)
    6478
    6579  ps_criterion = 0.15
    66   CALL getin('ps_criterion',ps_criterion)
     80  call getin('ps_criterion',ps_criterion)
    6781
    6882  dt_pem=1
    69   CALL getin('dt_pem', dt_pem)
     83  call getin('dt_pem', dt_pem)
    7084
    71 ! #---------- Subsurface parameters --------------#
    72 
     85! #---------- Subsurface parameters
    7386  soil_pem=.true.
    74   CALL getin('soil_pem', soil_pem)
     87  call getin('soil_pem', soil_pem)
    7588
    7689  adsorption_pem = .true.
    77   CALL getin('adsorption_pem',adsorption_pem)
     90  call getin('adsorption_pem',adsorption_pem)
    7891
    7992  co2glaciersflow = .true.
    80   CALL getin('co2glaciersflow', co2glaciersflow)
     93  call getin('co2glaciersflow', co2glaciersflow)
    8194
    8295  h2oglaciersflow = .true.
    83   CALL getin('h2oglaciersflow', h2oglaciersflow)
     96  call getin('h2oglaciersflow', h2oglaciersflow)
    8497
    8598  reg_thprop_dependp = .false.
    86   CALL getin('reg_thprop_dependp',reg_thprop_dependp)
    87   print*, 'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
     99  call getin('reg_thprop_dependp',reg_thprop_dependp)
     100  write(*,*) 'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
    88101
    89 ! #---------- Layering parameters --------------#
    90 
     102! #---------- Layering parameters
    91103  fluxgeo = 0.
    92   CALL getin('fluxgeo',fluxgeo)
    93   print*,'Flux Geothermal is set to',fluxgeo
     104  call getin('fluxgeo',fluxgeo)
     105  write(*,*) 'Flux Geothermal is set to',fluxgeo
    94106   
    95107  depth_breccia   = 10.
    96   CALL getin('depth_breccia',depth_breccia)
    97   print*,'Depth of breccia is set to',depth_breccia
     108  call getin('depth_breccia',depth_breccia)
     109  write(*,*) 'Depth of breccia is set to',depth_breccia
    98110
    99111  depth_bedrock   = 1000.
    100   CALL getin('depth_bedrock',depth_bedrock)
    101   print*,'Depth of bedrock is set to',depth_bedrock
     112  call getin('depth_bedrock',depth_bedrock)
     113  write(*,*) 'Depth of bedrock is set to',depth_bedrock
    102114
    103115   icetable_equilibrium = .true.
    104    CALL getin('icetable_equilibrium',icetable_equilibrium)
    105    print*, 'Do we compute the ice table at equilibrium?', icetable_equilibrium
     116   call getin('icetable_equilibrium',icetable_equilibrium)
     117   write(*,*)  'Is the ice table computed at equilibrium?', icetable_equilibrium
    106118
    107119   icetable_dynamic = .false.
    108    CALL getin('icetable_dynamic',icetable_dynamic)
    109    print*, 'Do we compute the ice table with the dynamic method?', icetable_dynamic
     120   call getin('icetable_dynamic',icetable_dynamic)
     121   write(*,*)  'Is the ice table computed with the dynamic method?', icetable_dynamic
    110122  if ((.not.soil_pem).and.((icetable_equilibrium).or.(icetable_dynamic))) then
    111        print*,'Ice table  must be used when soil_pem = T'
     123       write(*,*) 'Ice table  must be used when soil_pem = T'
    112124       call abort_physic(modname,"Ice table  must be used when soil_pem = T",1)
    113125  endif
    114126
    115127  if ((.not.soil_pem).and.adsorption_pem) then
    116        print*,'Adsorption must be used when soil_pem = T'
     128       write(*,*) 'Adsorption must be used when soil_pem = T'
    117129       call abort_physic(modname,"Adsorption must be used when soil_pem = T",1)
    118130  endif
    119131 
    120132  if ((.not.soil_pem).and.(fluxgeo.gt.0.)) then
    121        print*,'Soil is not activated but Flux Geo > 0.'
     133       write(*,*) 'Soil is not activated but Flux Geo > 0.'
    122134       call abort_physic(modname,"Soil is not activated but Flux Geo > 0.",1)
    123135  endif
    124136 
    125137  if ((.not.soil_pem).and.reg_thprop_dependp) then
    126      print*,'Regolith properties vary with Ps only when soil is set to true'
     138     write(*,*) 'Regolith properties vary with Ps only when soil is set to true'
    127139     call abort_physic(modname,'Regolith properties vary with Ps only when soil is set to true',1)
    128140  endif
    129141
    130   if (evol_orbit_pem.and.year_bp_ini.eq.0.) then
    131      print*,'You want to follow the file ob_ex_lsp.asc for changing orb parameters,'
    132      print*,'but you did not specify from which year to start.'
     142  if (evol_orbit_pem .and. year_bp_ini == 0.) then
     143     write(*,*)  'You want to follow the file obl_ecc_lsp.asc for changing orb parameters,'
     144     write(*,*)  'but you did not specify from which year to start.'
    133145     call abort_physic(modname,"evol_orbit_pem=.true. but year_bp_ini=0",1)
    134146  endif
    135147
    136148  water_reservoir_nom = 1e4
    137   CALL getin('water_reservoir_nom',water_reservoir_nom)
     149  call getin('water_reservoir_nom',water_reservoir_nom)
    138150
    139151  END SUBROUTINE conf_pem

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

    r2980 r3039  
    1515SUBROUTINE criterion_waterice_stop(cell_area,ini_surf,qsurf,STOPPING,ngrid,initial_h2o_ice)
    1616
    17   USE temps_mod_evol, ONLY: water_ice_criterion
    18   use comslope_mod, ONLY: subslope_dist,nslope
     17  use time_evol_mod, only: water_ice_criterion
     18  use comslope_mod,  only: subslope_dist,nslope
    1919
    2020      IMPLICIT NONE
    … …  
    7979SUBROUTINE criterion_co2_stop(cell_area,ini_surf,qsurf,STOPPING_ice,STOPPING_ps,ngrid,initial_co2_ice,global_ave_press_GCM,global_ave_press_new,nslope)
    8080
    81   USE temps_mod_evol, ONLY: co2_ice_criterion,ps_criterion
    82   use comslope_mod, ONLY: subslope_dist
     81  use time_evol_mod, only: co2_ice_criterion,ps_criterion
     82  use comslope_mod,  only: subslope_dist
    8383
    8484      IMPLICIT NONE

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

    r2980 r3039  
    88                             iim_input,jjm_input,ngrid,cell_area,nslope)
    99
    10   USE temps_mod_evol, ONLY: dt_pem
     10  use time_evol_mod, only: dt_pem
    1111
    1212      IMPLICIT NONE

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

    r3031 r3039  
    77SUBROUTINE evol_h2o_ice_s(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,qsurf,tendencies_h2o_ice_phys,STOPPING)
    88
    9   use temps_mod_evol, only: dt_pem
    10   use comslope_mod, only: subslope_dist,def_slope_mean
     9  use time_evol_mod,          only: dt_pem
     10  use comslope_mod,           only: subslope_dist, def_slope_mean
    1111  use criterion_pem_stop_mod, only: criterion_waterice_stop
    12   use comconst_mod,only: pi
     12  use comconst_mod,           only: pi
    1313
    1414  IMPLICIT NONE

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

    r2980 r3039  
    1 SUBROUTINE info_run_PEM(year_iter, criterion_stop)
     1SUBROUTINE info_run_PEM(year_iter,criterion_stop,i_myear,n_myear)
    22
    33!=======================================================================
    44!
    55! Purpose: Write in a file called info_run_PEM.txt the reason why the PEM did stop and the number of extrapolation year done
     6!          Update the file tmp_PEMyears.txt to count the number of simulated Martian years
    67!
    7 ! Author: RV 
     8! Author: RV, JBC
    89!=======================================================================
    910
    10    IMPLICIT NONE
     11  use time_evol_mod, only: convert_years
    1112
    12   INTEGER, intent(in) :: year_iter, criterion_stop ! # of year and reason to stop
     13  IMPLICIT NONE
    1314
    14   logical :: exist
     15  integer, intent(in) :: year_iter, criterion_stop ! # of year and reason to stop
     16  integer, intent(in) :: i_myear, n_myear ! Current simulated Martian year and maximum number of Martian years to be simulated
    1517
    16   inquire(file='info_run_PEM.txt', exist=exist)
    17   if (exist) then
    18     open(12, file='info_run_PEM.txt', status="old", position="append", action="write")
    19   else
    20     open(12, file='info_run_PEM.txt', status="new", action="write")
    21   end if
    22   write(12, *) year_iter, criterion_stop
    23   close(12)
     18  logical :: ok
     19
     20inquire(file = 'info_run_PEM.txt', exist = ok)
     21if (ok) then
     22    open(12,file = 'info_run_PEM.txt',status = "old",position = "append",action = "write")
     23else
     24    open(12,file = 'info_run_PEM.txt',status = "new",action = "write")
     25endif
     26write(12,*) year_iter, i_myear, criterion_stop
     27close(12)
     28
     29open(100,file = 'tmp_PEMyears.txt',status = 'replace')
     30write(100,*) i_myear, n_myear, convert_years
     31close(100)
    2432
    2533END SUBROUTINE info_run_PEM

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

    r2980 r3039  
    44!!! Purpose: Define parameters from Laskar et al., 2004 evolution
    55!!!
    6 !!! Author: RV
     6!!! Author: RV, JBC
    77!!!
    88!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    1010implicit none
    1111
    12   real,save,allocatable :: yearlask(:)          ! year before present from Laskar et al. Tab
    13   real,save,allocatable :: oblask(:)            ! obliquity    [deg]
    14   real,save,allocatable :: exlask(:)            ! excentricity [deg]
    15   real,save,allocatable :: lsplask(:)           ! ls perihelie [deg]
    16   integer, save :: last_ilask                   ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
     12  real,save,allocatable :: yearlask(:) ! year before present from Laskar et al. Tab
     13  real,save,allocatable :: obllask(:)  ! obliquity    [deg]
     14  real,save,allocatable :: ecclask(:)  ! excentricity [deg]
     15  real,save,allocatable :: lsplask(:)  ! ls perihelie [deg]
     16  integer, save         :: last_ilask  ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
    1717
    1818contains
    1919
    20   subroutine ini_lask_param_mod(nlask)
     20  subroutine ini_lask_param_mod
    2121 
    2222  implicit none
    23   integer,intent(in) :: nlask ! number of line
    24  
     23
     24  integer :: nlask ! number of lines in Laskar files
     25  integer :: ierr
     26
     27  nlask = 0
     28  open(1,file = 'obl_ecc_lsp.asc')
     29  do
     30      read(1,*,iostat = ierr)
     31      if (ierr /= 0) exit
     32      nlask = nlask + 1
     33  enddo
     34  close(1)
    2535  allocate(yearlask(nlask))
    26   allocate(oblask(nlask))
    27   allocate(exlask(nlask))
     36  allocate(obllask(nlask))
     37  allocate(ecclask(nlask))
    2838  allocate(lsplask(nlask))
    2939 
    … …  
    3646
    3747  if (allocated(yearlask)) deallocate(yearlask)
    38   if (allocated(oblask)) deallocate(oblask)
    39   if (allocated(exlask)) deallocate(exlask)
     48  if (allocated(obllask)) deallocate(obllask)
     49  if (allocated(ecclask)) deallocate(ecclask)
    4050  if (allocated(lsplask)) deallocate(lsplask)
    4151

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

    r3026 r3039  
    1       MODULE orbit_param_criterion_mod
     1MODULE orbit_param_criterion_mod
    22
    33!=======================================================================
    … …  
    99!=======================================================================
    1010
    11       IMPLICIT NONE
    12 
    13       CONTAINS
    14 
    15       SUBROUTINE orbit_param_criterion(year_iter_max)
     11IMPLICIT NONE
     12
     13CONTAINS
     14
     15SUBROUTINE orbit_param_criterion(i_myear,year_iter_max)
     16
    1617#ifdef CPP_IOIPSL
    17   use IOIPSL, only: getin
     18    use IOIPSL,          only: getin
    1819#else
    19   ! if not using IOIPSL, we still need to use (a local version of) getin
    20   use ioipsl_getincom, only: getin
     20    ! if not using IOIPSL, we still need to use (a local version of) getin
     21    use ioipsl_getincom, only: getin
    2122#endif
    22 
    23       USE temps_mod_evol, ONLY: year_bp_ini, year_PEM, var_obl, var_ex, var_lsp
    2423#ifndef CPP_STD
    25       USE planete_h, ONLY: e_elips, obliquit, timeperi
     24    use planete_h,   only: e_elips, obliquit, timeperi
    2625#else
    27       use planete_mod, only: e_elips, obliquit, timeperi
     26    use planete_mod, only: e_elips, obliquit, timeperi
    2827#endif
    29       USE comconst_mod, ONLY: pi
    30       USE lask_param_mod, only: yearlask,oblask,exlask,lsplask, &
    31                                 ini_lask_param_mod,last_ilask
    32 
    33       IMPLICIT NONE
     28  use time_evol_mod,  only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
     29  use comconst_mod,   only: pi
     30  use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, ini_lask_param_mod, last_ilask
     31
     32  IMPLICIT NONE
    3433
    3534!--------------------------------------------------------
    3635! Input Variables
    3736!--------------------------------------------------------
     37  integer, intent(in) :: i_myear ! Martian year of the simulation
    3838
    3939!--------------------------------------------------------
    4040! Output Variables
    4141!--------------------------------------------------------
    42 
    43       integer,intent(out) :: year_iter_max      ! Maximum number of iteration of the PEM before orb changes too much
     42  integer, intent(out) :: year_iter_max ! Maximum number of iteration of the PEM before orb changes too much
    4443
    4544!--------------------------------------------------------
    4645! Local variables
    4746!--------------------------------------------------------
    48 
    49       real :: Year                      ! Year of the simulation
    50       real :: Lsp                       ! Year of the simulation
    51       integer nlask, ilask, iylask      ! Loop variables
    52       parameter (nlask = 20001)         ! Number of line in Laskar file
    53 
    54       real max_change_obl,max_change_ex,max_change_lsp ! Percentage of change that is considered to be acceptible
    55 
    56       real max_obl,max_ex,max_lsp       ! Maximum value of orbit param given the acceptable percentage
    57       real min_obl,min_ex,min_lsp       ! Maximum value of orbit param given the acceptable percentage
    58       real max_obl_iter,max_ex_iter,max_lsp_iter ! Maximum year iteration before reaching an unacceptable value
    59       real xa,xb,ya,yb,yc
    60 
    61       ! **********************************************************************
    62       ! 0. Initializations
    63       ! **********************************************************************
    64 
    65           Year = year_bp_ini + year_PEM      ! We compute the current year
    66           Lsp = 360. - timeperi*360./(2.*pi) ! We convert in degree
    67 
    68           call ini_lask_param_mod(nlask) ! Allocation of variables
    69 
    70           print*, "orbit_param_criterion, Year in pem.def=", year_bp_ini
    71           print*, "orbit_param_criterion, Year in the startpem.nc =", year_PEM
    72           print*, "orbit_param_criterion, Current year=", Year
    73           print*, "orbit_param_criterion, Current obl=", obliquit
    74           print*, "orbit_param_criterion, Current ex=", e_elips
    75           print*, "orbit_param_criterion, Current lsp=", Lsp
     47  integer :: Year                                           ! Year of the simulation
     48  real    :: Lsp                                            ! Ls perihelion
     49  real    :: max_change_obl, max_change_ecc, max_change_lsp ! Maximum admissible change
     50  real    :: max_obl, max_ecc, max_lsp                      ! Maximum value of orbit param given the admissible change
     51  real    :: min_obl, min_ecc, min_lsp                      ! Minimum value of orbit param given the admissible change
     52  integer :: max_obl_iter, max_ecc_iter, max_lsp_iter       ! Maximum year iteration before reaching an inadmissible value of orbit param
     53  integer :: ilask, iylask                                  ! Loop variables
     54  real    :: xa, xb, ya, yb, yc_obl, yc_ecc, yc_lsp         ! Variables for interpolation
     55
     56! **********************************************************************
     57! 0. Initializations
     58! **********************************************************************
     59Year = year_bp_ini + i_myear        ! We compute the current year
     60Lsp = modulo(-timeperi*180./pi,360.) ! We convert in degree
     61
     62call ini_lask_param_mod ! Allocation of variables
     63
     64write(*,*) 'orbit_param_criterion, Current year =', Year
     65write(*,*) 'orbit_param_criterion, Current obl. =', obliquit
     66write(*,*) 'orbit_param_criterion, Current ecc. =', e_elips
     67write(*,*) 'orbit_param_criterion, Current Lsp  =', Lsp
    7668
    7769! We read the file
    78 
    79           open(73,file='ob_ex_lsp.asc')
    80           do ilask = 1,nlask
    81             read(73,*) yearlask(ilask), oblask(ilask), exlask(ilask), lsplask(ilask)
    82             yearlask(ilask) = yearlask(ilask)*1000.
    83             if (yearlask(ilask) > Year) last_ilask = ilask + 1
    84           enddo
    85           close(73)
    86 
    87        print*, "Coresponding line in the ob_ex_lsp.asc file=", last_ilask
    88 
    89 !Constant max change case
    90 
    91         max_change_obl = 0.5
    92         max_change_ex = 0.1
    93         max_change_lsp = 40.
    94 
    95         CALL getin('max_change_obl', max_change_obl)
    96         max_obl = obliquit + max_change_obl
    97         min_obl = obliquit - max_change_obl
    98         print*, "Maximum obliquity accepted=", max_obl
    99         print*, "Minimum obliquity accepted=", min_obl
    100 
    101         CALL getin('max_change_ex', max_change_ex)
    102         max_ex = e_elips + max_change_ex
    103         min_ex = e_elips - max_change_ex
    104         print*, "Maximum excentricity accepted=", max_ex
    105         print*, "Minimum excentricity accepted=", min_ex
    106 
    107         CALL getin('max_change_lsp', max_change_lsp)
    108         max_lsp = Lsp + max_change_lsp
    109         min_lsp = Lsp - max_change_lsp
    110         print*, "Maximum lsp accepted=", max_lsp
    111         print*, "Minimum lsp accepted=", min_lsp
    112 
    113 !End Constant max change case
     70open(73,file = 'obl_ecc_lsp.asc')
     71do ilask = 1,size(yearlask,1)
     72    read(73,*) yearlask(ilask), obllask(ilask), ecclask(ilask), lsplask(ilask)
     73    yearlask(ilask) = yearlask(ilask)*1000./convert_years ! Conversion from Laskar's kilo Earth years to PEM Martian years
     74    if (yearlask(ilask) > Year) last_ilask = ilask + 1
     75enddo
     76close(73)
     77write(*,*) 'Corresponding line in the obl_ecc_lsp.asc file =', last_ilask
     78
     79! Constant max change case
     80max_change_obl = 0.6
     81max_change_ecc = 2.8e-3
     82max_change_lsp = 18.
     83
     84call getin('max_change_obl', max_change_obl)
     85max_obl = obliquit + max_change_obl
     86min_obl = obliquit - max_change_obl
     87write(*,*) 'Maximum obliquity accepted =', max_obl
     88write(*,*) 'Minimum obliquity accepted =', min_obl
     89
     90call getin('max_change_ecc', max_change_ecc)
     91max_ecc = min(e_elips + max_change_ecc,1. - 1.e-10) ! ex < 1.
     92min_ecc = max(e_elips - max_change_ecc,0.) ! ex >= 0.
     93write(*,*) 'Maximum eccentricity accepted =', max_ecc
     94write(*,*) 'Minimum eccentricity accepted =', min_ecc
     95
     96call getin('max_change_lsp', max_change_lsp)
     97max_lsp = Lsp + max_change_lsp
     98min_lsp = Lsp - max_change_lsp
     99write(*,*) 'Maximum Lsp accepted =', max_lsp
     100write(*,*) 'Minimum Lsp accepted =', min_lsp
     101! End Constant max change case
    114102
    115103! If we do not want some orb parameter to change, they should not be a stopping criterion,
    116104! So the number of iteration corresponding is set to maximum
    117         max_obl_iter = 999999
    118         max_ex_iter = 999999
    119         max_lsp_iter = 999999
     105max_obl_iter = 1000000
     106max_ecc_iter = 1000000
     107max_lsp_iter = 1000000
    120108
    121109        ! Tendency of the orbital parameter for the considered year gives the limitation between min and max
    122         do ilask = last_ilask,2,-1
    123           xa = yearlask(ilask)
    124           xb = yearlask(ilask - 1)
    125           if (xa <= Year .and. Year < xb) then
    126             ! Obliquity
    127             if (var_obl) then
    128               ya = oblask(ilask)
    129               yb = oblask(ilask - 1)
    130               if (ya < yb) then ! Increasing -> max is the limitation
    131                 yc = max_obl
    132               else ! Decreasing -> min is the limitation
    133                 yc = min_obl
    134               endif
    135             endif
    136 
    137             ! Excentricity
    138             if (var_ex) then
    139               ya = exlask(ilask)
    140               yb = exlask(ilask - 1)
    141               if (ya < yb) then ! Increasing -> max is the limitation
    142                 yc = max_ex
    143               else ! Decreasing -> min is the limitation
    144                 yc = min_ex
    145               endif
    146             endif
    147 
    148             ! Lsp
    149             if (var_lsp) then
    150               ya = lsplask(ilask)
    151               yb = lsplask(ilask - 1)
    152               if (ya < yb) then ! Increasing -> max is the limitation
    153                 if (yb - ya > 300.) then ! If modulo is "crossed" then it is the other way around
    154                   !yb = yb - 360.
    155                   yc = min_lsp
    156                 else
    157                   yc = max_lsp
    158                 endif
    159               else ! Decreasing -> min is the limitation
    160                 if (ya - yb > 300.) then ! If modulo is "crossed" then it is the other way around
    161                   !yb = yb + 360.
    162                   yc = max_lsp
    163                 else
    164                   yc = min_lsp
    165                 endif
    166               endif
    167             endif
    168             iylask = ilask
    169             exit ! The loop is left as soon as the right interval is found
    170           endif
    171         enddo
    172         if (ilask == 1) then
    173           write(*,*) 'The year does not match with Laskar data in the file ob_ex_lsp.asc.'
    174           stop
    175         endif
    176 
    177         ! Linear interpolation gives the maximum reachable year according to the limitation
     110do ilask = last_ilask,2,-1
     111    xa = yearlask(ilask)
     112    xb = yearlask(ilask - 1)
     113    if (xa <= Year .and. Year < xb) then
    178114        ! Obliquity
    179115        if (var_obl) then
    180           do ilask = iylask,2,-1
    181             ya = oblask(ilask)
    182             yb = oblask(ilask - 1)
    183             if (ya <= yc .and. yc < yb) then
    184               xa = yearlask(ilask)
    185               xb = yearlask(ilask - 1)
    186               max_obl_iter = int((max_obl - ya)*(xb - xa)/(yb - ya) + xa) - Year
    187               exit ! The loop is left as soon as the right interval is found
     116            ya = obllask(ilask)
     117            yb = obllask(ilask - 1)
     118            if (ya < yb) then ! Increasing -> max is the limitation
     119                yc_obl = max_obl
     120            else ! Decreasing -> min is the limitation
     121                yc_obl = min_obl
    188122            endif
    189           enddo
    190         endif
    191         ! Excentricity
    192         if (var_ex) then
    193           do ilask = iylask,2,-1
    194             ya = exlask(ilask)
    195             yb = exlask(ilask - 1)
    196             if (ya <= yc .and. yc < yb) then
    197               xa = yearlask(ilask)
    198               xb = yearlask(ilask - 1)
    199               max_ex_iter = int((max_ex - ya)*(xb - xa)/(yb - ya) + xa) - Year
    200               exit ! The loop is left as soon as the right interval is found
    201             endif           
    202           enddo
    203         endif
     123        endif
     124
     125        ! Eccentricity
     126        if (var_ecc) then
     127            ya = ecclask(ilask)
     128            yb = ecclask(ilask - 1)
     129            if (ya < yb) then ! Increasing -> max is the limitation
     130                yc_ecc = max_ecc
     131            else ! Decreasing -> min is the limitation
     132                yc_ecc = min_ecc
     133            endif
     134        endif
     135
    204136        ! Lsp
    205137        if (var_lsp) then
    206           do ilask = iylask,2,-1
    207138            ya = lsplask(ilask)
    208139            yb = lsplask(ilask - 1)
    209             if (ya <= yc .and. yc < yb) then
    210               xa = yearlask(ilask)
    211               xb = yearlask(ilask - 1)
    212               max_lsp_iter = int((max_lsp - ya)*(xb - xa)/(yb - ya) + xa) - Year
    213               exit ! The loop is left as soon as the right interval is found
     140             if (ya < yb) then ! Increasing -> max is the limitation
     141                if (yb - ya > 300.) then ! If modulo is "crossed" then it is the other way around
     142                    yc_lsp = min_lsp
     143                else
     144                    yc_lsp = max_lsp
     145                endif
     146            else ! Decreasing -> min is the limitation
     147                if (ya - yb > 300.) then ! If modulo is "crossed" then it is the other way around
     148                    yc_lsp = max_lsp
     149                else
     150                    yc_lsp = min_lsp
     151                endif
    214152            endif
    215           enddo
    216         endif
    217 
    218       print*, "Maximum number of iteration for the obl. parameter=", max_obl_iter
    219       print*, "Maximum number of iteration for the ex. parameter=", max_ex_iter
    220       print*, "Maximum number of iteration for the lsp. parameter=", max_lsp_iter
    221 
    222       year_iter_max = min(max_obl_iter,max_ex_iter,max_lsp_iter)
    223       if (year_iter_max > 0) then
    224         print*, "So the max. number of iteration (year) for the orbital parameter=", year_iter_max
    225       else
    226         write(*,*) 'The max. number of iteration (year) is not compatible with Laskar data in the file ob_ex_lsp.asc.'
    227         stop
    228       endif
    229 
    230         END SUBROUTINE orbit_param_criterion
     153        endif
     154        iylask = ilask
     155        exit ! The loop is left as soon as the right interval is found
     156    endif
     157enddo
     158
     159if (ilask == 1) then
     160    write(*,*) 'The year does not match with Laskar data in the file obl_ecc_lsp.asc.'
     161    stop
     162endif
     163
     164! Linear interpolation gives the maximum reachable year according to the limitation
     165! Obliquity
     166if (var_obl) then
     167    do ilask = iylask,2,-1
     168        ya = obllask(ilask)
     169        yb = obllask(ilask - 1)
     170        if (min(ya,yb) <= yc_obl .and. yc_obl < max(ya,yb)) then
     171            xa = yearlask(ilask)
     172            xb = yearlask(ilask - 1)
     173            max_obl_iter = floor((yc_obl - ya)*(xb - xa)/(yb - ya) + xa) - Year
     174            exit ! The loop is left as soon as the right interval is found
     175        endif
     176    enddo
     177endif
     178! Eccentricity
     179if (var_ecc) then
     180    do ilask = iylask,2,-1
     181        ya = ecclask(ilask)
     182        yb = ecclask(ilask - 1)
     183        if (min(ya,yb) <= yc_ecc .and. yc_ecc < max(ya,yb)) then
     184            xa = yearlask(ilask)
     185            xb = yearlask(ilask - 1)
     186            max_ecc_iter = floor((yc_ecc - ya)*(xb - xa)/(yb - ya) + xa) - Year
     187            exit ! The loop is left as soon as the right interval is found
     188        endif
     189    enddo
     190endif
     191! Lsp
     192if (var_lsp) then
     193    do ilask = iylask,2,-1
     194        ya = lsplask(ilask)
     195        yb = lsplask(ilask - 1)
     196        if (min(ya,yb) <= yc_lsp .and. yc_lsp < max(ya,yb)) then
     197            xa = yearlask(ilask)
     198            xb = yearlask(ilask - 1)
     199            max_lsp_iter = floor((yc_lsp - ya)*(xb - xa)/(yb - ya) + xa) - Year
     200            exit ! The loop is left as soon as the right interval is found
     201        endif
     202    enddo
     203endif
     204
     205write(*,*) 'Max. number of iterations for the obl. parameter =', max_obl_iter
     206write(*,*) 'Max. number of iterations for the ecc. parameter =', max_ecc_iter
     207write(*,*) 'Max. number of iterations for the Lsp parameter  =', max_lsp_iter
     208
     209year_iter_max = min(max_obl_iter,max_ecc_iter,max_lsp_iter)
     210if (year_iter_max > 0) then
     211    write(*,*) 'So the max. number of iterations for the orbital criteria =', year_iter_max
     212else
     213    write(*,*) 'The max. number of iterations is not compatible with Laskar data in the file obl_ecc_lsp.asc.'
     214    stop
     215endif
     216
     217END SUBROUTINE orbit_param_criterion
    231218
    232219!********************************************************************************   
    233      
    234       END MODULE orbit_param_criterion_mod
     220
     221END MODULE orbit_param_criterion_mod
     222

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

    r3032 r3039  
    5353                                      ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
    5454                                      co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
    55 use temps_mod_evol,             only: dt_pem, evol_orbit_pem, Max_iter_pem
     55use time_evol_mod,              only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
    5656use orbit_param_criterion_mod,  only: orbit_param_criterion
    5757use recomp_orb_param_mod,       only: recomp_orb_param
    … …  
    5959                                      ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
    6060use soil_thermalproperties_mod, only: update_soil_thermalproperties
     61use time_phylmdz_mod,           only: daysec, dtphys, day_end
    6162
    6263#ifndef CPP_STD
    … …  
    7980                                  PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, &
    8081                                  ALBEDO_CO2_ICE_SPECTV, phys_state_var_init
    81     use aerosol_mod,  only: iniaerosol
    82     use planete_mod,  only: apoastr, periastr, year_day, peri_day, obliquit
    83     use comcstfi_mod, only: r, mugaz
     82    use aerosol_mod,        only: iniaerosol
     83    use planete_mod,        only: apoastr, periastr, year_day, peri_day, obliquit
     84    use comcstfi_mod,       only: r, mugaz
    8485#endif
    8586
    8687#ifndef CPP_1D
    87     use iniphysiq_mod,            only: iniphysiq
    88     use control_mod,              only: iphysiq, day_step, nsplit_phys
     88    use iniphysiq_mod,    only: iniphysiq
     89    use control_mod,      only: iphysiq, day_step, nsplit_phys
    8990#else
    90     use time_phylmdz_mod,         only: iphysiq, day_step
     91    use time_phylmdz_mod, only: iphysiq, day_step
    9192#endif
    92 use time_phylmdz_mod, only: daysec,dtphys
     93
    9394#ifdef CPP_1D
    9495    use regular_lonlat_mod,       only: init_regular_lonlat
    … …  
    104105include "comgeom.h"
    105106include "iniprint.h"
     107include "callkeys.h"
    106108
    107109integer ngridmx
    … …  
    114116integer :: day_ini   ! First day of the simulation
    115117real    :: pday      ! Physical day
     118real    :: ptime     ! Physical time
    116119real    :: time_phys ! Same as GCM
    117120real    :: ptimestep ! Same as GCM
    … …  
    154157real                                :: global_ave_press_old  ! constant: Global average pressure of initial/previous time step
    155158real                                :: global_ave_press_new  ! constant: Global average pressure of current time step
    156 real, dimension(:,:),   allocatable ::  zplev_new            ! Physical x Atmospheric field : mass of the atmospheric  layers in the pem at current time step [kg/m^2]
    157 real, dimension(:,:),   allocatable ::  zplev_gcm            ! same but retrieved from the gcm [kg/m^2]
    158 real, dimension(:,:,:), allocatable ::  zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
     159real, dimension(:,:),   allocatable :: zplev_new             ! Physical x Atmospheric field : mass of the atmospheric  layers in the pem at current time step [kg/m^2]
     160real, dimension(:,:),   allocatable :: zplev_gcm             ! same but retrieved from the gcm [kg/m^2]
     161real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
    159162real, dimension(:,:,:), allocatable :: zplev_old_timeseries  ! same but with the time series, for oldest time step
    160163logical                             :: STOPPING_water        ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
    … …  
    185188real, dimension(:,:),   allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
    186189real, dimension(:,:),   allocatable :: tendencies_h2o_ice     ! physical pointx slope  field : Tendency of evolution of perenial h2o ice
    187 real, dimension(:,:),   allocatable :: flag_co2flow(:,:)      !(ngrid,nslope): Flag where there is a CO2 glacier flow
    188 real, dimension(:),     allocatable :: flag_co2flow_mesh(:)   !(ngrid)       : Flag where there is a CO2 glacier flow
    189 real, dimension(:,:),   allocatable :: flag_h2oflow(:,:)      !(ngrid,nslope): Flag where there is a H2O glacier flow
    190 real, dimension(:),     allocatable :: flag_h2oflow_mesh(:)   !(ngrid)       : Flag where there is a H2O glacier flow
     190real, dimension(:,:),   allocatable :: flag_co2flow(:,:)      ! (ngrid,nslope): Flag where there is a CO2 glacier flow
     191real, dimension(:),     allocatable :: flag_co2flow_mesh(:)   ! (ngrid)       : Flag where there is a CO2 glacier flow
     192real, dimension(:,:),   allocatable :: flag_h2oflow(:,:)      ! (ngrid,nslope): Flag where there is a H2O glacier flow
     193real, dimension(:),     allocatable :: flag_h2oflow_mesh(:)   ! (ngrid)       : Flag where there is a H2O glacier flow
    191194
    192195! Variables for surface and soil
    … …  
    194197real, allocatable :: tsoil_ave(:,:,:)            ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K]
    195198real, allocatable :: tsurf_GCM_timeseries(:,:,:) ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K]
    196 real, allocatable :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
    197 real, allocatable :: tsoil_GCM_timeseries(:,:,:,:)      !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
     199real, allocatable :: tsoil_phys_PEM_timeseries(:,:,:,:) ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
     200real, allocatable :: tsoil_GCM_timeseries(:,:,:,:)      ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
    198201real, allocatable :: tsurf_ave_yr1(:,:)                 ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
    199202real, allocatable :: TI_locslope(:,:)    ! Physic x Soil: Intermediate thermal inertia  to compute Tsoil [SI]
    … …  
    216219integer         :: year_iter     ! number of iteration
    217220integer         :: year_iter_max ! maximum number of iterations before stopping
     221integer         :: i_myear       ! Global number of Martian years of the chained simulations
     222integer         :: n_myear       ! Maximum number of Martian years of the chained simulations
    218223real            :: timestep      ! timestep [s]
    219224real            :: watercap_sum  ! total mass of water cap [kg/m^2]
    … …  
    247252
    248253! Loop variables
    249 integer :: i, j, ig0, l, ig, nnq, t, islope, ig_loop, islope_loop, iloop, isoil,icap
     254integer :: i, j, ig0, l, ig, nnq, t, islope, ig_loop, islope_loop, iloop, isoil, icap
    250255
    251256! Parallel variables
    … …  
    318323#endif
    319324
    320 call conf_pem
     325call conf_pem(i_myear,n_myear)
    321326
    322327!------------------------
    … …  
    449454    if (noms(nnq) == "co2") igcm_co2 = nnq
    450455enddo
    451 r = 8.314511E+0*1000.E+0/mugaz
     456r = 8.314511*1000./mugaz
    452457
    453458!------------------------
    … …  
    661666     call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
    662667#else
    663      call iniorbit(apoastr, periastr, year_day, peri_day,obliquit)
     668     call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
    664669#endif
    665670
    666671if (evol_orbit_pem) then
    667     call orbit_param_criterion(year_iter_max)
     672    call orbit_param_criterion(i_myear,year_iter_max)
    668673else
    669674    year_iter_max = Max_iter_pem
    … …  
    677682!------------------------
    678683year_iter = 0
    679 do while (year_iter < year_iter_max)
    680 
     684
     685do while (year_iter < year_iter_max .and. i_myear < n_myear)
    681686! II.a.1. Compute updated global pressure
    682687    write(*,*) "Recomputing the new pressure..."
    … …  
    911916
    912917        call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM, &
    913                                           delta_h2o_icetablesublim)        ! Mass of H2O exchange between the ssi and the atmosphere.
     918                                          delta_h2o_icetablesublim)        ! Mass of H2O exchange between the ssi and the atmosphere
    914919
    915920        write(*,*) "Update soil propreties"
    … …  
    962967
    963968    year_iter = year_iter + dt_pem
     969    i_myear = i_myear + dt_pem
    964970
    965971    write(*,*) "Checking all the stopping criterion."
    … …  
    984990        criterion_stop = 4
    985991    endif
     992    if (i_myear >= n_myear) then
     993        write(*,*) "STOPPING because maximum number of Martian years to be simulated reached"
     994        criterion_stop = 5
     995    endif
    986996
    987997    if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure)  then
    … …  
    989999    else
    9901000        write(*,*) "We continue!"
    991         write(*,*) "Number of iteration of the PEM : year_iter=", year_iter
     1001        write(*,*) "Number of iterations of the PEM: year_iter =", year_iter
     1002        write(*,*) "Number of simulated Martian years: i_myear =", i_myear
    9921003    endif
    9931004
    … …  
    11181129enddo
    11191130
    1120 if (evol_orbit_pem) call recomp_orb_param(year_iter)
     1131if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter)
    11211132
    11221133!------------------------
    … …  
    11261137! III_b.1 Write restart_evol.nc
    11271138ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys
     1139print*, ptimestep,dtphys,iphysiq*daysec/real(day_step),nsplit_phys
    11281140pday = day_ini
    1129 ztime_fin = 0.
     1141ptime = 0.
     1142ztime_fin = pday - day_end + ptime + ptimestep/(float(iphysiq)*daysec)
    11301143
    11311144allocate(p(ip1jmp1,nlayer + 1))
    … …  
    11331146    call pression (ip1jmp1,ap,bp,ps,p)
    11341147    call massdair(p,masse)
    1135     call dynredem0("restart_evol.nc", day_ini, phis)
    1136     call dynredem1("restart_evol.nc",   &
    1137     time_0,vcov,ucov,teta,q,masse,ps)
     1148    call dynredem0("restart_evol.nc",day_ini,phis)
     1149    call dynredem1("restart_evol.nc",time_0,vcov,ucov,teta,q,masse,ps)
    11381150    write(*,*) "restart_evol.nc has been written"
    11391151#else
    … …  
    11731185
    11741186
    1175 call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope ,          &
    1176              tsoil_PEM, TI_PEM, porefillingice_depth,porefillingice_thickness, &
     1187call pemdem1("restartfi_PEM.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
     1188             TI_PEM, porefillingice_depth,porefillingice_thickness,        &
    11771189             co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
    1178 
    1179 call info_run_PEM(year_iter,criterion_stop)
    1180 
    11811190write(*,*) "restartfi_PEM.nc has been written"
    1182 write(*,*) "The PEM had run for ", year_iter, " years."
    1183 write(*,*) "LL & RV : So far so good"
     1191
     1192call info_run_PEM(year_iter,criterion_stop,i_myear,n_myear)
     1193
     1194write(*,*) "The PEM has run for", year_iter, "Martian years."
     1195write(*,*) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
     1196write(*,*) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
     1197write(*,*) "LL & RV & JBC: so far, so good!"
    11841198
    11851199deallocate(vmr_co2_gcm,ps_timeseries,tsurf_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)

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

    r3019 r3039  
    1    SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness, &
     1SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness, &
    22                      tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
    33                      global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,  &
    44                      m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir)
    55   
    6    use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
    7    use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,fluxgeo,inertiedat_PEM,water_reservoir_nom,depth_breccia,depth_bedrock,index_breccia,index_bedrock
    8    use comsoil_h, only:  volcapa,inertiedat
    9    use adsorption_mod, only : regolith_adsorption,adsorption_pem
    10    USE temps_mod_evol, ONLY: year_PEM
    11    USE ice_table_mod, only: computeice_table_equilibrium
    12    USE constants_marspem_mod,only: alpha_clap_h2o,beta_clap_h2o, &
    13                                    TI_breccia,TI_bedrock
    14    use soil_thermalproperties_mod, only: update_soil_thermalproperties
    15    use tracer_mod, only: mmol,igcm_h2o_vap ! tracer names and molar masses
     6  use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
     7  use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, water_reservoir_nom, depth_breccia, depth_bedrock, index_breccia, index_bedrock
     8  use comsoil_h,                  only: volcapa,inertiedat
     9  use adsorption_mod,             only: regolith_adsorption, adsorption_pem
     10  use ice_table_mod,              only: computeice_table_equilibrium
     11  use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
     12  use soil_thermalproperties_mod, only: update_soil_thermalproperties
     13  use tracer_mod,                 only: mmol, igcm_h2o_vap ! tracer names and molar masses
    1614
    1715#ifndef CPP_STD
    18       USE comcstfi_h, only: r, mugaz
     16    use comcstfi_h,   only: r, mugaz
     17    use surfdat_h,    only: watercaptag
    1918#else
    20       USE comcstfi_mod, only: r, mugaz
     19    use comcstfi_mod, only: r, mugaz
    2120#endif
    2221
    23 
    24 #ifndef CPP_STD   
    25    use surfdat_h, only: watercaptag
    26 #endif
    27 
    2822  implicit none
    29  
     23
     24  include "callkeys.h"
     25
    3026  character(len=*), intent(in) :: filename                    ! name of the startfi_PEM.nc
    3127  integer,intent(in) :: ngrid                                 ! # of physical grid points
    … …  
    7369   real :: alph_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computation []
    7470   real :: beta_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computatio []                               
    75    real :: year_PEM_read                                              ! Year of the PEM previous run
    7671   LOGICAL :: startpem_file                                           ! boolean to check if we read the startfile or not
    7772#ifdef CPP_STD   
    … …  
    9792!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    9893
    99 
    10094!0. Check if the start_PEM exist.
    10195
    … …  
    110104   call open_startphy(filename)
    111105
    112    call get_var("Time",year_PEM_read,found)
    113    year_PEM=INT(year_PEM_read)
    114    if(.not.found) then
    115       write(*,*)'PEMetat0: failed loading year_PEM; take default=0'
    116    else
    117       write(*,*)'year_PEM of startpem=', year_PEM
    118    endif
    119 
    120     if(soil_pem) then
     106    if (soil_pem) then
    121107 
    122108!1. Thermal Inertia
    … …  
    161147ENDDO ! islope
    162148
    163 print *,'PEMETAT0: THERMAL INERTIA DONE'
     149write(*,*) 'PEMETAT0: THERMAL INERTIA DONE'
    164150
    165151! b. Special case for inertiedat, inertiedat_PEM
    … …  
    204190
    205191!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    206 
    207192!2. Soil Temperature
    208 
    209193DO islope=1,nslope
    210194  write(num,fmt='(i2.2)') islope
    … …  
    260244        enddo
    261245      enddo
    262 
    263    
    264246ENDDO ! islope
    265247
    266 print *,'PEMETAT0: SOIL TEMP  DONE'
    267 
    268 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    269 
     248write(*,*) 'PEMETAT0: SOIL TEMP  DONE'
     249
     250!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    270251!3. Ice Table
    271 
    272 
    273252  call get_field("ice_table",ice_table,found)
    274253   if(.not.found) then
    … …  
    282261   endif
    283262
    284 print *,'PEMETAT0: ICE TABLE  DONE'
    285 
    286 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    287 
     263write(*,*) 'PEMETAT0: ICE TABLE  DONE'
     264
     265!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    288266!4. CO2 & H2O Adsorption
    289 
    290267 if(adsorption_pem) then
    291268  DO islope=1,nslope
    … …  
    318295       deltam_h2o_regolith_phys(:) = 0. 
    319296    endif
    320 print *,'PEMETAT0: CO2 & H2O adsorption done  '
     297write(*,*) 'PEMETAT0: CO2 & H2O adsorption DONE'
    321298    endif
    322299endif ! soil_pem
    323300
    324 
    325 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    326 
     301!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    327302!. 5 water reservoir
    328 
    329303#ifndef CPP_STD   
    330304   call get_field("water_reservoir",water_reservoir,found)
    … …  
    342316#endif
    343317
    344 
    345318call close_startphy
    346319
    347320else !No startfi, let's build all by hand
    348321
    349     year_PEM=0
    350 
    351     if(soil_pem) then
     322    if (soil_pem) then
    352323
    353324!a) Thermal inertia
    … …  
    407378enddo
    408379
    409 
    410380!!! zone de transition
    411381delta = depth_bedrock
    … …  
    424394 enddo
    425395
    426 print *,'PEMETAT0: TI  DONE'
    427 
    428 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    429 
    430 
     396write(*,*) 'PEMETAT0: TI DONE'
     397
     398!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    431399!b) Soil temperature
    432400    do islope = 1,nslope
    … …  
    458426        enddo
    459427     enddo
    460  
    461       do isoil = nsoil_GCM+1,nsoil_PEM
    462         do ig = 1,ngrid
    463         watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
     428
     429        do isoil = nsoil_GCM+1,nsoil_PEM
     430          do ig = 1,ngrid
     431            watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
     432          enddo
    464433        enddo
    465       enddo
    466434enddo !islope
    467 print *,'PEMETAT0: TSOIL DONE  '
    468 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    469 !c) Ice table   
    470      call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
    471      call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
    472      do islope = 1,nslope
    473      call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    474      enddo
    475    
    476 
    477 
    478 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    479 
    480 print *,'PEMETAT0: Ice table DONE  '
    481 
    482 
     435write(*,*) 'PEMETAT0: TSOIL DONE'
     436
     437!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     438!c) Ice table
     439       call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
     440       call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
     441       do islope = 1,nslope
     442         call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
     443       enddo
     444       write(*,*) 'PEMETAT0: Ice table DONE'
    483445
    484446!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    485447!d) Regolith adsorbed
    486 
    487448 if(adsorption_pem) then
    488449    m_co2_regolith_phys(:,:,:) = 0.
    … …  
    496457 endif
    497458
    498 print *,'PEMETAT0: CO2 adsorption done  '
    499 
     459write(*,*) 'PEMETAT0: CO2 adsorption DONE'
    500460endif !soil_pem
    501461
    502 
    503 
    504 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    505    
    506 
    507 
     462!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    508463!. e) water reservoir
    509 
    510464#ifndef CPP_STD   
    511 
    512465      do ig=1,ngrid
    513466        if(watercaptag(ig)) then
    … …  
    517470        endif
    518471      enddo
    519 
    520472#endif
    521473
    … …  
    527479      DO islope = 1,nslope
    528480        DO iloop = 1,nsoil_PEM
    529            if(isnan(tsoil_PEM(ig,iloop,islope))) then
    530          call abort_pem("PEM - pemetat0","NAN detected in Tsoil",1)
    531            endif
     481          if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
    532482        ENDDO
    533483      ENDDO
    … …  
    535485endif!soil_pem
    536486
    537 
    538 
    539    END SUBROUTINE
     487END SUBROUTINE

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

    r2961 r3039  
    1313contains
    1414
    15 subroutine pemdem0(filename,lonfi,latfi,cell_area,nsoil_PEM,ngrid, &
    16                          day_ini,time,nslope,def_slope,subslope_dist)
     15subroutine pemdem0(filename,lonfi,latfi,cell_area,nsoil_PEM,ngrid,day_ini,time,nslope,def_slope,subslope_dist)
    1716! create physics restart file and write time-independent variables
    18   use comsoil_h_PEM, only: soil_pem,mlayer_PEM,fluxgeo,inertiedat_PEM
    19   use iostart_PEM, only : open_restartphy, close_restartphy, &
    20                       put_var, put_field, length
    21   use mod_grid_phy_lmdz, only : klon_glo
     17  use comsoil_h_PEM,     only: soil_pem,mlayer_PEM,fluxgeo,inertiedat_PEM
     18  use iostart_PEM,       only: open_restartphy, close_restartphy, put_var, put_field, length
     19  use mod_grid_phy_lmdz, only: klon_glo
     20  use time_phylmdz_mod,  only: daysec
    2221#ifndef CPP_STD
    23   use planete_h, only: aphelie, emin_turb, lmixmin, obliquit, &
    24                        peri_day, periheli, year_day
    25   use comcstfi_h, only: g, mugaz, omeg, rad, rcp
     22  use planete_h,    only: aphelie, emin_turb, lmixmin, obliquit, peri_day, periheli, year_day
     23  use comcstfi_h,   only: g, mugaz, omeg, rad, rcp
    2624#else
    27   use planete_mod, only: apoastr, emin_turb, lmixmin, obliquit, &
    28                        peri_day, periastr, year_day
    29   USE comconst_mod, ONLY: g, omeg, rad
     25  use planete_mod,  only: apoastr, emin_turb, lmixmin, obliquit, peri_day, periastr, year_day
     26  use comconst_mod, only: g, omeg, rad
    3027  use comcstfi_mod, only: mugaz, rcp
    3128#endif
    32   use time_phylmdz_mod, only: daysec
     29
    3330  implicit none
    3431 
    3532  character(len=*), intent(in) :: filename
    36   real,intent(in) :: lonfi(ngrid)
    37   real,intent(in) :: latfi(ngrid)
    38   integer,intent(in) :: nsoil_PEM
    39   integer,intent(in) :: ngrid
    40   real,intent(in) :: day_ini
    41   real,intent(in) :: time
    42   real, intent(in) :: cell_area(ngrid) !boundaries for bining of the slopes
    43   integer,intent(in) :: nslope
    44   real, intent(in) :: def_slope(nslope+1) !boundaries for bining of the slopes
    45   real, intent(in) :: subslope_dist(ngrid,nslope) !undermesh statistics
     33  real, intent(in)            :: lonfi(ngrid)
     34  real, intent(in)            :: latfi(ngrid)
     35  integer, intent(in)          :: nsoil_PEM
     36  integer, intent(in)          :: ngrid
     37  real, intent(in)            :: day_ini
     38  real, intent(in)            :: time
     39  real, intent(in)             :: cell_area(ngrid) !boundaries for bining of the slopes
     40  integer, intent(in)          :: nslope
     41  real, intent(in)             :: def_slope(nslope+1) !boundaries for bining of the slopes
     42  real, intent(in)             :: subslope_dist(ngrid,nslope) !undermesh statistics
    4643 
    4744  ! Create physics start file
    … …  
    7067end subroutine pemdem0
    7168
    72 subroutine pemdem1(filename,year_iter,nsoil_PEM,ngrid,nslope, &
     69subroutine pemdem1(filename,i_myear,nsoil_PEM,ngrid,nslope, &
    7370                    tsoil_slope_PEM,inertiesoil_slope_PEM,ice_table_depth,ice_table_thickness, &
    7471                    m_co2_regolith,m_h2o_regolith,water_reservoir)
    … …  
    7673
    7774  ! write time-dependent variable to restart file
    78   use iostart_PEM, only : open_restartphy, close_restartphy, &
    79                       put_var, put_field
    80 
    81   use comsoil_h_PEM, only: mlayer_PEM,fluxgeo,inertiedat_PEM,soil_pem
    82   USE temps_mod_evol, ONLY: year_PEM
     75  use iostart_PEM,   only: open_restartphy, close_restartphy, put_var, put_field
     76  use comsoil_h_PEM, only: mlayer_PEM,fluxgeo, inertiedat_PEM,soil_pem
     77  use time_evol_mod, only: year_bp_ini, convert_years
    8378               
    8479  implicit none
    … …  
    8883#endif
    8984 
    90   character(len=*),intent(in) :: filename
    91   integer,intent(in) :: nsoil_PEM
    92   integer,intent(in) :: ngrid
    93   real,intent(IN) :: tsoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
    94   real,intent(IN) :: inertiesoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
    95   real,intent(IN) :: ice_table_depth(ngrid,nslope) !under mesh bining according to slope
    96   real,intent(IN) :: ice_table_thickness(ngrid,nslope) !under mesh bining according to slope
    97   integer,intent(IN) :: nslope
    98   real,intent(in) :: m_co2_regolith(ngrid,nsoil_PEM,nslope)
    99   real,intent(in) :: m_h2o_regolith(ngrid,nsoil_PEM,nslope)
    100   real,intent(in) :: water_reservoir(ngrid)
    101   integer :: islope
    102   CHARACTER*2 :: num 
    103   integer, intent(IN) :: year_iter
     85  character(len=*), intent(in) :: filename
     86  integer, intent(in)          :: nsoil_PEM, ngrid, nslope, i_myear
     87  real, intent(in)             :: tsoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
     88  real, intent(in)             :: inertiesoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
     89  real, intent(in)             :: ice_table_depth(ngrid,nslope) !under mesh bining according to slope
     90  real, intent(in)             :: ice_table_thickness(ngrid,nslope) !under mesh bining according to slope
     91  real, intent(in)             :: m_co2_regolith(ngrid,nsoil_PEM,nslope)
     92  real, intent(in)             :: m_h2o_regolith(ngrid,nsoil_PEM,nslope)
     93  real, intent(in)             :: water_reservoir(ngrid)
    10494
    105   real :: year_tot
    106  
    107   integer :: iq
    108   character(len=30) :: txt ! to store some text
    109   ! indexes of water vapour & water ice tracers (if any):
    110 
    111   year_tot=real(year_iter+year_PEM)
    112   print *, "Year of simulation=", year_tot 
     95  integer           :: islope
     96  character*2       :: num
     97  integer           :: iq
     98  character(len=30) :: txt  ! To store some text
     99  real              :: Year ! Year of the simulation
    113100
    114101  ! Open file
    … …  
    117104  ! First variable to write must be "Time", in order to correctly
    118105  ! set time counter in file
    119   call put_var("Time","Year of simulation",year_tot)
    120   call put_field('water_reservoir','water_reservoir', &
    121            water_reservoir,year_tot)
     106  Year = (year_bp_ini + i_myear)*convert_years
     107  call put_var("Time","Year of simulation",Year)
     108  call put_field('water_reservoir','water_reservoir',water_reservoir,Year)
     109
    122110    if(soil_pem) then
    123111 
    124112  ! Multidimensionnal variables (undermesh slope statistics)
    125 DO islope=1,nslope
     113do islope = 1,nslope
    126114  write(num,fmt='(i2.2)') islope
    127115  call put_field("tsoil_PEM_slope"//num,"Soil temperature by slope type", &
    128                  tsoil_slope_PEM(:,:,islope),year_tot)
     116                 tsoil_slope_PEM(:,:,islope),Year)
    129117  call put_field("TI_PEM_slope"//num,"Soil Thermal Inertia by slope type", &
    130                  inertiesoil_slope_PEM(:,:,islope),year_tot)
     118                 inertiesoil_slope_PEM(:,:,islope),Year)
    131119
    132  call put_field("mco2_reg_ads_slope"//num, "Mass of co2 adsorbded in the regolith", &
    133                     m_co2_regolith(:,:,islope), year_tot)
     120  call put_field("mco2_reg_ads_slope"//num, "Mass of co2 adsorbded in the regolith", &
     121                m_co2_regolith(:,:,islope),Year)
    134122  call put_field("mh2o_reg_ads_slope"//num, "Mass of h2o adsorbded in the regolith", &
    135                     m_h2o_regolith(:,:,islope), year_tot)
     123                 m_h2o_regolith(:,:,islope),Year)
     124enddo
    136125
    137 ENDDO
     126  call put_field("ice_table_depth","Depth of ice table",ice_table_depth,Year)
    138127
    139   call put_field("ice_table_depth","Depth of ice table", &
    140                  ice_table_depth,year_tot)
     128  call put_field("ice_table_thickness","Depth of ice table",ice_table_thickness,Year)
    141129
    142   call put_field("ice_table_thickness","Depth of ice table", &
    143                  ice_table_thickness,year_tot)
    144 
    145   call put_field("inertiedat_PEM","Thermal inertie of PEM ", &
    146                  inertiedat_PEM,year_tot)
     130  call put_field("inertiedat_PEM","Thermal inertie of PEM ",inertiedat_PEM,Year)
    147131
    148132   endif !soil_pem

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

    r3023 r3039  
    1       MODULE recomp_orb_param_mod
     1MODULE recomp_orb_param_mod
    22!=======================================================================
    33!
    … …  
    66! Author: RV, JBC
    77!=======================================================================
    8       IMPLICIT NONE
     8IMPLICIT NONE
    99
    10       CONTAINS
     10CONTAINS
    1111
    12       SUBROUTINE recomp_orb_param(final_iter)
     12SUBROUTINE recomp_orb_param(i_myear,year_iter)
    1313
    14       USE temps_mod_evol, ONLY: year_bp_ini, year_PEM, var_obl, var_ex, var_lsp
     14  use time_evol_mod, only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
    1515#ifndef CPP_STD     
    16       USE comconst_mod, ONLY: pi
    17       USE planete_h, ONLY: e_elips, obliquit, timeperi, periheli,aphelie,p_elips,peri_day,year_day
     16    use comconst_mod, only: pi
     17    use planete_h,    only: e_elips, obliquit, timeperi, periheli, aphelie, p_elips, peri_day, year_day
    1818#else
    19       use planete_mod, only: e_elips, obliquit, timeperi,periastr,apoastr,p_elips,peri_day,year_day
    20       USE comcstfi_mod, only: pi
    21 
     19    use planete_mod,  only: e_elips, obliquit, timeperi, periastr, apoastr, p_elips, peri_day, year_day
     20    use comcstfi_mod, only: pi
    2221#endif
    2322
    24       USE lask_param_mod, only: yearlask,oblask,exlask,lsplask, &
    25                                 end_lask_param_mod,last_ilask
     23  use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, end_lask_param_mod, last_ilask
    2624
    27       IMPLICIT NONE
     25  IMPLICIT NONE
    2826
    2927!--------------------------------------------------------
    3028! Input Variables
    3129!--------------------------------------------------------
    32 
    33       integer,intent(in) :: final_iter ! Number of year iteration of the PEM
     30  integer, intent(in) :: i_myear   ! Number of simulated Martian years
     31  integer, intent(in) :: year_iter ! Number of iterations of the PEM
    3432
    3533!--------------------------------------------------------
    … …  
    4038! Local variables
    4139!--------------------------------------------------------
     40  integer :: Year           ! Year of the simulation
     41  real    :: Lsp            ! Ls perihelion
     42  integer :: ilask          ! Loop variable
     43  real    :: halfaxe        ! Million of km
     44  real    :: unitastr
     45  real    :: xa, xb, ya, yb ! Variables for interpolation
    4246
    43       real :: Year     ! Year of the simulation
    44       real :: Lsp      ! time of peri in ls
    45       integer ilask    ! Loop variable
    46       real :: halfaxe  ! Million of km
    47       real :: unitastr
     47! **********************************************************************
     48! 0. Initializations
     49! **********************************************************************
     50#ifdef CPP_STD
     51    real aphelie
     52    real periheli
    4853
    49       real xa,xb,ya,yb
    50 
    51 
    52       ! **********************************************************************
    53       ! 0. Initializations
    54       ! **********************************************************************
    55 #ifdef CPP_STD
    56       real aphelie
    57       real periheli
    58 
    59       aphelie=apoastr
    60       periheli=periastr
     54    aphelie=apoastr
     55    periheli=periastr
    6156#endif
    6257
    63           Year = year_bp_ini + year_PEM + final_iter ! We compute the current year
    64           Lsp = 360. - timeperi*360./(2.*pi)        ! We convert in degree
     58Year = year_bp_ini + i_myear        ! We compute the current year
     59Lsp = modulo(-timeperi*180./pi,360.) ! We convert in degree
    6560
    66 !          Year=-953200
    67           print*, "recomp_orb_param, Old year=", year_bp_ini+year_PEM
    68           print*, "recomp_orb_param, New year=", Year
    69           print*, "recomp_orb_param, Old obl=", obliquit
    70           print*, "recomp_orb_param, Old ex=", e_elips
    71           print*, "recomp_orb_param, Old lsp=", Lsp
    72           print*, "recomp_orb_param, Old timeperi=", timeperi
     61write(*,*) 'recomp_orb_param, Old year =', Year - year_iter
     62write(*,*) 'recomp_orb_param, Old obl. =', obliquit
     63write(*,*) 'recomp_orb_param, Old ecc. =', e_elips
     64write(*,*) 'recomp_orb_param, Old Lsp  =', Lsp
    7365
    74         do ilask = last_ilask,2,-1
    75           xa = yearlask(ilask)
    76           xb = yearlask(ilask - 1)
    77           if (xa <= Year .and. Year < xb) then
    78             ! Obliquity
    79             if (var_obl) then
    80               ya = oblask(ilask)
    81               yb = oblask(ilask - 1)
    82               obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya
     66do ilask = last_ilask,2,-1
     67    xa = yearlask(ilask)
     68    xb = yearlask(ilask - 1)
     69    if (xa <= Year .and. Year < xb) then
     70        ! Obliquity
     71        if (var_obl) then
     72            ya = obllask(ilask)
     73            yb = obllask(ilask - 1)
     74            obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya
     75        endif
     76
     77        ! Eccentricity
     78        if (var_ecc) then
     79            ya = ecclask(ilask)
     80            yb = ecclask(ilask - 1)
     81            e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya
     82        endif
     83
     84        ! Lsp
     85        if (var_lsp) then
     86            ya = lsplask(ilask)
     87            yb = lsplask(ilask - 1)
     88            if (yb - ya > 300.) then ! If modulo is "crossed"
     89                if (ya < yb) then ! Increasing
     90                    yb = yb - 360.
     91                else ! Decreasing
     92                    yb = yb + 360.
     93                endif
    8394            endif
     95            Lsp = (Year - xa)*(yb - ya)/(xb - xa) + ya
     96            Lsp = modulo(Lsp,360.)
     97        endif
     98        exit ! The loop is left as soon as the right interval is found
     99    endif
     100enddo
     101halfaxe = 227.94
     102timeperi = Lsp*pi/180.
     103periheli = halfaxe*(1 - e_elips)
     104aphelie =  halfaxe*(1 + e_elips)
     105unitastr = 149.597927
     106p_elips = 0.5*(periheli + aphelie)*(1. - e_elips*e_elips)/unitastr
    84107
    85             ! Excentricity
    86             if (var_ex) then
    87               ya = exlask(ilask)
    88               yb = exlask(ilask - 1)
    89               e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya
    90             endif
    91 
    92             ! Lsp
    93             if (var_lsp) then
    94               ya = lsplask(ilask)
    95               yb = lsplask(ilask - 1)
    96               if (yb - ya > 300.) then ! If modulo is "crossed"
    97                 if (ya < yb) then ! Increasing
    98                   yb = yb - 360.
    99                 else ! Decreasing
    100                   yb = yb + 360.
    101                 endif
    102               endif
    103               Lsp = (Year - xa)*(yb - ya)/(xb - xa) + ya
    104               Lsp = modulo(Lsp,360.)
    105             endif
    106             exit ! The loop is left as soon as the right interval is found
    107           endif
    108         enddo
    109         halfaxe = 227.94
    110         timeperi = Lsp*2.*pi/360.
    111         periheli = halfaxe*(1 - e_elips)
    112         aphelie =  halfaxe*(1 + e_elips)
    113         unitastr = 149.597927
    114         p_elips = 0.5*(periheli + aphelie)*(1. - e_elips*e_elips)/unitastr
    115108#ifndef CPP_STD 
    116         call call_dayperi(timeperi,e_elips,peri_day,year_day)
     109    call call_dayperi(timeperi,e_elips,peri_day,year_day)
    117110#endif
    118111
    119        print*, "recomp_orb_param, Final year of the PEM run:", Year
    120        print*, "recomp_orb_param, New obl=", obliquit
    121        print*, "recomp_orb_param, New ex=", e_elips
    122        print*, "recomp_orb_param, New timeperi=", timeperi
     112write(*,*) 'recomp_orb_param, New year =', Year
     113write(*,*) 'recomp_orb_param, New obl. =', obliquit
     114write(*,*) 'recomp_orb_param, New ecc. =', e_elips
     115write(*,*) 'recomp_orb_param, New Lsp  =', Lsp
    123116
    124       END SUBROUTINE recomp_orb_param
     117END SUBROUTINE recomp_orb_param
    125118
    126119!********************************************************************************   
    127120     
    128       END MODULE recomp_orb_param_mod
     121END MODULE recomp_orb_param_mod
     122

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

    r3038 r3039  
    1 MODULE temps_mod_evol
     1MODULE time_evol_mod
    22
    33IMPLICIT NONE 
    44
    5   INTEGER   year_bp_ini         !     year_bp_ini : Initial year of the simulation of the PEM (in evol.def)
    6   INTEGER   dt_pem              !     dt_pem  : in Planetary years, the time step used by the PEM   
    7   REAL      water_ice_criterion !     Percentage of change of the surface of water ice sublimating before stopping the PEM
    8   REAL      co2_ice_criterion   !     Percentage of change of the surface of co2 ice sublimating before stopping the PEM
    9   REAL      ps_criterion        !     Percentage of change of averaged surface pressure before stopping the PEM
    10   INTEGER   year_PEM            !     year written in startfiPEM.nc
    11   INTEGER  Max_iter_pem        !     Maximal number of iteration when converging to a steady state, read in evol.def
    12   LOGICAL   evol_orbit_pem      !     True if we want to follow the orbital parameters of ob_ex_lsp.asc, read in evol.def
    13   LOGICAL   var_obl             !     True if we want the PEM to follow ob_ex_lsp.asc parameters for obliquity
    14   LOGICAL   var_ex              !     True if we want the PEM to follow ob_ex_lsp.asc parameters for excenticity
    15   LOGICAL   var_lsp             !     True if we want the PEM to follow ob_ex_lsp.asc parameters for ls perihelie
     5  integer :: year_bp_ini         !     Initial year (in Planetary years) of the simulation of the PEM defined in run.def (in Earth years)
     6  integer :: dt_pem              !     Time step used by the PEM in Planetary years
     7  real    :: convert_years       !     Conversion ratio from Planetary years to Earth years
     8  real    :: water_ice_criterion !     Percentage of change of the surface of water ice sublimating before stopping the PEM
     9  real    :: co2_ice_criterion   !     Percentage of change of the surface of co2 ice sublimating before stopping the PEM
     10  real    :: ps_criterion        !     Percentage of change of averaged surface pressure before stopping the PEM
     11  integer :: Max_iter_pem        !     Maximal number of iteration when converging to a steady state, read in evol.def
     12  logical :: evol_orbit_pem      !     True if we want to follow the orbital parameters of obl_ecc_lsp.asc, read in evol.def
     13  logical :: var_obl             !     True if we want the PEM to follow obl_ecc_lsp.asc parameters for obliquity
     14  logical :: var_ecc             !     True if we want the PEM to follow obl_ecc_lsp.asc parameters for excenticity
     15  logical :: var_lsp             !     True if we want the PEM to follow obl_ecc_lsp.asc parameters for ls perihelie
    1616
    17 END MODULE temps_mod_evol
     17END MODULE time_evol_mod

  • trunk/LMDZ.MARS/changelog.txt

    r3038 r3039  
    41854185* Some changes in run_PEM.def and xml files.
    41864186JBC
     4187
     4188
     4189== 11/09/2023 == JBC
     4190* New management of time in PEM. While definitions in "run.def"/"launching script" and Laskar's data are in Earth years, the PEM works in Martian years. It follows several changes and the new variable 'convert_years' to make the conversion;
     4191* New parameter for the years to be simulated. The user does not ask anymore for a number of PEM iterations to do but for a number of Earth years to simulate. The GCM years are now counted! To do so, a temporary file "tmp_PEMyears.txt" gives few basic information between the runs;
     4192* New values for the maximal admissible change of orbital parameters. They are now coherent and allows the PEM to run for ~1000 Martian years;
     4193* Laskar's data interpolation has been cleaned further;
     4194* Some cleaning and renaming of PEM subroutines in the course of the previous modifications.
Note: See TracChangeset for help on using the changeset viewer.