Changeset 3663 for trunk/LMDZ.GENERIC/libf/phystd

Timestamp:

Feb 27, 2025, 7:03:20 PM (4 months ago)

Author:

gmilcareck

Message:

Generic PCM:
Cleaning up the code. The gas constant and Avogadro number values
was not the same between the files in the model.
We choose the value recommended by the 2019 revision of the SI.
R=8.314463 J.K-1.mol-1 and NA=6.022141e23 mol-1.
GM

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• aeropacity.F90 (modified) (7 diffs)
• blackl.F (modified) (2 diffs)
• calc_rayleigh.F90 (modified) (2 diffs)
• comcstfi_mod.F90 (modified) (1 diff)
• dyn1d/phasechange_kcm.F90 (modified) (1 diff)
• dyn1d/rcm1d.F (modified) (1 diff)
• inifis_mod.F90 (modified) (1 diff)
• interpolateH2O_self_foreign.F90 (modified) (1 diff)
• newsedim.F (modified) (1 diff)
• optci.F90 (modified) (1 diff)
• optcv.F90 (modified) (1 diff)
• params_h.F90 (modified) (1 diff)
• stokes.F90 (modified) (1 diff)
• tabfi_mod.F90 (modified) (5 diffs)
• thermcell_env.F90 (modified) (1 diff)
• thermcell_plume.F90 (modified) (1 diff)
• watercommon_h.F90 (modified) (1 diff)

trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90

@@ -724 +724 @@
 ! p_bot <=> zb ; p_top <=> zb+zc ; h_bot <=> Hlo ; h_top <=> Hup
 ! p<p_top: N=No*(p/p_top)**(h_lay/h_top)      h_lay=RT/g  (in m)
-! p>p_bot: N=No*(p_bot/p)**(h_lay/h_bot)      R=8.31/mu (mu in kg/mol)
+! p>p_bot: N=No*(p_bot/p)**(h_lay/h_bot)      R=8.314463/mu (mu in kg/mol)
 ! N is in m-3
 !
@@ -745 +745 @@
            DO l=1,nlayer-1
 
-             h_lay=8.31*pt(ig,l)/(g*0.044)
+             h_lay=8.314463*pt(ig,l)/(g*0.044)
 
              !! 1. below 2e5 Pa: no aerosols
@@ -791 +791 @@
            DO l=1,nlayer-1
 
-             h_lay=8.31*pt(ig,l)/(g*0.044)
+             h_lay=8.314463*pt(ig,l)/(g*0.044)
 
              !! 1. below 2e5 Pa: no aerosols
@@ -837 +837 @@
            DO l=1,nlayer-1
 
-             h_lay=8.31*pt(ig,l)/(g*0.044)
+             h_lay=8.314463*pt(ig,l)/(g*0.044)
 
              !! 1. below 2e5 Pa: no aerosols
@@ -883 +883 @@
            DO l=1,nlayer-1
  
-              h_lay=8.31*pt(ig,l)/(g*0.044)
+              h_lay=8.314463*pt(ig,l)/(g*0.044)
 
              !! 1. below 2e5 Pa: no aerosols
@@ -929 +929 @@
            DO l=1,nlayer-1
 
-             h_lay=8.31*pt(ig,l)/(g*0.044)
+             h_lay=8.314463*pt(ig,l)/(g*0.044)
 
              !! 1. below 7e4 Pa: no aerosols
@@ -975 +975 @@
 !     ig=10
 !      do l=1,nlayer
-!          print*,8.31*pt(ig,l)/(g*0.044),pplay(ig,l),aerosol(ig,l,1),aerosol(ig,l,2),aerosol(ig,l,3),aerosol(ig,l,4)
+!          print*,8.314463*pt(ig,l)/(g*0.044),pplay(ig,l),aerosol(ig,l,1),aerosol(ig,l,2),aerosol(ig,l,3),aerosol(ig,l,4)
 !         print*,l,pplay(ig,l),aerosol(ig,l,5)
 !      enddo

trunk/LMDZ.GENERIC/libf/phystd/blackl.F

@@ -4 +4 @@
 
       ! physical constants
-      sigma=5.67032D-8
+      sigma=5.670374D-8
       pi=datan(1.d0)*4.d0
-      c0=2.9979d+08
-      h=6.6262d-34
-      cbol=1.3806d-23
+      c0=2.997925d+08
+      h=6.62607d-34
+      cbol=1.380649d-23
       rind=1.d0
       c=c0/rind
@@ -26 +26 @@
 
       ! physical constants
-      sigma=5.67032D-8
+      sigma=5.670374D-8
       pi=datan(1.d0)*4.d0
-      c0=2.9979d+08
-      h=6.6262d-34
-      cbol=1.3806d-23
+      c0=2.997925d+08
+      h=6.62607d-34
+      cbol=1.380649d-23
       rind=1.d0
       c=c0/rind

trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90

@@ -188 +188 @@
                      P0(:) = pmid(:)*scalep
                      lr = 0.589
-                     rho(:) = mass_frac(igas,:)*muvar(:)/massmol(igas)*P0(:)/(8.314462*T0(:)/muvar(:)/1000.)
+                     rho(:) = mass_frac(igas,:)*muvar(:)/massmol(igas)*P0(:)/(8.314463*T0(:)/muvar(:)/1000.)
                      rhos(:) = rho(:)/rhor
                      ts(:) = T0(:)/Tr
@@ -258 +258 @@
                  endif
                  
-                 tauconsti(igas,:) = 24.*pi**3 *6.022140857E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380648813E-23*T0(:)))**2)
+                 tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
                  ! N=P/(kB*T)
                  ! pmid*scalep -> mbar to Pa

trunk/LMDZ.GENERIC/libf/phystd/comcstfi_mod.F90

@@ -5 +5 @@
       REAL,SAVE :: rad ! radius of the planet (m)
       REAL,SAVE :: g ! gravity (m/s2)
-      REAL,SAVE :: r ! reduced gas constant (r=8.314511/(mugaz/1000.0))
+      REAL,SAVE :: r ! reduced gas constant (r=8.314463/(mugaz/1000.0))
       REAL,SAVE :: cpp ! Cp of the atmosphere
       REAL,SAVE :: rcp ! r/cpp
       REAL,SAVE :: mugaz ! molar mass of the atmosphere (g/mol)
       REAL,SAVE :: omeg ! planet rotation rate (rad/s)
-      REAL,SAVE :: avocado ! something like 6.022e23
+      REAL,SAVE :: avocado !  6.022141e23
 !$OMP THREADPRIVATE(pi,rad,g,r,cpp,rcp,mugaz,omeg,avocado)
 

trunk/LMDZ.GENERIC/libf/phystd/dyn1d/phasechange_kcm.F90

@@ -1254 +1254 @@
 !
 !    The specific gas constant R is related to the universal gas
-!    constant R-bar = 8.31441 J/(mol degrees Kelvin) by the molar mass 
+!    constant R-bar = 8.314463 J/(mol degrees Kelvin) by the molar mass 
 !    M = 18.0152 g/mol:
 !

trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F

@@ -659 +659 @@
       PRINT *,"--> mugaz = ",mugaz
       PRINT *,"--> cpp = ",cpp
-      r = 8.314511E+0 * 1000.E+0 / mugaz
+      r = 8.314463E+0 * 1000.E+0 / mugaz
       rcp = r / cpp
 

trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90

@@ -92 +92 @@
   r=pr
   rcp=r/cpp
-  mugaz=8.314*1000./pr ! dummy init
+  mugaz=8.314463*1000./pr ! dummy init
   pi=2.*asin(1.)
-  avocado = 6.02214179e23   ! added by RW
+  avocado = 6.022141e23   ! added by RW
 
   ! Initialize some "temporal and calendar" related variables

trunk/LMDZ.GENERIC/libf/phystd/interpolateH2O_self_foreign.F90

@@ -34 +34 @@
 
       double precision kB
-      parameter(kB=1.3806488e-23)
+      parameter(kB=1.380649e-23)
 
       double precision amagatS, amagatF, abcoefS, abcoefF, Nmolec

trunk/LMDZ.GENERIC/libf/phystd/newsedim.F

@@ -146 +146 @@
 !       - Constant  to compute gas mean free path
 !        l= (T/P)*a, avec a = (  0.707*8.31/(4*pi*molrad**2 * avogadro))
-         a = 0.707*8.31/(4*3.1416* molrad**2  * 6.023e23)
+         a = 0.707*8.314463/(4*3.1416* molrad**2  * 6.022141e23)
 
 c       - Correction to account for non-spherical shape (Murphy et al.  1990)

trunk/LMDZ.GENERIC/libf/phystd/optci.F90

@@ -154 +154 @@
      ! if we have continuum opacities, we need dz
      if(kastprof)then
-        dz(k) = dpr(k)*(1000.0d0*8.3145d0/muvar(k))*TMID(K)/(g*PMID(K))
+        dz(k) = dpr(k)*(1000.0d0*8.314463d0/muvar(k))*TMID(K)/(g*PMID(K))
         U(k)  = Cmk*DPR(k)*mugaz/muvar(k) 
      else

trunk/LMDZ.GENERIC/libf/phystd/optcv.F90

@@ -140 +140 @@
      ! if we have continuum opacities, we need dz
      if(kastprof)then
-        dz(k) = dpr(k)*(1000.0d0*8.3145d0/muvar(k))*TMID(K)/(g*PMID(K))
+        dz(k) = dpr(k)*(1000.0d0*8.314463d0/muvar(k))*TMID(K)/(g*PMID(K))
         U(k)  = Cmk*DPR(k)*mugaz/muvar(k) 
      else

trunk/LMDZ.GENERIC/libf/phystd/params_h.F90

@@ -12 +12 @@
       save ideal_v
 
-      parameter (rc  = 8.314462d0) ! ideal gas constant [J mol^-1 K^-1]
+      parameter (rc  = 8.314463d0) ! ideal gas constant [J mol^-1 K^-1]
 
       end module params_h

trunk/LMDZ.GENERIC/libf/phystd/stokes.F90

@@ -49 +49 @@
 
 
-         a = 0.707*8.31/(4*pi*molrad**2 * avocado)
+         a = 0.707*8.314463/(4.*pi*molrad**2 * avocado)
          b = (2./9.) * rho_aer * g / visc 
         !print*,'molrad=',molrad

trunk/LMDZ.GENERIC/libf/phystd/tabfi_mod.F90

@@ -105 +105 @@
           call abort_physic(modname,"Missing value for omega in def files!",1)
         endif
-        mugaz=(8.3144621/r)*1.E3
+        mugaz=(8.314463/r)*1.E3
         ! daysec already set by inifis
         ! dtphys alread set by inifis
@@ -180 +180 @@
       mugaz = tab_cntrl(tab0+8)
       rcp = tab_cntrl(tab0+9)
-      cpp=(8.314511/(mugaz/1000.0))/rcp
+      cpp=(8.314463/(mugaz/1000.0))/rcp
       daysec = tab_cntrl(tab0+10)
       dtphys = tab_cntrl(tab0+11)
@@ -455 +455 @@
           write(*,*) 
           write(*,*) ' mugaz (new value):',mugaz
-          r=8.314511/(mugaz/1000.0)
+          r=8.314463/(mugaz/1000.0)
           write(*,*) ' R (new value):',r
 
@@ -465 +465 @@
           write(*,*) 
           write(*,*) ' rcp (new value):',rcp
-          r=8.314511/(mugaz/1000.0)
+          r=8.314463/(mugaz/1000.0)
           cpp=r/rcp
           write(*,*) ' cpp (new value):',cpp
@@ -477 +477 @@
           write(*,*) ' cpp (new value):',cpp
           write(*,*) ' mugaz (new value):',mugaz
-          r=8.314511/(mugaz/1000.0)
+          r=8.314463/(mugaz/1000.0)
           rcp=r/cpp
           write(*,*) ' rcp (new value):',rcp

trunk/LMDZ.GENERIC/libf/phystd/thermcell_env.F90

@@ -132 +132 @@
                ENDDO
 
-               RV_generic = (8.314511*1000.)/(epsi_generic*mugaz)
+               RV_generic = (8.314463*1000.)/(epsi_generic*mugaz)
                RETV_generic = RV_generic/r-1.
                RLvCp_generic = RLVTT_generic/cpp

trunk/LMDZ.GENERIC/libf/phystd/thermcell_plume.F90

@@ -144 +144 @@
          RLvCp_comp = RLvCp 
       ELSEIF (generic_condensation .AND. .NOT. water ) THEN
-         RV_generic = (8.314511*1000.)/(epsi_generic*mugaz)
+         RV_generic = (8.314463*1000.)/(epsi_generic*mugaz)
          RETV_comp = RV_generic/r-1.
          RLvCp_comp = RLVTT_generic/cpp

trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90

@@ -49 +49 @@
 
       !RW = 1000.*R/mH2O
-      RW = 1000.*8.314/mH2O ! caution! R is R/mugaz already!
+      RW = 1000.*8.314463/mH2O ! caution! R is R/mugaz already!
 
       RCPV   = 1.88e3 ! specific heat capacity of water vapor at 350K