Changeset 3989 for trunk/LMDZ.COMMON/libf/evolution/orbit.F90

Timestamp:

Dec 11, 2025, 12:56:05 PM (5 weeks ago)

Author:

jbclement

Message:

PEM:
Massive structural refactor of the PEM codebase for improved readability, consistency and maintainability. The goal is to modernize, standardize and consolidate the code while removing legacy complexity. In detail, this change:

• Performs large-scale cleanup removing unused code, obsolete routines, duplicated functionality and deprecated initialization logic;
• Removes "*_mod" wrappers;
• Replaces mixed naming conventions with clearer variable names, domain-based file/module names and purpose-based routine names;
• Adds native reading/writing capabilities to the PEM removing the cumbersome dependency on Mars PCM subroutines;
• Centralizes soil/ice/orbital/PEM configuration logic into dedicated modules;
• Simplifies routines and updates calls/interfaces to match the new structure.

This refactor significantly clarifies the codebase and provides a cleaner foundation for forthcoming developments.
JBC

File:

• trunk/LMDZ.COMMON/libf/evolution/orbit.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90) (12 diffs)

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

@@ -1 @@
-MODULE orbit_param_criterion_mod
+MODULE orbit
 
 !=======================================================================
@@ -11 +11 @@
 implicit none
 
+real, dimension(:), allocatable :: year_lask  ! year before present from Laskar et al.
+real, dimension(:), allocatable :: obl_lask   ! obliquity    [deg]
+real, dimension(:), allocatable :: ecc_lask   ! eccentricity
+real, dimension(:), allocatable :: lsp_lask   ! ls perihelie [deg]
+integer                         :: iyear_lask ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
+
 !=======================================================================
 contains
 !=======================================================================
 
-SUBROUTINE orbit_param_criterion(i_myear,n_myear_leg)
+SUBROUTINE read_orbitpm(Year)
+
+use evolution, only: convert_years
+
+implicit none
+
+real, intent(in) :: Year ! Martian year of the simulation
+
+integer :: n ! number of lines in Laskar files
+integer :: ierr, i
+
+n = 0
+open(1,file = 'obl_ecc_lsp.asc',status = 'old',action = 'read')
+do
+    read(1,*,iostat = ierr)
+    if (ierr /= 0) exit
+    n = n + 1
+enddo
+close(1)
+allocate(year_lask(n),obl_lask(n),ecc_lask(n),lsp_lask(n))
+
+iyear_lask = 0
+do i = 1,size(year_lask,1)
+    read(1,*) year_lask(i), obl_lask(i), ecc_lask(i), lsp_lask(i)
+    year_lask(i) = year_lask(i)*1000./convert_years ! Conversion from Laskar's kilo Earth years to PEM Martian years
+    if (year_lask(i) > Year) iyear_lask = i + 1
+enddo
+close(1)
+if (iyear_lask == 0 .or. iyear_lask == size(year_lask,1) + 1) then
+    error stop 'The current year could not be found in the file "obl_ecc_lsp.asc".'
+else
+    write(*,*) 'The current year in the "obl_ecc_lsp.asc" file is at line:', iyear_lask
+endif
+
+END SUBROUTINE read_orbitpm
+
+!=======================================================================
+
+SUBROUTINE end_orbit
+
+implicit none
+
+if (allocated(year_lask)) deallocate(year_lask)
+if (allocated(obl_lask)) deallocate(obl_lask)
+if (allocated(ecc_lask)) deallocate(ecc_lask)
+if (allocated(lsp_lask)) deallocate(lsp_lask)
+
+END SUBROUTINE end_orbit
+!=======================================================================
+
+SUBROUTINE orbit_param_criterion(i_myear,nyears_maxorb)
 
 #ifdef CPP_IOIPSL
@@ -25 +81 @@
 use planete_h,      only: e_elips, obliquit, lsperi
 use phys_constants, only: pi
-use time_evol_mod,  only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
-use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, ini_lask_param_mod, last_ilask
+use evolution,      only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
 
 implicit none
@@ -38 +93 @@
 ! Output Variables
 !--------------------------------------------------------
-real, intent(out) :: n_myear_leg ! Maximum number of iteration of the PEM before orb changes too much
+real, intent(out) :: nyears_maxorb ! Maximum number of iteration of the PEM before orb changes too much
 
 !--------------------------------------------------------
@@ -48 +103 @@
 real                            :: max_obl, max_ecc, max_lsp                      ! Maximum value of orbit param given the admissible change
 real                            :: min_obl, min_ecc, min_lsp                      ! Minimum value of orbit param given the admissible change
-real                            :: max_obl_iter, max_ecc_iter, max_lsp_iter       ! Maximum year iteration before reaching an inadmissible value of orbit param
-integer                         :: ilask                                          ! Loop variable
+real                            :: nyears_maxobl, nyears_maxecc, nyears_maxlsp       ! Maximum year iteration before reaching an inadmissible value of orbit param
+integer                         :: i                                          ! Loop variable
 real                            :: xa, xb, ya, yb                                 ! Variables for interpolation
 logical                         :: found_obl, found_ecc, found_lsp                ! Flag variables for orbital parameters
@@ -62 +117 @@
 Lsp = lsperi*180./pi         ! We convert in degree
 
-call ini_lask_param_mod() ! Allocation of variables
+call read_orbitpm(Year) ! Allocation of variables
 
 write(*,*) 'orbit_param_criterion, Current year =', Year
@@ -69 +124 @@
 write(*,*) 'orbit_param_criterion, Current Lsp  =', Lsp
 
-! We read the file
-open(73,file = 'obl_ecc_lsp.asc',status = 'old',action = 'read')
-last_ilask = 0
-do ilask = 1,size(yearlask,1)
-    read(73,*) yearlask(ilask), obllask(ilask), ecclask(ilask), lsplask(ilask)
-    yearlask(ilask) = yearlask(ilask)*1000./convert_years ! Conversion from Laskar's kilo Earth years to PEM Martian years
-    if (yearlask(ilask) > Year) last_ilask = ilask + 1
+! Unwrap the Lsp changes in Laskar's file
+allocate(lsplask_unwrap(iyear_lask))
+lsplask_unwrap(iyear_lask) = lsp_lask(iyear_lask)
+do i = iyear_lask,2,-1
+    delta = lsp_lask(i - 1) - lsp_lask(i)
+    if (delta > 180.) then
+        lsp_corr = lsp_lask(i - 1) - 360.
+    else if (delta < -180.) then
+        lsp_corr = lsp_lask(i - 1) + 360.
+    else
+        lsp_corr = lsp_lask(i - 1)
+    endif
+    lsplask_unwrap(i - 1) = lsplask_unwrap(i) + lsp_corr - lsp_lask(i)
 enddo
-close(73)
-if (last_ilask == 0 .or. last_ilask == size(yearlask,1) + 1) then
-    error stop 'The current year could not be found in the file "obl_ecc_lsp.asc".'
-else
-    write(*,*) 'The current year in the "obl_ecc_lsp.asc" file is at line:', last_ilask
-endif
-
-! Unwrap the Lsp changes in Laskar's file
-allocate(lsplask_unwrap(last_ilask))
-lsplask_unwrap(last_ilask) = lsplask(last_ilask)
-do ilask = last_ilask,2,-1
-    delta = lsplask(ilask - 1) - lsplask(ilask)
-    if (delta > 180.) then
-        lsp_corr = lsplask(ilask - 1) - 360.
-    else if (delta < -180.) then
-        lsp_corr = lsplask(ilask - 1) + 360.
-    else
-        lsp_corr = lsplask(ilask - 1)
-    endif
-    lsplask_unwrap(ilask - 1) = lsplask_unwrap(ilask) + lsp_corr - lsplask(ilask)
-enddo
 
 ! Correct the current Lsp according to the unwrapping
-delta = Lsp - lsplask_unwrap(last_ilask)
+delta = Lsp - lsplask_unwrap(iyear_lask)
 if (delta > 180.) then
     Lsp = Lsp - 360.
@@ -133 +173 @@
 ! If we do not want some orb parameter to change, they should not be a stopping criterion,
 ! So the number of iteration corresponding is set to maximum
-max_obl_iter = 1000000.
-max_ecc_iter = 1000000.
-max_lsp_iter = 1000000.
+nyears_maxobl = 1000000.
+nyears_maxecc = 1000000.
+nyears_maxlsp = 1000000.
 
 ! The maximum reachable year is found by interpolation according to the maximum admissible change of orbital parameters
@@ -145 +185 @@
 if (.not. var_lsp) found_lsp = .true.
 
-do ilask = last_ilask,2,-1
-    xa = yearlask(ilask)
-    xb = yearlask(ilask - 1)
+do i = iyear_lask,2,-1
+    xa = year_lask(i)
+    xb = year_lask(i - 1)
 
     ! Obliquity
     if (var_obl .and. .not. found_obl) then
-        ya = obllask(ilask)
-        yb = obllask(ilask - 1)
+        ya = obl_lask(i)
+        yb = obl_lask(i - 1)
         if (yb < min_obl) then ! The minimum accepted is overtaken
-            max_obl_iter = (min_obl - ya)*(xb - xa)/(yb - ya) + xa - Year
+            nyears_maxobl = (min_obl - ya)*(xb - xa)/(yb - ya) + xa - Year
             found_obl = .true.
         else if (max_obl < yb) then ! The maximum accepted is overtaken
-            max_obl_iter = (max_obl - ya)*(xb - xa)/(yb - ya) + xa - Year
+            nyears_maxobl = (max_obl - ya)*(xb - xa)/(yb - ya) + xa - Year
             found_obl = .true.
         endif
@@ -164 +204 @@
     ! Eccentricity
     if (var_ecc .and. .not. found_ecc) then
-        ya = ecclask(ilask)
-        yb = ecclask(ilask - 1)
+        ya = ecc_lask(i)
+        yb = ecc_lask(i - 1)
         if (yb < min_ecc) then ! The minimum accepted is overtaken
-            max_ecc_iter = (min_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year
+            nyears_maxecc = (min_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year
             found_ecc = .true.
         else if (max_ecc < yb) then ! The maximum accepted is overtaken
-            max_ecc_iter = (max_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year
+            nyears_maxecc = (max_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year
             found_ecc = .true.
         endif
@@ -177 +217 @@
     ! Lsp
     if (var_lsp .and. .not. found_lsp) then
-        ya = lsplask_unwrap(ilask)
-        yb = lsplask_unwrap(ilask - 1)
+        ya = lsplask_unwrap(i)
+        yb = lsplask_unwrap(i - 1)
         if (yb < min_lsp) then ! The minimum accepted is overtaken
-            max_lsp_iter = (min_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year
+            nyears_maxlsp = (min_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year
             found_lsp = .true.
         else if (max_lsp < yb) then ! The maximum accepted is overtaken
-            max_lsp_iter = (max_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year
+            nyears_maxlsp = (max_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year
             found_lsp = .true.
         endif
@@ -197 +237 @@
 if (.not. found_lsp) write(*,*) 'The maximum reachable year for Lsp could not be found in the file "obl_ecc_lsp.asc".'
 
-write(*,*) 'Max. number of iterations for the obl. parameter =', max_obl_iter
-write(*,*) 'Max. number of iterations for the ecc. parameter =', max_ecc_iter
-write(*,*) 'Max. number of iterations for the Lsp parameter  =', max_lsp_iter
-
-n_myear_leg = min(max_obl_iter,max_ecc_iter,max_lsp_iter)
-if (n_myear_leg < 1.) n_myear_leg = 1. ! n_myear_leg should be at least equal to 1
-write(*,*) 'So the max. number of iterations for the orbital criteria =', n_myear_leg
+write(*,*) 'Max. number of years for the obl. parameter =', nyears_maxobl
+write(*,*) 'Max. number of years for the ecc. parameter =', nyears_maxecc
+write(*,*) 'Max. number of years for the Lsp parameter  =', nyears_maxlsp
+
+nyears_maxorb = min(nyears_maxobl,nyears_maxecc,nyears_maxlsp)
+if (nyears_maxorb < 1.) nyears_maxorb = 1. ! nyears_maxorb should be at least equal to 1
+write(*,*) 'So the max. number of years for the orbital criteria =', nyears_maxorb
 
 END SUBROUTINE orbit_param_criterion
 
 !***********************************************************************
-
-END MODULE orbit_param_criterion_mod
-
+! Purpose: Recompute orbit parameters based on values by Laskar et al.
+SUBROUTINE set_new_orbitpm(i_myear,i_myear_leg)
+
+use evolution,        only: year_bp_ini, var_obl, var_ecc, var_lsp
+use phys_constants,   only: pi
+use planete_h,        only: e_elips, obliquit, lsperi, periheli, aphelie, p_elips, peri_day, year_day
+use call_dayperi_mod, only: call_dayperi
+
+implicit none
+
+!--------------------------------------------------------
+! Input Variables
+!--------------------------------------------------------
+real, intent(in) :: i_myear     ! Number of simulated Martian years
+real, intent(in) :: i_myear_leg ! Number of iterations of the PEM
+
+!--------------------------------------------------------
+! Output Variables
+!--------------------------------------------------------
+
+!--------------------------------------------------------
+! Local variables
+!--------------------------------------------------------
+real    :: Year           ! Year of the simulation
+real    :: Lsp            ! Ls perihelion
+integer :: i          ! Loop variable
+real    :: halfaxe        ! Million of km
+real    :: unitastr
+real    :: xa, xb, ya, yb ! Variables for interpolation
+logical :: found_year     ! Flag variable
+
+! **********************************************************************
+! 0. Initializations
+! **********************************************************************
+Year = year_bp_ini + i_myear ! We compute the current year
+Lsp = lsperi*180./pi         ! We convert in degree
+
+write(*,*) 'set_new_orbitpm, Old year =', Year - i_myear_leg
+write(*,*) 'set_new_orbitpm, Old obl. =', obliquit
+write(*,*) 'set_new_orbitpm, Old ecc. =', e_elips
+write(*,*) 'set_new_orbitpm, Old Lsp  =', Lsp
+write(*,*) 'set_new_orbitpm, New year =', Year
+
+found_year = .false.
+if (Year < 0.) then
+    do i = iyear_lask,2,-1
+        xa = year_lask(i)
+        xb = year_lask(i - 1)
+        if (xa <= Year .and. Year < xb) then
+            ! Obliquity
+            if (var_obl) then
+                ya = obl_lask(i)
+                yb = obl_lask(i - 1)
+                obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya
+            endif
+
+            ! Eccentricity
+            if (var_ecc) then
+                ya = ecc_lask(i)
+                yb = ecc_lask(i - 1)
+                e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya
+            endif
+
+            ! Lsp
+            if (var_lsp) then
+                ya = lsp_lask(i)
+                yb = lsp_lask(i - 1)
+                if (abs(yb - ya) > 300.) then ! If modulo is "crossed" through the interval
+                    if (ya < yb) then ! Lsp should be decreasing
+                        ya = ya + 360.
+                    else ! Lsp should be increasing
+                        yb = yb + 360.
+                    endif
+                endif
+                Lsp = modulo((Year - xa)*(yb - ya)/(xb - xa) + ya,360.)
+            endif
+            found_year = .true.
+            exit ! The loop is left as soon as the right interval is found
+        endif
+    enddo
+else if (Year >= 0 .and. Year < 100.) then ! For present orbital characteristics
+    if (var_obl) obliquit = 25.19
+    if (var_ecc) e_elips = 0.0934
+    if (var_lsp) Lsp = 251.
+    found_year = .true.
+endif
+
+if (.not. found_year) error stop 'The new year could not be found in the file "obl_ecc_lsp.asc".'
+
+write(*,*) 'set_new_orbitpm, New obl. =', obliquit
+write(*,*) 'set_new_orbitpm, New ecc. =', e_elips
+write(*,*) 'set_new_orbitpm, New Lsp  =', Lsp
+
+halfaxe = 227.94
+Lsp = Lsp*pi/180.
+periheli = halfaxe*(1 - e_elips)
+aphelie =  halfaxe*(1 + e_elips)
+unitastr = 149.597927
+p_elips = 0.5*(periheli + aphelie)*(1. - e_elips*e_elips)/unitastr
+
+call call_dayperi(Lsp,e_elips,peri_day,year_day)
+
+call end_orbit()
+
+END SUBROUTINE set_new_orbitpm
+
+END MODULE orbit
+