Index: /LMDZ6/trunk/libf/phylmd/macv2sp.f90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/macv2sp.f90 (revision 6127) +++ /LMDZ6/trunk/libf/phylmd/macv2sp.f90 (revision 6128) @@ -1,3 +1,3 @@ -SUBROUTINE MACv2SP(pphis,pplay,paprs,xlon,xlat,tau_allaer,piz_allaer,cg_allaer) +SUBROUTINE macv2sp(pphis,pplay,paprs,xlon,xlat,tau_allaer,piz_allaer,cg_allaer,dNovrN) ! !--routine to read the MACv2SP plume and compute optical properties @@ -13,10 +13,10 @@ !--dNovrN = enhancement factor for CDNC ! - USE mo_simple_plumes, ONLY: sp_aop_profile, sp_setup, sp_initialized, aerosols_SP_forcing_year + USE mo_simple_plumes, ONLY: sp_aod550_profile, sp_aop_profile, sp_setup ! subroutines + USE mo_simple_plumes, ONLY: nplumes, sp_initialized, aerosols_SP_forcing_year ! variables USE phys_cal_mod, ONLY : year_cur, year_len, days_elapsed USE dimphy USE aero_mod - USE phys_local_var_mod, ONLY: t_seri, od443aer, od550aer, od865aer, ec550aer, dryod550aer, od550lt1aer, od550SPaer, dNovrN - !!USE YOMCST, ONLY : RD, RG + USE phys_local_var_mod, ONLY: t_seri, od443aer, od550aer, od865aer, ec550aer, dryod550aer, od550lt1aer, od550SPaer ! USE yomcst_mod_h @@ -35,11 +35,12 @@ REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(INOUT) :: cg_allaer ! asymmetry parameter aerosol ! + REAL,DIMENSION(klon),INTENT(OUT) :: dNovrN !< anthropogenic increase in cloud drop number concentration (factor) + ! + REAL,DIMENSION(klon,klev,nplumes) :: aod550 REAL,DIMENSION(klon,klev) :: aod_prof, ssa_prof, asy_prof REAL,DIMENSION(klon,klev) :: z, dz REAL,DIMENSION(klon) :: oro, zrho, zt ! - INTEGER, PARAMETER :: nmon = 12 - ! - REAL, PARAMETER :: l443 = 443.0, l550 = 550.0, l865 = 865.0 !--wavelengths in nm + REAL, PARAMETER :: l443 = 443.0, l865 = 865.0 !--wavelengths in nm for diagnostics (other than 550 nm calculated by default) ! INTEGER, PARAMETER :: Nwvmax=25 @@ -58,5 +59,4 @@ ! REAL :: zlambda, zweight -! REAL :: year_fr ! also a dimension name in SP aerosol file ; renamed more properly "decimal_year" REAL :: decimal_year ! @@ -83,11 +83,8 @@ IF (.NOT.sp_initialized) CALL sp_setup - !--fractional year - ! - ! Original version, bugged : - ! year_fr = FLOAT(year_cur) + (FLOAT(day_cur)-0.5) / FLOAT(year_len) - ! Correction ASima & FH : -! year_fr= float(year_cur) + float(days_elapsed)/float(year_len) -! Choice between yearly vs fixed_year forcing, depending on 'aerosols_SP_forcing_year' + !--fractional year : + ! Original version, bugged : year_fr = FLOAT(year_cur) + (FLOAT(day_cur)-0.5) / FLOAT(year_len) + ! Correction 2026-03, FH & ASima : year_fr= float(year_cur) + float(days_elapsed)/float(year_len) + ! Name changed in decimal_year ; choice between yearly vs fixed_year forcing, depending on 'aerosols_SP_forcing_year' IF (aerosols_SP_forcing_year.EQ.-9999) THEN decimal_year= float(year_cur) + float(days_elapsed)/float(year_len) @@ -99,25 +96,18 @@ CALL abort_physic ('macv2sp','year not supported by plume model',1) ENDIF - ! - !--call to sp routine -- 443 nm - ! - CALL sp_aop_profile ( & - klev ,klon ,l443 ,oro ,xlon ,xlat , & - decimal_year ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , & - asy_prof ) - ! - !--AOD calculations for diagnostics - od443aer(:)= od443aer(:)+SUM(aod_prof(:,:),dim=2) - ! - !--call to sp routine -- 550 nm - ! - CALL sp_aop_profile ( & - klev ,klon ,l550 ,oro ,xlon ,xlat , & - decimal_year ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , & - asy_prof ) - ! - !--AOD calculations for diagnostics - ! (ASima : new diagnostic od550SPaer; corrected od550lt1aer and dryod550aer) - !--a/ AOD of SP aerosols = vertical sum of SP aod profile + + ! + !--call SUBROUTINE sp_aod550_profile once ; all the wavelength-dependent profiles and 2D diagnostics will use aod550 profile + ! + CALL sp_aod550_profile ( & + klev ,klon ,oro ,xlon ,xlat , & + decimal_year ,z ,dz ,dNovrN , aod550 ) + + !--AOD550 calculations for diagnostics + ! + ! aod550 profile = sum over the 'nplumes' dimension + aod_prof(:,:)=SUM(aod550(:,:,:),dim=3) + ! + !--a/ AOD of SP anthropic aerosols = vertical sum of SP aod profile od550SPaer(:)=SUM(aod_prof(:,:),dim=2) ! @@ -126,26 +116,32 @@ ! !--c/ fine-mode AOD = Inca1850(fine mode) + SP - ! original, bugged : includes od550aer of Inca1850 - !od550lt1aer(:)=od550lt1aer(:)+od550aer(:) od550lt1aer(:)=od550lt1aer(:)+od550SPaer(:) ! !--d/ dry AOD - ! original, bugged : includes od550aer of Inca1850 - !dryod550aer(:)=dryod550aer(:)+od550aer(:) dryod550aer(:)=dryod550aer(:)+od550SPaer(:) - ! ! !--extinction coefficient for diagnostic ec550aer(:,:)=ec550aer(:,:)+aod_prof(:,:)/dz(:,:) ! + + ! Diagnostic AOD calculations for other wavelengths : 443 and 865 nm + !--call to sp routine -- 443 nm + ! + CALL sp_aop_profile ( & + klev ,klon , l443 , & + aod550 ,aod_prof ,ssa_prof , & + asy_prof ) + + od443aer(:)= od443aer(:)+SUM(aod_prof(:,:),dim=2) + !--call to sp routine -- 865 nm ! CALL sp_aop_profile ( & - klev ,klon ,l865 ,oro ,xlon ,xlat , & - decimal_year ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , & + klev ,klon , l865 , & + aod550 ,aod_prof ,ssa_prof , & asy_prof ) - ! - !--AOD calculations for diagnostics + od865aer(:)=od865aer(:)+SUM(aod_prof(:,:),dim=2) + ! !--re-weighting of piz and cg arrays before adding the anthropogenic aerosols @@ -184,8 +180,10 @@ ENDIF ! - CALL sp_aop_profile ( & - klev ,klon ,zlambda ,oro ,xlon ,xlat , & - decimal_year ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , & - asy_prof ) + CALL sp_aop_profile ( & + klev ,klon , zlambda , & + aod550 ,aod_prof ,ssa_prof , & + asy_prof ) + + ! !--adding up the quantities tau, piz*tau and cg*piz*tau @@ -200,3 +198,3 @@ piz_allaer(:,:,2,:)=piz_allaer(:,:,2,:)/tau_allaer(:,:,2,:) ! -END SUBROUTINE MACv2SP +END SUBROUTINE macv2sp Index: /LMDZ6/trunk/libf/phylmd/mo_simple_plumes.f90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/mo_simple_plumes.f90 (revision 6127) +++ /LMDZ6/trunk/libf/phylmd/mo_simple_plumes.f90 (revision 6128) @@ -20,5 +20,12 @@ !! !! @par Copyright -!! +!! +!! 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 ; +!! sp_aop_profile : uses optical properties (aod, ssa, asy) at 550 to calculate their profiles for another wavelength ; +!! at each timestep is called in macv2sp for every required wavelength separately : +!! - aod diagnostics at 443 and 865um, +!! - optical properties (aod, ssa, asy) for the Nwvmax=25 wavelengths filling the 6 RRTM bands ! MODULE mo_simple_plumes @@ -28,6 +35,9 @@ IMPLICIT NONE + INTEGER, PARAMETER, PUBLIC :: nplumes = 9 !< Number of plumes + +! INTEGER, PARAMETER :: & +! nplumes = 9 ,& !< Number of plumes INTEGER, PARAMETER :: & - nplumes = 9 ,& !< Number of plumes nfeatures = 2 ,& !< Number of features per plume ntimes = 52 ,& !< Number of times resolved per year (52 => weekly resolution) @@ -67,5 +77,5 @@ time_weight_bg (nfeatures,nplumes) !< as time_weight but for natural background in Twomey effect - PUBLIC sp_aop_profile, sp_setup + PUBLIC sp_aod550_profile, sp_aop_profile, sp_setup CONTAINS @@ -94,5 +104,7 @@ CALL getin_p('aerosols_SP_coef_bN', aerosols_SP_coef_bN) CALL getin_p('aerosols_SP_aod_bg_gl', aerosols_SP_aod_bg_gl) - !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 + print *, 'aerosols_SP_forcing_year=',aerosols_SP_forcing_year + print *, 'aerosols_SP_coef_bN = ', aerosols_SP_coef_bN + print *, 'aerosols_SP_aod_bg_gl = ', aerosols_SP_aod_bg_gl !--only one processor reads the input data @@ -295,8 +307,8 @@ ! REAL, INTENT(IN) :: & - decimal_year !< Fractional Year (1850.0 - 2100.99) + decimal_year !< Decimal Year (1850.0 - 2100.99) INTEGER :: & - iyear ,& !< Integer year values between 1 and 156 (1850-2100) + iyear ,& !< Integer year values between 1 and 251 (1850-2100) ; in 2026 : data available for 1850-2023 (NaN after) iweek ,& !< Integer index (between 1 and ntimes); for ntimes=52 this corresponds to weeks (roughly) iplume ! plume number @@ -321,13 +333,14 @@ END SUBROUTINE set_time_weight ! - ! ------------------------------------------------------------------------------------------------------------------------ - ! SP_AOP_PROFILE: This subroutine calculates the simple plume aerosol and cloud active optical properties based on the - ! the simple plume fit to the MPI Aerosol Climatology (Version 2). It sums over nplumes to provide a profile of aerosol - ! optical properties on a host models vertical grid. + !------------------------------------------------------------------------------------------------------------------------ + ! sp_aod550_profile : This subroutine calculates Simple Plume aerosol aod550(nm) profile from MACv2SP (550um) data + ! (Stevens et al 2027 ; designed to fit the MPI Aerosol Climatology (Version 2) by Kinne, 2018), + ! It sums over nplumes to provide aod550 profile on a host model's vertical grid. + ! It also computes the dNovrN factor used to mimic Twomey (aci) effect by multiplying preind. cloud droplet number ! - SUBROUTINE sp_aop_profile ( & - nlevels ,ncol ,lambda ,oro ,lon ,lat , & - decimal_year ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , & - asy_prof ) + SUBROUTINE sp_aod550_profile ( & + nlevels ,ncol ,oro ,lon ,lat , & + decimal_year ,z ,dz ,dNovrN , aod550 ) + ! ! ---------- @@ -338,5 +351,4 @@ REAL, INTENT(IN) :: & - lambda, & !< wavelength decimal_year, & !< Fractional Year (1903.0 is the 0Z on the first of January 1903, Gregorian) oro(ncol), & !< orographic height (m) @@ -348,7 +360,5 @@ REAL, INTENT(OUT) :: & dNovrN(ncol) , & !< anthropogenic increase in cloud drop number concentration (factor) - aod_prof(ncol,nlevels) , & !< profile of aerosol optical depth - ssa_prof(ncol,nlevels) , & !< profile of single scattering albedo - asy_prof(ncol,nlevels) !< profile of asymmetry parameter + aod550(ncol,nlevels,nplumes) !< anthropogenic aod550 profiles by individual plumes INTEGER :: iplume, icol, k @@ -359,6 +369,4 @@ prof(ncol,nlevels), & !< scaled profile (by beta function) beta_sum(ncol), & !< vertical sum of beta function - ssa(ncol), & !< single scattering albedo - asy(ncol), & !< asymmetry parameter 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 @@ -379,15 +387,9 @@ f2, & !< contribution from feature 2 f3, & !< contribution from feature 1 in natural background of Twomey effect - f4, & !< contribution from feature 2 in natural background of Twomey effect - aod_550, & !< aerosol optical depth at 550nm - aod_lmd, & !< aerosol optical depth at input wavelength - lfactor !< factor to compute wavelength dependence of optical properties + f4 !< contribution from feature 2 in natural background of Twomey effect ! ! ---------- ! - ! initialize input data (by calling setup at first instance) - ! -! "CALL sp_setup" moved to macv2sp -! IF (.NOT.sp_initialized) CALL sp_setup + ! input data are initialized in macv2sp routine (by calling sp_setup at first instance) ! ! get time weights @@ -395,4 +397,131 @@ CALL set_time_weight(decimal_year) ! + ! initialize variables + ! + 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.)) + ENDDO + ENDDO + DO icol=1,ncol + dNovrN(icol) = 1.0 + caod_sp(icol) = 0.0 + caod_bg(icol) = aerosols_SP_aod_bg_gl + ENDDO + ! + ! sum contribution from plumes to construct composite profiles of aerosol optical properties + ! + DO iplume=1,nplumes + ! + ! 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) + ENDDO + ENDDO + ! + ! calculate plume weights + ! + DO icol=1,ncol + ! + ! get plume-center relative spatial parameters for specifying amplitude of plume at given lat and lon + ! + delta_lat = lat(icol) - plume_lat(iplume) + delta_lon = lon(icol) - plume_lon(iplume) + delta_lon_t = MERGE (260., 180., iplume == 1) + delta_lon = MERGE ( delta_lon-SIGN(360.,delta_lon) , delta_lon , ABS(delta_lon) > delta_lon_t) + + a_plume1 = 0.5 / (MERGE(sig_lon_E(1,iplume), sig_lon_W(1,iplume), delta_lon > 0)**2) + b_plume1 = 0.5 / (MERGE(sig_lat_E(1,iplume), sig_lat_W(1,iplume), delta_lon > 0)**2) + a_plume2 = 0.5 / (MERGE(sig_lon_E(2,iplume), sig_lon_W(2,iplume), delta_lon > 0)**2) + b_plume2 = 0.5 / (MERGE(sig_lat_E(2,iplume), sig_lat_W(2,iplume), delta_lon > 0)**2) + ! + ! adjust for a plume specific rotation which helps match plume state to climatology. + ! + lon1 = COS(theta(1,iplume))*(delta_lon) + SIN(theta(1,iplume))*(delta_lat) + lat1 = - SIN(theta(1,iplume))*(delta_lon) + COS(theta(1,iplume))*(delta_lat) + lon2 = COS(theta(2,iplume))*(delta_lon) + SIN(theta(2,iplume))*(delta_lat) + lat2 = - SIN(theta(2,iplume))*(delta_lon) + COS(theta(2,iplume))*(delta_lat) + ! + ! calculate contribution to plume from its different features, to get a column weight for the anthropogenic + ! (cw_an) and the fine-mode natural background aerosol (cw_bg) + ! + f1 = time_weight(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2)))) + f2 = time_weight(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2)))) + f3 = time_weight_bg(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2)))) + f4 = time_weight_bg(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2)))) + + cw_an(icol) = f1 * aod_spmx(iplume) + f2 * aod_spmx(iplume) + cw_bg(icol) = f3 * aod_fmbg(iplume) + f4 * aod_fmbg(iplume) + ! + ENDDO + ! + ! distribute plume optical depth at 550 over the vertical profile + ! + DO k=1,nlevels + DO icol = 1,ncol + aod550(icol,k,iplume) = prof(icol,k) * cw_an(icol) + caod_sp(icol) = caod_sp(icol) + aod550(icol,k,iplume) + caod_bg(icol) = caod_bg(icol) + prof(icol,k) * cw_bg(icol) + ENDDO ! icol + ENDDO ! k + ENDDO ! iplume + ! + ! + ! calculate effective radius normalization (divisor) factor + ! Stevens et al (2017), eq (15) + DO icol=1,ncol + 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) + ENDDO + + RETURN + END SUBROUTINE sp_aod550_profile + + ! ------------------------------------------------------------------------------------------------------------------------ + ! sp_aop_profile: This subroutine for Simple Plume aerosols, sums over nplumes to provide profiles of optical properties + ! on a host model's grid for a given wavelength "lambda". + ! It uses aod550(ncol,nlevels,nplumes) output by the subroutine sp_aod550_profile. + ! ------------------------------------------------------------------------------------------------------------------------ + + SUBROUTINE sp_aop_profile ( & + nlevels ,ncol ,lambda , & + aod550 ,aod_prof ,ssa_prof , & + asy_prof ) + ! + ! ---------- + ! + INTEGER, INTENT(IN) :: & + nlevels, & !< number of levels + ncol !< number of columns + + REAL, INTENT(IN) :: & + lambda, & !< wavelength + aod550(ncol,nlevels,nplumes) !< anthropogenic aod550 profiles by individual plumes + + REAL, INTENT(OUT) :: & + aod_prof(ncol,nlevels) , & !< profile of aerosol optical depth + ssa_prof(ncol,nlevels) , & !< profile of single scattering albedo + asy_prof(ncol,nlevels) !< profile of asymmetry parameter + + INTEGER :: iplume, icol, k + + REAL :: & + ssa, & !< single scattering albedo + asy, & !< asymmetry parameter + aod_lmdz, & !< aerosol optical depth at input wavelength + lfactor, & !< factor to compute wavelength dependence of optical properties + lextinct !< anthropogenic aerosol extinction (function of wavelenth and aerosol type/plume) + ! initialize variables, including output ! @@ -402,93 +531,40 @@ ssa_prof(icol,k) = 0.0 asy_prof(icol,k) = 0.0 - 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.)) ENDDO ENDDO - DO icol=1,ncol - dNovrN(icol) = 1.0 - caod_sp(icol) = 0.0 -! caod_bg(icol) = 0.02 - caod_bg(icol) = aerosols_SP_aod_bg_gl - ENDDO - ! + + ! lfactor, eq(12) de Stevens et al 2017 inversed (LAMBDA=max(1,lambda/700. therein) + lfactor = MIN(1.0,700.0/lambda) + ! sum contribution from plumes to construct composite profiles of aerosol optical properties ! DO iplume=1,nplumes ! - ! 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) - ENDDO - ENDDO - ! - ! calculate plume weights - ! - DO icol=1,ncol - ! - ! get plume-center relative spatial parameters for specifying amplitude of plume at given lat and lon - ! - delta_lat = lat(icol) - plume_lat(iplume) - delta_lon = lon(icol) - plume_lon(iplume) - delta_lon_t = MERGE (260., 180., iplume == 1) - delta_lon = MERGE ( delta_lon-SIGN(360.,delta_lon) , delta_lon , ABS(delta_lon) > delta_lon_t) - - a_plume1 = 0.5 / (MERGE(sig_lon_E(1,iplume), sig_lon_W(1,iplume), delta_lon > 0)**2) - b_plume1 = 0.5 / (MERGE(sig_lat_E(1,iplume), sig_lat_W(1,iplume), delta_lon > 0)**2) - a_plume2 = 0.5 / (MERGE(sig_lon_E(2,iplume), sig_lon_W(2,iplume), delta_lon > 0)**2) - b_plume2 = 0.5 / (MERGE(sig_lat_E(2,iplume), sig_lat_W(2,iplume), delta_lon > 0)**2) - ! - ! adjust for a plume specific rotation which helps match plume state to climatology. - ! - lon1 = COS(theta(1,iplume))*(delta_lon) + SIN(theta(1,iplume))*(delta_lat) - lat1 = - SIN(theta(1,iplume))*(delta_lon) + COS(theta(1,iplume))*(delta_lat) - lon2 = COS(theta(2,iplume))*(delta_lon) + SIN(theta(2,iplume))*(delta_lat) - lat2 = - SIN(theta(2,iplume))*(delta_lon) + COS(theta(2,iplume))*(delta_lat) - ! - ! calculate contribution to plume from its different features, to get a column weight for the anthropogenic - ! (cw_an) and the fine-mode natural background aerosol (cw_bg) - ! - f1 = time_weight(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2)))) - f2 = time_weight(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2)))) - f3 = time_weight_bg(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2)))) - f4 = time_weight_bg(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2)))) - - cw_an(icol) = f1 * aod_spmx(iplume) + f2 * aod_spmx(iplume) - cw_bg(icol) = f3 * aod_fmbg(iplume) + f4 * aod_fmbg(iplume) - ! - ! calculate wavelength-dependent scattering properties - ! - lfactor = MIN(1.0,700.0/lambda) - ssa(icol) = (ssa550(iplume) * lfactor**4) / ((ssa550(iplume) * lfactor**4) + ((1-ssa550(iplume)) * lfactor)) - asy(icol) = asy550(iplume) * SQRT(lfactor) - ENDDO + ! calculate wavelength-dependent scattering properties + ! Stevens et al 2017, eqs (11)&(13) +! NOTES ASima : ssa and asy only depend on iplume (via ssa550, asy550) and lfactor(lambda), don't need icol dimension. +! Also, why eq(11) was written in this more complex form than in the paper ? + !ssa(icol) = (ssa550(iplume) * lfactor**4) / ((ssa550(iplume) * lfactor**4) + ((1-ssa550(iplume)) * lfactor)) + ssa = (ssa550(iplume) * lfactor**4) / ((ssa550(iplume) * lfactor**4) + ((1-ssa550(iplume)) * lfactor)) + !ssa = ssa550(iplume) / (ssa550(iplume) + (1-ssa550(iplume)) / lfactor**3) + asy = asy550(iplume) * SQRT(lfactor) ! ! distribute plume optical properties across its vertical profile weighting by optical depth and scaling for ! wavelength using the angstrom parameter. - ! - lfactor = EXP(-angstrom(iplume) * LOG(lambda/550.0)) + ! + ! lextinct = factor for aerosol extiction, eq(10) in Stevens et al 2017 + ! Depends on 'iplume' via 'angstrom' (Note : angstrom=2. is prescribed for all plumes) + lextinct = EXP(-angstrom(iplume) * LOG(lambda/550.0)) + +! NOTE : lextinct, ssa, asy ne dependent ni de icol, ni de k ; peut-on ptimiser ? DO k=1,nlevels DO icol = 1,ncol - aod_550 = prof(icol,k) * cw_an(icol) - aod_lmd = aod_550 * lfactor - caod_sp(icol) = caod_sp(icol) + aod_550 - caod_bg(icol) = caod_bg(icol) + prof(icol,k) * cw_bg(icol) - asy_prof(icol,k) = asy_prof(icol,k) + aod_lmd * ssa(icol) * asy(icol) - ssa_prof(icol,k) = ssa_prof(icol,k) + aod_lmd * ssa(icol) - aod_prof(icol,k) = aod_prof(icol,k) + aod_lmd - ENDDO - ENDDO - ENDDO + aod_lmdz = aod550(icol,k,iplume) * lextinct + asy_prof(icol,k) = asy_prof(icol,k) + aod_lmdz * ssa * asy + ssa_prof(icol,k) = ssa_prof(icol,k) + aod_lmdz * ssa + aod_prof(icol,k) = aod_prof(icol,k) + aod_lmdz + ENDDO ! k (levels) + ENDDO ! icol + ENDDO ! iplume ! ! complete optical depth weighting @@ -500,14 +576,8 @@ ENDDO ENDDO - ! - ! calculate effective radius normalization (divisor) factor - ! Stevens et al (2017), eq (15) - DO icol=1,ncol -! dNovrN(icol) = LOG((1000.0 * (caod_sp(icol) + caod_bg(icol))) + 1.0)/LOG((1000.0 * caod_bg(icol)) + 1.0) - 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) - ENDDO + RETURN END SUBROUTINE sp_aop_profile - + END MODULE mo_simple_plumes Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 6127) +++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 6128) @@ -4105,6 +4105,6 @@ IF (flag_aerosol .EQ. 7) THEN - CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg, & - tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm) + CALL macv2sp(pphis,pplay,paprs,longitude_deg,latitude_deg, & + tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm,dNovrN) ENDIF Index: /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 (revision 6127) +++ /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 (revision 6128) @@ -5644,6 +5644,6 @@ IF (flag_aerosol .EQ. 7) THEN - CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg, & - tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm) + CALL macv2sp(pphis,pplay,paprs,longitude_deg,latitude_deg, & + tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm,dNovrN) ENDIF