Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/conf_phys_m.F90 =================================================================== --- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/conf_phys_m.F90 (revision 3317) +++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/conf_phys_m.F90 (revision 3318) @@ -414,4 +414,5 @@ ! - flag_aerosol=5 => dust only ! - flag_aerosol=6 => all aerosol + ! - flag_aerosol=7 => natural aerosol + MACv2SP flag_aerosol_omp = 0 @@ -2448,10 +2449,10 @@ ! Flag_aerosol cannot be to zero if we are in coupled mode for aerosol - IF (aerosol_couple .AND. flag_aerosol .eq. 0 ) THEN + IF (aerosol_couple .AND. flag_aerosol .EQ. 0 ) THEN CALL abort_physic('conf_phys', 'flag_aerosol cannot be to zero if aerosol_couple=y ', 1) ENDIF ! flag_aerosol need to be different to zero if ok_cdnc is activated - IF (ok_cdnc .AND. flag_aerosol .eq. 0) THEN + IF (ok_cdnc .AND. flag_aerosol .EQ. 0) THEN CALL abort_physic('conf_phys', 'flag_aerosol cannot be to zero if ok_cdnc is activated ', 1) ENDIF @@ -2460,4 +2461,12 @@ IF (ok_aie .AND. .NOT. ok_cdnc) THEN CALL abort_physic('conf_phys', 'ok_cdnc must be set to y if ok_aie is activated',1) + ENDIF + + ! flag_aerosol=7 => MACv2SP climatology + IF (flag_aerosol.EQ.7.AND. iflag_rrtm.NE.1) THEN + CALL abort_physic('conf_phys', 'flag_aerosol=7 (MACv2SP) can only be activated with RRTM',1) + ENDIF + IF (flag_aerosol.EQ.7.AND. NSW.NE.6) THEN + CALL abort_physic('conf_phys', 'flag_aerosol=7 (MACv2SP) can only be activated with NSW=6',1) ENDIF Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/macv2sp.F90 =================================================================== --- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/macv2sp.F90 (revision 3318) +++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/macv2sp.F90 (revision 3318) @@ -0,0 +1,174 @@ +SUBROUTINE MACv2SP(pphis,pplay,paprs,xlon,xlat,tau_allaer,piz_allaer,cg_allaer) + ! + !--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 + ! + USE mo_simple_plumes, ONLY: sp_aop_profile + USE phys_cal_mod, ONLY : year_cur, day_cur, year_len + USE dimphy + USE aero_mod + USE phys_local_var_mod, ONLY: t_seri, od443aer, od550aer, od865aer, ec550aer, dryod550aer, od550lt1aer, dNovrN + !!USE YOMCST, ONLY : RD, RG + ! + IMPLICIT NONE + ! + include "YOMCST.h" + ! + 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(OUT) :: tau_allaer ! epaisseur optique aerosol + REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(OUT) :: piz_allaer ! single scattering albedo aerosol + REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(OUT) :: cg_allaer ! asymmetry parameter aerosol + ! + 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 + ! + 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 :: year_fr + ! + INTEGER band, i, k, Nwv + ! + ! define the height and dheight arrays + ! + oro(:) = pphis(:)/RG ! surface height in m + ! + 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 + ! + !--fractional year + ! + year_fr = FLOAT(year_cur) + (FLOAT(day_cur)-0.5) / FLOAT(year_len) + print *,'year_fr=',year_fr + ! + !--call to sp routine -- 443 nm + ! + CALL sp_aop_profile ( & + klev ,klon ,l443 ,oro ,xlon ,xlat , & + year_fr ,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 , & + year_fr ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , & + asy_prof ) + ! + !--AOD calculations for diagnostics + od550aer(:)=od550aer(:)+SUM(aod_prof(:,:),dim=2) + ! + !--dry AOD calculation for diagnostics + dryod550aer(:)=dryod550aer(:)+od550aer(:) + ! + !--fine-mode AOD calculation for diagnostics + od550lt1aer(:)=od550lt1aer(:)+od550aer(:) + ! + !--extinction coefficient for diagnostic + ec550aer(:,:)=ec550aer(:,:)+aod_prof(:,:)/dz(:,:) + ! + !--call to sp routine -- 865 nm + ! + CALL sp_aop_profile ( & + klev ,klon ,l865 ,oro ,xlon ,xlat , & + year_fr ,z ,dz ,dNovrN ,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 + !--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 ,oro ,xlon ,xlat , & + year_fr ,z ,dz ,dNovrN ,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 + ! + !--renpomalizing 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 Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/mo_simple_plumes.F90 =================================================================== --- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/mo_simple_plumes.F90 (revision 3318) +++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/mo_simple_plumes.F90 (revision 3318) @@ -0,0 +1,413 @@ +!> +!! +!! @brief Module MO_SIMPLE_PLUMES: provides anthropogenic aerosol optical properties as a function of lat, lon +!! height, time, and wavelength +!! +!! @remarks +!! +!! @author Bjorn Stevens, Stephanie Fiedler and Karsten Peters MPI-Met, Hamburg (v1 release 2016-11-10) +!! +!! @change-log: +!! - 2016-12-05: beta release (BS, SF and KP, MPI-Met) +!! - 2016-09-28: revised representation of Twomey effect (SF, MPI-Met) +!! - 2015-09-28: bug fixes (SF, MPI-Met) +!! - 2016-10-12: revised maximum longitudinal extent of European plume (KP, SF, MPI-Met) +!! $ID: n/a$ +!! +!! @par Origin +!! Based on code originally developed at the MPI-Met by Karsten Peters, Bjorn Stevens, Stephanie Fiedler +!! and Stefan Kinne with input from Thorsten Mauritsen and Robert Pincus +!! +!! @par Copyright +!! +! +MODULE MO_SIMPLE_PLUMES + + USE netcdf + + IMPLICIT NONE + + 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) + nyears = 251 !< Number of years of available forcing + + LOGICAL, SAVE :: & + sp_initialized = .FALSE. !< parameter determining whether input needs to be read +!$OMP THREADPRIVATE(sp_initialized) + + REAL :: & + plume_lat (nplumes) ,& !< latitude of plume center (AOD maximum) + plume_lon (nplumes) ,& !< longitude of plume center (AOD maximum) + beta_a (nplumes) ,& !< parameter a for beta function vertical profile + beta_b (nplumes) ,& !< parameter b for beta function vertical profile + aod_spmx (nplumes) ,& !< anthropogenic AOD maximum at 550 for plumes + aod_fmbg (nplumes) ,& !< anthropogenic AOD at 550 for fine-mode natural background (idealized to mimic Twomey effect) + asy550 (nplumes) ,& !< asymmetry parameter at 550nm for plume + ssa550 (nplumes) ,& !< single scattering albedo at 550nm for plume + angstrom (nplumes) ,& !< Angstrom parameter for plume + sig_lon_E (nfeatures,nplumes) ,& !< Eastward extent of plume feature + sig_lon_W (nfeatures,nplumes) ,& !< Westward extent of plume feature + sig_lat_E (nfeatures,nplumes) ,& !< Southward extent of plume feature + sig_lat_W (nfeatures,nplumes) ,& !< Northward extent of plume feature + theta (nfeatures,nplumes) ,& !< Rotation angle of plume feature + ftr_weight (nfeatures,nplumes) ,& !< Feature weights + time_weight (nfeatures,nplumes) ,& !< Time weights + time_weight_bg (nfeatures,nplumes) ,& !< as time_weight but for natural background in Twomey effect + year_weight (nyears,nplumes) ,& !< Yearly weight for plume + ann_cycle (nfeatures,ntimes,nplumes) !< annual cycle for plume feature +!$OMP THREADPRIVATE(plume_lat,plume_lon,beta_a,beta_b,aod_spmx,aod_fmbg,asy550,ssa550,angstrom) +!$OMP THREADPRIVATE(sig_lon_E,sig_lon_W,sig_lat_E,sig_lat_W,theta,ftr_weight,year_weight,ann_cycle) + + PUBLIC sp_aop_profile + +CONTAINS + ! + ! ------------------------------------------------------------------------------------------------------------------------ + ! SP_SETUP: This subroutine should be called at initialization to read the netcdf data that describes the simple plume + ! climatology. The information needs to be either read by each processor or distributed to processors. + ! + SUBROUTINE sp_setup + ! + USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root + USE mod_phys_lmdz_omp_data, ONLY: is_omp_root + USE mod_phys_lmdz_transfert_para, ONLY: bcast + ! + ! ---------- + ! + INTEGER :: iret, ncid, DimID, VarID, xdmy + ! + ! ---------- + !--only one processor reads the input data + IF (is_mpi_root.AND.is_omp_root) THEN + ! + iret = nf90_open("MACv2.0-SP_v1.nc", NF90_NOWRITE, ncid) + IF (iret /= NF90_NOERR) STOP 'NetCDF File not opened' + ! + ! read dimensions and make sure file conforms to expected size + ! + iret = nf90_inq_dimid(ncid, "plume_number" , DimId) + iret = nf90_inquire_dimension(ncid, DimId, len = xdmy) + IF (xdmy /= nplumes) STOP 'NetCDF improperly dimensioned -- plume_number' + ! + iret = nf90_inq_dimid(ncid, "plume_feature", DimId) + iret = nf90_inquire_dimension(ncid, DimId, len = xdmy) + IF (xdmy /= nfeatures) STOP 'NetCDF improperly dimensioned -- plume_feature' + ! + iret = nf90_inq_dimid(ncid, "year_fr" , DimId) + iret = nf90_inquire_dimension(ncid, DimID, len = xdmy) + IF (xdmy /= ntimes) STOP 'NetCDF improperly dimensioned -- year_fr' + ! + iret = nf90_inq_dimid(ncid, "years" , DimId) + iret = nf90_inquire_dimension(ncid, DimID, len = xdmy) + IF (xdmy /= nyears) STOP 'NetCDF improperly dimensioned -- years' + ! + ! read variables that define the simple plume climatology + ! + iret = nf90_inq_varid(ncid, "plume_lat", VarId) + iret = nf90_get_var(ncid, VarID, plume_lat(:), start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lat' + iret = nf90_inq_varid(ncid, "plume_lon", VarId) + iret = nf90_get_var(ncid, VarID, plume_lon(:), start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lon' + iret = nf90_inq_varid(ncid, "beta_a" , VarId) + iret = nf90_get_var(ncid, VarID, beta_a(:) , start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading beta_a' + iret = nf90_inq_varid(ncid, "beta_b" , VarId) + iret = nf90_get_var(ncid, VarID, beta_b(:) , start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading beta_b' + iret = nf90_inq_varid(ncid, "aod_spmx" , VarId) + iret = nf90_get_var(ncid, VarID, aod_spmx(:) , start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading aod_spmx' + iret = nf90_inq_varid(ncid, "aod_fmbg" , VarId) + iret = nf90_get_var(ncid, VarID, aod_fmbg(:) , start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading aod_fmbg' + iret = nf90_inq_varid(ncid, "ssa550" , VarId) + iret = nf90_get_var(ncid, VarID, ssa550(:) , start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading ssa550' + iret = nf90_inq_varid(ncid, "asy550" , VarId) + iret = nf90_get_var(ncid, VarID, asy550(:) , start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading asy550' + iret = nf90_inq_varid(ncid, "angstrom" , VarId) + iret = nf90_get_var(ncid, VarID, angstrom(:), start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading angstrom' + ! + iret = nf90_inq_varid(ncid, "sig_lat_W" , VarId) + iret = nf90_get_var(ncid, VarID, sig_lat_W(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lat_W' + iret = nf90_inq_varid(ncid, "sig_lat_E" , VarId) + iret = nf90_get_var(ncid, VarID, sig_lat_E(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lat_E' + iret = nf90_inq_varid(ncid, "sig_lon_E" , VarId) + iret = nf90_get_var(ncid, VarID, sig_lon_E(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lon_E' + iret = nf90_inq_varid(ncid, "sig_lon_W" , VarId) + iret = nf90_get_var(ncid, VarID, sig_lon_W(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lon_W' + iret = nf90_inq_varid(ncid, "theta" , VarId) + iret = nf90_get_var(ncid, VarID, theta(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading theta' + iret = nf90_inq_varid(ncid, "ftr_weight" , VarId) + iret = nf90_get_var(ncid, VarID, ftr_weight(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lat' + iret = nf90_inq_varid(ncid, "year_weight" , VarId) + iret = nf90_get_var(ncid, VarID, year_weight(:,:) , start=(/1,1/),count=(/nyears,nplumes /)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading year_weight' + iret = nf90_inq_varid(ncid, "ann_cycle" , VarId) + iret = nf90_get_var(ncid, VarID, ann_cycle(:,:,:) , start=(/1,1,1/),count=(/nfeatures,ntimes,nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading ann_cycle' + ! + iret = nf90_close(ncid) + ! + ENDIF !--root processor + ! + CALL bcast(plume_lat) + CALL bcast(plume_lon) + CALL bcast(beta_a) + CALL bcast(beta_b) + CALL bcast(aod_spmx) + CALL bcast(aod_fmbg) + CALL bcast(asy550) + CALL bcast(ssa550) + CALL bcast(angstrom) + CALL bcast(sig_lon_E) + CALL bcast(sig_lon_W) + CALL bcast(sig_lat_E) + CALL bcast(sig_lat_W) + CALL bcast(theta) + CALL bcast(ftr_weight) + CALL bcast(year_weight) + CALL bcast(ann_cycle) + ! + sp_initialized = .TRUE. + ! + RETURN + ! + END SUBROUTINE sp_setup + ! + ! ------------------------------------------------------------------------------------------------------------------------ + ! SET_TIME_WEIGHT: The simple plume model assumes that meteorology constrains plume shape and that only source strength + ! influences the amplitude of a plume associated with a given source region. This routine retrieves the temporal weights + ! for the plumes. Each plume feature has its own temporal weights which varies yearly. The annual cycle is indexed by + ! week in the year and superimposed on the yearly mean value of the weight. + ! + SUBROUTINE set_time_weight(year_fr) + ! + ! ---------- + ! + REAL, INTENT(IN) :: & + year_fr !< Fractional Year (1850.0 - 2100.99) + + INTEGER :: & + iyear ,& !< Integer year values between 1 and 156 (1850-2100) + iweek ,& !< Integer index (between 1 and ntimes); for ntimes=52 this corresponds to weeks (roughly) + iplume ! plume number + ! + ! ---------- + ! + iyear = FLOOR(year_fr) - 1849 + iweek = FLOOR((year_fr - FLOOR(year_fr)) * ntimes) + 1 + + IF ((iweek > ntimes) .OR. (iweek < 1) .OR. (iyear > nyears) .OR. (iyear < 1)) THEN + CALL abort_physic('set_time_weight','Time out of bounds') + ENDIF + + DO iplume=1,nplumes + time_weight(1,iplume) = year_weight(iyear,iplume) * ann_cycle(1,iweek,iplume) + time_weight(2,iplume) = year_weight(iyear,iplume) * ann_cycle(2,iweek,iplume) + time_weight_bg(1,iplume) = ann_cycle(1,iweek,iplume) + time_weight_bg(2,iplume) = ann_cycle(2,iweek,iplume) + END DO + + RETURN + 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. + ! + SUBROUTINE sp_aop_profile ( & + nlevels ,ncol ,lambda ,oro ,lon ,lat , & + year_fr ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , & + asy_prof ) + ! + ! ---------- + ! + INTEGER, INTENT(IN) :: & + nlevels, & !< number of levels + ncol !< number of columns + + REAL, INTENT(IN) :: & + lambda, & !< wavelength + year_fr, & !< 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 + z (ncol,nlevels), & !< height above sea-level (m) + dz(ncol,nlevels) !< level thickness (difference between half levels) (m) + + 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 + + INTEGER :: iplume, icol, k + + REAL :: & + eta(ncol,nlevels), & !< normalized height (by 15 km) + z_beta(ncol,nlevels), & !< profile for scaling column optical depth + 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 + caod_sp(ncol), & !< column simple plume anthropogenic AOD at 550 nm + caod_bg(ncol), & !< column fine-mode natural background AOD at 550 nm + a_plume1, & !< gaussian longitude factor for feature 1 + a_plume2, & !< gaussian longitude factor for feature 2 + b_plume1, & !< gaussian latitude factor for feature 1 + b_plume2, & !< gaussian latitude factor for feature 2 + delta_lat, & !< latitude offset + delta_lon, & !< longitude offset + delta_lon_t, & !< threshold for maximum longitudinal plume extent used in transition from 360 to 0 degrees + lon1, & !< rotated longitude for feature 1 + lat1, & !< rotated latitude for feature 2 + lon2, & !< rotated longitude for feature 1 + lat2, & !< rotated latitude for feature 2 + f1, & !< contribution from feature 1 + 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 + ! + ! ---------- + ! + ! initialize input data (by calling setup at first instance) + ! + IF (.NOT.sp_initialized) CALL sp_setup + ! + ! get time weights + ! + CALL set_time_weight(year_fr) + ! + ! initialize variables, including output + ! + DO k=1,nlevels + DO icol=1,ncol + aod_prof(icol,k) = 0.0 + 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.)) + END DO + END DO + DO icol=1,ncol + dNovrN(icol) = 1.0 + caod_sp(icol) = 0.0 + caod_bg(icol) = 0.02 + END DO + ! + ! 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. + END DO + 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) + END DO + END DO + DO k=1,nlevels + DO icol=1,ncol + prof(icol,k) = ( prof(icol,k) / beta_sum(icol) ) * z_beta(icol,k) + END DO + END DO + ! + ! 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) + END DO + ! + ! 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)) + 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 + END DO + END DO + END DO + ! + ! complete optical depth weighting + ! + DO k=1,nlevels + DO icol = 1,ncol + asy_prof(icol,k) = MERGE(asy_prof(icol,k)/ssa_prof(icol,k), 0.0, ssa_prof(icol,k) > TINY(1.)) + ssa_prof(icol,k) = MERGE(ssa_prof(icol,k)/aod_prof(icol,k), 1.0, aod_prof(icol,k) > TINY(1.)) + END DO + END DO + ! + ! calculate effective radius normalization (divisor) factor + ! + 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) + END DO + + RETURN + END SUBROUTINE sp_aop_profile + +END MODULE MO_SIMPLE_PLUMES Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/newmicro.F90 =================================================================== --- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/newmicro.F90 (revision 3317) +++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/newmicro.F90 (revision 3318) @@ -3,5 +3,5 @@ -SUBROUTINE newmicro(ok_cdnc, bl95_b0, bl95_b1, paprs, pplay, t, pqlwp, pclc, & +SUBROUTINE newmicro(flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, paprs, pplay, t, pqlwp, pclc, & pcltau, pclemi, pch, pcl, pcm, pct, pctlwp, xflwp, xfiwp, xflwc, xfiwc, & mass_solu_aero, mass_solu_aero_pi, pcldtaupi, re, fl, reliq, reice, & @@ -10,5 +10,6 @@ USE dimphy USE phys_local_var_mod, ONLY: scdnc, cldncl, reffclwtop, lcc, reffclws, & - reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, zfice + reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, & + zfice, dNovrN USE phys_state_var_mod, ONLY: rnebcon, clwcon USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14) @@ -142,4 +143,5 @@ ! within the grid cell) + INTEGER flag_aerosol LOGICAL ok_cdnc REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula @@ -257,8 +259,4 @@ ! to be within [20, 1000 cm^3] - ! --present-day case - cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), & - 1.E-4))/log(10.))*1.E6 !-m-3 - cdnc(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc(i,k))) ! --pre-industrial case @@ -266,4 +264,18 @@ 1.E-4))/log(10.))*1.E6 !-m-3 cdnc_pi(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc_pi(i,k))) + + ! --present-day case + ! --flag_aerosol=7 => MACv2SP climatology + ! in this case there is an enhancement factor + IF (flag_aerosol .EQ. 7) THEN + cdnc(i, k) = cdnc_pi(i,k)*dNovrN(i) + ELSE + !--standard case, present day + cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), & + 1.E-4))/log(10.))*1.E6 !-m-3 + cdnc(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc(i,k))) + END IF + + ! --present-day case Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_local_var_mod.F90 =================================================================== --- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_local_var_mod.F90 (revision 3317) +++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_local_var_mod.F90 (revision 3318) @@ -151,4 +151,6 @@ REAL, SAVE, ALLOCATABLE :: scdnc(:,:) !$OMP THREADPRIVATE(scdnc) + REAL, SAVE, ALLOCATABLE :: dNovrN(:) + !$OMP THREADPRIVATE(dNovrN) REAL, SAVE, ALLOCATABLE :: cldncl(:) !$OMP THREADPRIVATE(cldncl) @@ -615,4 +617,5 @@ ALLOCATE(tau3d_aero(klon,klev,nwave,naero_tot)) ALLOCATE(scdnc(klon, klev)) + ALLOCATE(dNovrN(klon)) ALLOCATE(cldncl(klon)) ALLOCATE(reffclwtop(klon)) @@ -909,4 +912,5 @@ DEALLOCATE(tau3d_aero) DEALLOCATE(scdnc) + DEALLOCATE(dNovrN) DEALLOCATE(cldncl) DEALLOCATE(reffclwtop) Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90 =================================================================== --- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90 (revision 3317) +++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90 (revision 3318) @@ -3564,4 +3564,10 @@ tausum_aero, drytausum_aero, tau3d_aero) #endif + + 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) + ENDIF + ! ELSE IF (NSW.EQ.2) THEN @@ -3688,5 +3694,5 @@ #endif ENDIF - CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, & + CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, & paprs, pplay, t_seri, cldliq, cldfra, & cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, & Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90 =================================================================== --- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90 (revision 3317) +++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90 (revision 3318) @@ -381,5 +381,5 @@ ALLOCATE (aerosol_name(nb_aer)) aerosol_name(1) = id_CIDUSTM_phy - ELSEIF (flag_aerosol .EQ. 6) THEN + ELSEIF (flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN nb_aer = 13 ALLOCATE (aerosol_name(nb_aer)) Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90 =================================================================== --- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90 (revision 3317) +++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90 (revision 3318) @@ -507,5 +507,5 @@ ALLOCATE (aerosol_name(nb_aer)) aerosol_name(1) = id_CIDUSTM_phy - ELSEIF (flag_aerosol .EQ. 6) THEN + ELSEIF (flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN nb_aer = 13 ALLOCATE (aerosol_name(nb_aer)) Index: /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_lw_rrtm.F90 =================================================================== --- /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_lw_rrtm.F90 (revision 3317) +++ /LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_lw_rrtm.F90 (revision 3318) @@ -46,5 +46,5 @@ IF (ok_alw) THEN !--aerosol LW effects ! - IF (flag_aerosol.EQ.5.OR.flag_aerosol.EQ.6) THEN !-Dust + IF (flag_aerosol.EQ.5.OR.flag_aerosol.EQ.6.OR.flag_aerosol.EQ.7) THEN !-Dust ! zdh(:,:)=pdel(:,:)/(RG*zrho(:,:)) ! m