Changeset 4191

Timestamp:

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

Author:

dubos

Message:

simple_physics : cleanup astronomy

Location:

dynamico_lmdz/simple_physics/phyparam

Files:

• param/iniorbit.F (deleted)
• param/iniphyparam.F (modified) (1 diff)
• param/phyparam.F (modified) (1 diff)
• param/zenang.F (modified) (3 diffs)
• physics/astronomy.F90 (modified) (1 diff)

dynamico_lmdz/simple_physics/phyparam/param/iniphyparam.F

@@ -64 +64 @@
       INTEGER ig,ierr,offset
  
-      EXTERNAL inifrict,iniorbit,orbite
+      EXTERNAL inifrict,orbite
  
       print*,'INIPHYPARAM'

dynamico_lmdz/simple_physics/phyparam/param/phyparam.F

@@ -253 +253 @@
          print*,tsoil(ngrid/2+1,nsoilmx/2+2)
          print*,'TS ',tsurf(igout),tsoil(igout,5)
-         CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
+         CALL iniorbit
 
          if (.not.callrad) then

dynamico_lmdz/simple_physics/phyparam/param/zenang.F

@@ -1 +1 @@
       SUBROUTINE zenang(klon,longi,gmtime,pdtrad,lat,long,
      s                  pmu0,frac)
+      USE astronomy
       IMPLICIT none
 c=============================================================
@@ -24 +25 @@
 c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
 c================================================================
-#include "comorbit.h"
       integer klon
 c================================================================
@@ -155 +155 @@
       END
 c===================================================================
-#ifdef PASLA
-      SUBROUTINE zenith (klon,longi, gmtime, lat, long,
-     s                   pmu0, fract)
-      IMPLICIT none
-c
-c Auteur(s): Z.X. Li (LMD/ENS)
-c
-c Objet: calculer le cosinus de l'angle zenithal du soleil en
-c        connaissant la declinaison du soleil, la latitude et la
-c        longitude du point sur la terre, et le temps universel
-c
-c Arguments d'entree:
-c     longi  : declinaison du soleil (en degres)
-c     gmtime : temps universel en second qui varie entre 0 et 86400
-c     lat    : latitude en degres
-c     long   : longitude en degres
-c Arguments de sortie:
-c     pmu0   : cosinus de l'angle zenithal
-c
-c====================================================================
-#include "YOMCST.h"
-c====================================================================
-      INTETER klon
-      REAL longi, gmtime
-      REAL lat(klon), long(klon), pmu0(klon), fract(klon)
-c=====================================================================
-      INTEGER n
-      REAL zpi, zpir, omega, zgmtime
-      REAL incl, lat_sun, lon_sun
-c----------------------------------------------------------------------
-      zpi = 4.0*ATAN(1.0)
-      zpir = zpi / 180.0
-      zgmtime=gmtime*86400.
-c
-      incl=R_incl * zpir
-c
-      lon_sun = longi * zpir
-      lat_sun = ASIN (SIN(lon_sun)*SIN(incl) )
-c
-c--initialisation a la nuit
-c
-      DO n =1, klon
-        pmu0(n)=0.
-        fract(n)=0.0
-      ENDDO
-c
-c 1 degre en longitude = 240 secondes en temps
-c
-      DO n = 1, klon
-         omega = zgmtime + long(n)*86400.0/360.0
-         omega = omega / 86400.0 * 2.0 * zpi
-         omega = MOD(omega + 2.0 * zpi, 2.0 * zpi)
-         omega = omega - zpi
-         pmu0(n) = sin(lat(n)*zpir) * sin(lat_sun)
-     .           + cos(lat(n)*zpir) * cos(lat_sun)
-     .           * cos(omega)
-         pmu0(n) = MAX (pmu0(n), 0.0)
-         IF (pmu0(n).GT.1.E-6) fract(n)=1.0
-      ENDDO
-c
-      RETURN
-      END
-#endif

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
+