source: trunk/LMDZ.COMMON/libf/evolution/deftank/inipem_orbit.sh @ 3554

Last change on this file since 3554 was 3351, checked in by jbclement, 7 months ago

PEM:

  • Updates for several scripts in the deftank.
  • Addition of the bash script "clean.sh" to clean a folder after a PEM simulation.

JBC

  • Property svn:executable set to *
File size: 4.6 KB
RevLine 
[3181]1#!/bin/bash
2######################################################################
3### Script to modify the orbital parameters of a file "startfi.nc" ###
4### according to the date set in the file "run_PEM.def"            ###
5######################################################################
6
7# Name of the file to be modified
8name_file="startfi.nc"
9
10# Name of the file containing the orbital data
11orb_data="obl_ecc_lsp.asc"
12
[3210]13
14######################################################################
[3181]15# Check if files necessary for the script exist
16if [ ! -f "$orb_data" ]; then
17    echo "Error: file \"$orb_data\" not found!"
18    exit 1
19fi
20if [ ! -f "startfi.nc" ]; then
21    echo "Error: file \"$name_file\" not found!"
22    exit 1
23fi
24if [ ! -f "run_PEM.def" ]; then
25    echo "Error: file \"run_PEM.def\" not found!"
26    exit 1
27fi
28
29# Get the number of Earth years before present to start the PEM run from "run_PEM.def"
30eyears_bp_ini=$(grep 'year_earth_bp_ini=' run_PEM.def | cut -d '=' -f 2 | tr -d '[:space:]')
31echo "The starting date is $eyears_bp_ini Earth years according to \"run_PEM.def\"."
32
33# Get the new values for the orbital parameters
[3351]34yearlask=$(echo "$eyears_bp_ini/1000." | bc -l)
[3181]35found=false
36read -r y1 obl1 ecc1 lsp1 < "$orb_data"
37while read -r y2 obl2 ecc2 lsp2; do
38    if [ "$first_line" = true ]; then
39        first_line=false
40        continue # Skip the first line as it's already read outside the loop
41    fi
42    if [ "$(echo "$y1 >= $yearlask && $yearlask > $y2" | bc)" -eq 1 ]; then
43        found=true
44        # Calculate linear interpolations for each orbital parameter
45        new_obl=$(echo "($obl2 - $obl1)/($y2 - $y1)*($yearlask - $y1) + $obl1" | bc -l)
46        new_ecc=$(echo "($ecc2 - $ecc1)/($y2 - $y1)*($yearlask - $y1) + $ecc1" | bc -l)
47        if [ "$(echo "$lsp2 - $lsp1 > 300. || $lsp2 - $lsp1 < -300." | bc)" -eq 1 ]; then # If modulo is "crossed" through the interval
48            if [ "$(echo "$lsp1 < $lsp2" | bc -l)" -eq 1 ]; then # Lsp should be decreasing
49                lsp1=$(echo "$lsp1 + 360" | bc -l)
50            else # Lsp should be increasing
51                lsp2=$(echo "$lsp2 + 360" | bc -l)
52            fi
53        fi
54        new_Lsp=$(echo "($lsp2 - $lsp1)/($y2 - $y1)*($yearlask - $y1) + $lsp1" | bc -l)
55        echo "New obliquity    = $new_obl"
56        echo "New eccentricity = $new_ecc"
57        echo "New Lsp          = $new_Lsp"
58        echo "Done!"
59        break
60    fi
61    y1=$y2
62    obl1=$obl2
63    ecc1=$ecc2
64    lsp1=$lsp2
65done < "$orb_data"
66if [ "$found" = false ]; then
67    echo "Error: the current year could not be found in the file \"$orb_data\"!"
68    exit 1
69fi
70
71
72# Compute the modified values in the file
73if [ ! -f "$name_file" ]; then
74    echo "Error: file \"$name_file\" not found!"
75    exit 1
76fi
77year_day=669.
78halfaxe=227.94
79pi=$(echo "4.*a(1.)" | bc -l)
80degrad=$(echo "180./$pi" | bc -l)
81
82periheli=$(echo "$halfaxe*(1. - $new_ecc)" | bc -l)
83aphelie=$(echo "$halfaxe*(1. + $new_ecc)" | bc -l)
84
85tan=$(echo "s(0.5*$new_Lsp/$degrad)/c(0.5*$new_Lsp/$degrad)" | bc -l)
86zx0=$(echo "-2.*a($tan*sqrt((1. - $new_ecc)/(1. + $new_ecc)))" | bc -l)
87if [ $(echo "$zx0 <= 0." | bc -l) -eq 1 ]; then
88    zx0=$(echo "$zx0 + 2.*$pi" | bc -l)
89fi
90peri_day=$(echo "$year_day*(1. - ($zx0 - $new_ecc*s($zx0))/(2.*$pi))" | bc -l)
91
92ztheta=$(echo "($new_iniLs - $new_Lsp)/$degrad" | bc -l)
93tan=$(echo "s(0.5*$ztheta)/c(0.5*$ztheta)" | bc -l)
94zx0=$(echo "2.*a($tan*sqrt((1. - $new_ecc)/(1. + $new_ecc)))" | bc -l)
95xref=$(echo "$zx0 - $new_ecc*s($zx0)" | bc -l)
96new_inisol=$(echo "$peri_day + $xref*$year_day/(2.*$pi)" | bc -l)
97if [ $(echo "$new_inisol < 0." | bc -l) -eq 1 ]; then
98    new_inisol=$(echo "$new_inisol + $year_day" | bc -l)
99fi
100if [ $(echo "$new_inisol >= $year_day" | bc -l) -eq 1 ]; then
101    new_inisol=$(echo "$new_inisol - $year_day" | bc -l)
102fi
103
104# Update the netCDF file
105# controle(15) = periheli ! min. Sun-Mars distance (Mkm)  ̃206.66
106# controle(16) = aphelie  ! max. SUn-Mars distance (Mkm)  ̃249.22
107# controle(17) = peri_day ! date of perihelion (sols since N. spring)
108# controle(18) = obliquit ! Obliquity of the planet (deg)  ̃23.98
109# controle(3)  = day_ini + int(time) ! initial day
110ncap2 -O -s "controle(17)=$new_obl" \
111         -s "controle(14)=$periheli" \
112         -s "controle(15)=$aphelie" \
113         -s "controle(16)=$peri_day" \
114         -s "controle(2)=$new_inisol" \
115         $name_file $name_file
116
117echo "In \"$name_file\":"
118echo "New obliquit     = $new_obl"
119echo "New eccentricity = $new_ecc -> new periheli = $periheli"
120echo "                              -> new aphelie  = $aphelie"
121echo "New Lsp          = $new_Lsp -> new peri_day = $peri_day"
122echo "New initial Ls   = $new_iniLs -> New initial sol = $new_inisol"
123echo "Done!"
Note: See TracBrowser for help on using the repository browser.