Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, & aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky) use radinc_h, only : L_TAUMAX,naerkind use aerosol_mod USE comgeomfi_h USE tracer_h implicit none !================================================================== ! ! Purpose ! ------- ! Compute aerosol optical depth in each gridbox. ! ! Authors ! ------- ! F. Forget ! F. Montmessin (water ice scheme) ! update J.-B. Madeleine (2008) ! dust removal, simplification by Robin Wordsworth (2009) ! ! Input ! ----- ! ngrid Number of horizontal gridpoints ! nlayer Number of layers ! nq Number of tracers ! pplev Pressure (Pa) at each layer boundary ! pq Aerosol mixing ratio ! reffrad(ngrid,nlayer,naerkind) Aerosol effective radius ! QREFvis3d(ngrid,nlayermx,naerkind) \ 3d extinction coefficients ! QREFir3d(ngrid,nlayermx,naerkind) / at reference wavelengths ! ! Output ! ------ ! aerosol Aerosol optical depth in layer l, grid point ig ! tau_col Total column optical depth at grid point ig ! !======================================================================= #include "dimensions.h" #include "dimphys.h" #include "callkeys.h" #include "comcstfi.h" #include "comvert.h" INTEGER ngrid,nlayer,nq REAL pplay(ngrid,nlayer) REAL pplev(ngrid,nlayer+1) REAL pq(ngrid,nlayer,nq) REAL aerosol(ngrid,nlayer,naerkind) REAL reffrad(ngrid,nlayer,naerkind) REAL QREFvis3d(ngrid,nlayermx,naerkind) REAL QREFir3d(ngrid,nlayermx,naerkind) REAL tau_col(ngrid) ! REAL tauref(ngrid), tau_col(ngrid) real cloudfrac(ngrid,nlayermx) real aerosol0 INTEGER l,ig,iq,iaer LOGICAL firstcall DATA firstcall/.true./ SAVE firstcall REAL CBRT EXTERNAL CBRT INTEGER,SAVE :: i_co2ice=0 ! co2 ice INTEGER,SAVE :: i_h2oice=0 ! water ice CHARACTER(LEN=20) :: tracername ! to temporarily store text ! for fixed dust profiles real topdust, expfactor, zp REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling ! BENJAMIN MODIFS real CLFtot real totcloudfrac(ngrid) logical clearsky ! identify tracers IF (firstcall) THEN write(*,*) "Tracers found in aeropacity:" do iq=1,nq tracername=noms(iq) if (tracername.eq."co2_ice") then i_co2ice=iq write(*,*) "i_co2ice=",i_co2ice endif if (tracername.eq."h2o_ice") then i_h2oice=iq write(*,*) "i_h2oice=",i_h2oice endif enddo if (noaero) then print*, "No active aerosols found in aeropacity" else print*, "If you would like to use aerosols, make sure any old" print*, "start files are updated in newstart using the option" print*, "q=0" write(*,*) "Active aerosols found in aeropacity:" endif if ((iaero_co2.ne.0).and.(.not.noaero)) then print*, 'iaero_co2= ',iaero_co2 endif if (iaero_h2o.ne.0) then print*,'iaero_h2o= ',iaero_h2o endif if (iaero_dust.ne.0) then print*,'iaero_dust= ',iaero_dust endif if (iaero_h2so4.ne.0) then print*,'iaero_h2so4= ',iaero_h2so4 endif firstcall=.false. ENDIF ! of IF (firstcall) ! --------------------------------------------------------- !================================================================== ! CO2 ice aerosols !================================================================== if (iaero_co2.ne.0) then iaer=iaero_co2 ! 1. Initialization aerosol(1:ngrid,1:nlayermx,iaer)=0.0 ! 2. Opacity calculation if (noaero) then ! aerosol set to zero aerosol(1:ngrid,1:nlayermx,iaer)=0.0 elseif (aerofixco2.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed aerosol(1:ngrid,1:nlayermx,iaer)=1.e-9 !aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option else DO ig=1, ngrid DO l=1,nlayer-1 ! to stop the rad tran bug aerosol0 = & ( 0.75 * QREFvis3d(ig,l,iaer) / & ( rho_co2 * reffrad(ig,l,iaer) ) ) * & ( pq(ig,l,i_co2ice) + 1.E-9 ) * & ( pplev(ig,l) - pplev(ig,l+1) ) / g aerosol0 = max(aerosol0,1.e-9) aerosol0 = min(aerosol0,L_TAUMAX) aerosol(ig,l,iaer) = aerosol0 ! aerosol(ig,l,iaer) = 0.0 ! print*, aerosol(ig,l,iaer) ! using cloud fraction ! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF)) ! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX) ENDDO ENDDO end if ! if fixed or varying end if ! if CO2 aerosols !================================================================== ! Water ice / liquid !================================================================== if (iaero_h2o.ne.0) then iaer=iaero_h2o ! 1. Initialization aerosol(1:ngrid,1:nlayermx,iaer)=0.0 ! 2. Opacity calculation if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then aerosol(1:ngrid,1:nlayermx,iaer) =1.e-9 ! put cloud at cloudlvl if(kastprof.and.(cloudlvl.ne.0.0))then ig=1 do l=1,nlayer if(int(cloudlvl).eq.l)then !if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then print*,'Inserting cloud at level ',l !aerosol(ig,l,iaer)=10.0 rho_ice=920.0 ! the Kasting approximation aerosol(ig,l,iaer) = & ( 0.75 * QREFvis3d(ig,l,iaer) / & ( rho_ice * reffrad(ig,l,iaer) ) ) * & !( pq(ig,l,i_h2oice) + 1.E-9 ) * & ( 4.0e-4 + 1.E-9 ) * & ( pplev(ig,l) - pplev(ig,l+1) ) / g open(115,file='clouds.out',form='formatted') write(115,*) l,aerosol(ig,l,iaer) close(115) return endif end do call abort endif else do ig=1, ngrid do l=1,nlayer-1 ! to stop the rad tran bug aerosol(ig,l,iaer) = & !modification by BC ( 0.75 * QREFvis3d(ig,l,iaer) / & ( rho_ice * reffrad(ig,l,iaer) ) ) * & ! pq(ig,l,i_h2oice) * & !JL I dropped the +1e-9 here to have the same !( pplev(ig,l) - pplev(ig,l+1) ) / g ! opacity in the clearsky=true and the ! clear=false/pq=0 case ( pq(ig,l,i_h2oice) + 1.E-9 ) * & ! Doing this makes the code unstable, so I have restored it (RW) ( pplev(ig,l) - pplev(ig,l+1) ) / g enddo enddo if(CLFvarying)then call totalcloudfrac(ngrid,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1:ngrid,1:nlayermx,iaer)) do ig=1, ngrid do l=1,nlayer-1 ! to stop the rad tran bug CLFtot = max(totcloudfrac(ig),0.01) aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9) enddo enddo else do ig=1, ngrid do l=1,nlayer-1 ! to stop the rad tran bug CLFtot = CLFfixval aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9) enddo enddo end if!(CLFvarying) endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky) end if ! End if h2o aerosol !================================================================== ! Dust !================================================================== if (iaero_dust.ne.0) then iaer=iaero_dust ! 1. Initialization aerosol(1:ngrid,1:nlayermx,iaer)=0.0 topdust=30.0 ! km (used to be 10.0 km) LK ! 2. Opacity calculation ! expfactor=0. DO l=1,nlayer-1 DO ig=1,ngrid ! Typical mixing ratio profile zp=(preff/pplay(ig,l))**(70./topdust) expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3) ! Vertical scaling function aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) & *expfactor ENDDO ENDDO ! Rescaling each layer to reproduce the choosen (or assimilated) ! dust extinction opacity at visible reference wavelength, which ! is scaled to the "preff" reference surface pressure available in comvert.h ! and stored in startfi.nc taudusttmp(1:ngrid)=0. DO l=1,nlayer DO ig=1,ngrid taudusttmp(ig) = taudusttmp(ig) & + aerosol(ig,l,iaer) ENDDO ENDDO DO l=1,nlayer-1 DO ig=1,ngrid aerosol(ig,l,iaer) = max(1E-20, & dusttau & * pplev(ig,1) / preff & * aerosol(ig,l,iaer) & / taudusttmp(ig)) ENDDO ENDDO end if ! If dust aerosol !================================================================== ! H2SO4 !================================================================== ! added by LK if (iaero_h2so4.ne.0) then iaer=iaero_h2so4 ! 1. Initialization aerosol(1:ngrid,1:nlayermx,iaer)=0.0 ! 2. Opacity calculation ! expfactor=0. DO l=1,nlayer-1 DO ig=1,ngrid ! Typical mixing ratio profile zp=(preff/pplay(ig,l))**(70./30) !emulating topdust expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3) ! Vertical scaling function aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor ENDDO ENDDO tauh2so4tmp(1:ngrid)=0. DO l=1,nlayer DO ig=1,ngrid tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer) ENDDO ENDDO DO l=1,nlayer-1 DO ig=1,ngrid aerosol(ig,l,iaer) = max(1E-20, & 1 & * pplev(ig,1) / preff & * aerosol(ig,l,iaer) & / tauh2so4tmp(ig)) ENDDO ENDDO ! 1/700. is assuming a "sulfurtau" of 1 ! Sulfur aerosol routine to be improved. ! aerosol0 = & ! ( 0.75 * QREFvis3d(ig,l,iaer) / & ! ( rho_h2so4 * reffrad(ig,l,iaer) ) ) * & ! ( pq(ig,l,i_h2so4) + 1.E-9 ) * & ! ( pplev(ig,l) - pplev(ig,l+1) ) / g ! aerosol0 = max(aerosol0,1.e-9) ! aerosol0 = min(aerosol0,L_TAUMAX) ! aerosol(ig,l,iaer) = aerosol0 ! ENDDO ! ENDDO end if ! -------------------------------------------------------------------------- ! Column integrated visible optical depth in each point (used for diagnostic) tau_col(:)=0.0 do iaer = 1, naerkind do l=1,nlayer do ig=1,ngrid tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer) end do end do end do do ig=1,ngrid do l=1,nlayer do iaer = 1, naerkind if(aerosol(ig,l,iaer).gt.1.e3)then print*,'WARNING: aerosol=',aerosol(ig,l,iaer) print*,'at ig=',ig,', l=',l,', iaer=',iaer print*,'QREFvis3d=',QREFvis3d(ig,l,iaer) print*,'reffrad=',reffrad(ig,l,iaer) endif end do end do end do do ig=1,ngrid if(tau_col(ig).gt.1.e3)then print*,'WARNING: tau_col=',tau_col(ig) print*,'at ig=',ig print*,'aerosol=',aerosol(ig,:,:) print*,'QREFvis3d=',QREFvis3d(ig,:,:) print*,'reffrad=',reffrad(ig,:,:) endif end do return end subroutine aeropacity