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

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

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

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

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r3991 r4065  
    1313!-----------------------------------------------------------------------
    1414
    15 ! DECLARATION
    16 ! -----------
    17 implicit none
    18 
    19 ! MODULE VARIABLES
    20 ! ----------------
    21 real, dimension(:), allocatable :: year_lask  ! year before present from Laskar et al.
    22 real, dimension(:), allocatable :: obl_lask   ! obliquity [deg]
    23 real, dimension(:), allocatable :: ecc_lask   ! eccentricity
    24 real, dimension(:), allocatable :: lsp_lask   ! ls perihelie [deg]
    25 integer                         :: iyear_lask ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
     15! DEPENDENCIES
     16! ------------
     17use numerics, only: dp, di, k4, largest_nb, minieps
     18
     19! DECLARATION
     20! -----------
     21implicit none
     22
     23! PARAMETERS
     24! ----------
     25! File-related:
     26character(15),                       parameter :: orbitfile_name = 'obl_ecc_lsp.asc'
     27integer(di),                         protected :: iyear_lask! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
     28real(dp), dimension(:), allocatable, protected :: year_lask ! Current_year before present from Laskar et al.
     29real(dp), dimension(:), allocatable, protected :: obl_lask  ! Obliquity in Laskar's data [deg]
     30real(dp), dimension(:), allocatable, protected :: ecc_lask  ! Eccentricity in Laskar's data
     31real(dp), dimension(:), allocatable, protected :: lsp_lask  ! Ls perihelie in Laskar's data [deg]
     32! Planet-related:
     33real(dp),                            protected :: obliquity    ! Planet obliquity [deg]
     34real(dp),                            protected :: eccentricity ! Orbit eccentricity
     35real(dp),                            protected :: ls_peri      ! Solar longitude of the perihelion [deg]
     36real(dp),                            protected :: perihelion   ! Perihelion [Mkm]
     37real(dp),                            protected :: semimaj_axis ! Semi-major axis [Mkm]
     38real(dp),                            protected :: aphelion     ! Aphelion [Mkm]
     39real(dp),                            protected :: date_peri    ! Date of perihelion [sol]
     40real(dp),                            protected :: sol_len_s    ! Length of a sol [s]
     41real(dp),                            protected :: yr_len_sol   ! Length of a year [sol]
     42! Evolution-related:
     43logical(k4), protected :: evol_orbit     ! Flag to follow the orbital parameters
     44logical(k4), protected :: evol_obl       ! Flag the PEM to follow obliquity
     45logical(k4), protected :: evol_ecc       ! Flag the PEM to follow eccenticity
     46logical(k4), protected :: evol_lsp       ! Flag the PEM to follow Ls perihelie
     47real(dp),    protected :: max_change_obl ! Maximum admissible change for obliquity
     48real(dp),    protected :: max_change_ecc ! Maximum admissible change for eccentricity
     49real(dp),    protected :: max_change_lsp ! Maximum admissible change for Lsp
    2650
    2751contains
    … …  
    2953
    3054!=======================================================================
    31 SUBROUTINE read_orbitpm(Year)
    32 !-----------------------------------------------------------------------
    33 ! NAME
    34 !     read_orbitpm
    35 !
    36 ! DESCRIPTION
    37 !     Read Laskar orbital file and prepare arrays; locate index for
    38 !     closest lower year relative to current simulation year.
     55SUBROUTINE set_orbit_config(evol_orbit_in,evol_obl_in,evol_ecc_in,evol_lsp_in,max_change_obl_in,max_change_ecc_in,max_change_lsp_in)
     56!-----------------------------------------------------------------------
     57! NAME
     58!     set_orbit_config
     59!
     60! DESCRIPTION
     61!     Setter for 'orbit' configuration parameters.
     62!
     63! AUTHORS & DATE
     64!     JB Clement, 02/2026
     65!
     66! NOTES
     67!
     68!-----------------------------------------------------------------------
     69
     70! DEPENDENCIES
     71! ------------
     72use utility, only: real2str, bool2str
     73use display, only: print_msg
     74
     75! DECLARATION
     76! -----------
     77implicit none
     78
     79! ARGUMENTS
     80! ---------
     81logical(k4), intent(in) :: evol_orbit_in, evol_obl_in, evol_ecc_in, evol_lsp_in
     82real(dp),    intent(in) :: max_change_obl_in, max_change_ecc_in, max_change_lsp_in
     83
     84! CODE
     85! ----
     86evol_orbit = evol_orbit_in
     87evol_obl = evol_obl_in
     88evol_ecc = evol_ecc_in
     89evol_lsp = evol_lsp_in
     90max_change_obl = max_change_obl_in
     91max_change_ecc = max_change_ecc_in
     92max_change_lsp = max_change_lsp_in
     93call print_msg('evol_orbit     = '//bool2str(evol_orbit))
     94call print_msg('evol_obl       = '//bool2str(evol_obl))
     95call print_msg('evol_ecc       = '//bool2str(evol_ecc))
     96call print_msg('evol_lsp       = '//bool2str(evol_lsp))
     97call print_msg('max_change_obl = '//real2str(max_change_obl))
     98call print_msg('max_change_ecc = '//real2str(max_change_ecc))
     99call print_msg('max_change_lsp = '//real2str(max_change_lsp))
     100
     101END SUBROUTINE set_orbit_config
     102!=======================================================================
     103
     104!=======================================================================
     105SUBROUTINE init_orbit(obliquity_in,perihelion_in,aphelion_in,date_peri_in,yr_len_sol_in,sol_len_s_in)
     106!-----------------------------------------------------------------------
     107! NAME
     108!     init_orbit
     109!
     110! DESCRIPTION
     111!     Initialize orbital characteristics of the planet.
     112!
     113! AUTHORS & DATE
     114!     JB Clement, 12/2025
     115!
     116! NOTES
     117!     Taken from 'iniorbit' in the Mars PCM.
     118!-----------------------------------------------------------------------
     119
     120! DEPENDENCIES
     121! ------------
     122use maths,   only: pi
     123use display, only: print_msg
     124use utility, only: real2str
     125
     126! DECLARATION
     127! -----------
     128implicit none
     129
     130! ARGUMENTS
     131! ---------
     132real(dp), intent(in) :: obliquity_in, perihelion_in, aphelion_in, date_peri_in, yr_len_sol_in, sol_len_s_in
     133
     134! LOCAL VARIABLES
     135! ---------------
     136real(dp)    :: zxref, zanom, zz, zx0, zdx
     137real(dp)    :: lsp ! Solar longitude of the perihelion [rad]
     138integer(di) :: iter
     139
     140! CODE
     141! ----
     142! Setters
     143obliquity = obliquity_in
     144perihelion = perihelion_in
     145aphelion = aphelion_in
     146date_peri = date_peri_in
     147sol_len_s = sol_len_s_in
     148yr_len_sol = yr_len_sol_in
     149
     150! Compute eccentricity and semi-major axis
     151eccentricity = (aphelion - perihelion)/(aphelion + perihelion)
     152semimaj_axis = 0.5_dp*(aphelion + perihelion)
     153
     154! Compute mean anomaly zanom
     155zz = (yr_len_sol - date_peri)/yr_len_sol
     156zanom = 2._dp*pi*(zz - nint(zz))
     157zxref = abs(zanom)
     158
     159! Solve equation zx0 - e*sin(zx0) = zxref for eccentric anomaly zx0 using Newton method
     160zx0 = zxref + eccentricity*sin(zxref)
     161do iter = 1,100
     162    zdx = -(zx0 - eccentricity*sin(zx0) - zxref)/(1._dp - eccentricity*cos(zx0))
     163    if (abs(zdx) <= 1.e-12_dp) exit
     164    zx0 = zx0 + zdx
     165end do
     166
     167zx0 = zx0 + zdx
     168if (zanom < 0.) zx0 = -zx0
     169
     170! Compute solar longitude of the perihelion
     171lsp = -2._dp*atan(sqrt((1._dp + eccentricity)/(1._dp - eccentricity))*tan(zx0/2._dp))
     172lsp = modulo(lsp,2._dp*pi)
     173ls_peri = lsp*180._dp/pi ! Conversion in degree
     174
     175call print_msg('Obliquity          [deg] = '//real2str(obliquity))
     176call print_msg('Eccentricity       [-]   = '//real2str(eccentricity))
     177call print_msg('Ls of perihelion   [deg] = '//real2str(ls_peri))
     178call print_msg('Perihelion         [Mkm] = '//real2str(perihelion))
     179call print_msg('Aphelion           [Mkm] = '//real2str(aphelion))
     180call print_msg('Semi-major axis    [Mkm] = '//real2str(semimaj_axis))
     181call print_msg('Date of perihelion [sol] = '//real2str(date_peri))
     182call print_msg('Length of a sol    [s]   = '//real2str(sol_len_s))
     183call print_msg('Length of a year   [sol] = '//real2str(yr_len_sol))
     184
     185END SUBROUTINE init_orbit
     186!=======================================================================
     187
     188!=======================================================================
     189SUBROUTINE lsp2datep(lsp,datep)
     190!-----------------------------------------------------------------------
     191! NAME
     192!     lsp2datep
     193!
     194! DESCRIPTION
     195!     Initialize orbital data of the planet.
     196!
     197! AUTHORS & DATE
     198!     JB Clement, 12/2025
     199!
     200! NOTES
     201!     Taken from 'lsp2solp' in the Mars PCM.
     202!-----------------------------------------------------------------------
     203
     204! DEPENDENCIES
     205! ------------
     206use maths, only: pi
     207
     208! DECLARATION
     209! -----------
     210implicit none
     211
     212! ARGUMENTS
     213! ---------
     214real(dp), intent(in)  :: lsp   ! Ls at perihelion [deg]
     215real(dp), intent(out) :: datep ! Date at perihelion [sol]
     216
     217! LOCAL VARIABLES
     218! ---------------
     219real(dp) :: zx0 ! Eccentric anomaly at Ls = 0
     220
     221! CODE
     222! ----
     223! Compute orbit geometrical parameters
     224zx0 = -2._dp*atan(tan(0.5_dp*lsp*pi/180._dp)*sqrt((1._dp - eccentricity)/(1._dp + eccentricity)))
     225if (zx0 <= 0._dp) zx0 = zx0 + 2._dp*pi
     226
     227! Compute sol at perihelion
     228datep = yr_len_sol*(1._dp - (zx0 - eccentricity*sin(zx0))/(2._dp*pi))
     229
     230END SUBROUTINE lsp2datep
     231!=======================================================================
     232
     233!=======================================================================
     234SUBROUTINE ini_orbit()
     235!-----------------------------------------------------------------------
     236! NAME
     237!     ini_orbit
     238!
     239! DESCRIPTION
     240!     Initialize the parameters of module 'orbit'.
    39241!
    40242! AUTHORS & DATE
    … …  
    43245!
    44246! NOTES
    45 !     Converts kilo Earth years to Martian years using 'convert_years'.
     247!
    46248!-----------------------------------------------------------------------
    47249
    48250! DEPENDENCIES
    49251! ------------
    50 use evolution, only: convert_years
     252use stoppage, only: stop_clean
     253
     254! DECLARATION
     255! -----------
     256implicit none
     257
     258! LOCAL VARIABLES
     259! ---------------
     260integer(di) :: n ! Number of lines in the file
     261integer(di) :: ierr, funit
     262logical(k4) :: here
     263
     264! CODE
     265! ----
     266inquire(file = orbitfile_name,exist = here)
     267if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//orbitfile_name//'"!',1)
     268
     269! Get the number of lines for Laskar's orbital file
     270open(newunit = funit,file = orbitfile_name,status = 'old',action = 'read',iostat = ierr)
     271if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'error opening file "'//orbitfile_name//'"!',ierr)
     272n = 0
     273do
     274    read(funit,*,iostat = ierr)
     275    if (ierr /= 0) exit
     276    n = n + 1
     277end do
     278close(funit)
     279
     280! Allocation of module arrays
     281if (.not. allocated(year_lask)) allocate(year_lask(n))
     282if (.not. allocated(obl_lask)) allocate(obl_lask(n))
     283if (.not. allocated(ecc_lask)) allocate(ecc_lask(n))
     284if (.not. allocated(lsp_lask)) allocate(lsp_lask(n))
     285
     286END SUBROUTINE ini_orbit
     287!=======================================================================
     288
     289!=======================================================================
     290SUBROUTINE end_orbit
     291!-----------------------------------------------------------------------
     292! NAME
     293!     end_orbit
     294!
     295! DESCRIPTION
     296!     Deallocate orbital data arrays.
     297!
     298! AUTHORS & DATE
     299!     R. Vandemeulebrouck
     300!     JB Clement, 2023-2025
     301!
     302! NOTES
     303!
     304!-----------------------------------------------------------------------
     305
     306! DECLARATION
     307! -----------
     308implicit none
     309
     310! CODE
     311! ----
     312if (allocated(year_lask)) deallocate(year_lask)
     313if (allocated(obl_lask)) deallocate(obl_lask)
     314if (allocated(ecc_lask)) deallocate(ecc_lask)
     315if (allocated(lsp_lask)) deallocate(lsp_lask)
     316
     317END SUBROUTINE end_orbit
     318!=======================================================================
     319
     320!=======================================================================
     321SUBROUTINE read_orbitpm(n_yr_sim)
     322!-----------------------------------------------------------------------
     323! NAME
     324!     read_orbitpm
     325!
     326! DESCRIPTION
     327!     Read Laskar's orbital file and locate index for
     328!     closest lower year relative to current simulation year.
     329!
     330! AUTHORS & DATE
     331!     JB Clement, 2023-2025
     332!
     333! NOTES
     334!
     335!-----------------------------------------------------------------------
     336
     337! DEPENDENCIES
     338! ------------
     339use evolution, only: pem_ini_date, r_plnt2earth_yr
     340use stoppage,  only: stop_clean
     341use display,   only: print_msg
     342use utility,   only: real2str, int2str
    51343
    52344! DECLARATION
    … …  
    56348! ARGUMENTS
    57349! ---------
    58 real, intent(in) :: Year ! Martian year of the simulation
     350real(dp), intent(in) :: n_yr_sim
    59351
    60352! LOCAL VARIABLES
    61353! ---------------
    62 integer :: n ! Number of lines in Laskar files
    63 integer :: ierr, i
    64 
    65 ! CODE
    66 ! ----
    67 n = 0
    68 open(1,file = 'obl_ecc_lsp.asc',status = 'old',action = 'read')
    69 do
    70     read(1,*,iostat = ierr)
    71     if (ierr /= 0) exit
    72     n = n + 1
    73 enddo
    74 close(1)
    75 allocate(year_lask(n),obl_lask(n),ecc_lask(n),lsp_lask(n))
    76 
    77 iyear_lask = 0
    78 do i = 1,size(year_lask,1)
    79     read(1,*) year_lask(i), obl_lask(i), ecc_lask(i), lsp_lask(i)
    80     year_lask(i) = year_lask(i)*1000./convert_years ! Conversion from Laskar's kilo Earth years to PEM Martian years
    81     if (year_lask(i) > Year) iyear_lask = i + 1
    82 enddo
    83 close(1)
    84 if (iyear_lask == 0 .or. iyear_lask == size(year_lask,1) + 1) then
    85     error stop 'The current year could not be found in the file "obl_ecc_lsp.asc".'
     354integer(di) :: ierr, funit, i, n
     355logical(k4) :: here
     356real(dp)    :: current_year
     357
     358! CODE
     359! ----
     360call print_msg('> Reading "'//orbitfile_name//'"')
     361! Compute the current year
     362current_year = pem_ini_date + n_yr_sim
     363call print_msg('Current year = '//real2str(current_year))
     364
     365inquire(file = orbitfile_name,exist = here)
     366if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//orbitfile_name//'"!',1)
     367
     368! Read the Laskars's orbital file
     369open(newunit = funit,file = orbitfile_name,status = 'old',action = 'read',iostat = ierr)
     370if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'error opening file "'//orbitfile_name//'"!',ierr)
     371n = size(year_lask,1) ! Number of lines in the file
     372iyear_lask = 0 ! Line number for closest lower year relative to current simulation year
     373do i = 1,n
     374    read(funit,*) year_lask(i), obl_lask(i), ecc_lask(i), lsp_lask(i)
     375    year_lask(i) = year_lask(i)*1000._dp/r_plnt2earth_yr ! Conversion from Laskar's kilo Earth years to planetary years
     376    if (year_lask(i) > current_year) iyear_lask = i + 1
     377end do
     378close(funit)
     379if (iyear_lask == 0 .or. iyear_lask == n + 1) then
     380    call stop_clean(__FILE__,__LINE__,'the current year could not be found in the file "'//orbitfile_name//'".',1)
    86381else
    87     write(*,*) 'The current year in the "obl_ecc_lsp.asc" file is at line:', iyear_lask
    88 endif
     382    call print_msg('The current year in "'//orbitfile_name//'" is at line '//int2str(iyear_lask)//'.')
     383end if
    89384
    90385END SUBROUTINE read_orbitpm
    … …  
    92387
    93388!=======================================================================
    94 SUBROUTINE end_orbit
    95 !-----------------------------------------------------------------------
    96 ! NAME
    97 !     end_orbit
    98 !
    99 ! DESCRIPTION
    100 !     Deallocate orbital data arrays.
     389SUBROUTINE compute_maxyr_orbit(n_yr_sim,nmax_yr_runorb)
     390!-----------------------------------------------------------------------
     391! NAME
     392!     compute_maxyr_orbit
     393!
     394! DESCRIPTION
     395!     Determine maximum number of years before orbital params change
     396!     beyond admissible thresholds.
    101397!
    102398! AUTHORS & DATE
    … …  
    105401!
    106402! NOTES
    107 !
    108 !-----------------------------------------------------------------------
    109 
    110 ! DECLARATION
    111 ! -----------
    112 implicit none
    113 
    114 ! CODE
    115 ! ----
    116 if (allocated(year_lask)) deallocate(year_lask)
    117 if (allocated(obl_lask))  deallocate(obl_lask)
    118 if (allocated(ecc_lask))  deallocate(ecc_lask)
    119 if (allocated(lsp_lask))  deallocate(lsp_lask)
    120 
    121 END SUBROUTINE end_orbit
    122 !=======================================================================
    123 
    124 !=======================================================================
    125 SUBROUTINE orbit_param_criterion(i_myear,nyears_maxorb)
    126 !-----------------------------------------------------------------------
    127 ! NAME
    128 !     orbit_param_criterion
    129 !
    130 ! DESCRIPTION
    131 !     Determine maximum PEM iterations before orbital params change
    132 !     beyond admissible thresholds.
    133 !
    134 ! AUTHORS & DATE
    135 !     R. Vandemeulebrouck
    136 !     JB Clement, 2023-2025
    137 !
    138 ! NOTES
    139 !
    140 !-----------------------------------------------------------------------
    141 #ifdef CPP_IOIPSL
    142     use IOIPSL,          only: getin
    143 #else
    144     ! if not using IOIPSL, we still need to use (a local version of) getin
    145     use ioipsl_getincom, only: getin
    146 #endif
    147 use planete_h,      only: e_elips, obliquit, lsperi
    148 use phys_constants, only: pi
    149 use evolution,      only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
     403!     Handles Lsp modulo crossings.
     404!-----------------------------------------------------------------------
     405
     406! DEPENDENCIES
     407! ------------
     408use maths,     only: pi
     409use evolution, only: pem_ini_date
     410use display,   only: print_msg
     411use utility,   only: real2str
    150412
    151413! DECLARATION
    … …  
    155417! ARGUMENTS
    156418! ---------
    157 real, intent(in)  :: i_myear       ! Martian year of the simulation
    158 real, intent(out) :: nyears_maxorb ! Maximum number of iteration of the PEM before orb changes too much
     419real(dp), intent(in)  :: n_yr_sim       ! Current year of the simulation
     420real(dp), intent(out) :: nmax_yr_runorb ! Maximum number of years of the PEM before orbital parameters change too much
    159421
    160422! LOCAL VARIABLES
    161423! ---------------
    162 real                            :: Year                                           ! Year of the simulation
    163 real                            :: Lsp                                            ! Ls perihelion
    164 real                            :: max_change_obl, max_change_ecc, max_change_lsp ! Maximum admissible change
    165 real                            :: max_obl, max_ecc, max_lsp                      ! Maximum value of orbit param given the admissible change
    166 real                            :: min_obl, min_ecc, min_lsp                      ! Minimum value of orbit param given the admissible change
    167 real                            :: nyears_maxobl, nyears_maxecc, nyears_maxlsp    ! Maximum year iteration before reaching an inadmissible value of orbit param
    168 real                            :: xa, xb, ya, yb                  ! Variables for interpolation
    169 logical                         :: found_obl, found_ecc, found_lsp ! Flag variables for orbital parameters
    170 real                            :: lsp_corr                        ! Lsp corrected if the "modulo is crossed"
    171 real                            :: delta                           ! Intermediate variable
    172 real, dimension(:), allocatable :: lsplask_unwrap                  ! Unwrapped sequence of Lsp from Laskar's file
    173 integer                         :: i
    174 
    175 ! CODE
    176 ! ----
    177 ! **********************************************************************
    178 ! 0. Initializations
    179 ! **********************************************************************
    180 Year = year_bp_ini + i_myear ! We compute the current year
    181 Lsp = lsperi*180./pi         ! We convert in degree
    182 
    183 call read_orbitpm(Year) ! Allocation of variables
    184 
    185 write(*,*) 'orbit_param_criterion, Current year =', Year
    186 write(*,*) 'orbit_param_criterion, Current obl. =', obliquit
    187 write(*,*) 'orbit_param_criterion, Current ecc. =', e_elips
    188 write(*,*) 'orbit_param_criterion, Current Lsp  =', Lsp
     424real(dp)                            :: current_year
     425real(dp)                            :: max_obl, max_ecc, max_lsp
     426real(dp)                            :: min_obl, min_ecc, min_lsp
     427real(dp)                            :: nmax_yr_runobl, nmax_yr_runecc, nmax_yr_runlsp ! Maximum year iteration before reaching an inadmissible value of orbit param
     428real(dp)                            :: xa, xb, ya, yb
     429logical(k4)                         :: found_obl, found_ecc, found_lsp
     430real(dp)                            :: lsp_corr ! Lsp corrected if the "modulo is crossed"
     431real(dp)                            :: delta
     432real(dp), dimension(:), allocatable :: lsplask_unwrap
     433integer(di)                         :: i
     434
     435! CODE
     436! ----
     437if (.not. evol_orbit) then
     438    nmax_yr_runorb = largest_nb ! Default huge value
     439    return
     440end if
     441call print_msg('> Computing the accepted maximum number of years due to orbital parameters')
     442current_year = pem_ini_date + n_yr_sim ! We compute the current year
    189443
    190444! Unwrap the Lsp changes in Laskar's file
    … …  
    193447do i = iyear_lask,2,-1
    194448    delta = lsp_lask(i - 1) - lsp_lask(i)
    195     if (delta > 180.) then
    196         lsp_corr = lsp_lask(i - 1) - 360.
    197     else if (delta < -180.) then
    198         lsp_corr = lsp_lask(i - 1) + 360.
     449    if (delta > 180._dp) then
     450        lsp_corr = lsp_lask(i - 1) - 360._dp
     451    else if (delta < -180._dp) then
     452        lsp_corr = lsp_lask(i - 1) + 360._dp
    199453    else
    200454        lsp_corr = lsp_lask(i - 1)
    201     endif
     455    end if
    202456    lsplask_unwrap(i - 1) = lsplask_unwrap(i) + lsp_corr - lsp_lask(i)
    203 enddo
     457end do
    204458
    205459! Correct the current Lsp according to the unwrapping
    206 delta = Lsp - lsplask_unwrap(iyear_lask)
    207 if (delta > 180.) then
    208     Lsp = Lsp - 360.
    209 else if (delta < -180.) then
    210     Lsp = Lsp + 360.
    211 endif
    212 
    213 ! Constant max change case
    214 max_change_obl = 1.
    215 max_change_ecc = 5.e-3
    216 max_change_lsp = 20.
    217 
    218 call getin('max_change_obl',max_change_obl)
    219 max_obl = obliquit + max_change_obl
    220 min_obl = obliquit - max_change_obl
    221 write(*,*) 'Maximum obliquity accepted =', max_obl
    222 write(*,*) 'Minimum obliquity accepted =', min_obl
    223 
    224 call getin('max_change_ecc',max_change_ecc)
    225 max_ecc = min(e_elips + max_change_ecc,1. - 1.e-10) ! ecc < 1.
    226 min_ecc = max(e_elips - max_change_ecc,0.) ! ecc >= 0.
    227 write(*,*) 'Maximum eccentricity accepted =', max_ecc
    228 write(*,*) 'Minimum eccentricity accepted =', min_ecc
    229 
    230 call getin('max_change_lsp',max_change_lsp)
    231 max_lsp = Lsp + max_change_lsp
    232 min_lsp = Lsp - max_change_lsp
    233 write(*,*) 'Maximum Lsp accepted =', max_lsp
    234 write(*,*) 'Minimum Lsp accepted =', min_lsp
    235 ! End Constant max change case
    236 
    237 ! If we do not want some orb parameter to change, they should not be a stopping criterion,
    238 ! So the number of iteration corresponding is set to maximum
    239 nyears_maxobl = 1000000.
    240 nyears_maxecc = 1000000.
    241 nyears_maxlsp = 1000000.
    242 
    243 ! The maximum reachable year is found by interpolation according to the maximum admissible change of orbital parameters
     460delta = ls_peri - lsplask_unwrap(iyear_lask)
     461if (delta > 180._dp) then
     462    ls_peri = ls_peri - 360._dp
     463else if (delta < -180._dp) then
     464    ls_peri = ls_peri + 360._dp
     465end if
     466
     467! Maximum change accepted for orbital parameters
     468max_obl = obliquity + max_change_obl
     469min_obl = obliquity - max_change_obl
     470call print_msg('Obl. (current|accepted min|accepted max): '//real2str(obliquity)//' | '//real2str(min_obl)//real2str(max_obl))
     471
     472max_ecc = min(eccentricity + max_change_ecc,1. - minieps) ! ecc < 1.
     473min_ecc = max(eccentricity - max_change_ecc,0._dp) ! ecc >= 0.
     474call print_msg('Ecc. (current|accepted min|accepted max): '//real2str(eccentricity)//' | '//real2str(min_ecc)//' | '//real2str(max_ecc))
     475
     476max_lsp = ls_peri + max_change_lsp
     477min_lsp = ls_peri - max_change_lsp
     478call print_msg('Lsp  (current|accepted min|accepted max): '//real2str(ls_peri)//' | '//real2str(min_lsp)//' | '//real2str(max_lsp))
     479
     480! Initialization
     481nmax_yr_runobl = largest_nb ! Default huge value
     482nmax_yr_runecc = largest_nb ! Default huge value
     483nmax_yr_runlsp =largest_nb ! Default huge value
    244484found_obl = .false.
    245485found_ecc = .false.
    246486found_lsp = .false.
    247 if (.not. var_obl) found_obl = .true.
    248 if (.not. var_ecc) found_ecc = .true.
    249 if (.not. var_lsp) found_lsp = .true.
    250 
     487if (.not. evol_obl) found_obl = .true.
     488if (.not. evol_ecc) found_ecc = .true.
     489if (.not. evol_lsp) found_lsp = .true.
     490
     491! The maximum reachable year is found by interpolation according to the maximum admissible change of orbital parameters
    251492do i = iyear_lask,2,-1
    252493    xa = year_lask(i)
    … …  
    254495
    255496    ! Obliquity
    256     if (var_obl .and. .not. found_obl) then
     497    if (evol_obl .and. .not. found_obl) then
    257498        ya = obl_lask(i)
    258499        yb = obl_lask(i - 1)
    259500        if (yb < min_obl) then ! The minimum accepted is overtaken
    260             nyears_maxobl = (min_obl - ya)*(xb - xa)/(yb - ya) + xa - Year
     501            nmax_yr_runobl = (min_obl - ya)*(xb - xa)/(yb - ya) + xa - current_year
    261502            found_obl = .true.
    262503        else if (max_obl < yb) then ! The maximum accepted is overtaken
    263             nyears_maxobl = (max_obl - ya)*(xb - xa)/(yb - ya) + xa - Year
     504            nmax_yr_runobl = (max_obl - ya)*(xb - xa)/(yb - ya) + xa - current_year
    264505            found_obl = .true.
    265         endif
    266     endif
     506        end if
     507    end if
    267508
    268509    ! Eccentricity
    269     if (var_ecc .and. .not. found_ecc) then
     510    if (evol_ecc .and. .not. found_ecc) then
    270511        ya = ecc_lask(i)
    271512        yb = ecc_lask(i - 1)
    272513        if (yb < min_ecc) then ! The minimum accepted is overtaken
    273             nyears_maxecc = (min_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year
     514            nmax_yr_runecc = (min_ecc - ya)*(xb - xa)/(yb - ya) + xa - current_year
    274515            found_ecc = .true.
    275516        else if (max_ecc < yb) then ! The maximum accepted is overtaken
    276             nyears_maxecc = (max_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year
     517            nmax_yr_runecc = (max_ecc - ya)*(xb - xa)/(yb - ya) + xa - current_year
    277518            found_ecc = .true.
    278         endif
    279     endif
     519        end if
     520    end if
    280521
    281522    ! Lsp
    282     if (var_lsp .and. .not. found_lsp) then
     523    if (evol_lsp .and. .not. found_lsp) then
    283524        ya = lsplask_unwrap(i)
    284525        yb = lsplask_unwrap(i - 1)
    285526        if (yb < min_lsp) then ! The minimum accepted is overtaken
    286             nyears_maxlsp = (min_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year
     527            nmax_yr_runlsp = (min_lsp - ya)*(xb - xa)/(yb - ya) + xa - current_year
    287528            found_lsp = .true.
    288529        else if (max_lsp < yb) then ! The maximum accepted is overtaken
    289             nyears_maxlsp = (max_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year
     530            nmax_yr_runlsp = (max_lsp - ya)*(xb - xa)/(yb - ya) + xa - current_year
    290531            found_lsp = .true.
    291         endif
    292     endif
     532        end if
     533    end if
    293534
    294535    if (found_obl .and. found_ecc .and. found_lsp) exit ! The loop is left as soon as everything is found
    295 enddo
     536end do
    296537
    297538deallocate(lsplask_unwrap)
    298539
    299 if (.not. found_obl) write(*,*) 'The maximum reachable year for obl. could not be found in the file "obl_ecc_lsp.asc".'
    300 if (.not. found_ecc) write(*,*) 'The maximum reachable year for ecc. could not be found in the file "obl_ecc_lsp.asc".'
    301 if (.not. found_lsp) write(*,*) 'The maximum reachable year for Lsp could not be found in the file "obl_ecc_lsp.asc".'
    302 
    303 write(*,*) 'Max. number of years for the obl. parameter =', nyears_maxobl
    304 write(*,*) 'Max. number of years for the ecc. parameter =', nyears_maxecc
    305 write(*,*) 'Max. number of years for the Lsp parameter  =', nyears_maxlsp
    306 
    307 nyears_maxorb = min(nyears_maxobl,nyears_maxecc,nyears_maxlsp)
    308 if (nyears_maxorb < 1.) nyears_maxorb = 1. ! nyears_maxorb should be at least equal to 1
    309 write(*,*) 'So the max. number of years for the orbital criteria =', nyears_maxorb
    310 
    311 END SUBROUTINE orbit_param_criterion
    312 !=======================================================================
    313 
    314 !=======================================================================
    315 SUBROUTINE set_new_orbitpm(i_myear,i_myear_leg)
    316 !-----------------------------------------------------------------------
    317 ! NAME
    318 !     set_new_orbitpm
    319 !
    320 ! DESCRIPTION
    321 !     Recompute orbit parameters from Laskar values at a new year and
    322 !     update planetary constants; handles Lsp modulo crossings.
     540if (.not. found_obl) call print_msg('The maximum reachable year for obl. could not be found in "'//orbitfile_name//'".')
     541if (.not. found_ecc) call print_msg('The maximum reachable year for ecc. could not be found in "'//orbitfile_name//'".')
     542if (.not. found_lsp) call print_msg('The maximum reachable year for Lsp could not be found in "'//orbitfile_name//'".')
     543
     544call print_msg('Number of years accepted for obl. = '//real2str(nmax_yr_runobl))
     545call print_msg('Number of years accepted for ecc. = '//real2str(nmax_yr_runecc))
     546call print_msg('Number of years accepted for Lsp  = '//real2str(nmax_yr_runlsp))
     547
     548nmax_yr_runorb = min(nmax_yr_runobl,nmax_yr_runecc,nmax_yr_runlsp)
     549if (nmax_yr_runorb < 1._dp) nmax_yr_runorb = 1._dp ! nmax_yr_runorb should be at least equal to 1
     550
     551call print_msg('So the max. number of years for the orbital criterion = '//real2str(nmax_yr_runorb))
     552
     553END SUBROUTINE compute_maxyr_orbit
     554!=======================================================================
     555
     556!=======================================================================
     557SUBROUTINE update_orbit(n_yr_sim,n_yr_run)
     558!-----------------------------------------------------------------------
     559! NAME
     560!     update_orbit
     561!
     562! DESCRIPTION
     563!     Update orbital parameters from Laskar values at a new year.
    323564!
    324565! AUTHORS & DATE
    … …  
    330571!-----------------------------------------------------------------------
    331572
    332 use evolution,        only: year_bp_ini, var_obl, var_ecc, var_lsp
    333 use phys_constants,   only: pi
    334 use planete_h,        only: e_elips, obliquit, lsperi, periheli, aphelie, p_elips, peri_day, year_day
     573use evolution,        only: pem_ini_date
     574use maths,            only: pi
    335575use call_dayperi_mod, only: call_dayperi
     576use stoppage,         only: stop_clean
     577use display,          only: print_msg
     578use utility,          only: real2str
    336579
    337580! DECLARATION
    … …  
    341584! ARGUMENTS
    342585! ---------
    343 real, intent(in) :: i_myear     ! Number of simulated Martian years
    344 real, intent(in) :: i_myear_leg ! Number of iterations of the PEM
     586real(dp), intent(in) :: n_yr_sim ! Number of simulated years
     587real(dp), intent(in) :: n_yr_run ! Number of iterations of the PEM
    345588
    346589! LOCAL VARIABLES
    347590! ---------------
    348 real    :: Year           ! Year of the simulation
    349 real    :: Lsp            ! Ls perihelion
    350 real    :: halfaxe        ! Million of km
    351 integer :: i
    352 real    :: unitastr
    353 real    :: xa, xb, ya, yb ! Variables for interpolation
    354 logical :: found_year     ! Flag variable
    355 
    356 ! CODE
    357 ! ----
    358 ! **********************************************************************
    359 ! 0. Initializations
    360 ! **********************************************************************
    361 Year = year_bp_ini + i_myear ! We compute the current year
    362 Lsp = lsperi*180./pi         ! We convert in degree
    363 
    364 write(*,*) 'set_new_orbitpm, Old year =', Year - i_myear_leg
    365 write(*,*) 'set_new_orbitpm, Old obl. =', obliquit
    366 write(*,*) 'set_new_orbitpm, Old ecc. =', e_elips
    367 write(*,*) 'set_new_orbitpm, Old Lsp  =', Lsp
    368 write(*,*) 'set_new_orbitpm, New year =', Year
     591integer(di) :: i
     592real(dp)    :: current_year, old_year, obl_old, ecc_old, lsp_old
     593real(dp)    :: xa, xb, ya, yb ! Variables for interpolation
     594logical(k4) :: found_year
     595
     596! CODE
     597! ----
     598call print_msg('> Updating the orbital parameters for the new year')
     599current_year = pem_ini_date + n_yr_sim ! We compute the current year
     600old_year = current_year - n_yr_run
     601obl_old = obliquity
     602ecc_old = eccentricity
     603lsp_old = ls_peri
    369604
    370605found_year = .false.
    371 if (Year < 0.) then
     606if (current_year < 0._dp) then
    372607    do i = iyear_lask,2,-1
    373608        xa = year_lask(i)
    374609        xb = year_lask(i - 1)
    375         if (xa <= Year .and. Year < xb) then
     610        if (xa <= current_year .and. current_year < xb) then
    376611            ! Obliquity
    377             if (var_obl) then
     612            if (evol_obl) then
    378613                ya = obl_lask(i)
    379614                yb = obl_lask(i - 1)
    380                 obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya
    381             endif
     615                obliquity = (current_year - xa)*(yb - ya)/(xb - xa) + ya
     616            end if
    382617
    383618            ! Eccentricity
    384             if (var_ecc) then
     619            if (evol_ecc) then
    385620                ya = ecc_lask(i)
    386621                yb = ecc_lask(i - 1)
    387                 e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya
    388             endif
     622                eccentricity = (current_year - xa)*(yb - ya)/(xb - xa) + ya
     623            end if
    389624
    390625            ! Lsp
    391             if (var_lsp) then
     626            if (evol_lsp) then
    392627                ya = lsp_lask(i)
    393628                yb = lsp_lask(i - 1)
    394                 if (abs(yb - ya) > 300.) then ! If modulo is "crossed" through the interval
     629                if (abs(yb - ya) > 300._dp) then ! If modulo is "crossed" through the interval
    395630                    if (ya < yb) then ! Lsp should be decreasing
    396                         ya = ya + 360.
     631                        ya = ya + 360._dp
    397632                    else ! Lsp should be increasing
    398                         yb = yb + 360.
    399                     endif
    400                 endif
    401                 Lsp = modulo((Year - xa)*(yb - ya)/(xb - xa) + ya,360.)
    402             endif
     633                        yb = yb + 360._dp
     634                    end if
     635                end if
     636                ls_peri = modulo((current_year - xa)*(yb - ya)/(xb - xa) + ya,360._dp)
     637            end if
    403638            found_year = .true.
    404639            exit ! The loop is left as soon as the right interval is found
    405         endif
    406     enddo
    407 else if (Year >= 0 .and. Year < 100.) then ! For present orbital characteristics
    408     if (var_obl) obliquit = 25.19
    409     if (var_ecc) e_elips = 0.0934
    410     if (var_lsp) Lsp = 251.
     640        end if
     641    end do
     642else if (current_year >= 0._dp .and. current_year < 100._dp) then ! For present orbital characteristics
     643    if (evol_obl) obliquity = obl_lask(1)
     644    if (evol_ecc) eccentricity = ecc_lask(1)
     645    if (evol_lsp) ls_peri = lsp_lask(1)
    411646    found_year = .true.
    412 endif
    413 
    414 if (.not. found_year) error stop 'The new year could not be found in the file "obl_ecc_lsp.asc".'
    415 
    416 write(*,*) 'set_new_orbitpm, New obl. =', obliquit
    417 write(*,*) 'set_new_orbitpm, New ecc. =', e_elips
    418 write(*,*) 'set_new_orbitpm, New Lsp  =', Lsp
    419 
    420 halfaxe = 227.94
    421 Lsp = Lsp*pi/180.
    422 periheli = halfaxe*(1 - e_elips)
    423 aphelie =  halfaxe*(1 + e_elips)
    424 unitastr = 149.597927
    425 p_elips = 0.5*(periheli + aphelie)*(1. - e_elips*e_elips)/unitastr
    426 
    427 call call_dayperi(Lsp,e_elips,peri_day,year_day)
    428 
    429 call end_orbit()
    430 
    431 END SUBROUTINE set_new_orbitpm
     647end if
     648
     649if (.not. found_year) call stop_clean(__FILE__,__LINE__,'the new year could not be found in the file "'//orbitfile_name//'".',1)
     650
     651call lsp2datep(ls_peri,date_peri)
     652perihelion = semimaj_axis*(1._dp - eccentricity)
     653aphelion = semimaj_axis*(1._dp + eccentricity)
     654
     655call print_msg('Year (ini|now) = '//real2str(old_year)//' | '//real2str(current_year))
     656call print_msg('Obl. (ini|now) = '//real2str(obl_old)//' | '//real2str(obliquity))
     657call print_msg('Ecc. (ini|now) = '//real2str(ecc_old)//' | '//real2str(eccentricity))
     658call print_msg('Lsp  (ini|now) = '//real2str(lsp_old)//' | '//real2str(ls_peri))
     659
     660END SUBROUTINE update_orbit
    432661!=======================================================================
    433662
    434663END MODULE orbit
    435 
Note: See TracChangeset for help on using the changeset viewer.