Changeset 3850 for trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90

Timestamp:

Jul 15, 2025, 2:43:41 PM (3 days ago)

Author:

jbclement

Message:

PEM:

• Handling the case where the new date is positive (present-day climate) in "recomp_orb_param_mod.F90", especially because of rounding errors.
• Correction in "visu_evol_layering.py" to get the proper extent of the stratication + replacing 'imshow' by 'pcolormesh' to better process the data.
• Addition of a graphic to plot the changes in obliquity directly onto the evolution of the stratification over time for better analysis. Obliquity curve is colored according to stratification gain/loss.
• Correction of numbering in "lib_launchPEM.sh" when requesting a relaunch.
• Update of "PEMrun.job" due to new parsing of arguments introduced in r3836/r3837.

JBC

File:

• trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90 (modified) (4 diffs)

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

@@ -31 +31 @@
 ! Input Variables
 !--------------------------------------------------------
-real, intent(in) :: i_myear   ! Number of simulated Martian years
+real, intent(in) :: i_myear     ! Number of simulated Martian years
 real, intent(in) :: i_myear_leg ! Number of iterations of the PEM
 
@@ -62 +62 @@
 Lsp = lsperi*180./pi         ! We convert in degree
 
+print*, year_bp_ini,i_myear,i_myear_leg
+
 write(*,*) 'recomp_orb_param, Old year =', Year - i_myear_leg
 write(*,*) 'recomp_orb_param, Old obl. =', obliquit
@@ -69 +71 @@
 
 found_year = .false.
-do ilask = last_ilask,2,-1
-    xa = yearlask(ilask)
-    xb = yearlask(ilask - 1)
-    if (xa <= Year .and. Year < xb) then
-        ! Obliquity
-        if (var_obl) then
-            ya = obllask(ilask)
-            yb = obllask(ilask - 1)
-            obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya
+if (Year < 0.) then
+    do ilask = last_ilask,2,-1
+        xa = yearlask(ilask)
+        xb = yearlask(ilask - 1)
+        if (xa <= Year .and. Year < xb) then
+            ! Obliquity
+            if (var_obl) then
+                ya = obllask(ilask)
+                yb = obllask(ilask - 1)
+                obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya
+            endif
+
+            ! Eccentricity
+            if (var_ecc) then
+                ya = ecclask(ilask)
+                yb = ecclask(ilask - 1)
+                e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya
+            endif
+
+            ! Lsp
+            if (var_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
+                        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
-
-        ! Eccentricity
-        if (var_ecc) then
-            ya = ecclask(ilask)
-            yb = ecclask(ilask - 1)
-            e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya
-        endif
-
-        ! Lsp
-        if (var_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
-                    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
+    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(*,*) 'recomp_orb_param, New obl. =', obliquit
+write(*,*) 'recomp_orb_param, New ecc. =', e_elips
+write(*,*) 'recomp_orb_param, New Lsp  =', Lsp
 
 halfaxe = 227.94
@@ -118 +131 @@
 #endif
 
-write(*,*) 'recomp_orb_param, New obl. =', obliquit
-write(*,*) 'recomp_orb_param, New ecc. =', e_elips
-write(*,*) 'recomp_orb_param, New Lsp  =', Lsp
-
 END SUBROUTINE recomp_orb_param