Changeset 3197 for trunk/LMDZ.PLUTO

Timestamp:

Feb 1, 2024, 11:17:06 AM (10 months ago)

Author:

afalco

Message:

Pluto PCM:
Fixed some flags for aerosols
AF

Location:

trunk/LMDZ.PLUTO

Files:

• deftank/callphys.def (modified) (1 diff)
• deftank/traceur.def (modified) (1 diff)
• libf/phypluto/callcorrk.F90 (modified) (4 diffs)
• libf/phypluto/inifis_mod.F90 (modified) (28 diffs)
• libf/phypluto/physiq_mod.F90 (modified) (1 diff)
• libf/phypluto/turbdiff_mod.F90 (modified) (20 diffs)

trunk/LMDZ.PLUTO/deftank/callphys.def

@@ -87 +87 @@
 sedimentation = .true.
 
-naerkind = 0
+naerkind = 1
 hazeprop = optprop_tholins_fractal100nm.dat
 

trunk/LMDZ.PLUTO/deftank/traceur.def

@@ -1 @@
-0
+1
+n2

trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90

@@ -486 +486 @@
       end do
 
-      ! Get 3D aerosol optical properties.
-      call aeroptproperties(ngrid,nlayer,reffrad,nueffrad,         &
-           QVISsQREF3d,omegaVIS3d,gVIS3d,                          &
-           QIRsQREF3d,omegaIR3d,gIR3d,                             &
-           QREFvis3d,QREFir3d)
-
-      ! Get aerosol optical depths.
-      call aeropacity(ngrid,nlayer,nq,pplay,pplev, pt,pq,aerosol,      &
-           reffrad,nueffrad,QREFvis3d,QREFir3d,                             &
-           tau_col,cloudfrac,totcloudfrac,clearsky)
+      if (aerohaze) then
+      !    if (haze_radproffix) then
+      !       call haze_reffrad_fix(ngrid,nlayer,zzlay,&
+      !           reffrad,nueffrad)
+      !     endif
+
+         ! Get 3D aerosol optical properties.
+         call aeroptproperties(ngrid,nlayer,reffrad,nueffrad,         &
+            QVISsQREF3d,omegaVIS3d,gVIS3d,                          &
+            QIRsQREF3d,omegaIR3d,gIR3d,                             &
+            QREFvis3d,QREFir3d)
+
+         ! Get aerosol optical depths.
+         call aeropacity(ngrid,nlayer,nq,pplay,pplev, pt,pq,aerosol,      &
+            reffrad,nueffrad,QREFvis3d,QREFir3d,                             &
+            tau_col,cloudfrac,totcloudfrac,clearsky)
+      endif
 
 !-----------------------------------------------------------------------
@@ -513 +520 @@
 
 
+         if (aerohaze) then
             do iaer=1,naerkind
                ! Shortwave.
@@ -620 +628 @@
             end do ! naerkind
 
+         endif ! aerohaze
 !-----------------------------------------------------------------------
 !     Aerosol optical depths
 !-----------------------------------------------------------------------
-
+        IF (aerohaze) THEN
          do iaer=1,naerkind     ! a bug was here
             do k=0,nlayer-1
@@ -642 +651 @@
 
          end do ! naerkind
+        ELSE
+         tauaero(:,:)=0
+        ENDIF ! aerohaze
 
          ! Albedo and Emissivity.

trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90

@@ -34 +34 @@
 !   -------
 !
-!   Initialisation for the physical parametrisations of the LMD 
+!   Initialisation for the physical parametrisations of the LMD
 !   Generic Model.
 !
@@ -78 +78 @@
   INTEGER ig
   CHARACTER(len=20) :: rname="inifis" ! routine name, for messages
- 
+
   EXTERNAL iniorbit,orbite
   EXTERNAL SSUM
   REAL SSUM
- 
+
   ! Initialize flags lunout, prt_level, debug (in print_control_mod)
   CALL init_print_control
@@ -110 +110 @@
   ! do we read a startphy.nc file? (default: .true.)
   call getin_p("startphy_file",startphy_file)
-  
+
 ! --------------------------------------------------------------
 !  Reading the "callphys.def" file controlling some key options
 ! --------------------------------------------------------------
-     
+
   IF (is_master) THEN
     ! check if 'callphys.def' file is around
@@ -139 +139 @@
      if (is_master) write(*,*) "Haze optical properties datafile"
      hazeprop="../../../datadir"  ! default path
-     call getin_p("hazeprop",hazeprop) 
+     call getin_p("hazeprop",hazeprop)
      if (is_master) write(*,*) trim(rname)//" hazeprop = ",trim(hazeprop)
 
@@ -175 +175 @@
      call getin_p("season",season)
      if (is_master) write(*,*) trim(rname)//": season = ",season
-     
+
      if (is_master) write(*,*) trim(rname)//&
        ": No seasonal cycle: initial day to lock the run during restart"
@@ -235 +235 @@
      call getin_p("tplanckmin",tplanckmin)
      if (is_master) write(*,*) trim(rname)//": tplanckmin = ",tplanckmin
- 
+
      if (is_master) write(*,*) trim(rname)//&
        ": Maximum atmospheric temperature for Planck function integration ?"
@@ -241 +241 @@
      call getin_p("tplanckmax",tplanckmax)
      if (is_master) write(*,*) trim(rname)//": tplanckmax = ",tplanckmax
- 
+
      if (is_master) write(*,*) trim(rname)//&
        ": Temperature step for Planck function integration ?"
@@ -247 +247 @@
      call getin_p("dtplanck",dtplanck)
      if (is_master) write(*,*) trim(rname)//": dtplanck = ",dtplanck
- 
+
      if (is_master) write(*,*) trim(rname)//&
        ": call gaseous absorption in the visible bands?"//&
@@ -254 +254 @@
      call getin_p("callgasvis",callgasvis)
      if (is_master) write(*,*) trim(rname)//": callgasvis = ",callgasvis
-        
+
      if (is_master) write(*,*) trim(rname)//&
        ": call continuum opacities in radiative transfer ?"//&
@@ -290 +290 @@
      call getin_p("callsoil",callsoil)
      if (is_master) write(*,*) trim(rname)//": callsoil = ",callsoil
-         
+
      if (is_master) write(*,*)trim(rname)//&
        ": Rad transfer is computed every iradia", &
@@ -297 +297 @@
      call getin_p("iradia",iradia)
      if (is_master) write(*,*)trim(rname)//": iradia = ",iradia
-       
+
      if (is_master) write(*,*)trim(rname)//": Rayleigh scattering ?"
      rayleigh=.false.
@@ -347 +347 @@
      call getin_p("nsoilmx",nsoilmx)
      if (is_master) write(*,*)trim(rname)//": nsoilmx=",nsoilmx
-     
+
      if (is_master) write(*,*)trim(rname)//&
        ": Thickness of topmost soil layer (m)?"
@@ -353 +353 @@
      call getin_p("lay1_soil",lay1_soil)
      if (is_master) write(*,*)trim(rname)//": lay1_soil = ",lay1_soil
-     
+
      if (is_master) write(*,*)trim(rname)//&
        ": Coefficient for soil layer thickness distribution?"
@@ -380 +380 @@
       call getin_p("global1d",global1d)
       write(*,*) "global1d = ",global1d
-      
+
       ! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
       if (global1d.and.diurnal) then
@@ -464 +464 @@
      if(kastprof)then
         radfixed=.true.
-     endif  
+     endif
 
       if (is_master) write(*,*)trim(rname)//&
@@ -472 +472 @@
      if (is_master) write(*,*)trim(rname)//": naerkind = ",naerkind
 
-     
+
 !***************************************************************
          !! TRACERS options
@@ -478 +478 @@
      if (is_master)write(*,*)trim(rname)//&
       "call N2 condensation ?"
-     N2cond=.true. ! default value
-     call getin_p("N2cond",N2cond)
-     if (is_master)write(*,*)trim(rname)//&
-      " N2cond = ",N2cond
+     n2cond=.true. ! default value
+     call getin_p("n2cond",n2cond)
+     if (is_master)write(*,*)trim(rname)//&
+      " n2cond = ",n2cond
 
      if (is_master)write(*,*)trim(rname)//&
@@ -655 +655 @@
      if (is_master)write(*,*)trim(rname)//&
      "nb_monomer = ",nb_monomer
-     
+
      if (is_master)write(*,*)trim(rname)//&
      "fixed radius profile from txt file ?"
@@ -750 +750 @@
      if (is_master)write(*,*)trim(rname)//&
       " thres_ch4ice = ",thres_ch4ice
-     fdch4_latn=-1. 
-     fdch4_lats=0. 
-     fdch4_lone=0. 
-     fdch4_lonw=-1. 
-     fdch4_depalb=0.5 
-     fdch4_finalb=0.95 
-     fdch4_maxalb=0.99 
-     fdch4_ampl=200. 
-     fdch4_maxice=100. 
+     fdch4_latn=-1.
+     fdch4_lats=0.
+     fdch4_lone=0.
+     fdch4_lonw=-1.
+     fdch4_depalb=0.5
+     fdch4_finalb=0.95
+     fdch4_maxalb=0.99
+     fdch4_ampl=200.
+     fdch4_maxice=100.
      call getin_p("  fdch4_latn",fdch4_latn)
      call getin_p("  fdch4_lats",fdch4_lats)
@@ -903 +903 @@
      if (is_master)write(*,*)trim(rname)//&
       " tholateq = ",tholateq
-     tholatn=-1. 
-     tholats=0. 
-     tholone=0. 
-     tholonw=-1. 
+     tholatn=-1.
+     tholats=0.
+     tholone=0.
+     tholonw=-1.
      alb_tho_spe=0.1 ! default value
      emis_tho_spe=1. ! default value
@@ -1039 +1039 @@
      if (is_master) write(*,*)trim(rname)//&
        ": optprop_back2lay_ir = ",trim(optprop_back2lay_ir)
-     
+
      if (is_master) write(*,*)trim(rname)//&
        ": TWOLAY AEROSOL: pres_bottom_tropo? in pa"
@@ -1132 +1132 @@
 
      ! This is necessary, we just set the number of aerosol layers
-     IF(.NOT.ALLOCATED(aeronlay_tauref))      ALLOCATE(aeronlay_tauref(nlayaero))     
-     IF(.NOT.ALLOCATED(aeronlay_lamref))      ALLOCATE(aeronlay_lamref(nlayaero))     
-     IF(.NOT.ALLOCATED(aeronlay_choice))      ALLOCATE(aeronlay_choice(nlayaero))     
-     IF(.NOT.ALLOCATED(aeronlay_pbot))        ALLOCATE(aeronlay_pbot(nlayaero))     
-     IF(.NOT.ALLOCATED(aeronlay_ptop))        ALLOCATE(aeronlay_ptop(nlayaero))     
-     IF(.NOT.ALLOCATED(aeronlay_sclhght))     ALLOCATE(aeronlay_sclhght(nlayaero))     
-     IF(.NOT.ALLOCATED(aeronlay_size))        ALLOCATE(aeronlay_size(nlayaero))     
-     IF(.NOT.ALLOCATED(aeronlay_nueff))       ALLOCATE(aeronlay_nueff(nlayaero))     
-     IF(.NOT.ALLOCATED(optprop_aeronlay_ir))  ALLOCATE(optprop_aeronlay_ir(nlayaero))     
-     IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))     
+     IF(.NOT.ALLOCATED(aeronlay_tauref))      ALLOCATE(aeronlay_tauref(nlayaero))
+     IF(.NOT.ALLOCATED(aeronlay_lamref))      ALLOCATE(aeronlay_lamref(nlayaero))
+     IF(.NOT.ALLOCATED(aeronlay_choice))      ALLOCATE(aeronlay_choice(nlayaero))
+     IF(.NOT.ALLOCATED(aeronlay_pbot))        ALLOCATE(aeronlay_pbot(nlayaero))
+     IF(.NOT.ALLOCATED(aeronlay_ptop))        ALLOCATE(aeronlay_ptop(nlayaero))
+     IF(.NOT.ALLOCATED(aeronlay_sclhght))     ALLOCATE(aeronlay_sclhght(nlayaero))
+     IF(.NOT.ALLOCATED(aeronlay_size))        ALLOCATE(aeronlay_size(nlayaero))
+     IF(.NOT.ALLOCATED(aeronlay_nueff))       ALLOCATE(aeronlay_nueff(nlayaero))
+     IF(.NOT.ALLOCATED(optprop_aeronlay_ir))  ALLOCATE(optprop_aeronlay_ir(nlayaero))
+     IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
 
      if (is_master) write(*,*)trim(rname)//&
@@ -1161 +1161 @@
        ": Generic n-layer aerosols: Vertical profile choice : "
        write(*,*)trim(rname)//&
-       "             (1) Tau btwn ptop and pbot follows atm. scale height" 
+       "             (1) Tau btwn ptop and pbot follows atm. scale height"
        write(*,*)trim(rname)//&
        "         or  (2) Tau above pbot follows its own scale height"
@@ -1214 +1214 @@
      if (is_master) write(*,*)trim(rname)//&
        ": optprop_aeronlay_ir = ",optprop_aeronlay_ir
-     
+
 
 !=================================
@@ -1250 +1250 @@
 
      if (is_master) write(*,*)trim(rname)//": Generic Condensation of tracers ?"
-     generic_condensation=.false. !default value 
+     generic_condensation=.false. !default value
      call getin_p("generic_condensation",generic_condensation)
      if (is_master) write(*,*)trim(rname)//": generic_condensation = ",generic_condensation
-     
+
 ! Test of incompatibility:
      if (is_master) write(*,*)trim(rname)//": Spectral Dependant albedo ?"
@@ -1269 +1269 @@
      call getin_p("albedosnow",albedosnow)
      if (is_master) write(*,*)trim(rname)//": albedosnow = ",albedosnow
-             
+
      if (is_master) write(*,*)trim(rname)//": N2 ice albedo ?"
      albedon2ice=0.5       ! default value
@@ -1331 +1331 @@
        if (is_master) write(*,*)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
        if (is_master) write(*,*)'New Mg from callfis is ',mugaz,'amu'
-  
+
      endif ! of if (cpp_mugaz_mode == 1)
      call su_gases
@@ -1382 +1382 @@
   ! initialize variables and allocate arrays in radcommon_h
   call ini_radcommon_h(naerkind)
-   
+
   ! allocate "comsoil_h" arrays
   call ini_comsoil_h(ngrid)
-    
+
   END SUBROUTINE inifis
 

trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90

@@ -19 +19 @@
       use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind, corrkdir, banddir
       use generic_cloud_common_h, only : epsi_generic, Psat_generic
-      use thermcell_mod, only: init_thermcell_mod
       use gases_h, only: gnom, gfrac
       use radcommon_h, only: sigma, glat, grav, BWNV, WNOI, DWNI, DWNV, WNOV

trunk/LMDZ.PLUTO/libf/phypluto/turbdiff_mod.F90

@@ -1 +1 @@
 module turbdiff_mod
-      
+
 implicit none
-      
+
 contains
-      
+
       subroutine turbdiff(ngrid,nlay,nq,          &
-          ptimestep,pcapcal,                    &   
+          ptimestep,pcapcal,                    &
           pplay,pplev,pzlay,pzlev,pz0,                 &
           pu,pv,pt,ppopsk,pq,ptsrf,pemis,pqsurf,       &
@@ -24 +24 @@
 
 !==================================================================
-!     
+!
 !     Purpose
 !     -------
 !     Turbulent diffusion (mixing) for pot. T, U, V and tracers
-!     
+!
 !     Implicit scheme
 !     We start by adding to variables x the physical tendencies
@@ -34 +34 @@
 !
 !     x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
-!     
+!
 !     Authors
-!     ------- 
+!     -------
 !     F. Hourdin, F. Forget, R. Fournier (199X)
 !     R. Wordsworth, B. Charnay (2010)
@@ -43 +43 @@
 !               by accounting for dissipation of turbulent kinetic energy.
 !         - Accounting for lost mean flow kinetic energy should come soon.
-!     
+!
 !==================================================================
 
@@ -75 +75 @@
 !     Tracers
 !     --------
-      integer,intent(in) :: nq 
+      integer,intent(in) :: nq
       real,intent(in) :: pqsurf(ngrid,nq)
-      real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) 
-      real,intent(out) :: pdqdif(ngrid,nlay,nq) 
-      real,intent(out) :: pdqsdif(ngrid,nq) 
-      
+      real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
+      real,intent(out) :: pdqdif(ngrid,nlay,nq)
+      real,intent(out) :: pdqsdif(ngrid,nq)
+
 !     local
 !     -----
@@ -110 +110 @@
       LOGICAL,SAVE :: firstcall=.true.
 !$OMP THREADPRIVATE(firstcall)
-      
+
 !     Tracers
 !     -------
@@ -130 +130 @@
       real cd0, roughratio
 
-      real dqsdif_total(ngrid) 
-      real zq0(ngrid) 
+      real dqsdif_total(ngrid)
+      real zq0(ngrid)
 
 
@@ -199 +199 @@
 !
 !     Source of turbulent kinetic energy at the surface
-!     ------------------------------------------------- 
+!     -------------------------------------------------
 !     Formula is Cd_0 = (karman / log[1+z1/z0])^2
 
@@ -210 +210 @@
        if(nosurf)then
          zcdv_true(ig)=0.D+0 !JL12 disable atm/surface momentum flux
-         zcdh_true(ig)=0.D+0 !JL12 disable sensible heat flux 
+         zcdh_true(ig)=0.D+0 !JL12 disable sensible heat flux
        endif
       ENDDO
@@ -221 +221 @@
 
 !     Turbulent diffusion coefficients in the boundary layer
-!     ------------------------------------------------------ 
-
-      call vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay,pu,pv,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature 
-      
+!     ------------------------------------------------------
+
+      call vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay,pu,pv,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature
+
 !     Adding eddy mixing to mimic 3D general circulation in 1D
 !     R. Wordsworth & F. Forget (2010)
@@ -246 +246 @@
 !            do ig=1,ngrid
 !               zkh(ig,ilev) = max(kzz_eddy(ilev),zkh(ig,ilev))
-!               zkv(ig,ilev) = max(kzz_eddy(ilev),zkv(ig,ilev)) 
+!               zkv(ig,ilev) = max(kzz_eddy(ilev),zkv(ig,ilev))
 !            end do
 
             do ig=1,ngrid
                zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
-               zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev)) 
+               zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
             end do
          end do
@@ -263 +263 @@
 !JL12 calculate the flux coefficients (tables multiplied elements by elements)
       zfluxv(1:ngrid,1:nlay)=zkv(1:ngrid,1:nlay)*zb0(1:ngrid,1:nlay)
-      
+
 !-----------------------------------------------------------------------
 !     4. Implicit inversion of u
@@ -312 +312 @@
 
       DO ig=1,ngrid
-         z1(ig)=1./(zmass(ig,nlay)+zfluxv(ig,nlay)) 
+         z1(ig)=1./(zmass(ig,nlay)+zfluxv(ig,nlay))
          zcv(ig,nlay)=zmass(ig,nlay)*zv(ig,nlay)*z1(ig)
          zdv(ig,nlay)=zfluxv(ig,nlay)*z1(ig)
@@ -412 +412 @@
          DO ig=1,ngrid
             z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1) &
-                + pfluxsrf(ig)*ptimestep + zdplanck(ig)*ptsrf(ig) 
-            z2(ig) = pcapcal(ig)+zdplanck(ig)+cpp*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1)) 
+                + pfluxsrf(ig)*ptimestep + zdplanck(ig)*ptsrf(ig)
+            z2(ig) = pcapcal(ig)+zdplanck(ig)+cpp*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1))
             ztsrf(ig) = z1(ig) / z2(ig)
             pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep
             zt(ig,1)   = zct(ig,1) + zdt(ig,1)*ztsrf(ig)
          ENDDO
-! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme 
+! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
 
 
@@ -448 +448 @@
 !     Implicit inversion of q
 !     -----------------------
-         do iq=1,nq 
+         do iq=1,nq
 
             if (iq.ne.ivap) then
@@ -456 +456 @@
                   zcq(ig,nlay)=zmass(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
                   zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
-               ENDDO 
-            
+               ENDDO
+
                DO ilay=nlay-1,2,-1
                   DO ig=1,ngrid
@@ -466 +466 @@
                ENDDO
 
+               do ig=1,ngrid
+                  z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
+                  zcq(ig,1)=(zmass(ig,1)*zq(ig,1,iq)+zfluxq(ig,2)*zcq(ig,2)+(-pdqsdif(ig,iq))*ptimestep)*z1(ig)
+                       ! tracer flux from surface
+                       ! currently pdqsdif always zero here,
+                       ! so last line is superfluous
+               enddo
+
 !     Starting upward calculations for simple tracer mixing (e.g., dust)
                do ig=1,ngrid
@@ -483 +491 @@
 !     and assuming an infinite source of water on the ground
 !     ------------------------------------------------------------------
-        
+
       endif                     ! tracer
 
@@ -498 +506 @@
          enddo
       enddo
-      
+
       DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
          sensibFlux(ig)=cpp*zfluxt(ig,1)/ptimestep*(zt(ig,1)*zovExner(ig,1)-ztsrf(ig))
@@ -511 +519 @@
             enddo
          enddo
-      endif 
+      endif
    end subroutine turbdiff
-      
+
 end module turbdiff_mod