Changeset 3047 for trunk/LMDZ.COMMON/libf

Timestamp:

Sep 21, 2023, 11:56:12 AM (21 months ago)

Author:

jbclement

Message:

Mars PEM:
Correction of a case where maximum admissible change of orbital parameters could not be found, in particular for Lsp because of modulo + some improvements.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• orbit_param_criterion_mod.F90 (modified) (8 diffs)
• recomp_orb_param_mod.F90 (modified) (3 diffs)

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

@@ -52 +52 @@
   integer :: max_obl_iter, max_ecc_iter, max_lsp_iter       ! Maximum year iteration before reaching an inadmissible value of orbit param
   integer :: ilask, iylask                                  ! Loop variables
-  real    :: xa, xb, ya, yb, yc_obl, yc_ecc, yc_lsp         ! Variables for interpolation
+  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
 
 ! **********************************************************************
@@ -69 +71 @@
 ! We read the file
 open(73,file = 'obl_ecc_lsp.asc')
+last_ilask = 0
 do ilask = 1,size(yearlask,1)
     read(73,*) yearlask(ilask), obllask(ilask), ecclask(ilask), lsplask(ilask)
@@ -75 +78 @@
 enddo
 close(73)
-write(*,*) 'Corresponding line in the obl_ecc_lsp.asc file =', last_ilask
+if (last_ilask == 0 .or. last_ilask == size(yearlask,1) + 1) then
+    write(*,*) 'The current year could not be find in the file "obl_ecc_lsp.asc".'
+    stop
+else
+    write(*,*) 'The current year in the "obl_ecc_lsp.asc" file is at line:', last_ilask
+endif
 
 ! Constant max change case 
@@ -82 +90 @@
 max_change_lsp = 18.
 
-call getin('max_change_obl', max_change_obl)
+call getin('max_change_obl',max_change_obl)
 max_obl = obliquit + max_change_obl
 min_obl = obliquit - max_change_obl
@@ -88 +96 @@
 write(*,*) 'Minimum obliquity accepted =', min_obl
 
-call getin('max_change_ecc', max_change_ecc)
+call getin('max_change_ecc',max_change_ecc)
 max_ecc = min(e_elips + max_change_ecc,1. - 1.e-10) ! ex < 1.
 min_ecc = max(e_elips - max_change_ecc,0.) ! ex >= 0.
@@ -94 +102 @@
 write(*,*) 'Minimum eccentricity accepted =', min_ecc
 
-call getin('max_change_lsp', max_change_lsp)
+call getin('max_change_lsp',max_change_lsp)
 max_lsp = Lsp + max_change_lsp
 min_lsp = Lsp - max_change_lsp
@@ -107 +115 @@
 max_lsp_iter = 1000000
 
-        ! Tendency of the orbital parameter for the considered year gives the limitation between min and max
+! The maximum reachable year is found by interpolation according to the maximum admissible change of orbital parameters
+found_obl = .false.
+found_ecc = .false.
+found_lsp = .false.
+if (.not. var_obl) found_obl = .true.
+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
     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)
-            if (ya < yb) then ! Increasing -> max is the limitation
-                yc_obl = max_obl
-            else ! Decreasing -> min is the limitation
-                yc_obl = min_obl
-            endif
-        endif
-
-        ! Eccentricity
-        if (var_ecc) then
-            ya = ecclask(ilask)
-            yb = ecclask(ilask - 1)
-            if (ya < yb) then ! Increasing -> max is the limitation
-                yc_ecc = max_ecc
-            else ! Decreasing -> min is the limitation
-                yc_ecc = min_ecc
-            endif
-        endif
-
-        ! Lsp
-        if (var_lsp) then
-            ya = lsplask(ilask)
-            yb = lsplask(ilask - 1)
-             if (ya < yb) then ! Increasing -> max is the limitation
-                if (yb - ya > 300.) then ! If modulo is "crossed" then it is the other way around
-                    yc_lsp = min_lsp
+
+    ! Obliquity
+    if (var_obl .and. .not. found_obl) then
+        ya = obllask(ilask)
+        yb = obllask(ilask - 1)
+        if (yb < min_obl) then ! The minimum accepted is overtaken
+            max_obl_iter = floor((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 = floor((max_obl - ya)*(xb - xa)/(yb - ya) + xa) - Year
+            found_obl = .true.
+        endif
+    endif
+
+    ! Eccentricity
+    if (var_ecc .and. .not. found_ecc) then
+        ya = ecclask(ilask)
+        yb = ecclask(ilask - 1)
+        if (yb < min_ecc) then ! The minimum accepted is overtaken
+            max_ecc_iter = floor((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 = floor((max_ecc - ya)*(xb - xa)/(yb - ya) + xa) - Year
+            found_ecc = .true.
+        endif
+    endif
+
+    ! 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
-                    yc_lsp = max_lsp
+                    crossed_modulo_i = .false.
                 endif
-            else ! Decreasing -> min is the limitation
-                if (ya - yb > 300.) then ! If modulo is "crossed" then it is the other way around
-                    yc_lsp = max_lsp
+            else ! Lsp should be increasing
+                if (.not. crossed_modulo_d) then
+                    yb = yb + 360.
+                    if (max_lsp > 360.) crossed_modulo_i = .true.
                 else
-                    yc_lsp = min_lsp
+                    crossed_modulo_d = .false.
                 endif
             endif
-        endif
-        iylask = ilask
-        exit ! The loop is left as soon as the right interval is found
+        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
+        if (yb < min_lsp) then ! The minimum accepted is overtaken
+            max_lsp_iter = floor((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 = floor((max_lsp - ya)*(xb - xa)/(yb - ya) + xa) - Year
+            found_lsp = .true.
+        endif
     endif
+
+    if (found_obl .and. found_ecc .and. found_lsp) exit ! The loop is left as soon as everything is found
 enddo
 
-if (ilask == 1) then
-    write(*,*) 'The year does not match with Laskar data in the file obl_ecc_lsp.asc.'
-    stop
-endif
-
-! Linear interpolation gives the maximum reachable year according to the limitation
-! Obliquity
-if (var_obl) then
-    do ilask = iylask,2,-1
-        ya = obllask(ilask)
-        yb = obllask(ilask - 1)
-        if (min(ya,yb) <= yc_obl .and. yc_obl < max(ya,yb)) then
-            xa = yearlask(ilask)
-            xb = yearlask(ilask - 1)
-            max_obl_iter = floor((yc_obl - ya)*(xb - xa)/(yb - ya) + xa) - Year
-            exit ! The loop is left as soon as the right interval is found
-        endif
-    enddo
-endif
-! Eccentricity
-if (var_ecc) then
-    do ilask = iylask,2,-1
-        ya = ecclask(ilask)
-        yb = ecclask(ilask - 1)
-        if (min(ya,yb) <= yc_ecc .and. yc_ecc < max(ya,yb)) then
-            xa = yearlask(ilask)
-            xb = yearlask(ilask - 1)
-            max_ecc_iter = floor((yc_ecc - ya)*(xb - xa)/(yb - ya) + xa) - Year
-            exit ! The loop is left as soon as the right interval is found
-        endif
-    enddo
-endif
-! Lsp
-if (var_lsp) then
-    do ilask = iylask,2,-1
-        ya = lsplask(ilask)
-        yb = lsplask(ilask - 1)
-        if (min(ya,yb) <= yc_lsp .and. yc_lsp < max(ya,yb)) then
-            xa = yearlask(ilask)
-            xb = yearlask(ilask - 1)
-            max_lsp_iter = floor((yc_lsp - ya)*(xb - xa)/(yb - ya) + xa) - Year
-            exit ! The loop is left as soon as the right interval is found
-        endif
-    enddo
-endif
+if (.not. found_obl) write(*,*) 'The maximum reachable year for obl. could not be find in the file "obl_ecc_lsp.asc".'
+if (.not. found_ecc) write(*,*) 'The maximum reachable year for ecc. could not be find in the file "obl_ecc_lsp.asc".'
+if (.not. found_lsp) write(*,*) 'The maximum reachable year for Lsp could not be find in the file "obl_ecc_lsp.asc".'
 
 write(*,*) 'Max. number of iterations for the obl. parameter =', max_obl_iter
@@ -208 +205 @@
 
 year_iter_max = min(max_obl_iter,max_ecc_iter,max_lsp_iter)
-if (year_iter_max > 0) then
-    write(*,*) 'So the max. number of iterations for the orbital criteria =', year_iter_max
-else
-    write(*,*) 'The max. number of iterations is not compatible with Laskar data in the file obl_ecc_lsp.asc.'
-    stop
-endif
+write(*,*) 'So the max. number of iterations for the orbital criteria =', year_iter_max
 
 END SUBROUTINE orbit_param_criterion

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

@@ -44 +44 @@
   real    :: unitastr
   real    :: xa, xb, ya, yb ! Variables for interpolation
+  logical :: found_year     ! Flag variable
 
 ! **********************************************************************
@@ -64 +65 @@
 write(*,*) 'recomp_orb_param, Old Lsp  =', Lsp
 
+found_year = .false.
 do ilask = last_ilask,2,-1
     xa = yearlask(ilask)
@@ -86 +88 @@
             ya = lsplask(ilask)
             yb = lsplask(ilask - 1)
-            if (yb - ya > 300.) then ! If modulo is "crossed"
-                if (ya < yb) then ! Increasing
-                    yb = yb - 360.
-                else ! Decreasing
+            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 = (Year - xa)*(yb - ya)/(xb - xa) + ya
-            Lsp = modulo(Lsp,360.)
+            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
+
+if (.not. found_year) then
+    write(*,*) 'The new year could not be found in the file "obl_ecc_lsp.asc".'
+    stop
+endif
+
 halfaxe = 227.94
 Lsp = Lsp*pi/180.