Changeset 1382

Timestamp:

Mar 4, 2015, 10:03:37 AM (10 years ago)

Author:

emillour

Message:

Mars GCM:
Some cleanup in iniorbit.F
FF+EM

Location:

trunk/LMDZ.MARS

Files:

• README (modified) (1 diff)
• libf/phymars/datareadnc.F (modified) (1 diff)
• libf/phymars/iniorbit.F (modified) (2 diffs)
• libf/phymars/planete_h.F90 (modified) (1 diff)

trunk/LMDZ.MARS/README

@@ -2158 +2158 @@
   http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir
   or /u/lmdz/WWW/planets/mars/datadir for local users.
+
+== 04/03/2015 == FF+EM
+- Some cleanup in iniorbit.F

trunk/LMDZ.MARS/libf/phymars/datareadnc.F

@@ -113 +113 @@
       write(*,*) 'datareadnc: opening file surface.nc'
 
-      datafile="/u/forget/WWW/datagcm/datafile" ! default path to surface.nc
+      datafile="/u/lmdz/WWW/planets/mars/datadir" ! default path to surface.nc
       call getin("datadir",datafile) ! but users may specify another path
       

trunk/LMDZ.MARS/libf/phymars/iniorbit.F

@@ -1 +1 @@
       SUBROUTINE iniorbit
      $     (paphelie,pperiheli,pyear_day,pperi_day,pobliq)
-      use planete_h
-      USE comcstfi_h
+      use planete_h, only: aphelie, periheli, year_day, peri_day,
+     &                     obliquit, unitastr, e_elips, p_elips,
+     &                     timeperi
+      use comcstfi_h, only: pi
       IMPLICIT NONE
 
-c=======================================================================
-c
-c   Auteur:
-c   -------
-c     Frederic Hourdin      22 Fevrier 1991
-c
-c   Objet:
-c   ------
-c    Initialisation du sous programme orbite qui calcule
-c    a une date donnee de l'annee de duree year_day commencant
-c    a l'equinoxe de printemps et dont le perihelie se situe
-c    a la date peri_day, la distance au soleil et la declinaison.
-c
-c   Interface:
-c   ----------
-c   - Doit etre appele avant d'utiliser orbite.
-c   - initialise une partie du common planete.h
-c
-c   Arguments:
-c   ----------
-c
-c   Input:
-c   ------
-c   aphelie       \   aphelie et perihelie de l'orbite
-c   periheli      /   en millions de kilometres.
-c
-c=======================================================================
+!=======================================================================
+! Initialisation of orbital parameters (stored in planete_h module)
+!=======================================================================
 
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
+!   Arguments:
+!   ----------
 
-c   Arguments:
-c   ----------
+      REAL,INTENT(IN) :: paphelie,pperiheli,pyear_day,pperi_day,pobliq
 
-      REAL paphelie,pperiheli,pyear_day,pperi_day,pobliq
-
-c   Local:
-c   ------
+!   Local:
+!   ------
 
       REAL zxref,zanom,zz,zx0,zdx
       INTEGER iter
 
-c-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
 
       pi=2.*asin(1.)
@@ -58 +32 @@
       peri_day=pperi_day
 
-      PRINT*,'Perihelie en Mkm  ',periheli
-      PRINT*,'Aphelise  en Mkm  ',aphelie 
-      PRINT*,'obliquite en degres  :',obliquit
-      unitastr=149.597927
+      PRINT*,'iniorbit: Perihelion in Mkm  ',periheli
+      PRINT*,'iniorbit: Aphelion  in Mkm  ',aphelie 
+      PRINT*,'iniorbit: Obliquity in degrees  :',obliquit
+      unitastr=149.597927 ! 1 UA, in Mkm
       e_elips=(aphelie-periheli)/(periheli+aphelie)
       p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr
 
-      print*,'e_elips',e_elips
-      print*,'p_elips',p_elips
+      print*,'iniorbit: e_elips',e_elips
+      print*,'iniorbit: p_elips',p_elips
 
-c-----------------------------------------------------------------------
-c calcul de l'angle polaire et de la distance au soleil :
-c -------------------------------------------------------
+!-----------------------------------------------------------------------
+! compute polar angle and distance to the Sun:
+! -------------------------------------------------------
 
-c  calcul de l'zanomalie moyenne
+!  compute mean anomaly zanom
 
       zz=(year_day-pperi_day)/year_day
       zanom=2.*pi*(zz-nint(zz))
       zxref=abs(zanom)
-      PRINT*,'zanom  ',zanom
+      PRINT*,'iniorbit: zanom  ',zanom
 
-c  resolution de l'equation horaire  zx0 - e * sin (zx0) = zxref
-c  methode de Newton
+!  solve equation  zx0 - e * sin (zx0) = zxref for eccentric anomaly zx0
+!  using Newton method
 
       zx0=zxref+e_elips*sin(zxref)
-      DO 110 iter=1,100
+      DO iter=1,100
          zdx=-(zx0-e_elips*sin(zx0)-zxref)/(1.-e_elips*cos(zx0))
-         if(abs(zdx).le.(1.e-12)) goto 120
+         if(abs(zdx).le.(1.e-12)) exit
          zx0=zx0+zdx
-110   continue
-120   continue
+      ENDDO
+
       zx0=zx0+zdx
       if(zanom.lt.0.) zx0=-zx0
-      PRINT*,'zx0   ',zx0
-
-c zteta est la longitude solaire
+      PRINT*,'iniorbit: zx0   ',zx0
 
       timeperi=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
-      PRINT*,'longitude solaire du perihelie timeperi = ',timeperi
+      PRINT*,'iniorbit: Perihelion solar long. Ls (deg)=',
+     &       360.-timeperi*180./pi
 
-      RETURN
       END

trunk/LMDZ.MARS/libf/phymars/planete_h.F90

@@ -2 +2 @@
       IMPLICIT NONE
 
-      REAL aphelie
-      REAL periheli
-      REAL year_day
-      REAL peri_day
-      REAL obliquit
+      REAL aphelie ! Aphelion, in Mkm
+      REAL periheli ! Perihelion, in Mkm
+      REAL year_day ! number of days in the year
+      REAL peri_day ! date of perihelion, in days
+      REAL obliquit ! obliquity of the planet, in degrees
       REAL lmixmin
       REAL emin_turb
       REAL coefvis
       REAL coefir
-      REAL timeperi
-      REAL e_elips
+      REAL timeperi ! angle of the perihelion, in rad
+      REAL e_elips ! orbit eccentricity
       REAL p_elips
-      REAL unitastr
+      REAL unitastr ! Astronomical unit AU, in Mkm
 
       END MODULE planete_h