Changeset 4190


Ignore:
Timestamp:
Dec 20, 2019, 12:07:05 AM (5 years ago)
Author:
dubos
Message:

simple_physics : cleanup astronomy

Location:
dynamico_lmdz/simple_physics/phyparam
Files:
7 edited
1 moved

Legend:

Unmodified
Added
Removed

  • dynamico_lmdz/simple_physics/phyparam/param/iniorbit.F

    r4176 r4190  
    11      SUBROUTINE iniorbit
    22     $     (paphelie,pperiheli,pyear_day,pperi_day,pobliq)
     3      USE astronomy
    34      IMPLICIT NONE
    45
    … …  
    3233
    3334c-----------------------------------------------------------------------
    34 c   Declarations:
    35 c   -------------
    36 
    37 #include "comorbit.h"
    3835
    3936c   Arguments:
    … …  
    5653      peri_day=pperi_day
    5754
    58       pi=2.*asin(1.)
    5955      PRINT*,'Perihelie en Mkm  ',periheli
    6056      PRINT*,'Aphelise  en Mkm  ',aphelie
    6157      PRINT*,'obliquite en degres  :',obliquit
    62       unitastr=149.597927
     58
    6359      e_elips=(aphelie-periheli)/(periheli+aphelie)
    6460      p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr

  • dynamico_lmdz/simple_physics/phyparam/param/iniphyparam.F

    r4189 r4190  
    1111      use comsaison
    1212      USE geometry_mod, ONLY : longitude,latitude,cell_area
    13       USE phys_const, ONLY : pi,rad,g,r,cpp,rcp,dtphys,unjours,mugaz
     13      USE phys_const, ONLY : rad,g,r,cpp,rcp,dtphys,unjours,mugaz
    1414      USE planet, ONLY : coefir, coefvis
    15       USE planet, ONLY : aphelie, periheli, peri_day, obliquit, year_day
     15      USE astronomy
    1616      USE vdif_mod, ONLY : lmixmin, emin_turb
    1717      IMPLICIT NONE
    … …  
    197197      coslon(:)=cos(long(:))
    198198
    199       pi=2.*asin(1.)
    200199      prad=rad
    201200      pg=g

  • dynamico_lmdz/simple_physics/phyparam/param/orbite.F

    r4176 r4190  
    11      SUBROUTINE orbite(pls,pdist_sol,pdecli)
     2      USE astronomy
    23      IMPLICIT NONE
    34
    … …  
    3233c   -------------
    3334
    34 #include "comorbit.h"
    35 
    3635c arguments:
    3736c ----------

  • dynamico_lmdz/simple_physics/phyparam/param/paramdef.F

    r4189 r4190  
    33      USE vdif_mod, ONLY : lmixmin, emin_turb
    44      USE planet, ONLY : coefir, coefvis
    5       USE planet, ONLY : aphelie, periheli, obliquit, peri_day, year_day
     5      USE astronomy
    66      IMPLICIT NONE
    77c-----------------------------------------------------------------------
    … …  
    2020      character str1*1
    2121      INTEGER ig,ierr
    22       real zz,z1,z2,pi
     22      real zz,z1,z2
    2323
    2424      REAL I_mer,I_ter,z0_mer,z0_ter
    … …  
    2929c-----------------------------------------------------------------------
    3030
    31       pi=2.*asin(1.)
    3231      year_day=360.
    3332      periheli=173.

  • dynamico_lmdz/simple_physics/phyparam/param/phyparam.F

    r4188 r4190  
    1212      USE phys_const
    1313      USE planet
     14      USE astronomy
    1415      USE vdif_mod, ONLY : vdif
    15 c
     16c     
    1617      IMPLICIT NONE
    1718c=======================================================================

  • dynamico_lmdz/simple_physics/phyparam/physics/astronomy.F90

    r4178 r4190  
    1       SUBROUTINE solarlong(pday,psollong)
    2       IMPLICIT NONE
     1MODULE astronomy
     2  IMPLICIT NONE
     3  SAVE
    34
    4 c=======================================================================
    5 c
    6 c   Objet:
    7 c   ------
    8 c
    9 c      Calcul de la distance soleil-planete et de la declinaison
    10 c   en fonction du jour de l'annee.
    11 c
    12 c
    13 c   Methode:
    14 c   --------
    15 c
    16 c      Calcul complet de l'elipse
    17 c
    18 c   Interface:
    19 c   ----------
    20 c
    21 c      Uncommon comprenant les parametres orbitaux.
    22 c
    23 c   Arguments:
    24 c   ----------
    25 c
    26 c   Input:
    27 c   ------
    28 c   pday          jour de l'annee (le jour 0 correspondant a l'equinoxe)
    29 c   lwrite        clef logique pour sorties de controle
    30 c
    31 c   Output:
    32 c   -------
    33 c   pdist_sol     distance entre le soleil et la planete
    34 c                 ( en unite astronomique pour utiliser la constante
    35 c                  solaire terrestre 1370 Wm-2 )
    36 c   pdecli        declinaison ( en radians )
    37 c
    38 c=======================================================================
    39 c-----------------------------------------------------------------------
    40 c   Declarations:
    41 c   -------------
     5  REAL :: aphelie, periheli, year_day, peri_day, obliquit, &
     6       timeperi, e_elips,p_elips
     7  REAL, PARAMETER :: unitastr=149.597927, & ! millions of km
     8       pi=2.*ASIN(1.)
     9CONTAINS
     10 
     11  SUBROUTINE solarlong(pday,psollong)
    4212
    43 #include "comorbit.h"
     13!=======================================================================
     14!      Calcul de la distance soleil-planete et de la declinaison
     15!   en fonction du jour de l'annee.
     16!
     17!   Methode:
     18!   --------
     19!      Calcul complet de l'ellipse
     20!
     21!   Input:
     22!   ------
     23!   pday          jour de l'annee (le jour 0 correspondant a l'equinoxe)
     24!   lwrite        clef logique pour sorties de controle
     25!
     26!   Output:
     27!   -------
     28!   pdist_sol     distance entre le soleil et la planete
     29!                 ( en unite astronomique pour utiliser la constante
     30!                  solaire terrestre 1370 Wm-2 )
     31!   pdecli        declinaison ( en radians )
     32!
     33!=======================================================================
    4434
    45 c arguments:
    46 c ----------
     35    USE planet
     36    REAL, INTENT(IN) :: pday
     37    REAL, INTENT(OUT) :: psollong
     38    LOGICAL, PARAMETER ::  lwrite=.TRUE.
     39    REAL pdist_sol, pdecli
    4740
    48       REAL pday,pdist_sol,pdecli,psollong
    49       LOGICAL lwrite
     41! Local:
     42! ------
     43    REAL zanom,xref,zx0,zdx,zteta,zz
     44    INTEGER iter
    5045
    51 c Local:
    52 c ------
     46!--------------------------------------------------------
     47! calcul de l'angle polaire et de la distance au soleil :
     48! -------------------------------------------------------
    5349
    54       REAL zanom,xref,zx0,zdx,zteta,zz
    55       INTEGER iter
     50!  calcul de l'zanomalie moyenne
    5651
     52    zz=(pday-peri_day)/year_day
     53    zanom=2.*pi*(zz-nint(zz))
     54    xref=abs(zanom)
    5755
    58 c-----------------------------------------------------------------------
    59 c calcul de l'angle polaire et de la distance au soleil :
    60 c -------------------------------------------------------
    61 
    62 c  calcul de l'zanomalie moyenne
    63 
    64       zz=(pday-peri_day)/year_day
    65       zanom=2.*pi*(zz-nint(zz))
    66       xref=abs(zanom)
    67 
    68 c  resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
    69 c  methode de Newton
    70 
    71       zx0=xref+e_elips*sin(xref)
    72       DO 110 iter=1,10
    73          zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
    74          if(abs(zdx).le.(1.e-7)) goto 120
    75          zx0=zx0+zdx
    76 110   continue
    77 120   continue
    78       zx0=zx0+zdx
    79       if(zanom.lt.0.) zx0=-zx0
    80 
    81 c zteta est la longitude solaire
    82 
    83       zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
    84 
    85       psollong=zteta-timeperi
    86 
    87       IF(psollong.LT.0.) psollong=psollong+2.*pi
    88       IF(psollong.GT.2.*pi) psollong=psollong-2.*pi
    89 c-----------------------------------------------------------------------
    90 c   sorties eventuelles:
    91 c   ---------------------
    92 
    93       IF (lwrite) THEN
    94          PRINT*,'jour de l"annee   :',pday
    95          PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol
    96          PRINT*,'declinaison (en degres) :',pdecli*180./pi
    97       ENDIF
    98 
    99       RETURN
    100       END
     56!  resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
     57!  methode de Newton
     58   
     59    zx0=xref+e_elips*sin(xref)
     60    DO iter=1,10
     61       zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
     62       zx0=zx0+zdx
     63    END DO
     64    zx0=zx0+zdx
     65    if(zanom.lt.0.) zx0=-zx0
     66   
     67    ! zteta est la longitude solaire
     68   
     69    zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
     70   
     71    psollong=zteta-timeperi
     72   
     73    IF(psollong.LT.0.) psollong=psollong+2.*pi
     74    IF(psollong.GT.2.*pi) psollong=psollong-2.*pi
     75    !-----------------------------------------------------------------------
     76    !   sorties eventuelles:
     77    !   ---------------------
     78   
     79    IF (lwrite) THEN
     80       PRINT*,'jour de l"annee   :',pday
     81       PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol
     82       PRINT*,'declinaison (en degres) :',pdecli*180./pi
     83    ENDIF
     84   
     85  END SUBROUTINE solarlong
     86 
     87END MODULE astronomy

  • dynamico_lmdz/simple_physics/phyparam/physics/phys_const.F90

    r4183 r4190  
    11MODULE phys_const
    22 
    3   REAL pi,rad,g,r,cpp,rcp,dtphys,unjours,mugaz,omeg
     3  REAL rad,g,r,cpp,rcp,dtphys,unjours,mugaz,omeg
    44
    55END MODULE phys_const

  • dynamico_lmdz/simple_physics/phyparam/physics/planet.F90

    r4189 r4190  
    33  SAVE
    44
    5   REAL :: aphelie, periheli, year_day, peri_day, obliquit, &
    6        coefvis, coefir
     5  REAL :: coefvis, coefir
    76 
    87END MODULE planet
Note: See TracChangeset for help on using the changeset viewer.