Changeset 3038 for trunk/LMDZ.MARS/deftank/pem/modify_startfi_orbit.sh

Timestamp:

Sep 8, 2023, 7:26:52 PM (16 months ago)

Author:

jbclement

Message:

Mars PEM:

• New version of scripts (launch_pem.sh and exeGCM.sh) to launch the chained simulation of GCM and PEM runs. It should be simpler and clearer for users;
• Update of README accordingly;
• ob_ex_lsp.asc is now in Earth years like original Laskar's data + rename obl_ecc_lsp.asc;
• Some improvements for the script modify_startfi_orbit.sh;
• Some changes in run_PEM.def and xml files.

JBC

File:

• trunk/LMDZ.MARS/deftank/pem/modify_startfi_orbit.sh (modified) (3 diffs)

trunk/LMDZ.MARS/deftank/pem/modify_startfi_orbit.sh

@@ -4 +4 @@
 name_file="startfi.nc"
 
-new_obl=30.0286
-new_exc=0.08212066
-new_Lsp=347.443947402
+new_obl=25.18941689
+new_ecc=0.09340902
+new_Lsp=251.06881714
+new_iniLs=0.
 
 
@@ -12 +13 @@
 year_day=669.
 halfaxe=227.94
-pi=$(echo "4*a(1)" | bc -l)
-degrad=$(echo "360./(2.*$pi)" | bc -l)
+pi=$(echo "4.*a(1.)" | bc -l)
+degrad=$(echo "180./$pi" | bc -l)
+
+periheli=$(echo "$halfaxe*(1. - $new_ecc)" | bc -l)
+aphelie=$(echo "$halfaxe*(1. + $new_ecc)" | bc -l)
+
 tan=$(echo "s(0.5*$new_Lsp/$degrad)/c(0.5*$new_Lsp/$degrad)" | bc -l)
-zx0=$(echo "-2.*a($tan*sqrt((1.-$new_exc)/(1.+$new_exc)))" | bc -l)
-if [ $(echo "$zx0 <= 0" | bc -l) -eq 1 ]; then
+zx0=$(echo "-2.*a($tan*sqrt((1. - $new_ecc)/(1. + $new_ecc)))" | bc -l)
+if [ $(echo "$zx0 <= 0." | bc -l) -eq 1 ]; then
   zx0=$(echo "$zx0 + 2.*$pi" | bc -l)
 fi
-solp=$(echo "$year_day*(1.-($zx0 - $new_exc*s($zx0))/(2*$pi))" | bc -l)
+peri_day=$(echo "$year_day*(1. - ($zx0 - $new_ecc*s($zx0))/(2.*$pi))" | bc -l)
+
+ztheta=$(echo "($new_iniLs - $new_Lsp)/$degrad" | bc -l)
+tan=$(echo "s(0.5*$ztheta)/c(0.5*$ztheta)" | bc -l)
+zx0=$(echo "2.*a($tan*sqrt((1. - $new_ecc)/(1. + $new_ecc)))" | bc -l)
+xref=$(echo "$zx0 - $new_ecc*s($zx0)" | bc -l)
+new_inisol=$(echo "$peri_day + $xref*$year_day/(2.*$pi)" | bc -l)
+if [ $(echo "$new_inisol < 0." | bc -l) -eq 1 ]; then
+  new_inisol=$(echo "$new_inisol + $year_day" | bc -l)
+fi
+if [ $(echo "$new_inisol >= $year_day" | bc -l) -eq 1 ]; then
+  new_inisol=$(echo "$new_inisol - $year_day" | bc -l)
+fi
 
 # Update the netCDF file
@@ -26 +43 @@
 # controle(17) = peri_day ! date of perihelion (sols since N. spring)
 # controle(18) = obliquit ! Obliquity of the planet (deg)  ̃23.98
+# controle(3)  = day_ini + int(time) ! initial day
 ncap2 -O -s "controle(17)=$new_obl" \
-      -s "controle(14)=$halfaxe*(1-$new_exc)" \
-      -s "controle(15)=$halfaxe*(1+$new_exc)" \
-      -s "controle(16)=$solp" \
-      $name_file $name_file.temp
+         -s "controle(14)=$periheli" \
+         -s "controle(15)=$aphelie" \
+         -s "controle(16)=$peri_day" \
+         -s "controle(2)=$new_inisol" \
+         $name_file $name_file.temp
 
 # Rename the temporary file back to the original filename
 mv $name_file.temp $name_file
 
-echo New obliquit = $new_obl
-echo New periheli = $(echo "$halfaxe*(1-$new_exc)" | bc -l)
-echo New aphelie = $(echo "$halfaxe*(1+$new_exc)" | bc -l)
-echo New peri_day = $solp
+echo "In $name_file:"
+echo "New obliquit     = $new_obl"
+echo "New eccentricity = $new_ecc -> new periheli = $periheli"
+echo "                              -> new aphelie  = $aphelie"
+echo "New Lsp          = $new_Lsp -> new peri_day = $peri_day"
+echo "New initial Ls   = $new_iniLs -> New initial sol = $new_inisol"
+echo "Done!"