Changeset 3791

Timestamp:

Jun 4, 2025, 6:27:23 PM (4 days ago)

Author:

jbclement

Message:

PEM:
Correction and simplification to compute the maximum number of years according to the changes in Lsp + update of "run_PEM.def".
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• deftank/run_PEM.def (modified) (2 diffs)
• orbit_param_criterion_mod.F90 (modified) (6 diffs)

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

@@ -697 +697 @@
 
 == 03/06/2025 == JBC
-Correction of a small bug for the launching script when restarting the simulation from an initial PCM run + renaming of variables in the visualization scripts for the layering structure
+Correction of a small bug for the launching script when restarting the simulation from an initial PCM run + renaming of variables in the visualization scripts for the layering structure.
+
+== 04/06/2025 == JBC
+Correction and simplification to compute the maximum number of years according to the changes in Lsp + update of "run_PEM.def".

trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def

@@ -26 +26 @@
 # var_lsp=.false.
 
+# Time step length of the PEM in Martian years? Default = 1.
+# dt=1
+
 #---------- Stopping criteria parameters ----------#
 # If evol_orbit_pem=.false., maximal number of iterations if no stopping criterion is reached? Default=100000000
@@ -39 +42 @@
 # ps_criterion = 0.15
 
-# Time step length of the PEM in Martian years? Default = 1.
-# dt=1
+# Acceptance of change for obliquity? Default = 1.
+# max_change_obl=1.
+
+# Acceptance of change for eccentricity? Default = 5.e-3
+# max_change_ecc=5.e-3
+
+# Acceptance of change for Lsp? Default = 20.
+# max_change_lsp=20.
 
 #---------- Subsurface parameters ----------#

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

@@ -48 +48 @@
 ! Local variables
 !--------------------------------------------------------
-real    :: Year                                           ! Year of the simulation
-real    :: Lsp                                            ! Ls perihelion
-real    :: max_change_obl, max_change_ecc, max_change_lsp ! Maximum admissible change
-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    :: xa, xb, ya, yb                                 ! Variables for interpolation
-logical :: found_obl, found_ecc, found_lsp                ! Flag variables for orbital parameters
-logical :: crossed_modulo_d, crossed_modulo_i             ! Flag variables for modulo of Lsp
+real                            :: Year                                           ! Year of the simulation
+real                            :: Lsp                                            ! Ls perihelion
+real                            :: max_change_obl, max_change_ecc, max_change_lsp ! Maximum admissible change
+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                            :: xa, xb, ya, yb                                 ! Variables for interpolation
+logical                         :: found_obl, found_ecc, found_lsp                ! Flag variables for orbital parameters
+real                            :: lsp_corr                                       ! Lsp corrected if the "modulo is crossed"
+real                            :: delta                                          ! Intermediate variable
+real, dimension(:), allocatable :: lsplask_unwrap                                 ! Unwrapped sequence of Lsp from Laskar's file
 
 ! **********************************************************************
@@ -65 +67 @@
 Lsp = lsperi*180./pi         ! We convert in degree
 
-call ini_lask_param_mod ! Allocation of variables
+call ini_lask_param_mod() ! Allocation of variables
 
 write(*,*) 'orbit_param_criterion, Current year =', Year
@@ -87 +89 @@
 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)
+if (delta > 180.) then
+    Lsp = Lsp - 360.
+else if (delta < -180.) then
+    Lsp = Lsp + 360.
+endif
+
 ! Constant max change case
-max_change_obl = 0.6
-max_change_ecc = 2.8e-3
-max_change_lsp = 18.
+max_change_obl = 1.
+max_change_ecc = 5.e-3
+max_change_lsp = 20.
 
 call getin('max_change_obl',max_change_obl)
@@ -124 +149 @@
 if (.not. var_ecc) found_ecc = .true.
 if (.not. var_lsp) found_lsp = .true.
-crossed_modulo_d = .false.
-crossed_modulo_i = .false.
 
 do ilask = last_ilask,2,-1
@@ -159 +182 @@
     ! Lsp
     if (var_lsp .and. .not. found_lsp) then
-        ya = lsplask(ilask)
-        yb = lsplask(ilask - 1)
-        if (abs(yb - ya) > 300.) then ! If modulo is "crossed" through the interval
-            if (ya < yb) then ! Lsp should be decreasing
-                if (.not. crossed_modulo_i) then
-                    yb = yb - 360.
-                    if (min_lsp < 0.) crossed_modulo_d = .true.
-                else
-                    crossed_modulo_i = .false.
-                endif
-            else ! Lsp should be increasing
-                if (.not. crossed_modulo_d) then
-                    yb = yb + 360.
-                    if (max_lsp > 360.) crossed_modulo_i = .true.
-                else
-                    crossed_modulo_d = .false.
-                endif
-            endif
-        else
-            if (crossed_modulo_d) then ! If decreasing Lsp "crossed" the modulo but min_lsp has not been met yet
-                ya = ya - 360.
-                yb = yb - 360.
-            else if (crossed_modulo_i) then ! If increasing Lsp "crossed" the modulo but max_lsp has not been met yet
-                ya = ya + 360.
-                yb = yb + 360.
-            endif
-        endif
+        ya = lsplask_unwrap(ilask)
+        yb = lsplask_unwrap(ilask - 1)
         if (yb < min_lsp) then ! The minimum accepted is overtaken
             max_lsp_iter = (min_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year
@@ -198 +196 @@
 enddo
 
+deallocate(lsplask_unwrap)
+
 if (.not. found_obl) write(*,*) 'The maximum reachable year for obl. could not be found in the file "obl_ecc_lsp.asc".'
 if (.not. found_ecc) write(*,*) 'The maximum reachable year for ecc. could not be found in the file "obl_ecc_lsp.asc".'