Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, & aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky) use radinc_h, only : naerkind, L_TAUMAX 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(ngridmx,nlayermx,naerkind) \ 3d extinction coefficients ! QREFir3d(ngridmx,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 "comgeomfi.h" #include "tracer.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(ngridmx,nlayermx,naerkind) REAL QREFir3d(ngridmx,nlayermx,naerkind) REAL tau_col(ngrid) ! REAL tauref(ngrid), tau_col(ngrid) real cloudfrac(ngridmx,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(ngridmx) ! Temporary dust opacity used before scaling ! BENJAMIN MODIFS real CLFtot,CLF1,CLF2 real totcloudfrac(ngridmx) logical clearsky ! identify tracers IF (firstcall) THEN ! are these tests of any real use ? IF(ngrid.NE.ngridmx) THEN PRINT*,'STOP in aeropacity!' PRINT*,'problem of dimensions:' PRINT*,'ngrid =',ngrid PRINT*,'ngridmx =',ngridmx STOP ENDIF if (nq.gt.nqmx) then write(*,*) 'STOP in aeropacity: (nq .gt. nqmx)!' write(*,*) 'nq=',nq,' nqmx=',nqmx stop endif do iq=1,nqmx tracername=noms(iq) if (tracername.eq."co2_ice") then i_co2ice=iq endif if (tracername.eq."h2o_ice") then i_h2oice=iq endif enddo write(*,*) "aeropacity: i_co2ice=",i_co2ice write(*,*) " i_h2oice=",i_h2oice if(watercond.and.(.not.aerofixed))then if(naerkind.lt.2)then print*,'Cannot have active H2O clouds with naerkind less than 2!' call abort endif endif if(dusttau.gt.0.0)then if(naerkind.lt.3)then print*,'Cannot have active dust with naerkind less than 3!' call abort endif endif firstcall=.false. ENDIF ! of IF (firstcall) DO iaer = 1, naerkind ! Loop on aerosol kind (NOT tracer) ! --------------------------------------------------------- aerkind: SELECT CASE (iaer) !================================================================== CASE(1) aerkind ! CO2 ice (iaer=1) !================================================================== ! 1. Initialization aerosol(:,:,iaer)=0.0 ! 2. Opacity calculation if (aerofixed.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed do ig=1, ngrid do l=1,nlayer aerosol(ig,l,iaer)=1.e-9 end do !aerosol(ig,12,iaer)=4.0 ! single cloud layer option end do 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 ! 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 !================================================================== CASE(2) aerkind ! Water ice / liquid (iaer=2) !================================================================== ! 1. Initialization aerosol(:,:,iaer)=0.0 ! 2. Opacity calculation if (aerofixed.or.(i_h2oice.eq.0).or.clearsky) then do ig=1, ngrid ! to stop the rad tran bug do l=1,nlayer aerosol(ig,l,iaer) =1.e-9 end do end do else do ig=1, ngrid do l=1,nlayer-1 ! to stop the rad tran bug if(CLFvarying)then CLF1 = max(cloudfrac(ig,l),0.01) CLFtot = max(totcloudfrac(ig),0.01) CLF2 = CLF1/CLFtot CLF2 = min(1.,CLF2) else CLF1=1.0 CLF2=CLFfixval endif aerosol0 = & ( 0.75 * QREFvis3d(ig,l,iaer) / & ( rho_ice * reffrad(ig,l,iaer) ) ) * & ( pq(ig,l,i_h2oice) + 1.E-9 ) * & ( pplev(ig,l) - pplev(ig,l+1) ) / g / & CLF1 aerosol(ig,l,iaer) = -log(1 - CLF2 + CLF2*exp(-aerosol0)) ! why no division in exponential? ! because its already done in aerosol0 aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9) aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),aerosol0) enddo enddo end if !================================================================== CASE(3) aerkind ! Dust (iaer=3) !================================================================== ! 1. Initialization aerosol(:,:,iaer)=0.0 topdust=10.0 ! km print*,'WARNING, dust is experimental for Early Mars' ! 2. Opacity calculation do l=1,nlayer do ig=1,ngrid-1 ! to stop the rad tran bug zp=700./pplay(ig,l) aerosol(ig,l,1)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) & *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 ) & *QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer) ! takes into account particle variation enddo enddo !================================================================== END SELECT aerkind ENDDO ! iaer (loop on aerosol kind) ! -------------------------------------------------------------------------- ! 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