SUBROUTINE macv2sp(pphis,pplay,paprs,xlon,xlat,tau_allaer,piz_allaer,cg_allaer,dNovrN) ! !--routine to read the MACv2SP plume and compute optical properties !--requires flag_aerosol = 7 !--feeds into aerosol optical properties and newmicro cloud droplet size if ok_cdnc activated !--for this one needs to feed natural (pre-industrial) aerosols twice for nat and 1980 files !--pre-ind aerosols (index=1) are not changed, present-day aerosols (index=2) are incremented !--uses model year so year_cur needs to be correct in the model simulation ! !--aod_prof = AOD per layer !--ssa_prof = SSA !--asy_prof = asymetry parameter !--dNovrN = enhancement factor for CDNC ! !ASima : *_allaer arguments and non-SP "aer" diagnostics arrive from rrtm/aeropt_5wv_rrtm.F90 with values for inca1850 aerosols USE mo_simple_plumes, ONLY: sp_aod550_profile, sp_aop_profile, sp_aod_diag, 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 ! SP-specific diagnostics (an SP extinction coefficient could be also defined: ec550SPaer = od550SPaer_prof/dz) USE phys_local_var_mod, ONLY: od550SPaer, od550SPaer_prof ! USE yomcst_mod_h IMPLICIT NONE ! ! REAL,DIMENSION(klon),INTENT(IN) :: pphis ! Geopotentiel de surface REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour les interfaces de chaque couche (en Pa) REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point ! REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(INOUT) :: tau_allaer ! epaisseur optique aerosol REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(INOUT) :: piz_allaer ! single scattering albedo aerosol 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,nplumes) :: aod550_pl REAL,DIMENSION(klon,klev) :: aod_prof, ssa_prof, asy_prof REAL,DIMENSION(klon) :: aod_diag REAL,DIMENSION(klon,klev) :: z, dz REAL,DIMENSION(klon) :: oro, zrho, zt ! REAL, PARAMETER :: l443 = 443.0, l865 = 865.0 !--wavelengths in nm for diagnostics (other than 550 nm calculated by default) ! INTEGER, PARAMETER :: Nwvmax=25 REAL, DIMENSION(0:Nwvmax), PARAMETER :: lambda=(/ 240.0, & !--this one is for band 1 280.0, 300.0, 330.0, 360.0, 400.0, & !--these are bounds of Streamer bands 440.0, 480.0, 520.0, 570.0, 640.0, & 690.0, 750.0, 780.0, 870.0, 1000.0, & 1100.0, 1190.0, 1280.0, 1530.0, 1640.0, & 2130.0, 2380.0, 2910.0, 3420.0, 4000.0 /) ! REAL, DIMENSION(1:Nwvmax-1), PARAMETER :: weight =(/ & !--and the weights to be given to the bands 0.01, 4.05, 9.51, 15.99, 26.07, 33.10, & !--corresponding to a typical solar spectrum 33.07, 39.91, 52.67, 27.89, 43.60, 13.67, & 42.22, 40.12, 32.70, 14.44, 19.48, 14.23, & 13.43, 16.42, 8.33, 0.95, 0.65, 2.76 /) ! REAL :: zlambda, zweight REAL :: decimal_year ! INTEGER band, k, Nwv ! ! define the height and dheight arrays ! oro(:) = pphis(:)/RG ! surface height in m ! ASima : the following calculations of z and dz are also done in phys_output_write_mod.F90 ! for wind100m, under "IF (vars_defined)" ; "z" therein does not include "oro" ! Layers' top and middle altitudes zhalf et zfull are also calculated for diagnostics ! in phys_output_write_mod.F90, under "IF (vars_defined)", using the geopotential "phi" ... DO k = 1, klev zrho(:) = pplay(:,k)/t_seri(:,k)/RD ! air density in kg/m3 dz(:,k) = (paprs(:,k)-paprs(:,k+1))/zrho(:)/RG ! layer thickness in m IF (k==1) THEN z(:,1) = oro(:) + (paprs(:,1)-pplay(:,1))/zrho(:)/RG ! altitude middle of first layer in m zt(:) = oro(:) + dz(:,1) ! altitude top of first layer in m ELSE z(:,k) = zt(:) + (paprs(:,k)-pplay(:,k))/zrho(:)/RG ! altitude middle of layer k in m zt(:) = zt(:) + dz(:,k) ! altitude top of layer k in m ENDIF ENDDO ! ! Call sp_setup to retrieve 'aerosol_SP_forcing_year' value used in decimal_year calculation below 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 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) ELSE decimal_year= float(aerosols_SP_forcing_year) + float(days_elapsed)/float(year_len) ENDIF IF (decimal_year.LT.1850.0.OR.decimal_year.GE.2024.0) THEN CALL abort_physic ('macv2sp','year not supported by plume model',1) ENDIF ! !--call SUBROUTINE sp_aod550_profile once ; all the wavelength-dependent profiles and 2D diagnostics will use aod550 profile ! CALL sp_aod550_profile ( & klev ,klon ,xlon ,xlat , & decimal_year ,z ,dz ,dNovrN , aod550 ) !--AOD550 calculations for diagnostics ! ! AOD profile = sum over the 'nplumes' dimension od550SPaer_prof(:,:)=SUM(aod550(:,:,:),dim=3) ! !--a/ AOD of SP anthropic aerosols = vertical sum of SP aod profile od550SPaer(:)=SUM(od550SPaer_prof(:,:),dim=2) ! !--b/ total AOD = Inca1850 + SP od550aer(:)=od550aer(:)+od550SPaer(:) ! !--c/ fine-mode AOD = Inca1850(fine mode) + SP od550lt1aer(:)=od550lt1aer(:)+od550SPaer(:) ! !--d/ dry AOD dryod550aer(:)=dryod550aer(:)+od550SPaer(:) ! !--extinction coefficient for diagnostic ec550aer(:,:)=ec550aer(:,:)+od550SPaer_prof(:,:)/dz(:,:) ! ! ! AOD calculations for diagnostics at other wavelengths : 443 and 865 nm ! using AOD550(icol,iplume) : aod550_pl(:,:)=SUM(aod550(:,:,:),dim=2) ! !--call to sp routine -- 443 nm ! CALL sp_aod_diag ( & klon , l443 , & aod550_pl ,aod_diag ) od443aer(:)= od443aer(:)+aod_diag(:) ! !--call to sp routine -- 865 nm ! CALL sp_aod_diag ( & klon , l865 , & aod550_pl ,aod_diag ) od865aer(:)= od865aer(:)+aod_diag(:) ! !--re-weighting of piz and cg arrays before adding the anthropogenic aerosols !--index 2 = all natural + anthropogenic aerosols piz_allaer(:,:,2,:)=piz_allaer(:,:,2,:)*tau_allaer(:,:,2,:) cg_allaer(:,:,2,:) =cg_allaer(:,:,2,:)*piz_allaer(:,:,2,:) ! !--now computing the same at many wavelengths to fill the model bands ! DO Nwv=0,Nwvmax-1 IF (Nwv.EQ.0) THEN !--RRTM spectral band 1 zlambda=lambda(Nwv) zweight=1.0 band=1 ELSEIF (Nwv.LE.5) THEN !--RRTM spectral band 2 zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1)) zweight=weight(Nwv)/SUM(weight(1:5)) band=2 ELSEIF (Nwv.LE.10) THEN !--RRTM spectral band 3 zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1)) zweight=weight(Nwv)/SUM(weight(6:10)) band=3 ELSEIF (Nwv.LE.16) THEN !--RRTM spectral band 4 zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1)) zweight=weight(Nwv)/SUM(weight(11:16)) band=4 ELSEIF (Nwv.LE.21) THEN !--RRTM spectral band 5 zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1)) zweight=weight(Nwv)/SUM(weight(17:21)) band=5 ELSE !--RRTM spectral band 6 zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1)) zweight=weight(Nwv)/SUM(weight(22:Nwvmax-1)) band=6 ENDIF ! 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 tau_allaer(:,:,2,band)=tau_allaer(:,:,2,band)+zweight*MAX(aod_prof(:,:),1.e-15) piz_allaer(:,:,2,band)=piz_allaer(:,:,2,band)+zweight*MAX(aod_prof(:,:),1.e-15)*ssa_prof(:,:) cg_allaer(:,:,2,band) =cg_allaer(:,:,2,band) +zweight*MAX(aod_prof(:,:),1.e-15)*ssa_prof(:,:)*asy_prof(:,:) ! ENDDO ! !--renormalizing cg and piz now that MACv2SP increments have been added cg_allaer(:,:,2,:) =cg_allaer(:,:,2,:) /piz_allaer(:,:,2,:) piz_allaer(:,:,2,:)=piz_allaer(:,:,2,:)/tau_allaer(:,:,2,:) ! END SUBROUTINE macv2sp