Changeset 4190 for dynamico_lmdz/simple_physics/phyparam/physics

Timestamp:

Dec 20, 2019, 12:07:05 AM (6 years ago)

Author:

dubos

Message:

simple_physics : cleanup astronomy

Location:

dynamico_lmdz/simple_physics/phyparam/physics

Files:

• astronomy.F90 (moved) (moved from dynamico_lmdz/simple_physics/phyparam/param/solarlong.F) (1 diff)
• phys_const.F90 (modified) (1 diff)
• planet.F90 (modified) (1 diff)

dynamico_lmdz/simple_physics/phyparam/physics/astronomy.F90

@@ -1 @@
-      SUBROUTINE solarlong(pday,psollong)
-      IMPLICIT NONE
+MODULE astronomy
+  IMPLICIT NONE
+  SAVE
 
-c=======================================================================
-c
-c   Objet:
-c   ------
-c
-c      Calcul de la distance soleil-planete et de la declinaison
-c   en fonction du jour de l'annee.
-c
-c
-c   Methode:
-c   --------
-c
-c      Calcul complet de l'elipse
-c
-c   Interface:
-c   ----------
-c
-c      Uncommon comprenant les parametres orbitaux.
-c
-c   Arguments:
-c   ----------
-c
-c   Input:
-c   ------
-c   pday          jour de l'annee (le jour 0 correspondant a l'equinoxe)
-c   lwrite        clef logique pour sorties de controle
-c
-c   Output:
-c   -------
-c   pdist_sol     distance entre le soleil et la planete
-c                 ( en unite astronomique pour utiliser la constante 
-c                  solaire terrestre 1370 Wm-2 )
-c   pdecli        declinaison ( en radians )
-c
-c=======================================================================
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
+  REAL :: aphelie, periheli, year_day, peri_day, obliquit, &
+       timeperi, e_elips,p_elips
+  REAL, PARAMETER :: unitastr=149.597927, & ! millions of km
+       pi=2.*ASIN(1.)
+CONTAINS
+ 
+  SUBROUTINE solarlong(pday,psollong)
 
-#include "comorbit.h"
+!=======================================================================
+!      Calcul de la distance soleil-planete et de la declinaison
+!   en fonction du jour de l'annee.
+!
+!   Methode:
+!   --------
+!      Calcul complet de l'ellipse
+!
+!   Input:
+!   ------
+!   pday          jour de l'annee (le jour 0 correspondant a l'equinoxe)
+!   lwrite        clef logique pour sorties de controle
+!
+!   Output:
+!   -------
+!   pdist_sol     distance entre le soleil et la planete
+!                 ( en unite astronomique pour utiliser la constante 
+!                  solaire terrestre 1370 Wm-2 )
+!   pdecli        declinaison ( en radians )
+!
+!=======================================================================
 
-c arguments:
-c ----------
+    USE planet
+    REAL, INTENT(IN) :: pday
+    REAL, INTENT(OUT) :: psollong
+    LOGICAL, PARAMETER ::  lwrite=.TRUE.
+    REAL pdist_sol, pdecli
 
-      REAL pday,pdist_sol,pdecli,psollong
-      LOGICAL lwrite
+! Local:
+! ------
+    REAL zanom,xref,zx0,zdx,zteta,zz
+    INTEGER iter
 
-c Local:
-c ------
+!--------------------------------------------------------
+! calcul de l'angle polaire et de la distance au soleil :
+! -------------------------------------------------------
 
-      REAL zanom,xref,zx0,zdx,zteta,zz
-      INTEGER iter
+!  calcul de l'zanomalie moyenne
 
+    zz=(pday-peri_day)/year_day
+    zanom=2.*pi*(zz-nint(zz))
+    xref=abs(zanom)
 
-c-----------------------------------------------------------------------
-c calcul de l'angle polaire et de la distance au soleil :
-c -------------------------------------------------------
-
-c  calcul de l'zanomalie moyenne
-
-      zz=(pday-peri_day)/year_day
-      zanom=2.*pi*(zz-nint(zz))
-      xref=abs(zanom)
-
-c  resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
-c  methode de Newton
-
-      zx0=xref+e_elips*sin(xref)
-      DO 110 iter=1,10
-         zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
-         if(abs(zdx).le.(1.e-7)) goto 120
-         zx0=zx0+zdx
-110   continue
-120   continue
-      zx0=zx0+zdx
-      if(zanom.lt.0.) zx0=-zx0
-
-c zteta est la longitude solaire
-
-      zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
-
-      psollong=zteta-timeperi
-
-      IF(psollong.LT.0.) psollong=psollong+2.*pi
-      IF(psollong.GT.2.*pi) psollong=psollong-2.*pi
-c-----------------------------------------------------------------------
-c   sorties eventuelles:
-c   ---------------------
-
-      IF (lwrite) THEN
-         PRINT*,'jour de l"annee   :',pday
-         PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol
-         PRINT*,'declinaison (en degres) :',pdecli*180./pi
-      ENDIF
-
-      RETURN
-      END
+!  resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
+!  methode de Newton
+    
+    zx0=xref+e_elips*sin(xref)
+    DO iter=1,10
+       zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
+       zx0=zx0+zdx
+    END DO
+    zx0=zx0+zdx
+    if(zanom.lt.0.) zx0=-zx0
+    
+    ! zteta est la longitude solaire
+    
+    zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
+    
+    psollong=zteta-timeperi
+    
+    IF(psollong.LT.0.) psollong=psollong+2.*pi
+    IF(psollong.GT.2.*pi) psollong=psollong-2.*pi
+    !-----------------------------------------------------------------------
+    !   sorties eventuelles:
+    !   ---------------------
+    
+    IF (lwrite) THEN
+       PRINT*,'jour de l"annee   :',pday
+       PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol
+       PRINT*,'declinaison (en degres) :',pdecli*180./pi
+    ENDIF
+    
+  END SUBROUTINE solarlong
+  
+END MODULE astronomy

dynamico_lmdz/simple_physics/phyparam/physics/phys_const.F90

@@ -1 +1 @@
 MODULE phys_const
   
-  REAL pi,rad,g,r,cpp,rcp,dtphys,unjours,mugaz,omeg
+  REAL rad,g,r,cpp,rcp,dtphys,unjours,mugaz,omeg
 
 END MODULE phys_const

dynamico_lmdz/simple_physics/phyparam/physics/planet.F90

@@ -3 +3 @@
   SAVE
 
-  REAL :: aphelie, periheli, year_day, peri_day, obliquit, &
-       coefvis, coefir
+  REAL :: coefvis, coefir
   
 END MODULE planet