Changeset 4191 for dynamico_lmdz/simple_physics/phyparam/physics

Timestamp:

Dec 20, 2019, 12:41:10 AM (6 years ago)

Author:

dubos

Message:

simple_physics : cleanup astronomy

File:

• dynamico_lmdz/simple_physics/phyparam/physics/astronomy.F90 (modified) (1 diff)

dynamico_lmdz/simple_physics/phyparam/physics/astronomy.F90

@@ -85 +85 @@
   END SUBROUTINE solarlong
   
+  SUBROUTINE iniorbit
+    !=======================================================================
+    !
+    !   Auteur:
+    !   -------
+    !     Frederic Hourdin      22 Fevrier 1991
+    !
+    !   Objet:
+    !   ------
+    !    Initialisation du sous programme orbite qui calcule
+    !    a une date donnee de l'annee de duree year_day commencant
+    !    a l'equinoxe de printemps et dont le perihelie se situe
+    !    a la date peri_day, la distance au soleil et la declinaison.
+    !
+    !   Interface:
+    !   ----------
+    !   - initialise certaines variables de ce module
+    !   - Doit etre appele avant d'utiliser orbite.
+    !
+    !   Arguments:
+    !   ----------
+    !
+    !   Input:
+    !   ------
+    !   aphelie       \   aphelie et perihelie de l'orbite
+    !   periheli      /   en millions de kilometres.
+    !
+    !=======================================================================
+    
+    !-----------------------------------------------------------------------
+    
+    !   Local:
+    !   ------
+    REAL zxref,zanom,zz,zx0,zdx
+    INTEGER iter
+    
+    !-----------------------------------------------------------------------
+    
+    PRINT*,'Perihelie en Mkm  ',periheli
+    PRINT*,'Aphelise  en Mkm  ',aphelie 
+    PRINT*,'obliquite en degres  :',obliquit
+    
+    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
+    
+    !-----------------------------------------------------------------------
+    ! calcul de l'angle polaire et de la distance au soleil :
+    ! -------------------------------------------------------
+    
+    !  calcul de l'zanomalie moyenne
+    
+    zz=(year_day-peri_day)/year_day
+    zanom=2.*pi*(zz-nint(zz))
+    zxref=abs(zanom)
+    PRINT*,'zanom  ',zanom
+    
+    !  resolution de l'equation horaire  zx0 - e * sin (zx0) = zxref
+    !  methode de Newton
+    
+    zx0=zxref+e_elips*sin(zxref)
+    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
+       zx0=zx0+zdx
+    END DO
+    
+    zx0=zx0+zdx
+    if(zanom.lt.0.) zx0=-zx0
+    PRINT*,'zx0   ',zx0
+    
+    ! zteta est la longitude solaire
+    
+    timeperi=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
+    PRINT*,'longitude solaire du perihelie timeperi = ',timeperi
+    
+    
+  END SUBROUTINE iniorbit
+
 END MODULE astronomy
+