Changeset 3185 for trunk/LMDZ.MARS

Timestamp:

Jan 26, 2024, 10:10:05 AM (10 months ago)

Author:

csegonne

Message:

MARS PCM:

• conc_mod.F : Added new subroutines init_r_cp_mu and update_r_cp_mu_ak to replace initialization of rnew, cpnew, mmean and akknew as constants and their update if callthermos or photochem (same update as in concentrations.F) that were done in physic_mod.F.
• rnew, cpnew, mmean and akknew are now protected meaning they cannot be modified anywhere in the model other than in conc_mod.F .
• concentrations.F has been removed.
• Knowing where rnew, cpnew, mmean come from is now more explicit. These modifications respond to the need expressed in ticket #90.

Location:

trunk/LMDZ.MARS

Files:

• changelog.txt (modified) (4 diffs)
• libf/phymars/conc_mod.F90 (modified) (2 diffs)
• libf/phymars/physiq_mod.F (modified) (4 diffs)

trunk/LMDZ.MARS/changelog.txt

@@ -2466 +2466 @@
 microphysics (just like water clouds). The amplitude is spantCO2, also read from callphys.def
 - The previous logical microphysco2 has been removed.
-- Cloud opacity at 1µm is now computed in the co2cloud.F routine
+- Cloud opacity at 1µm is now computed in the co2cloud.F routine
 - Most of the co2 ice clouds scheme writediagfi are now in co2clouds.F
 - Some cleaning of the CO2 ice clouds routine has been done. Not perfect yet!
@@ -2482 +2482 @@
   if a gravity wave propagates to the mesosphere or not (wheter the wave saturates on its way or not). 
   The sub-grid temperature distribution is cancelled if the saturation index value of the column (between 12 and 80 km) is > 0.1, and is applied (+/- 3K) otherwise. if the keyword satindexco2=.true. in callphys.def (only applies if CLFvaryingCO2=.true. anyway). If set to .false., deactivate this filter and the sub grid T distribution applies everywhere and everytime. See comments in co2clouds.F
-- All that you need to launch a GCM run with co2 clouds is described in co2clouds.F comments. Ehouarn Millour and Anni Määttänen have the files. 
+- All that you need to launch a GCM run with co2 clouds is described in co2clouds.F comments. Ehouarn Millour and Anni Määttänen have the files. 
  
 == 15/11/2017 == EM
@@ -4291 +4291 @@
 d_coef=X for the diffusion coeficent (4e-4 is the defualt)
 <h2o_ice_depth=X> if its zero or lower there is no subsurface ice (default -1).
-outputs: zdq_ssi, momentarily zdqsdif after the interaction with the SSIzdq_subtimestep,  momentarily flux after the interaction with the SSIzdq_ssi_frost momentarily zdqsdif after the interaction of the frost with the SSI
+outputs: zdq_ssi, momentarily zdqsdif after the interaction with the SSIzdq_subtimestep,  momentarily flux after the interaction with the SSIzdq_ssi_frost momentarily zdqsdif after the interaction of the frost with the SSI
 q1, atmospheric first layer water vapor quatity
 qeq, SSI water vapor quantity
@@ -4445 +4445 @@
 == 25/01/2024 == JBC
 Some changes to prepare the introduction of slopes in 1D: in particular, the subroutine "getslopes.F90" and "param_slope.F90" are now inside the module "slope_mod.F90" + Few small cleanings throughout the code.
+
+== 26/01/2024 == CS
+- conc_mod.F : Added new subroutines init_r_cp_mu and update_r_cp_mu_ak to replace initialization of rnew, cpnew, mmean and akknew as constants and their update if callthermos or photochem (same update as in concentrations.F) that were done in physic_mod.F. 
+- rnew, cpnew, mmean and akknew are now protected meaning they cannot be modified anywhere in the model other than in conc_mod.F .
+- concentrations.F has been removed.
+- Knowing where rnew, cpnew, mmean come from is now more explicit. These modifications respond to the need expressed in ticket #90.

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

@@ -3 +3 @@
 implicit none
 
-  real,save,allocatable :: mmean(:,:)  ! mean molecular mass of the atmosphere
-  real,save,allocatable :: Akknew(:,:) ! thermal conductivity cofficient
-  real,save,allocatable :: cpnew(:,:)  ! specicic heat
-  real,save,allocatable :: rnew(:,:)   ! specific gas constant
+ real,save,allocatable,protected :: mmean(:,:)  ! mean molecular mass of the atmosphere
+ real,save,allocatable,protected :: Akknew(:,:) ! thermal conductivity coefficient
+ real,save,allocatable,protected :: cpnew(:,:)  ! specific heat
+ real,save,allocatable,protected :: rnew(:,:)   ! specific gas constant
 
 !$OMP THREADPRIVATE(mmean,Akknew,cpnew,rnew)
@@ -12 +12 @@
 contains
 
-  subroutine ini_conc_mod(ngrid,nlayer)
+ subroutine ini_conc_mod(ngrid,nlayer)
   
-  implicit none
-  integer,intent(in) :: ngrid ! number of atmospheric columns
-  integer,intent(in) :: nlayer ! number of atmospheric levels
+ implicit none
+ integer,intent(in) :: ngrid ! number of atmospheric columns
+ integer,intent(in) :: nlayer ! number of atmospheric levels
   
-    allocate(mmean(ngrid,nlayer))
-    allocate(Akknew(ngrid,nlayer))
-    allocate(cpnew(ngrid,nlayer))
-    allocate(rnew(ngrid,nlayer))
+  allocate(mmean(ngrid,nlayer))
+  allocate(Akknew(ngrid,nlayer))
+  allocate(cpnew(ngrid,nlayer))
+  allocate(rnew(ngrid,nlayer))
     
-  end subroutine ini_conc_mod
-
-
-  subroutine end_conc_mod
-
-  implicit none
-
-    if (allocated(mmean)) deallocate(mmean)
-    if (allocated(Akknew)) deallocate(Akknew)
-    if (allocated(cpnew)) deallocate(cpnew)
-    if (allocated(rnew)) deallocate(rnew)
-
-  end subroutine end_conc_mod
+ end subroutine ini_conc_mod
+
+
+ subroutine end_conc_mod
+
+ implicit none
+
+  if (allocated(mmean)) deallocate(mmean)
+  if (allocated(Akknew)) deallocate(Akknew)
+  if (allocated(cpnew)) deallocate(cpnew)
+  if (allocated(rnew)) deallocate(rnew)
+
+ end subroutine end_conc_mod
+
+ subroutine init_r_cp_mu(ngrid,nlayer)
+
+ USE comcstfi_h, only: r, cpp, mugaz
+ implicit none
+! Input  
+ integer,intent(in) :: ngrid ! number of atmospheric columns
+ integer,intent(in) :: nlayer ! number of atmospheric layers
+
+! local variables
+ integer l, ig
+
+! Initialize R, Cp and mu as constant
+!        do l=1,nlayer
+!          do ig=1,ngrid
+!            rnew(ig,l)=r
+!            cpnew(ig,l)=cpp
+ !           mmean(ig,l)=mugaz
+!          enddo
+ !       enddo
+         rnew(:,:)=r
+         cpnew(:,:)=cpp
+         mmean(:,:)=mugaz
+     return           
+ end subroutine init_r_cp_mu
+
+
+ subroutine update_r_cp_mu_ak(ngrid,nlayer,nq,pplay,pt,pdt,pq,pdq,ptimestep)
+
+ use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,             &
+&                      igcm_o2, igcm_o3, igcm_h, igcm_h2,               &
+&                      igcm_oh, igcm_ho2, igcm_n2, igcm_ar,             &
+&                      igcm_h2o_vap, igcm_n, igcm_no, igcm_no2, &
+&                      igcm_n2d, igcm_co2plus, igcm_oplus,              &
+&                      igcm_o2plus, igcm_coplus, igcm_cplus,            &
+&                      igcm_nplus, igcm_noplus, igcm_n2plus,            &
+&                      igcm_hplus, igcm_hco2plus, mmol,         &
+&                      igcm_he, igcm_elec 
+ implicit none 
+ !=======================================================================
+! CALCULATION OF MEAN MOLECULAR MASS, Cp, Akk and R
+!
+! mmean(ngrid,nlayer)   amu
+! cpnew(ngrid,nlayer)   J/kg/K
+! rnew(ngrid,nlayer)    J/kg/K
+! akknew(ngrid,nlayer)  coefficient of thermal conduction
+!
+! version: April 2012 - Franck Lefevre
+!=======================================================================
+
+!     declarations
+ 
+    include "callkeys.h"
+
+!     input/output
+
+    integer,intent(in) :: ngrid ! number of atmospheric columns
+    integer,intent(in) :: nlayer ! number of atmospheric layers
+    integer,intent(in) :: nq ! number of tracers
+    real,intent(in) :: pplay(ngrid,nlayer)
+    real,intent(in) :: pt(ngrid,nlayer)
+    real,intent(in) :: pdt(ngrid,nlayer)
+    real,intent(in) :: pq(ngrid,nlayer,nq)
+    real,intent(in) :: pdq(ngrid,nlayer,nq)
+    real,intent(in) :: ptimestep
+
+!     local variables
+
+    integer       :: i, l, ig, iq
+    integer, save :: nbq
+    integer,allocatable,save :: niq(:)
+    real          :: ni(nq), ntot
+    real          :: zq(ngrid, nlayer, nq)
+    real          :: zt(ngrid, nlayer)
+    real,allocatable,save    :: aki(:)
+    real,allocatable,save    :: cpi(:)
+
+    logical, save :: firstcall = .true.
+
+!$OMP THREADPRIVATE(nbq,niq,aki,cpi,firstcall)
+
+    if (firstcall) then
+
+      ! allocate local saved arrays:
+       allocate(aki(nq))
+       allocate(cpi(nq))
+       allocate(niq(nq))
+!        find index of chemical tracers to use
+!        initialize thermal conductivity and specific heat coefficients
+!        !? values are estimated
+
+       nbq = 0 ! to count number of tracers used in this subroutine
+
+       if (igcm_co2 /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_co2
+          aki(nbq) = 3.072e-4
+          cpi(nbq) = 0.834e3
+       end if
+       if (igcm_co /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_co
+          aki(nbq) = 4.87e-4
+          cpi(nbq) = 1.034e3
+       end if
+       if (igcm_o /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_o
+          aki(nbq) = 7.59e-4
+          cpi(nbq) = 1.3e3
+       end if
+       if (igcm_o1d /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_o1d
+          aki(nbq) = 7.59e-4  !?
+          cpi(nbq) = 1.3e3    !?
+       end if
+       if (igcm_o2 /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_o2
+          aki(nbq) = 5.68e-4
+          cpi(nbq) = 0.9194e3
+       end if
+       if (igcm_o3 /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_o3
+          aki(nbq) = 3.00e-4  !?
+          cpi(nbq) = 0.800e3  !?
+       end if
+       if (igcm_h /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_h
+          aki(nbq) = 0.0
+          cpi(nbq) = 20.780e3
+       end if
+       if (igcm_h2 /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_h2
+          aki(nbq) = 36.314e-4
+          cpi(nbq) = 14.266e3
+       end if
+       if (igcm_oh /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_oh
+          aki(nbq)  = 7.00e-4 !?
+          cpi(nbq)  = 1.045e3
+       end if
+       if (igcm_ho2 /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_ho2
+          aki(nbq) = 0.0
+          cpi(nbq) = 1.065e3  !?
+       end if
+       if (igcm_n2 /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_n2
+          aki(nbq) = 5.6e-4
+          cpi(nbq) = 1.034e3
+       end if
+       if (igcm_ar /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_ar
+          aki(nbq) = 3.4e-4
+          cpi(nbq) = 5.2e2
+         ! aki/cpi from Leslie A. Guildner,
+         ! JOURNAL OF RESEARCH of the National Bureau of Standards- 
+         ! A. Physics and Chemistry
+         ! Vol. 79A, No.2, March-April 1975
+       end if
+       if (igcm_h2o_vap /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_h2o_vap
+          aki(nbq) = 0.0
+          cpi(nbq) = 1.870e3
+       end if
+       if (igcm_n /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_n
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+       if(igcm_no /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_no
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+       if(igcm_no2 /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_no2
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+       if(igcm_n2d /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_n2d
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+       if(igcm_he /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_he
+          aki(nbq) = 30.e-4
+          cpi(nbq) = 5.2e3
+         ! aki/cpi from Leslie A. Guildner,
+         ! JOURNAL OF RESEARCH of the National Bureau of Standards- 
+         ! A. Physics and Chemistry
+         ! Vol. 79A, No.2, March-April 1975
+       endif
+       if(igcm_co2plus /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_co2plus
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+       if(igcm_oplus /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_oplus
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+       if(igcm_o2plus /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_o2plus
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+       if(igcm_coplus /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_coplus
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+       if(igcm_cplus /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_cplus
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+       if(igcm_nplus /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_nplus
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+       if(igcm_noplus /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_noplus
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+       if(igcm_n2plus /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_n2plus
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+       if(igcm_hplus /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_hplus
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+       if(igcm_hco2plus /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_hco2plus
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+       if(igcm_elec /= 0) then
+          nbq = nbq + 1
+          niq(nbq) = igcm_elec
+          aki(nbq) = 0.0
+          cpi(nbq) = 0.0
+       endif
+
+      ! tell the world about it:
+       write(*,*) "concentrations: firstcall, nbq=",nbq
+       write(*,*) "  niq(1:nbq)=",niq(1:nbq)
+       write(*,*) "  aki(1:nbq)=",aki(1:nbq)
+       write(*,*) "  cpi(1:nbq)=",cpi(1:nbq)
+
+      firstcall = .false.
+
+    end if ! if (firstcall)
+
+!     update temperature
+
+    do l = 1,nlayer
+       do ig = 1,ngrid
+          zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
+       end do
+    end do
+
+!     update tracers
+
+    do l = 1,nlayer
+       do ig = 1,ngrid
+          do i = 1,nbq
+             iq = niq(i) 
+             zq(ig,l,iq) = max(1.e-30, pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep)
+         end do
+      end do
+    end do
+
+!     mmean : mean molecular mass
+!     rnew  : specific gas constant
+
+    mmean(:,:)  = 0.
+
+    do l = 1,nlayer
+       do ig = 1,ngrid
+          do i = 1,nbq
+             iq = niq(i) 
+             mmean(ig,l) = mmean(ig,l) + zq(ig,l,iq)/mmol(iq)
+          end do
+          mmean(ig,l) = 1./mmean(ig,l)
+          rnew(ig,l) = 8.314/mmean(ig,l)*1.e3     ! J/kg/K           
+       end do
+    end do
+
+!     cpnew  : specific heat
+!     akknew : thermal conductivity coefficient
+
+    cpnew(:,:)  = 0.
+    akknew(:,:) = 0.
+
+    do l = 1,nlayer
+       do ig = 1,ngrid
+          ntot = pplay(ig,l)/(1.381e-23*zt(ig,l))*1.e-6  ! in #/cm3
+          do i = 1,nbq
+             iq = niq(i) 
+             ni(iq) = ntot*zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
+             cpnew(ig,l) = cpnew(ig,l) + ni(iq)*cpi(i)
+             akknew(ig,l) = akknew(ig,l) + ni(iq)*aki(i)
+          end do 
+          cpnew(ig,l) = cpnew(ig,l)/ntot
+          akknew(ig,l) = akknew(ig,l)/ntot
+       end do
+    end do
+
+    return 
+ end subroutine update_r_cp_mu_ak
 
 end module conc_mod

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

@@ -66 +66 @@
       use nirco2abs_mod, only: nirco2abs
       use slope_mod, only: theta_sl, psi_sl, getslopes, param_slope
-      use conc_mod, only: rnew, cpnew, mmean
+      use conc_mod, only: init_r_cp_mu, update_r_cp_mu_ak, rnew, 
+     &                    cpnew, mmean
       use time_phylmdz_mod, only: iphysiq, day_step, ecritstart, daysec
       use dimradmars_mod, only: aerosol, totcloudfrac,
@@ -757 +758 @@
          endif
 #endif
-c        Initialize R and Cp as constant
-
-         if (.not.callthermos .and. .not.photochem) then
-                 do l=1,nlayer
-                  do ig=1,ngrid
-                   rnew(ig,l)=r
-                   cpnew(ig,l)=cpp
-                   mmean(ig,l)=mugaz
-                   enddo
-                  enddo  
-         endif         
+
+c        Initialize rnew cpnew and mmean as constant
+c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         call init_r_cp_mu(ngrid,nlayer)      
 
          if(callnlte.and.nltemodel.eq.2) call nlte_setup
@@ -779 +773 @@
      .                  albedo_h2o_cap
         ENDIF
+
+
 
 #ifndef MESOSCALE
@@ -894 +890 @@
 c-----------------------------------------------------------------------
 c    1.2.1 Compute mean mass, cp, and R
-c    concentrations outputs rnew(ngrid,nlayer), cpnew(ngrid,nlayer)
+c    update_r_cp_mu_ak outputs rnew(ngrid,nlayer), cpnew(ngrid,nlayer)
 c       , mmean(ngrid,nlayer) and Akknew(ngrid,nlayer)
 c    --------------------------------
 
       if(photochem.or.callthermos) then
-         call concentrations(ngrid,nlayer,nq,
+         call update_r_cp_mu_ak(ngrid,nlayer,nq,
      &                       zplay,pt,pdt,pq,pdq,ptimestep)
       endif