source: trunk/UTIL/PYTHON/time.F @ 578

Last change on this file since 578 was 403, checked in by acolaitis, 13 years ago

PYTHON.

############################################
1/------------------------------------------
############################################
# First version of a new command line routine
# that takes in input a diagfi and a MCS datafile
# as formatted by Luca Montabone
# and find the values in the diagfi temp field
# corresponding to daytime and nightime localtimes in MCS file
# (at each longitude and latitude, which is especially important at the poles
# where the spacecraft covers multiple local times for each bin)
# The selected diagfi data is then binned at the Ls of the MCS file
# and written as a diagfi_MCS.nc file.
#
# Example: (you have to be connected on Penn, as the file is on the local disk)

/san0/acolmd/SVN/trunk/UTIL/PYTHON/mcs.py -f /home/local_home/colaitis/Polar_simulations/Polar_All_Th/diagfi16.nc -m /san0/acolmd/DATA/MCS/MCS_processeddata/MCSdata_binned_MY29.nc

# It is then easy to compare the resulting gcm re-arranged data to the MCS dataset using
# your favorite pp command. Note that in this first version of the routine, only the temperature field is
# recomputed. A more complex version will come out later, with variable selection and conservation of
# fields needed to then use zrecast or hrecast.

############################################
2/------------------------------------------
############################################

# Added python versions of ls2sol and sol2ls
# which can work either with single values or arrays
# These routines are in times.py so that you can add
# in your script:

from times import ls2sol,sol2ls

############################################
2/------------------------------------------
############################################

# Added new routine in mymath, which finds the index of the
# element nearest to "value" in "array". The array
# can be masked or not, the output is an array.
# If the value is a NaN, the routine returns a NaN.
# One can specify the option "strict=True" so that
# the routine will return a NaN if the value asked is not
# between the min and max of the array.
File size: 3.0 KB
Line 
1cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2      real function sol2ls(sol)
3
4c  convert a given martian day number (sol)
5c  into corresponding solar longitude, Ls (in degr.),
6c  where sol=0=Ls=0 is the
7c  northern hemisphere spring equinox.
8
9      implicit none
10
11c     input
12      real    sol
13c     output
14      real    ls
15
16c     local variables
17      double precision    year_day,peri_day,timeperi,e_elips
18      double precision    pi,radtodeg
19c number of martian days (sols) in a martian year
20      parameter (year_day=668.6)
21c perihelion date (in sols)
22      parameter (peri_day=485.35)
23c orbital eccentricity
24      parameter (e_elips=0.09340)
25      parameter (pi=3.14159265358979)
26c  radtodeg: 180/pi
27      parameter (radtodeg=57.2957795130823)
28c  timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
29      parameter (timeperi=1.90258341759902)
30
31      double precision    zanom,xref,zx0,zdx,zteta,zz
32c  xref: mean anomaly, zx0: eccentric anomaly, zteta: true anomaly     
33      integer iter
34
35      zz=(sol-peri_day)/year_day
36      zanom=2.*pi*(zz-nint(zz))
37      xref=dabs(zanom)
38
39c  The equation zx0 - e * sin (zx0) = xref, solved by Newton
40      zx0=xref+e_elips*dsin(xref)
41      do 110 iter=1,10
42         zdx=-(zx0-e_elips*dsin(zx0)-xref)/(1.-e_elips*dcos(zx0))
43         if(dabs(zdx).le.(1.d-7)) then
44           goto 120
45         endif
46         zx0=zx0+zdx
47  110 continue
48  120 continue
49      zx0=zx0+zdx
50      if(zanom.lt.0.) zx0=-zx0
51
52c compute true anomaly zteta, now that eccentric anomaly zx0 is known
53      zteta=2.*datan(dsqrt((1.+e_elips)/(1.-e_elips))*dtan(zx0/2.))
54
55c compute Ls
56      ls=real(zteta-timeperi)
57      if(ls.lt.0.) ls=ls+2.*real(pi)
58      if(ls.gt.2.*pi) ls=ls-2.*real(pi)
59c convert Ls in deg.
60      ls=real(radtodeg)*ls
61
62      sol2ls = ls
63      return
64      end
65
66cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
67      real function ls2sol(ls)
68
69c  Returns solar longitude, Ls (in deg.), from day number (in sol),
70c  where sol=0=Ls=0 at the northern hemisphere spring equinox
71
72      implicit none
73
74c  Arguments:
75      real ls
76
77c  Local:
78      double precision xref,zx0,zteta,zz
79      double precision peri_day,timeperi,e_elips
80      double precision pi,degrad
81      double precision,parameter :: year_day=668.6
82      parameter (peri_day=485.35)
83      parameter (timeperi=1.90258341759902)
84      parameter (e_elips=0.0934) 
85      parameter (pi=3.14159265358979)
86      parameter (degrad=57.2957795130823)
87
88      if (abs(ls).lt.1.0e-5) then
89         if (ls.ge.0.0) then
90            ls2sol = 0.0
91         else
92            ls2sol = real(year_day)
93         end if
94         return
95      end if
96
97      zteta = ls/degrad + timeperi
98      zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips)))
99      xref = zx0-e_elips*dsin(zx0)
100      zz = xref/(2.*pi)
101      ls2sol = real(zz*year_day + peri_day)
102      if (ls2sol.lt.0.0) ls2sol = ls2sol + real(year_day)
103      if (ls2sol.ge.year_day) ls2sol = ls2sol - real(year_day)
104
105      return
106      end
Note: See TracBrowser for help on using the repository browser.