Ignore:
Timestamp:
Mar 26, 2026, 7:09:02 PM (2 weeks ago)
Author:
asima
Message:

Major changes (don't change results):
Subroutine sp_aop_profile in the mo_simple_plume.f90 module is split in 2 subroutines:

  • sp_aod550_profile, for calculations at 550um ; only called in macv2sp once per timestep ;
  • sp_aop_profile only keeps the calculations of profiles of optical properties at other wavelengths ;

it is called in macv2sp for aod diagnostics at 443 and 865um,
and for each of the 25 wavelengths filling the 6 RRTM bands.

Minor :

  • subroutine name MACv2SP changed in macv2sp for coherence with file name ;
  • dNovrN becomes an argument INTENT(OUT) of macv2sp (instead of making USE from phys_local_var_mod) ;
  • lfactor variable name was used with 2 different meanings : eq 10 and eq 12 from Stevens et al. (2017); celui from eq(10) was renamed as 'lextinct'
  • parameter nmon = 12 : removed because not used ;
  • variable aod_lmd : renamed as aod_lmdz.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/mo_simple_plumes.f90

    r6102 r6128  
    2020!!
    2121!! @par Copyright
    22 !!
     22!!
     23!! 2026-03, A.Sima (LMD): for calculation optimisation, the main subroutine sp_aop_profile is split in 2 subroutines :
     24!!      sp_aod550_profile : calculates aod550um profile from macv2sp data (which are for 550um themselves), and dNovrN factor   
     25!!           at each timestep is only called in macv2sp once ;
     26!!      sp_aop_profile : uses optical properties (aod, ssa, asy) at 550 to calculate their profiles for another wavelength ;
     27!!           at each timestep is called in macv2sp for every required wavelength separately :
     28!!            - aod diagnostics at 443 and 865um,
     29!!            - optical properties (aod, ssa, asy) for the Nwvmax=25 wavelengths filling the 6 RRTM bands   
    2330!
    2431MODULE mo_simple_plumes
    … …  
    2835  IMPLICIT NONE
    2936
     37  INTEGER, PARAMETER, PUBLIC :: nplumes   = 9       !< Number of plumes
     38
     39!  INTEGER, PARAMETER ::                        &
     40!       nplumes   = 9                          ,& !< Number of plumes
    3041  INTEGER, PARAMETER ::                        &
    31        nplumes   = 9                          ,& !< Number of plumes
    3242       nfeatures = 2                          ,& !< Number of features per plume
    3343       ntimes    = 52                         ,& !< Number of times resolved per year (52 => weekly resolution)
    … …  
    6777       time_weight_bg (nfeatures,nplumes)        !< as time_weight but for natural background in Twomey effect
    6878
    69   PUBLIC sp_aop_profile, sp_setup
     79  PUBLIC sp_aod550_profile, sp_aop_profile, sp_setup
    7080
    7181CONTAINS
    … …  
    94104    CALL getin_p('aerosols_SP_coef_bN', aerosols_SP_coef_bN)
    95105    CALL getin_p('aerosols_SP_aod_bg_gl', aerosols_SP_aod_bg_gl)
    96     !print *,'aerosols_SP_forcing_year=',aerosols_SP_forcing_year, 'aerosols_SP_coef_bN = ', aerosols_SP_coef_bN, 'aerosols_SP_aod_bg_gl = ', aerosols_SP_aod_bg_gl
     106    print *, 'aerosols_SP_forcing_year=',aerosols_SP_forcing_year
     107    print *, 'aerosols_SP_coef_bN = ', aerosols_SP_coef_bN
     108    print *, 'aerosols_SP_aod_bg_gl = ', aerosols_SP_aod_bg_gl
    97109 
    98110    !--only one processor reads the input data
    … …  
    295307    !
    296308    REAL, INTENT(IN) ::  &
    297          decimal_year      !< Fractional Year (1850.0 - 2100.99)
     309         decimal_year      !< Decimal Year (1850.0 - 2100.99)
    298310
    299311    INTEGER          ::  &
    300          iyear          ,& !< Integer year values between 1 and 156 (1850-2100)
     312         iyear          ,& !< Integer year values between 1 and 251 (1850-2100) ; in 2026 : data available for 1850-2023 (NaN after)
    301313         iweek          ,& !< Integer index (between 1 and ntimes); for ntimes=52 this corresponds to weeks (roughly)
    302314         iplume            ! plume number
    … …  
    321333  END SUBROUTINE set_time_weight
    322334  !
    323   ! ------------------------------------------------------------------------------------------------------------------------
    324   ! SP_AOP_PROFILE:  This subroutine calculates the simple plume aerosol and cloud active optical properties based on the
    325   ! the simple plume fit to the MPI Aerosol Climatology (Version 2).  It sums over nplumes to provide a profile of aerosol
    326   ! optical properties on a host models vertical grid.
     335  !------------------------------------------------------------------------------------------------------------------------
     336  ! sp_aod550_profile : This subroutine calculates Simple Plume aerosol aod550(nm) profile from MACv2SP (550um) data
     337  !             (Stevens et al 2027 ; designed to fit the MPI Aerosol Climatology (Version 2) by Kinne, 2018),
     338  !              It sums over nplumes to provide aod550 profile on a host model's vertical grid.
     339  !              It also computes the dNovrN factor used to mimic Twomey (aci) effect by multiplying preind. cloud droplet number
    327340  !
    328   SUBROUTINE sp_aop_profile                                                                           ( &
    329        nlevels        ,ncol           ,lambda         ,oro            ,lon            ,lat            , &
    330        decimal_year        ,z              ,dz             ,dNovrN         ,aod_prof       ,ssa_prof       , &
    331        asy_prof       )
     341  SUBROUTINE sp_aod550_profile                                                                           ( &
     342       nlevels        ,ncol           ,oro            ,lon            ,lat            , &
     343       decimal_year   ,z              ,dz             ,dNovrN         , aod550   )
     344
    332345    !
    333346    ! ----------
    … …  
    338351
    339352    REAL, INTENT(IN)           :: &
    340          lambda,                  & !< wavelength
    341353         decimal_year,            & !< Fractional Year (1903.0 is the 0Z on the first of January 1903, Gregorian)
    342354         oro(ncol),               & !< orographic height (m)
    … …  
    348360    REAL, INTENT(OUT)          :: &
    349361         dNovrN(ncol)           , & !< anthropogenic increase in cloud drop number concentration (factor)
    350          aod_prof(ncol,nlevels) , & !< profile of aerosol optical depth
    351          ssa_prof(ncol,nlevels) , & !< profile of single scattering albedo
    352          asy_prof(ncol,nlevels)     !< profile of asymmetry parameter
     362         aod550(ncol,nlevels,nplumes) !< anthropogenic aod550 profiles by individual plumes
    353363
    354364    INTEGER                    :: iplume, icol, k
    … …  
    359369         prof(ncol,nlevels),       & !< scaled profile (by beta function)
    360370         beta_sum(ncol),           & !< vertical sum of beta function
    361          ssa(ncol),                & !< single scattering albedo
    362          asy(ncol),                & !< asymmetry parameter
    363371         cw_an(ncol),              & !< column weight for simple plume (anthropogenic) AOD at 550 nm
    364372         cw_bg(ncol),              & !< column weight for fine-mode natural background AOD at 550 nm
    … …  
    379387         f2,                       & !< contribution from feature 2
    380388         f3,                       & !< contribution from feature 1 in natural background of Twomey effect
    381          f4,                       & !< contribution from feature 2 in natural background of Twomey effect
    382          aod_550,                  & !< aerosol optical depth at 550nm
    383          aod_lmd,                  & !< aerosol optical depth at input wavelength
    384          lfactor                     !< factor to compute wavelength dependence of optical properties
     389         f4                          !< contribution from feature 2 in natural background of Twomey effect
    385390    !
    386391    ! ----------
    387392    !
    388     ! initialize input data (by calling setup at first instance)
    389     !
    390 ! "CALL sp_setup" moved to macv2sp
    391 !    IF (.NOT.sp_initialized) CALL sp_setup
     393    ! input data are initialized in macv2sp routine (by calling sp_setup at first instance)
    392394    !
    393395    ! get time weights
    … …  
    395397    CALL set_time_weight(decimal_year)
    396398    !
     399    ! initialize variables
     400    !
     401    DO k=1,nlevels
     402      DO icol=1,ncol
     403        z_beta(icol,k)   = MERGE(1.0, 0.0, z(icol,k) >= oro(icol))
     404        eta(icol,k)      = MAX(0.0,MIN(1.0,z(icol,k)/15000.))
     405      ENDDO
     406    ENDDO
     407    DO icol=1,ncol
     408      dNovrN(icol)   = 1.0
     409      caod_sp(icol)  = 0.0
     410      caod_bg(icol)  = aerosols_SP_aod_bg_gl
     411    ENDDO
     412    !
     413    ! sum contribution from plumes to construct composite profiles of aerosol optical properties
     414    !
     415    DO iplume=1,nplumes
     416      !
     417      ! calculate vertical distribution function from parameters of beta distribution
     418      !
     419      DO icol=1,ncol
     420        beta_sum(icol) = 0.
     421      ENDDO
     422      DO k=1,nlevels
     423        DO icol=1,ncol
     424          prof(icol,k)   = (eta(icol,k)**(beta_a(iplume)-1.) * (1.-eta(icol,k))**(beta_b(iplume)-1.)) * dz(icol,k)
     425          beta_sum(icol) = beta_sum(icol) + prof(icol,k)
     426        ENDDO
     427      ENDDO
     428      DO k=1,nlevels
     429        DO icol=1,ncol
     430          prof(icol,k)   = ( prof(icol,k) / beta_sum(icol) ) * z_beta(icol,k)
     431        ENDDO
     432      ENDDO
     433      !
     434      ! calculate plume weights
     435      !
     436      DO icol=1,ncol
     437        !
     438        ! get plume-center relative spatial parameters for specifying amplitude of plume at given lat and lon
     439        !
     440        delta_lat   = lat(icol) - plume_lat(iplume)
     441        delta_lon   = lon(icol) - plume_lon(iplume)
     442        delta_lon_t = MERGE (260., 180., iplume == 1)
     443        delta_lon   = MERGE ( delta_lon-SIGN(360.,delta_lon) , delta_lon , ABS(delta_lon) > delta_lon_t)
     444
     445        a_plume1  = 0.5 / (MERGE(sig_lon_E(1,iplume), sig_lon_W(1,iplume), delta_lon > 0)**2)
     446        b_plume1  = 0.5 / (MERGE(sig_lat_E(1,iplume), sig_lat_W(1,iplume), delta_lon > 0)**2)
     447        a_plume2  = 0.5 / (MERGE(sig_lon_E(2,iplume), sig_lon_W(2,iplume), delta_lon > 0)**2)
     448        b_plume2  = 0.5 / (MERGE(sig_lat_E(2,iplume), sig_lat_W(2,iplume), delta_lon > 0)**2)
     449        !
     450        ! adjust for a plume specific rotation which helps match plume state to climatology.
     451        !
     452        lon1 =   COS(theta(1,iplume))*(delta_lon) + SIN(theta(1,iplume))*(delta_lat)
     453        lat1 = - SIN(theta(1,iplume))*(delta_lon) + COS(theta(1,iplume))*(delta_lat)
     454        lon2 =   COS(theta(2,iplume))*(delta_lon) + SIN(theta(2,iplume))*(delta_lat)
     455        lat2 = - SIN(theta(2,iplume))*(delta_lon) + COS(theta(2,iplume))*(delta_lat)
     456        !
     457        ! calculate contribution to plume from its different features, to get a column weight for the anthropogenic
     458        ! (cw_an) and the fine-mode natural background aerosol (cw_bg)
     459        !
     460        f1 = time_weight(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2))))
     461        f2 = time_weight(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2))))
     462        f3 = time_weight_bg(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2))))
     463        f4 = time_weight_bg(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2))))
     464
     465        cw_an(icol) = f1 * aod_spmx(iplume) + f2 * aod_spmx(iplume) 
     466        cw_bg(icol) = f3 * aod_fmbg(iplume) + f4 * aod_fmbg(iplume)
     467        !
     468      ENDDO
     469      !
     470      ! distribute plume optical depth at 550 over the vertical profile
     471      !
     472      DO k=1,nlevels
     473        DO icol = 1,ncol
     474          aod550(icol,k,iplume) = prof(icol,k)     * cw_an(icol)
     475          caod_sp(icol)    = caod_sp(icol)    + aod550(icol,k,iplume)
     476          caod_bg(icol)    = caod_bg(icol)    + prof(icol,k) * cw_bg(icol)
     477        ENDDO  ! icol
     478      ENDDO  ! k
     479    ENDDO  ! iplume
     480    !
     481    !
     482    ! calculate effective radius normalization (divisor) factor
     483    ! Stevens et al (2017), eq (15)
     484    DO icol=1,ncol
     485      dNovrN(icol) = LOG((aerosols_SP_coef_bN * (caod_sp(icol) + caod_bg(icol))) + 1.0)/LOG((aerosols_SP_coef_bN * caod_bg(icol)) + 1.0)
     486    ENDDO
     487
     488    RETURN
     489  END SUBROUTINE sp_aod550_profile
     490
     491  ! ------------------------------------------------------------------------------------------------------------------------
     492  ! sp_aop_profile:  This subroutine for Simple Plume aerosols, sums over nplumes to provide profiles of optical properties
     493  !    on a host model's grid for a given wavelength "lambda".
     494  !    It uses aod550(ncol,nlevels,nplumes) output by the subroutine sp_aod550_profile.
     495  ! ------------------------------------------------------------------------------------------------------------------------
     496
     497  SUBROUTINE sp_aop_profile                                                                        ( &
     498       nlevels        ,ncol           ,lambda                     , &
     499       aod550         ,aod_prof       ,ssa_prof       , &
     500       asy_prof       )
     501    !
     502    ! ----------
     503    !
     504    INTEGER, INTENT(IN)        :: &
     505         nlevels,                 & !< number of levels
     506         ncol                       !< number of columns
     507
     508    REAL, INTENT(IN)           :: &
     509         lambda,                  & !< wavelength
     510         aod550(ncol,nlevels,nplumes) !< anthropogenic aod550 profiles by individual plumes
     511
     512    REAL, INTENT(OUT)          :: &
     513         aod_prof(ncol,nlevels) , & !< profile of aerosol optical depth
     514         ssa_prof(ncol,nlevels) , & !< profile of single scattering albedo
     515         asy_prof(ncol,nlevels)     !< profile of asymmetry parameter
     516
     517    INTEGER                    :: iplume, icol, k
     518
     519    REAL                       ::  &
     520         ssa,                      & !< single scattering albedo
     521         asy,                      & !< asymmetry parameter
     522         aod_lmdz,                  & !< aerosol optical depth at input wavelength
     523         lfactor,                  & !< factor to compute wavelength dependence of optical properties
     524         lextinct                    !< anthropogenic aerosol extinction (function of wavelenth and aerosol type/plume)
     525
    397526    ! initialize variables, including output
    398527    !
    … …  
    402531        ssa_prof(icol,k) = 0.0
    403532        asy_prof(icol,k) = 0.0
    404         z_beta(icol,k)   = MERGE(1.0, 0.0, z(icol,k) >= oro(icol))
    405         eta(icol,k)      = MAX(0.0,MIN(1.0,z(icol,k)/15000.))
    406533      ENDDO
    407534    ENDDO
    408     DO icol=1,ncol
    409       dNovrN(icol)   = 1.0
    410       caod_sp(icol)  = 0.0
    411 !      caod_bg(icol)  = 0.02
    412       caod_bg(icol)  = aerosols_SP_aod_bg_gl
    413     ENDDO
    414     !
     535
     536    ! lfactor, eq(12) de Stevens et al 2017 inversed (LAMBDA=max(1,lambda/700. therein)
     537    lfactor   = MIN(1.0,700.0/lambda)
     538
    415539    ! sum contribution from plumes to construct composite profiles of aerosol optical properties
    416540    !
    417541    DO iplume=1,nplumes
    418542      !
    419       ! calculate vertical distribution function from parameters of beta distribution
    420       !
    421       DO icol=1,ncol
    422         beta_sum(icol) = 0.
    423       ENDDO
    424       DO k=1,nlevels
    425         DO icol=1,ncol
    426           prof(icol,k)   = (eta(icol,k)**(beta_a(iplume)-1.) * (1.-eta(icol,k))**(beta_b(iplume)-1.)) * dz(icol,k)
    427           beta_sum(icol) = beta_sum(icol) + prof(icol,k)
    428         ENDDO
    429       ENDDO
    430       DO k=1,nlevels
    431         DO icol=1,ncol
    432           prof(icol,k)   = ( prof(icol,k) / beta_sum(icol) ) * z_beta(icol,k)
    433         ENDDO
    434       ENDDO
    435       !
    436       ! calculate plume weights
    437       !
    438       DO icol=1,ncol
    439         !
    440         ! get plume-center relative spatial parameters for specifying amplitude of plume at given lat and lon
    441         !
    442         delta_lat   = lat(icol) - plume_lat(iplume)
    443         delta_lon   = lon(icol) - plume_lon(iplume)
    444         delta_lon_t = MERGE (260., 180., iplume == 1)
    445         delta_lon   = MERGE ( delta_lon-SIGN(360.,delta_lon) , delta_lon , ABS(delta_lon) > delta_lon_t)
    446 
    447         a_plume1  = 0.5 / (MERGE(sig_lon_E(1,iplume), sig_lon_W(1,iplume), delta_lon > 0)**2)
    448         b_plume1  = 0.5 / (MERGE(sig_lat_E(1,iplume), sig_lat_W(1,iplume), delta_lon > 0)**2)
    449         a_plume2  = 0.5 / (MERGE(sig_lon_E(2,iplume), sig_lon_W(2,iplume), delta_lon > 0)**2)
    450         b_plume2  = 0.5 / (MERGE(sig_lat_E(2,iplume), sig_lat_W(2,iplume), delta_lon > 0)**2)
    451         !
    452         ! adjust for a plume specific rotation which helps match plume state to climatology.
    453         !
    454         lon1 =   COS(theta(1,iplume))*(delta_lon) + SIN(theta(1,iplume))*(delta_lat)
    455         lat1 = - SIN(theta(1,iplume))*(delta_lon) + COS(theta(1,iplume))*(delta_lat)
    456         lon2 =   COS(theta(2,iplume))*(delta_lon) + SIN(theta(2,iplume))*(delta_lat)
    457         lat2 = - SIN(theta(2,iplume))*(delta_lon) + COS(theta(2,iplume))*(delta_lat)
    458         !
    459         ! calculate contribution to plume from its different features, to get a column weight for the anthropogenic
    460         ! (cw_an) and the fine-mode natural background aerosol (cw_bg)
    461         !
    462         f1 = time_weight(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2))))
    463         f2 = time_weight(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2))))
    464         f3 = time_weight_bg(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2))))
    465         f4 = time_weight_bg(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2))))
    466 
    467         cw_an(icol) = f1 * aod_spmx(iplume) + f2 * aod_spmx(iplume) 
    468         cw_bg(icol) = f3 * aod_fmbg(iplume) + f4 * aod_fmbg(iplume)
    469         !
    470         ! calculate wavelength-dependent scattering properties
    471         !
    472         lfactor   = MIN(1.0,700.0/lambda)
    473         ssa(icol) = (ssa550(iplume) * lfactor**4) / ((ssa550(iplume) * lfactor**4) + ((1-ssa550(iplume)) * lfactor))
    474         asy(icol) =  asy550(iplume) * SQRT(lfactor)
    475       ENDDO
     543      ! calculate wavelength-dependent scattering properties
     544      !  Stevens et al 2017, eqs (11)&(13)
     545!    NOTES ASima : ssa and asy only depend on iplume (via ssa550, asy550) and lfactor(lambda), don't need icol dimension.
     546!      Also, why eq(11) was written in this more complex form than in the paper ?
     547      !ssa(icol) = (ssa550(iplume) * lfactor**4) / ((ssa550(iplume) * lfactor**4) + ((1-ssa550(iplume)) * lfactor))
     548      ssa = (ssa550(iplume) * lfactor**4) / ((ssa550(iplume) * lfactor**4) + ((1-ssa550(iplume)) * lfactor))
     549      !ssa = ssa550(iplume) / (ssa550(iplume) + (1-ssa550(iplume)) / lfactor**3)
     550      asy = asy550(iplume) * SQRT(lfactor)
    476551      !
    477552      ! distribute plume optical properties across its vertical profile weighting by optical depth and scaling for
    478553      ! wavelength using the angstrom parameter.
    479       !     
    480       lfactor = EXP(-angstrom(iplume) * LOG(lambda/550.0))
     554      !
     555      !   lextinct = factor for aerosol extiction, eq(10) in Stevens et al 2017
     556      !    Depends on 'iplume' via 'angstrom' (Note : angstrom=2. is prescribed for all plumes)
     557      lextinct = EXP(-angstrom(iplume) * LOG(lambda/550.0))
     558
     559! NOTE : lextinct, ssa, asy ne dependent ni de icol, ni de k ; peut-on ptimiser ?
    481560      DO k=1,nlevels
    482561        DO icol = 1,ncol
    483           aod_550          = prof(icol,k)     * cw_an(icol)
    484           aod_lmd          = aod_550          * lfactor
    485           caod_sp(icol)    = caod_sp(icol)    + aod_550
    486           caod_bg(icol)    = caod_bg(icol)    + prof(icol,k) * cw_bg(icol)
    487           asy_prof(icol,k) = asy_prof(icol,k) + aod_lmd * ssa(icol) * asy(icol)
    488           ssa_prof(icol,k) = ssa_prof(icol,k) + aod_lmd * ssa(icol)
    489           aod_prof(icol,k) = aod_prof(icol,k) + aod_lmd
    490         ENDDO
    491       ENDDO
    492     ENDDO
     562           aod_lmdz          = aod550(icol,k,iplume)     * lextinct
     563           asy_prof(icol,k) = asy_prof(icol,k) + aod_lmdz * ssa  * asy
     564           ssa_prof(icol,k) = ssa_prof(icol,k) + aod_lmdz * ssa
     565           aod_prof(icol,k) = aod_prof(icol,k) + aod_lmdz
     566        ENDDO ! k (levels)
     567      ENDDO ! icol
     568    ENDDO ! iplume
    493569    !
    494570    ! complete optical depth weighting
    … …  
    500576      ENDDO
    501577    ENDDO
    502     !
    503     ! calculate effective radius normalization (divisor) factor
    504     ! Stevens et al (2017), eq (15)
    505     DO icol=1,ncol
    506 !      dNovrN(icol) = LOG((1000.0 * (caod_sp(icol) + caod_bg(icol))) + 1.0)/LOG((1000.0 * caod_bg(icol)) + 1.0)
    507       dNovrN(icol) = LOG((aerosols_SP_coef_bN * (caod_sp(icol) + caod_bg(icol))) + 1.0)/LOG((aerosols_SP_coef_bN * caod_bg(icol)) + 1.0)
    508     ENDDO
     578
    509579
    510580    RETURN
    511581  END SUBROUTINE sp_aop_profile
    512   
     582 
    513583END MODULE mo_simple_plumes
Note: See TracChangeset for help on using the changeset viewer.