Changeset 1384 for trunk/LMDZ.GENERIC/libf/phystd/iniorbit.F

Timestamp:

Mar 4, 2015, 5:00:33 PM (10 years ago)

Author:

emillour

Message:

Generic GCM:

• Some code cleanup: turning comcstfi.h into module comcstfi_mod.F90

EM

File:

• trunk/LMDZ.GENERIC/libf/phystd/iniorbit.F (modified) (3 diffs)

trunk/LMDZ.GENERIC/libf/phystd/iniorbit.F

@@ -4 +4 @@
       USE planete_mod, only: apoastr, periastr, year_day, obliquit,
      &                       peri_day, e_elips, p_elips, timeperi
+      use comcstfi_mod, 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 periastre 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_mod
-c
-c   Arguments:
-c   ----------
-c
-c   Input:
-c   ------
-c   apoastr       \   apoastron and periastron of the orbit in AU
-c   periastr      /   
-c
-c=======================================================================
-
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
-
-!#include "planete.h"
-#include "comcstfi.h"
+!=======================================================================
+! Initialisation of orbital parameters (stored in planete_h module)
+!=======================================================================
 
 c   Arguments:
 c   ----------
 
-      REAL papoastr,pperiastr,pyear_day,pperi_day,pobliq
+      REAL,INTENT(IN) :: papoastr,pperiastr,pyear_day,pperi_day,pobliq
 
 c   Local:
@@ -80 +50 @@
 
  
-      PRINT*,'Periastron in AU  ',periastr
-      PRINT*,'Apoastron in AU  ',apoastr 
-      PRINT*,'Obliquity in degrees  :',obliquit
+      PRINT*,'iniorbit: Periastron in AU  ',periastr
+      PRINT*,'iniorbit: Apoastron in AU  ',apoastr 
+      PRINT*,'iniorbit: Obliquity in degrees  :',obliquit
 
 
@@ -88 +58 @@
       p_elips=0.5*(periastr+apoastr)*(1-e_elips*e_elips)
 
-      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*,'Solar longitude of periastron timeperi = ',timeperi
+      PRINT*,'iniorbit: Perihelion solar long. Ls (deg)=',
+     &       360.-timeperi*180./pi
 
-      RETURN
       END