Changeset 3974 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Nov 21, 2025, 4:28:06 PM (2 months ago)
Author:
jbclement
Message:

Mars PCM:
Simplification of the orbital parameters initialization (resolve ticket #109): removing redundant actions in subroutine 'iniorbit' in regard of what is done in "tabfi.F" and moving the subroutine 'lsp2solp' from "tabfi.F" to "planete_h.F90" + some cleanings.
JBC

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90

    r3964 r3974  
    3535use comcstfi_h,          only: pi
    3636use geometry_mod,        only: latitude
    37 use soil_settings_mod, only: soil_settings
    38 use tabfi_mod, only: tabfi
    39 use callkeys_mod, only: startphy_file, rdstorm, hdo
    40 use callkeys_mod, only: CLFvarying, CLFfixval
     37use soil_settings_mod,   only: soil_settings
     38use tabfi_mod,           only: tabfi
     39use callkeys_mod,        only: startphy_file, rdstorm, hdo, CLFvarying, CLFfixval
    4140
    4241implicit none

  • trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

    r3918 r3974  
    8585      use turb_mod, only: q2, wstar, ustar, sensibFlux,
    8686     &                    zmax_th, hfmax_th, turb_resolved
    87       use planete_h, only: aphelie, periheli, year_day, peri_day,
    88      &                     obliquit
    89       use planete_h, only: iniorbit
     87      use planete_h, only: obliquit
    9088      use orbite_mod, only: orbite
    9189      USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad
    … …  
    718716c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    719717         CALL surfini(ngrid,nslope,qsurf)
    720          CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
     718
    721719c        initialize soil
    722720c        ~~~~~~~~~~~~~~~

  • trunk/LMDZ.MARS/libf/phymars/planete_h.F90

    r3095 r3974  
    11MODULE planete_h
    2 IMPLICIT NONE
     2
     3implicit none
    34
    45REAL,SAVE :: aphelie   ! Aphelion, in Mkm
    … …  
    2324
    2425CONTAINS
     26!-----------------------------------------------------------------------
     27! Gives sol at perihelion for ls at perihelion (for precession cycles)
     28  SUBROUTINE lsp2solp(lsp,solp)
     29
     30  use comcstfi_h, only: pi
     31
     32  implicit none
     33
     34  ! Arguments:
     35  real, intent(in) :: lsp ! ls at perihelion
     36  real, intent(out) :: solp ! sol at perihelion
     37
     38  ! Locals:
     39  real :: zx0 ! eccentric anomaly at Ls=0
     40  real, parameter :: degrad = 180.d0/pi
     41
     42  ! Compute orbit geometrical parameters
     43  e_elips = (aphelie - periheli)/(aphelie + periheli)
     44  zx0 = -2.0*datan(dtan(0.5*lsp/degrad)*dsqrt((1. - e_elips)/(1. + e_elips)))
     45  if (zx0 <= 0.) zx0 = zx0 + 2.*pi
     46
     47  ! Compute sol at perihelion
     48  solp = year_day*(1. - (zx0 - e_elips*dsin(zx0))/(2.*pi))
     49
     50  END SUBROUTINE lsp2solp
     51
     52!-----------------------------------------------------------------------
     53! Initialize orbital parameters
     54  SUBROUTINE iniorbit()
     55
     56  use comcstfi_h, only: pi
     57
     58  implicit none
     59
     60  ! Locals:
     61  character(8), parameter :: myname = "iniorbit"
     62  real    :: zxref, zanom, zz, zx0, zdx
     63  integer :: iter
     64
     65  ! Info about orbital values
     66  write(*,*)myname,': Aphelion in Mkm                 =',aphelie
     67  write(*,*)myname,': Perihelion in Mkm               =',periheli
     68  write(*,*)myname,': Number of days in the year      =',year_day
     69  write(*,*)myname,': Date of perihelion in days      =',peri_day
     70  write(*,*)myname,': Obliquity in degrees            =',obliquit
     71
     72  ! Compute orbit geometrical parameters
     73  e_elips = (aphelie - periheli)/(aphelie + periheli)
     74  p_elips = 0.5*(periheli + aphelie)*(1. - e_elips*e_elips)/unitastr
     75  write(*,*)myname,': Orbit eccentricity              =',e_elips
     76  write(*,*)myname,': Ellipse semi-latus rectum       =',p_elips
    2577     
    26   SUBROUTINE iniorbit(paphelie,pperiheli,pyear_day,pperi_day,pobliq)
    27   use comcstfi_h, only: pi
    28   IMPLICIT NONE
    29   ! initialize orbital parameters
    30  
    31   REAL,INTENT(IN) :: paphelie,pperiheli,pyear_day,pperi_day,pobliq
     78  ! Compute mean anomaly zanom
     79  zz = (year_day - peri_day)/year_day
     80  zanom = 2.*pi*(zz - nint(zz))
     81  zxref = abs(zanom)
     82  write(*,*)myname,': zanom                           =',zanom
    3283
    33   CHARACTER(LEN=8),PARAMETER :: myname="iniorbit"
    34   REAL :: zxref,zanom,zz,zx0,zdx
    35   INTEGER :: iter
    36 
    37   ! copy over input values
    38   aphelie =paphelie
    39   periheli=pperiheli
    40   year_day=pyear_day
    41   obliquit=pobliq
    42   peri_day=pperi_day
    43  
    44   write(*,*)myname,': Perihelion in Mkm  ',periheli
    45   write(*,*)myname,': Aphelion  in Mkm  ',aphelie
    46   write(*,*)myname,': Obliquity in degrees  :',obliquit
    47      
    48   ! compute orbit geometrical parameters
    49   e_elips=(aphelie-periheli)/(periheli+aphelie)
    50   p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr
    51      
    52   write(*,*)myname,': e_elips',e_elips
    53   write(*,*)myname,': p_elips',p_elips
    54      
    55   !  compute mean anomaly zanom
    56 
    57   zz=(year_day-pperi_day)/year_day
    58   zanom=2.*pi*(zz-nint(zz))
    59   zxref=abs(zanom)
    60   write(*,*)myname,': zanom  ',zanom
    61 
    62   ! solve equation  zx0 - e * sin (zx0) = zxref for eccentric anomaly zx0
     84  ! Solve equation  zx0 - e * sin (zx0) = zxref for eccentric anomaly zx0
    6385  ! using Newton method
    6486
    65   zx0=zxref+e_elips*sin(zxref)
    66   DO iter=1,100
    67     zdx=-(zx0-e_elips*sin(zx0)-zxref)/(1.-e_elips*cos(zx0))
    68     if(abs(zdx).le.(1.e-12)) exit
    69     zx0=zx0+zdx
     87  zx0 = zxref + e_elips*sin(zxref)
     88  DO iter = 1,100
     89    zdx = -(zx0 - e_elips*sin(zx0) - zxref)/(1. - e_elips*cos(zx0))
     90    if (abs(zdx).le.(1.e-12)) exit
     91    zx0 = zx0 + zdx
    7092  ENDDO
    7193
    72   zx0=zx0+zdx
    73   if(zanom.lt.0.) zx0=-zx0
    74   write(*,*)myname,': zx0   ',zx0
     94  zx0 = zx0 + zdx
     95  if (zanom < 0.) zx0 = -zx0
     96  write(*,*)myname,': Eccentric anomaly zx0           =',zx0
    7597
    76   lsperi=-2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
     98  lsperi = -2.*atan(sqrt((1. + e_elips)/(1. - e_elips))*tan(zx0/2.))
    7799  lsperi = modulo(lsperi,2.*pi)
    78   write(*,*)myname,': Perihelion solar long. Ls (deg)=', &
    79                           lsperi*180./pi
     100  write(*,*)myname,': Perihelion solar long. Ls (deg) =', lsperi*180./pi
    80101
    81102  END SUBROUTINE iniorbit

  • trunk/LMDZ.MARS/libf/phymars/tabfi.F

    r3902 r3974  
    6161      use time_phylmdz_mod, only: daysec, dtphys
    6262      use planete_h, only: aphelie, emin_turb, lmixmin, obliquit,
    63      &                     peri_day, periheli, year_day
     63     &                     peri_day, periheli, year_day,
     64     &                     lsp2solp, iniorbit
     65
    6466      implicit none
    6567 
    … …  
    460462          write(*,*) 'perihelion =',periheli
    461463          write(*,*) 'and',year_day,'sols/year'
    462           call lsp2solp(peri_ls,peri_day,aphelie,periheli,year_day)
     464          call lsp2solp(peri_ls,peri_day)
    463465          write(*,*) 'peri_day (new value):',peri_day
    464466
    … …  
    543545
    544546c-----------------------------------------------------------------------
     547c       Initialization of orbital parameters
     548c-----------------------------------------------------------------------
     549      call iniorbit()
     550
     551c-----------------------------------------------------------------------
    545552c       Case when using a start file from before March 1996 (without iceradius...
    546553c-----------------------------------------------------------------------
    … …  
    555562       end if
    556563
    557 c-----------------------------------------------------------------------
    558564      END SUBROUTINE tabfi
    559565
    560 
    561 
    562      
    563 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    564 ! gives sol at perihelion for ls at perihelion (for precession cycles)
    565       subroutine lsp2solp(lsp,solp,aphelie,periheli,year_day)
    566 
    567       implicit none
    568 !  Arguments:
    569       real lsp     ! Input: ls at perihelion
    570       real solp    ! Output: sol at perihelion
    571       real aphelie,periheli,year_day ! Input: parameters
    572  
    573 !  Local:
    574       double precision zx0 ! eccentric anomaly at Ls=0
    575       double precision e_elips
    576       double precision pi,degrad
    577      
    578       parameter (pi=4.d0*atan(1.d0))
    579       parameter (degrad=180.d0/pi)
    580 
    581       e_elips=(aphelie-periheli)/(aphelie+periheli)     
    582       zx0 = -2.0*datan(dtan(0.5*lsp/degrad)
    583      .          *dsqrt((1.-e_elips)/(1.+e_elips)))
    584       if (zx0 .le. 0.) zx0 = zx0 + 2.*pi
    585      
    586       solp  = year_day*(1.-(zx0-e_elips*dsin(zx0))/(2.*pi))
    587 
    588 
    589       end subroutine lsp2solp
    590 
    591 
    592566      END MODULE tabfi_mod
    593 
Note: See TracChangeset for help on using the changeset viewer.