Changeset 3095

Timestamp:

Oct 22, 2023, 3:50:26 PM (14 months ago)

Author:

emillour

Message:

Mars PCM:
Some code tidying:
Made pi in module comcstfi_h a parameter (and not a variable recomputed at
various points by various routines) and added module routine init_comcstfi_h
to cleanly initialize module variables.
Moved iniorbit.F to be a module routine of planete_h since it initializes
(some of ) the module variables it contains.
EM

Location:

trunk/LMDZ.MARS

Files:

• changelog.txt (modified) (1 diff)
• libf/phymars/comcstfi_h.F90 (modified) (1 diff)
• libf/phymars/dyn1d/init_testphys1d_mod.F90 (modified) (1 diff)
• libf/phymars/iniorbit.F (deleted)
• libf/phymars/phys_state_var_init_mod.F90 (modified) (2 diffs)
• libf/phymars/physiq_mod.F (modified) (1 diff)
• libf/phymars/planete_h.F90 (modified) (1 diff)
• libf/phymars/solarlong.F (modified) (1 diff)

trunk/LMDZ.MARS/changelog.txt

@@ -4272 +4272 @@
 clauses in saved module variables in comcstfi_h.F90 and planete_h.F90.
 Prettyfied solarlong.F and made it a module. Likewise for conf_phys.F
+
+Some code tidying:
+Made pi in module comcstfi_h a parameter (and not a variable recomputed at
+various points by various routines) and added module routine init_comcstfi_h
+to cleanly initialize module variables.
+Moved iniorbit.F to be a module routine of planete_h since it initializes
+(some of ) the module variables it contains.

trunk/LMDZ.MARS/libf/phymars/comcstfi_h.F90

@@ -1 @@
-      MODULE comcstfi_h
-      IMPLICIT NONE
+MODULE comcstfi_h
+IMPLICIT NONE
 
-      REAL,SAVE :: pi ! something like 3.14159
-      REAL,SAVE :: rad ! radius of the planet (m)
-      REAL,SAVE :: g ! gravity (m/s2)
-!$OMP THREADPRIVATE(pi,rad,g)
-      REAL,SAVE :: r ! reduced gas constant (r=8.314511/(mugaz/1000.0))
-      REAL,SAVE :: cpp ! Cp of the atmosphere
-      REAL,SAVE :: rcp ! r/cpp
+REAL,PARAMETER :: pi=acos(-1.d0) ! something like 3.14159
+REAL,SAVE :: rad ! radius of the planet (m)
+REAL,SAVE :: g ! gravity (m/s2)
+!$OMP THREADPRIVATE(rad,g)
+REAL,SAVE :: r ! reduced gas constant (r=8.314511/(mugaz/1000.0))
+REAL,SAVE :: cpp ! Cp of the atmosphere
+REAL,SAVE :: rcp ! r/cpp
 !$OMP THREADPRIVATE(r,cpp,rcp)
-      REAL,SAVE :: mugaz ! molar mass of the atmosphere (g/mol)
-      REAL,SAVE :: omeg ! planet rotation rate (rad/s)
+REAL,SAVE :: mugaz ! molar mass of the atmosphere (g/mol)
+REAL,SAVE :: omeg ! planet rotation rate (rad/s)
 !$OMP THREADPRIVATE(mugaz,omeg)
-      
-      END MODULE comcstfi_h
+
+! NB: Ideally all module variables should be "protected"...
+
+CONTAINS
+
+SUBROUTINE init_comcstfi_h(prad,pcpp,pg,pr)
+IMPLICIT NONE
+
+REAL,INTENT(IN) :: prad
+REAL,INTENT(IN) :: pcpp
+REAL,INTENT(IN) :: pg
+REAL,INTENT(IN) :: pr
+
+  rad=prad
+  cpp=pcpp
+  g=pg
+  r=pr
+  rcp=r/cpp
+
+END SUBROUTINE init_comcstfi_h
+
+END MODULE comcstfi_h

trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90

@@ -169 +169 @@
 ! Prescribed constants to be set here
 !------------------------------------------------------
-pi = 2.*asin(1.)
 
 ! Mars planetary constants

trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90

@@ -45 +45 @@
       use conc_mod, only: ini_conc_mod,end_conc_mod
       use turb_mod, only: ini_turb_mod,end_turb_mod
-      use comcstfi_h, only: pi,rad,cpp,g,r,rcp
+      use comcstfi_h, only: init_comcstfi_h
       use tracer_mod, only: ini_tracer_mod,end_tracer_mod
       use time_phylmdz_mod, only: init_time
@@ -109 +109 @@
       
       ! set parameters in comcstfi_h
-      pi=2.*asin(1.) 
-      rad=prad
-      cpp=pcpp
-      g=pg
-      r=pr
-      rcp=r/cpp
+      call init_comcstfi_h(prad,pcpp,pg,pr)
 
       ! Initialize some "temporal and calendar" related variables

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

@@ -80 +80 @@
       use planete_h, only: aphelie, periheli, year_day, peri_day,
      &                     obliquit
+      use planete_h, only: iniorbit
       USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad 
       USE calldrag_noro_mod, ONLY: calldrag_noro

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

@@ -1 @@
-      MODULE planete_h
-      IMPLICIT NONE
+MODULE planete_h
+IMPLICIT NONE
 
-      REAL,SAVE :: aphelie   ! Aphelion, in Mkm
-      REAL,SAVE :: periheli  ! Perihelion, in Mkm
-      REAL,SAVE :: year_day  ! Number of days in the year
+REAL,SAVE :: aphelie   ! Aphelion, in Mkm
+REAL,SAVE :: periheli  ! Perihelion, in Mkm
+REAL,SAVE :: year_day  ! Number of days in the year
 !$OMP THREADPRIVATE(aphelie,periheli,year_day)
-      REAL,SAVE :: peri_day  ! Date of perihelion, in days
-      REAL,SAVE :: obliquit  ! Obliquity of the planet, in degrees
-      REAL,SAVE :: lmixmin
+REAL,SAVE :: peri_day  ! Date of perihelion, in days
+REAL,SAVE :: obliquit  ! Obliquity of the planet, in degrees
+REAL,SAVE :: lmixmin
 !$OMP THREADPRIVATE(peri_day,obliquit,lmixmin)
-      REAL,SAVE :: emin_turb
-      REAL,SAVE :: coefvis
-      REAL,SAVE :: coefir
+REAL,SAVE :: emin_turb
+REAL,SAVE :: coefvis
+REAL,SAVE :: coefir
 !$OMP THREADPRIVATE(emin_turb,coefvis,coefir)
-      REAL,SAVE :: lsperi    ! Solar longitude of the perihelion, angle in rad
-      REAL,SAVE :: e_elips   ! Orbit eccentricity
-      REAL,SAVE :: p_elips   ! Ellipse semi-latus rectum
-      REAL,SAVE :: unitastr  ! Astronomical unit AU, in Mkm
-!$OMP THREADPRIVATE(lsperi,e_elips,p_elips,unitastr)
+REAL,SAVE :: lsperi    ! Solar longitude of the perihelion, angle in rad
+REAL,SAVE :: e_elips   ! Orbit eccentricity
+REAL,SAVE :: p_elips   ! Ellipse semi-latus rectum
+!$OMP THREADPRIVATE(lsperi,e_elips,p_elips)
+REAL,PARAMETER :: unitastr=149.597927 ! Astronomical unit AU, in Mkm
 
-      END MODULE planete_h
+! NB: Ideally all module variables should be "protected"...
+
+CONTAINS
+      
+  SUBROUTINE iniorbit(paphelie,pperiheli,pyear_day,pperi_day,pobliq)
+  use comcstfi_h, only: pi
+  IMPLICIT NONE
+  ! initialize orbital parameters
+  
+  REAL,INTENT(IN) :: paphelie,pperiheli,pyear_day,pperi_day,pobliq
+
+  CHARACTER(LEN=8),PARAMETER :: myname="iniorbit" 
+  REAL :: zxref,zanom,zz,zx0,zdx
+  INTEGER :: iter
+
+  ! copy over input values
+  aphelie =paphelie
+  periheli=pperiheli
+  year_day=pyear_day
+  obliquit=pobliq
+  peri_day=pperi_day
+  
+  write(*,*)myname,': Perihelion in Mkm  ',periheli
+  write(*,*)myname,': Aphelion  in Mkm  ',aphelie
+  write(*,*)myname,': Obliquity in degrees  :',obliquit
+      
+  ! compute orbit geometrical parameters
+  e_elips=(aphelie-periheli)/(periheli+aphelie)
+  p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr
+      
+  write(*,*)myname,': e_elips',e_elips
+  write(*,*)myname,': p_elips',p_elips
+      
+  !  compute mean anomaly zanom
+
+  zz=(year_day-pperi_day)/year_day
+  zanom=2.*pi*(zz-nint(zz))
+  zxref=abs(zanom)
+  write(*,*)myname,': zanom  ',zanom
+
+  ! solve equation  zx0 - e * sin (zx0) = zxref for eccentric anomaly zx0
+  ! using Newton method
+
+  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)) exit
+    zx0=zx0+zdx
+  ENDDO
+
+  zx0=zx0+zdx
+  if(zanom.lt.0.) zx0=-zx0
+  write(*,*)myname,': zx0   ',zx0
+
+  lsperi=-2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
+  lsperi = modulo(lsperi,2.*pi)
+  write(*,*)myname,': Perihelion solar long. Ls (deg)=', &
+                          lsperi*180./pi
+
+  END SUBROUTINE iniorbit
+
+END MODULE planete_h

trunk/LMDZ.MARS/libf/phymars/solarlong.F

@@ -35 +35 @@
 
       zz=(pday-peri_day)/year_day
-      pi=2.*asin(1.)
       zanom=2.*pi*(zz-nint(zz))
       xref=abs(zanom)