Changeset 3094 for trunk/LMDZ.MARS/libf/phymars

Timestamp:

Oct 22, 2023, 1:42:30 PM (22 months ago)

Author:

emillour

Message:

Mars PCM:
Fix issues with dates when in parallel with OpenMP (missing copyin
statement when opening parallel section in iniphysiq). Added some threadprivate
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
EM

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• comcstfi_h.F90 (modified) (1 diff)
• conf_phys.F (modified) (4 diffs)
• dyn1d/init_testphys1d_mod.F90 (modified) (1 diff)
• physiq_mod.F (modified) (2 diffs)
• planete_h.F90 (modified) (1 diff)
• solarlong.F (modified) (2 diffs)

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

@@ -2 +2 @@
       IMPLICIT NONE
 
-      REAL pi ! something like 3.14159
-      REAL rad ! radius of the planet (m)
-      REAL g ! gravity (m/s2)
-      REAL r ! reduced gas constant (r=8.314511/(mugaz/1000.0))
-      REAL cpp ! Cp of the atmosphere
-      REAL rcp ! r/cpp
-      REAL mugaz ! molar mass of the atmosphere (g/mol)
-      REAL omeg ! planet rotation rate (rad/s)
-
+      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
+!$OMP THREADPRIVATE(r,cpp,rcp)
+      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

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

@@ +1 @@
+      MODULE conf_phys_mod
+      
+      IMPLICIT NONE
+      
+      CONTAINS
+
       SUBROUTINE conf_phys(ngrid,nlayer,nq)
  
@@ -6 +12 @@
 !   -------
 !
-!   Initialisation for the physical parametrisations of the LMD 
-!   martian atmospheric general circulation modele.
-!
-!   author: Frederic Hourdin 15 / 10 /93
-!   -------
-!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
-!             Ehouarn Millour (oct. 2008) tracers are now identified
-!              by their names and may not be contiguously
-!              stored in the q(:,:,:,:) array
-!             E.M. (june 2009) use getin routine to load parameters
-!             adapted to the mesoscale use - Aymeric Spiga - 01/2007-07/2011
-!             separated inifis into conf_phys and phys_state_var_init (A. Spiga)
-!
-!
-!   arguments:
-!   ----------
-!
-!   input:
-!   ------
-!
-!    nq                    Number of tracers
-!
-!=======================================================================
-!
+!   Initialisation for the physical parametrisations 
+!   flags (i.e. run-time options) of the Mars PCM 
 !-----------------------------------------------------------------------
-!   declarations:
-!   -------------
+
       USE ioipsl_getin_p_mod, ONLY : getin_p
       use tracer_mod, only : nuice_sed, ccn_factor, nuiceco2_sed,
@@ -61 +44 @@
       include "callkeys.h"
 
-      INTEGER,INTENT(IN) :: ngrid,nlayer,nq
+      INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
+      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
+      INTEGER,INTENT(IN) :: nq ! number of tracers
 
       INTEGER ierr,j
@@ -1226 +1211 @@
       call bcast(semi)
 
-      END
+      END SUBROUTINE conf_phys
+
+      END MODULE conf_phys_mod

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

@@ -39 +39 @@
 use turb_mod,                 only: q2
 use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd
+use conf_phys_mod, only: conf_phys
 ! Mostly for XIOS outputs:
 use mod_const_mpi,            only: COMM_LMDZ

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

@@ -61 +61 @@
       use comsaison_h, only: dist_sol, declin, zls, 
      &                       mu0, fract, local_time
+      use solarlong_mod, only: solarlong
       use nirdata_mod, only: NIR_leedat
       use nirco2abs_mod, only: nirco2abs
@@ -2605 +2606 @@
            ELSE ! IF LMDZ
 
-          if (ecritstart.GT.0) then !IF MULTIPLE RESTARTS nothing change
-          ztime_fin = pday - day_ini + ptime  
-     .               + ptimestep/(float(iphysiq)*daysec)
-          else !IF ONE RESTART final time in top of day_end
-          ztime_fin = pday - day_ini-(day_end-day_ini)
-     .               + ptime  + ptimestep/(float(iphysiq)*daysec)
-          endif
+            if (ecritstart.GT.0) then !IF MULTIPLE RESTARTS nothing change
+              ztime_fin = pday - day_ini + ptime  
+     &                    + ptimestep/(float(iphysiq)*daysec)
+            else !IF ONE RESTART final time in top of day_end
+              ztime_fin = pday - day_ini-(day_end-day_ini)
+     &                    + ptime  + ptimestep/(float(iphysiq)*daysec)
+            endif
 
            ENDIF ! of IF (grid_type==unstructured)

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

@@ -2 +2 @@
       IMPLICIT NONE
 
-      REAL aphelie   ! Aphelion, in Mkm
-      REAL periheli  ! Perihelion, in Mkm
-      REAL year_day  ! Number of days in the year
-      REAL peri_day  ! Date of perihelion, in days
-      REAL obliquit  ! Obliquity of the planet, in degrees
-      REAL lmixmin
-      REAL emin_turb
-      REAL coefvis
-      REAL coefir
-      REAL lsperi    ! Solar longitude of the perihelion, angle in rad
-      REAL e_elips   ! Orbit eccentricity
-      REAL p_elips   ! Ellipse semi-latus rectum
-      REAL unitastr  ! Astronomical unit AU, in Mkm
+      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
+!$OMP THREADPRIVATE(peri_day,obliquit,lmixmin)
+      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)
 
       END MODULE planete_h

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

@@ +1 @@
+      MODULE solarlong_mod
+
+      IMPLICIT NONE
+
+      CONTAINS
+      
       SUBROUTINE solarlong(pday,psollong)
-      use planete_h
-      USE comcstfi_h
+      use planete_h, only: lsperi, peri_day, year_day, e_elips
+      use comcstfi_h, only: pi
       IMPLICIT NONE
 
 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   Compute solar longitude psollong (in radians) for a given Mars date
+c   pday (where pday==0 at northern Spring Equinox and pday in sols 
+c   and fractions thereof)
 c=======================================================================
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
 
 c arguments:
 c ----------
 
-      REAL pday,pdist_sol,pdecli,psollong
-      LOGICAL lwrite
+      REAL,INTENT(IN) :: pday
+      REAL,INTENT(OUT) :: psollong
+      
+c Local variables:
+c ----------------
 
-c Local:
-c ------
-
-      REAL zanom,xref,zx0,zdx,zteta,zz
-      INTEGER iter
+      REAL :: xref ! mean anomaly
+      REAL :: zx0 ! eccentric anomaly
+      REAL :: zteta ! true anomaly
+      REAL :: zanom,zdx,zz
+      INTEGER :: iter
 
 
-c-----------------------------------------------------------------------
-c calcul de l'angle polaire et de la distance au soleil :
-c -------------------------------------------------------
-
-c  calcul de l'zanomalie moyenne
+c  compute mean anomaly
 
       zz=(pday-peri_day)/year_day
@@ -67 +39 @@
       xref=abs(zanom)
 
-c  resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
-c  methode de Newton
+c  solve equation  zx0 - e * sin (zx0) = xref
+c  using Newton method
 
       zx0=xref+e_elips*sin(xref)
-      DO 110 iter=1,10
+      DO iter=1,10 ! 10 is overkill, typically converges in 1-3 iterations
          zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
-         if(abs(zdx).le.(1.e-7)) goto 120
+         if(abs(zdx).le.(1.e-7)) exit
          zx0=zx0+zdx
-110   continue
-120   continue
+      ENDDO
+
       zx0=zx0+zdx
       if(zanom.lt.0.) zx0=-zx0
 
-c zteta est la longitude solaire
+c compute true anomaly zteta, now that eccentric anomaly zx0 is known
 
       zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
 
+c compute solar longitude from zteta
       psollong=zteta+lsperi
 
+c handle limit cases where psollong lands outside [0:2*pi]
       IF(psollong.LT.0.) psollong=psollong+2.*pi
       IF(psollong.GT.2.*pi) psollong=psollong-2.*pi
-c-----------------------------------------------------------------------
-c   sorties eventuelles:
-c   ---------------------
 
-c     IF (lwrite) THEN
-c        PRINT*,'jour de l"annee   :',pday
-c        PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol
-c        PRINT*,'declinaison (en degres) :',pdecli*180./pi
-c     ENDIF
+      END SUBROUTINE solarlong
 
-      RETURN
-      END
+      END MODULE solarlong_mod