Index: LMDZ6/trunk/libf/phylmd/conf_phys_m.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/conf_phys_m.F90 (revision 3273) +++ LMDZ6/trunk/libf/phylmd/conf_phys_m.F90 (revision 3274) @@ -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/trunk/libf/phylmd/macv2sp.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/macv2sp.F90 (revision 3274) +++ LMDZ6/trunk/libf/phylmd/macv2sp.F90 (revision 3274) @@ -0,0 +1,167 @@ +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 YOERAD, ONLY : NLW + USE YOMCST, ONLY : RD, RG + ! + 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(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 + ! + 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 + ! + z(:,1) = 0. ! altitude of surface taken as 0 + DO k = 1, klev-1 + oro(:) = pphis(:)/RG ! surface height in m + 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 + z(:,k+1) = z(:,k) + dz(:,k) ! height of interfaces in m, starting from 0 at surface + 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/trunk/libf/phylmd/mo_simple_plumes_v1.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/mo_simple_plumes_v1.F90 (revision 3274) +++ LMDZ6/trunk/libf/phylmd/mo_simple_plumes_v1.F90 (revision 3274) @@ -0,0 +1,383 @@ +!> +!! +!! @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 + + 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 + + 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 + ! + ! ---------- + ! + INTEGER :: iret, ncid, DimID, VarID, xdmy + ! + ! ---------- + ! + 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) + + 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/trunk/libf/phylmd/newmicro.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/newmicro.F90 (revision 3273) +++ LMDZ6/trunk/libf/phylmd/newmicro.F90 (revision 3274) @@ -1,5 +1,5 @@ ! $Id$ -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, & @@ -8,5 +8,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) @@ -141,4 +142,5 @@ ! within the grid cell) + INTEGER flag_aerosol LOGICAL ok_cdnc REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula @@ -216,6 +218,6 @@ xflwc(i, k) = (1.-zfice(i,k))*pqlwp(i, k) xfiwc(i, k) = zfice(i, k)*pqlwp(i, k) - END DO - END DO + ENDDO + ENDDO ELSE ! of IF (iflag_t_glace.EQ.0) DO k = 1, klev @@ -234,6 +236,6 @@ xflwc(i, k) = (1.-zfice(i,k))*pqlwp(i, k) xfiwc(i, k) = zfice(i, k)*pqlwp(i, k) - END DO - END DO + ENDDO + ENDDO ENDIF @@ -244,13 +246,7 @@ DO k = 1, klev DO i = 1, klon - ! Formula "D" of Boucher and Lohmann, Tellus, 1995 ! Cloud droplet number concentration (CDNC) is restricted ! 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(1000.E6, max(cdnc_min_m3,cdnc(i,k))) ! --pre-industrial case @@ -258,4 +254,42 @@ 1.E-4))/log(10.))*1.E6 !-m-3 cdnc_pi(i, k) = min(1000.E6, max(cdnc_min_m3,cdnc_pi(i,k))) + + ENDDO + ENDDO + + !--flag_aerosol=7 => MACv2SP climatology + !--in this case there is an enhancement factor + IF (flag_aerosol .EQ. 7) THEN + + !--present-day + DO k = 1, klev + DO i = 1, klon + cdnc(i, k) = cdnc_pi(i,k)*dNovrN(i) + ENDDO + ENDDO + + !--standard case + ELSE + + DO k = 1, klev + DO i = 1, klon + + ! Formula "D" of Boucher and Lohmann, Tellus, 1995 + ! Cloud droplet number concentration (CDNC) is restricted + ! 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(1000.E6, max(cdnc_min_m3,cdnc(i,k))) + + ENDDO + ENDDO + + ENDIF !--flag_aerosol + + !--computing cloud droplet size + DO k = 1, klev + DO i = 1, klon ! --present-day case @@ -292,8 +326,8 @@ zfiwp_var*(3.448E-03+2.431/rei) - END IF - - END DO - END DO + ENDIF + + ENDDO + ENDDO ELSE !--not ok_cdnc @@ -305,14 +339,14 @@ rad_chaud(i, k) = rad_chau2 rad_chaud_pi(i, k) = rad_chau2 - END DO - END DO + ENDDO + ENDDO DO k = min(3, klev) + 1, klev DO i = 1, klon rad_chaud(i, k) = rad_chau1 rad_chaud_pi(i, k) = rad_chau1 - END DO - END DO - - END IF !--ok_cdnc + ENDDO + ENDDO + + ENDIF !--ok_cdnc ! --computation of cloud optical depth and emissivity @@ -389,5 +423,5 @@ pclemi(i, k) = 1.0 - exp(-coef_chau*zflwp_var-df*k_ice*zfiwp_var) - END IF + ENDIF reice(i, k) = rei @@ -396,6 +430,6 @@ xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k) - END DO - END DO + ENDDO + ENDDO ! --if cloud droplet radius is fixed, then pcldtaupi=pcltau @@ -406,7 +440,7 @@ pcldtaupi(i, k) = pcltau(i, k) reice_pi(i, k) = reice(i, k) - END DO - END DO - END IF + ENDDO + ENDDO + ENDIF DO k = 1, klev @@ -415,6 +449,6 @@ reliq_pi(i, k) = rad_chaud_pi(i, k) reice_pi(i, k) = reice(i, k) - END DO - END DO + ENDDO + ENDDO ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS @@ -432,5 +466,5 @@ pcl(i) = 1.0 pctlwp(i) = 0.0 - END DO + ENDDO ! --calculation of liquid water path @@ -439,6 +473,6 @@ DO i = 1, klon pctlwp(i) = pctlwp(i) + pqlwp(i, k)*rhodz(i, k) - END DO - END DO + ENDDO + ENDDO ! --calculation of cloud properties with cloud overlap @@ -462,8 +496,8 @@ (i),kind=8),1.-zepsec)) zcloudl(i) = pclc(i, k) - END IF + ENDIF zcloud(i) = pclc(i, k) - END DO - END DO + ENDDO + ENDDO ELSE IF (novlp==2) THEN DO k = klev, 1, -1 @@ -477,7 +511,7 @@ ELSE IF (paprs(i,k)>=prlmc) THEN pcl(i) = min(pclc(i,k), pcl(i)) - END IF - END DO - END DO + ENDIF + ENDDO + ENDDO ELSE IF (novlp==3) THEN DO k = klev, 1, -1 @@ -491,8 +525,8 @@ ELSE IF (paprs(i,k)>=prlmc) THEN pcl(i) = pcl(i)*(1.0-pclc(i,k)) - END IF - END DO - END DO - END IF + ENDIF + ENDDO + ENDDO + ENDIF DO i = 1, klon @@ -500,5 +534,5 @@ pcm(i) = 1. - pcm(i) pcl(i) = 1. - pcl(i) - END DO + ENDDO ! ======================================================== @@ -521,8 +555,8 @@ ELSE lcc3d(i, k) = pclc(i, k)*phase3d(i, k) - END IF + ENDIF scdnc(i, k) = lcc3d(i, k)*cdnc(i, k) ! m-3 - END DO - END DO + ENDDO + ENDDO DO i = 1, klon @@ -532,5 +566,5 @@ IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1. IF (novlp.EQ.2) tcc(i) = 0. - END DO + ENDDO DO i = 1, klon @@ -546,8 +580,8 @@ WRITE (*, *) 'Hypothese de recouvrement: MAXIMUM' first = .FALSE. - END IF + ENDIF flag_max = -1. ftmp(i) = max(tcc(i), pclc(i,k)) - END IF + ENDIF IF (novlp.EQ.3) THEN @@ -555,8 +589,8 @@ WRITE (*, *) 'Hypothese de recouvrement: RANDOM' first = .FALSE. - END IF + ENDIF flag_max = 1. ftmp(i) = tcc(i)*(1-pclc(i,k)) - END IF + ENDIF IF (novlp.EQ.1) THEN @@ -566,9 +600,9 @@ & RANDOM' first = .FALSE. - END IF + ENDIF flag_max = 1. ftmp(i) = tcc(i)*(1.-max(pclc(i,k),pclc(i,k+1)))/(1.-min(pclc(i, & k+1),1.-thres_neb)) - END IF + ENDIF ! Effective radius of cloud droplet at top of cloud (m) reffclwtop(i) = reffclwtop(i) + rad_chaud(i, k)*1.0E-06*phase3d(i, & @@ -582,10 +616,10 @@ tcc(i) = ftmp(i) - END IF ! is there a visible, not-too-small cloud? - END DO ! loop over k + ENDIF ! is there a visible, not-too-small cloud? + ENDDO ! loop over k IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1. - tcc(i) - END DO ! loop over i + ENDDO ! loop over i ! ! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC @@ -615,6 +649,6 @@ reffclws(i, k) = radius reffclws(i, k) = reffclws(i, k)*lcc3dstra(i, k) - END DO !klev - END DO !klon + ENDDO !klev + ENDDO !klon ! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D @@ -628,5 +662,5 @@ lcc_integrat(i) = lcc_integrat(i) + lcc3d(i, k)*dh(i, k) height(i) = height(i) + dh(i, k) - END DO ! klev + ENDDO ! klev lcc_integrat(i) = lcc_integrat(i)/height(i) IF (lcc_integrat(i)<=1.0E-03) THEN @@ -634,6 +668,6 @@ ELSE cldnvi(i) = cldnvi(i)*lcc(i)/lcc_integrat(i) - END IF - END DO ! klon + ENDIF + ENDDO ! klon DO i = 1, klon @@ -649,12 +683,12 @@ IF (icc3dstra(i,k)<=0.0) icc3dstra(i, k) = 0.0 !FC (CAUSES) - END DO + ENDDO IF (reffclwtop(i)<=0.0) reffclwtop(i) = 0.0 IF (cldncl(i)<=0.0) cldncl(i) = 0.0 IF (cldnvi(i)<=0.0) cldnvi(i) = 0.0 IF (lcc(i)<=0.0) lcc(i) = 0.0 - END DO - - END IF !ok_cdnc + ENDDO + + ENDIF !ok_cdnc first=.false. !to be sure Index: LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 3273) +++ LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 3274) @@ -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) @@ -608,4 +610,5 @@ ALLOCATE(tau3d_aero(klon,klev,nwave,naero_tot)) ALLOCATE(scdnc(klon, klev)) + ALLOCATE(dNovrN(klon)) ALLOCATE(cldncl(klon)) ALLOCATE(reffclwtop(klon)) @@ -899,4 +902,5 @@ DEALLOCATE(tau3d_aero) DEALLOCATE(scdnc) + DEALLOCATE(dNovrN) DEALLOCATE(cldncl) DEALLOCATE(reffclwtop) Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 3273) +++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 3274) @@ -3582,4 +3582,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 @@ -3706,5 +3712,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, &