Changeset 6149 for LMDZ6/trunk/libf/phylmd/mo_simple_plumes.f90

Timestamp:

Apr 1, 2026, 6:36:00 PM (9 days ago)

Author:

asima

Message:

Bug fix for Simple Plume aerosols:
the orography was incorrectly accounted for in the aod vertical profile (since the original CMIP6 implementation)
ASima, FHourdin

File:

• LMDZ6/trunk/libf/phylmd/mo_simple_plumes.f90 (modified) (14 diffs)

LMDZ6/trunk/libf/phylmd/mo_simple_plumes.f90

@@ -21 +21 @@
 !! @par Copyright
 !!
-!! 2026-03, A.Sima (LMD): for calculation optimisation, the main subroutine sp_aop_profile is split in 2 subroutines :
+!! 2026-03 
+!!  A.Sima (LMD): for calculation optimisation, the main subroutine sp_aop_profile is split in 2 subroutines :
 !!      sp_aod550_profile : calculates aod550um profile from macv2sp data (which are for 550um themselves), and dNovrN factor   
 !!           at each timestep is only called in macv2sp once ;
@@ -28 +29 @@
 !!            - aod diagnostics at 443 and 865um,
 !!            - optical properties (aod, ssa, asy) for the Nwvmax=25 wavelengths filling the 6 RRTM bands   
+!!  ASima, FH : bug corrected in aod vertical profile calculation, concerning the impact of orography
 !
 MODULE mo_simple_plumes
@@ -37 +39 @@
   INTEGER, PARAMETER, PUBLIC :: nplumes   = 9       !< Number of plumes
 
-!  INTEGER, PARAMETER ::                        &
-!       nplumes   = 9                          ,& !< Number of plumes
   INTEGER, PARAMETER ::                        &
        nfeatures = 2                          ,& !< Number of features per plume
@@ -73 +73 @@
 !$OMP THREADPRIVATE(sig_lon_E,sig_lon_W,sig_lat_E,sig_lat_W,theta,ftr_weight,year_weight,ann_cycle)
 
+  REAL, SAVE ::  beta_sum       (nplumes)                ! numerically calculated integral of beta function (ASima,FH)
+!$OMP THREADPRIVATE(beta_sum)
+
   REAL ::                                      &
        time_weight    (nfeatures,nplumes)     ,& !< Time weights 
@@ -85 +88 @@
   ! climatology.  The information needs to be either read by each processor or distributed to processors.
   !
+  ! (ASima,FH : added call to sp_beta_sup to compute the integral of beta profile by plume, only once at initialisation
+  ! 
   SUBROUTINE sp_setup
     !
@@ -290 +295 @@
     CALL bcast(ann_cycle)
     !
+    !(ASima,FH) Calculate beta function integral for the 9 plumes
+    CALL sp_beta_sum
+    !
     sp_initialized = .TRUE.
     !
@@ -295 +303 @@
     !
   END SUBROUTINE sp_setup
+  !
+  ! ------------------------------------------------------------------------------------------------------------------------
+  !  (ASima, FH) sp_beta_sum : numerically calculates the integral of beta function for each of the 9 plumes : 
+  !     B(p_i, q_i) in eq.(8) of Stevens et al (2017)
+  !  To be only called once in the beginning of an integration period, same as (and just after) sp_setup.
+  !  The beta-function vertical profiles of optical properties apply from 0 (sea level) to zmax=15 km.
+  !  They only depend of plume caracteristics read in SP aerosol file.
+  !  
+  SUBROUTINE sp_beta_sum
+    !
+    ! ---------- 
+    !
+    ! 
+  INTEGER, PARAMETER :: nmax=1000       ! nb. of layers for discretisation of 0 to 1 interval
+
+!REAL, DIMENSION(nplumes) :: beta_sum   ! numerically calculated integral of beta function
+
+  INTEGER :: iplume, k
+
+  REAL :: eta, d_eta                    ! normalised height and thickness of a discretisation layer 
+  REAL :: beta_eta                      ! beta function value at relative height eta
+
+
+  d_eta = 1/float(nmax)   !  = layer thickness in discretisation of 0 to 1 interval
+
+  ! For each plume, sum over beta function values at the middle of each d_eta layer
+  DO iplume=1,nplumes  
+   
+      beta_sum(iplume) = 0.  !  initialisation
+
+      DO k=1,nmax
+          eta = (k - 0.5)/float(nmax)
+          beta_eta = eta**(beta_a(iplume)-1.) * (1.-eta)**(beta_b(iplume)-1.)
+          beta_sum(iplume) = beta_sum(iplume) + beta_eta * d_eta
+      ENDDO
+   ENDDO
+!   DO iplume=1,nplumes
+!    print*,'beta_sum(',iplume,') = ', beta_sum(iplume)
+!   ENDDO
+   !
+   RETURN
+   !
+  END SUBROUTINE sp_beta_sum
   !
   ! ------------------------------------------------------------------------------------------------------------------------
@@ -339 +390 @@
   !              It also computes the dNovrN factor used to mimic Twomey (aci) effect by multiplying preind. cloud droplet number
   !
-  SUBROUTINE sp_aod550_profile                                                                           ( &
-       nlevels        ,ncol           ,oro            ,lon            ,lat            , &
+  SUBROUTINE sp_aod550_profile                                                 ( &
+       nlevels        ,ncol           ,lon            ,lat            , &
        decimal_year   ,z              ,dz             ,dNovrN         , aod550   )
 
@@ -352 +403 @@
     REAL, INTENT(IN)           :: &
          decimal_year,            & !< Fractional Year (1903.0 is the 0Z on the first of January 1903, Gregorian)
-         oro(ncol),               & !< orographic height (m)
          lon(ncol),               & !< longitude 
          lat(ncol),               & !< latitude
@@ -363 +413 @@
 
     INTEGER                    :: iplume, icol, k
+    
+    REAL, PARAMETER            :: zmax = 15000 ! max altitude up to which beta-function profiles apply
 
     REAL                       ::  &
-         eta(ncol,nlevels),        & !< normalized height (by 15 km)
-         z_beta(ncol,nlevels),     & !< profile for scaling column optical depth
+         eta_lmdz(ncol,nlevels),   & !< normalized height (by zmax = 15 km)  
+         d_eta_lmdz(ncol,nlevels), & !< normalized height (by zmax = 15 km)
          prof(ncol,nlevels),       & !< scaled profile (by beta function)
-         beta_sum(ncol),           & !< vertical sum of beta function
          cw_an(ncol),              & !< column weight for simple plume (anthropogenic) AOD at 550 nm
          cw_bg(ncol),              & !< column weight for fine-mode natural background AOD at 550 nm
@@ -401 +452 @@
     DO k=1,nlevels
       DO icol=1,ncol
-        z_beta(icol,k)   = MERGE(1.0, 0.0, z(icol,k) >= oro(icol))
-        eta(icol,k)      = MAX(0.0,MIN(1.0,z(icol,k)/15000.))
+         eta_lmdz(icol,k)      = MIN(1.0,z(icol,k)/15000.)
+         d_eta_lmdz(icol,k)      = dz(icol,k)/15000.
       ENDDO
     ENDDO
@@ -411 +462 @@
     ENDDO
     !
+!    DO iplume=1,nplumes
+!     print*,'IN sp_aod550 : beta_sum(',iplume,') = ', beta_sum(iplume)
+!    ENDDO
+
     ! sum contribution from plumes to construct composite profiles of aerosol optical properties
     !
@@ -417 +472 @@
       ! calculate vertical distribution function from parameters of beta distribution
       !
-      DO icol=1,ncol
-        beta_sum(icol) = 0.
-      ENDDO
       DO k=1,nlevels
         DO icol=1,ncol
-          prof(icol,k)   = (eta(icol,k)**(beta_a(iplume)-1.) * (1.-eta(icol,k))**(beta_b(iplume)-1.)) * dz(icol,k)
-          beta_sum(icol) = beta_sum(icol) + prof(icol,k)
-        ENDDO
-      ENDDO
-      DO k=1,nlevels
-        DO icol=1,ncol
-          prof(icol,k)   = ( prof(icol,k) / beta_sum(icol) ) * z_beta(icol,k)
+           prof(icol,k)   = (eta_lmdz(icol,k)**(beta_a(iplume)-1.) * (1.-eta_lmdz(icol,k))**(beta_b(iplume)-1.)) * d_eta_lmdz(icol,k)/beta_sum(iplume)
         ENDDO
       ENDDO
@@ -495 +541 @@
   ! ------------------------------------------------------------------------------------------------------------------------ 
 
-  SUBROUTINE sp_aop_profile                                                                        ( &
-       nlevels        ,ncol           ,lambda                     , &
-       aod550         ,aod_prof       ,ssa_prof       , &
-       asy_prof       )
+  SUBROUTINE sp_aop_profile                                    ( &
+       nlevels        ,ncol           ,lambda                  , &
+       aod550         ,aod_prof       ,ssa_prof       ,asy_prof  )
     !
     ! ----------