Changeset 3549 for trunk/LMDZ.MARS

Timestamp:

Dec 12, 2024, 2:10:37 PM (11 days ago)

Author:

jbclement

Message:

Mars PCM:
In the 1D model, merging under the 'paleomars' flag of the modifications for orbital parameters taken from "callphys.def".
JBC

Location:

trunk/LMDZ.MARS

Files:

• changelog.txt (modified) (1 diff)
• libf/phymars/dyn1d/init_testphys1d_mod.F90 (modified) (3 diffs)

trunk/LMDZ.MARS/changelog.txt

@@ -4742 +4742 @@
 == 06/12/2024 == JBC
 Correction of the initial time for the calendar in the "diagsoil.nc" file which caused an abnormal closing of ncview window when displaying variables.
+
+== 12/12/2024 == JBC
+In the 1D model, merging under the 'paleomars' flag of the modifications for orbital parameters taken from "callphys.def".

trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90

@@ -335 +335 @@
 ! ------------------
 tauvis = 0.2 ! Default value for tauvis (dust opacity)
-write(*,'("Reference dust opacity at ",f4.0," Pa?")')odpref
+write(*,'("Reference dust opacity at ",f4.0," Pa?")') odpref
 call getin("tauvis",tauvis)
-write(*,*) " tauvis = ",tauvis
+write(*,*) " tauvis =",tauvis
 
 ! Orbital parameters
@@ -344 +344 @@
     paleomars = .false. ! Default: no water ice reservoir
     call getin("paleomars",paleomars)
+    write(*,*) "paleomars=", paleomars
     if (paleomars) then
-        write(*,*) "paleomars=", paleomars
-        write(*,*) "Orbital parameters from callphys.def"
-        write(*,*) "Enter eccentricity & Lsperi"
-        write(*,*) 'Martian eccentricity (0<e<1)?'
-        call getin('eccentric ',eccentric)
-        write(*,*)"eccentric =",eccentric
-        write(*,*) 'Solar longitude of perihelion (0<Ls<360)?'
-        call getin('Lsperi',Lsperi )
-        write(*,*)"Lsperi=",Lsperi
-        Lsperi = Lsperi*pi/180.0 ! Put it in rad for peri_day
-        periheli = halfaxe*(1 - eccentric)
-        aphelie = halfaxe*(1 + eccentric)
-        call call_dayperi(Lsperi,eccentric,peri_day,year_day)
-        write(*,*) "Corresponding orbital params for GCM"
-        write(*,*) " periheli = ",periheli
-        write(*,*) " aphelie = ",aphelie
-        write(*,*) "date of perihelion (sol)",peri_day
-    else
-        write(*,*) "paleomars=", paleomars
-        write(*,*) "Default present-day orbital parameters"
-        write(*,*) "Unless specified otherwise"
-        write(*,*)'Min. distance Sun-Mars (Mkm)?'
-        call getin("periheli",periheli)
-        write(*,*) " periheli = ",periheli
-
-        write(*,*)'Max. distance Sun-Mars (Mkm)?'
-        call getin("aphelie",aphelie)
-        write(*,*) " aphelie = ",aphelie
-
-        write(*,*)'Day of perihelion?'
-        call getin("periday",peri_day)
-        write(*,*) " periday = ",peri_day
+        write(*,*) "Orbital parameters taken from callphys.def"
 
         write(*,*)'Obliquity?'
         call getin("obliquit",obliquit)
-        write(*,*) " obliquit = ",obliquit
+        write(*,*) " obliquit =",obliquit
+
+        write(*,*) 'Eccentricity (0<e<1)?'
+        call getin('eccentric',eccentric)
+        write(*,*)"eccentric =",eccentric
+        if (eccentric < 0. .or. eccentric > 1.) then
+            write(*,*) 'Value for ''eccentric'' is not between 0 and 1.'
+            error stop 'Please, specify a correct value!'
+        endif
+
+        call getin('Lsperi',Lsperi)
+        write(*,*) 'Solar longitude of perihelion (0<Ls<360)?'
+        write(*,*)"Lsperi =",Lsperi
+        if (eccentric < 0. .or. eccentric > 360.) then
+            write(*,*) 'Value for ''Lsperi'' is not between 0 and 360.'
+            error stop 'Please, specify a correct value!'
+        endif
+        
+        Lsperi = Lsperi*pi/180. ! Convert to rad for peri_day
+        periheli = halfaxe*(1. - eccentric)
+        aphelie = halfaxe*(1. + eccentric)
+        call call_dayperi(Lsperi,eccentric,peri_day,year_day)
+        write(*,*) "Corresponding orbital params for PCM"
+        write(*,*) "   periheli =",periheli
+        write(*,*) "   aphelie  =",aphelie
+        write(*,*) "   peri_day =",peri_day
+    else
+        write(*,*) "Default present-day orbital parameters"
      endif
+else
+    write(*,*) "Orbital parameters taken from startfi.nc"
 endif !(.not. therestartfi)
 
@@ -387 +386 @@
 ! ------------------
 latitude = 0. ! default value for latitude
-write(*,*)'latitude (in degrees)?'
+write(*,*) 'Latitude (in degrees)?'
 call getin("latitude",latitude(1))
-write(*,*) " latitude = ",latitude
+write(*,*) " latitude =",latitude
 latitude = latitude*pi/180.
 longitude = 0.