Subroutine aeropacity(ngrid,nlayer,pplay,pplev,pt, & aerosol,reffrad,nueffrad,QREFvis3d,tau_col) use radinc_h, only : naerkind use aerosol_mod implicit none #include "YOMCST.h" !================================================================== ! ! 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) ! Generic n-layer aerosol - J. Vatant d'Ollone (2020) ! ! Input ! ----- ! ngrid Number of horizontal gridpoints ! nlayer Number of layers ! pplev Pressure (Pa) at each layer boundary ! reffrad(ngrid,nlayer,naerkind) Aerosol effective radius ! QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients ! QREFir3d(ngrid,nlayer,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 ! !======================================================================= INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa) REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa) REAL,INTENT(IN) :: pt(ngrid,nlayer) ! mid-layer temperatre (K) REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind) ! aerosol effective variance REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth INTEGER l,ig,iaer,ll LOGICAL,SAVE :: firstcall=.true. !$OMP THREADPRIVATE(firstcall) character*80 abort_message character*20 modname ! for venus clouds real,ALLOCATABLE,SAVE :: pbot2(:,:),pbot(:,:),ptop(:,:),ptop2(:,:) real,ALLOCATABLE,SAVE :: hbot2(:,:),hbot(:,:),htop(:,:),htop2(:,:) real,ALLOCATABLE,SAVE :: nmode(:,:) !$OMP THREADPRIVATE(pbot2,pbot,ptop,ptop2) !$OMP THREADPRIVATE(hbot2,hbot,htop,htop2) !$OMP THREADPRIVATE(nmode) real :: mode_dens,hlay ! --------------------------------------------------------- ! --------------------------------------------------------- ! First call only IF (firstcall) THEN modname = 'aeropacity' ! verify tracers write(*,*) "Active aerosols found in aeropacity, designed for Venus" if (iaero_venus1.eq.1) then print*,'iaero_venus1= ',iaero_venus1 else abort_message='iaero_venus1 is not 1...' call abort_physic(modname,abort_message,1) endif if (iaero_venus2.eq.2) then print*,'iaero_venus2= ',iaero_venus2 else abort_message='iaero_venus2 is not 2...' call abort_physic(modname,abort_message,1) endif if (iaero_venus2p.eq.3) then print*,'iaero_venus2p= ',iaero_venus2p else abort_message='iaero_venus2p is not 3...' call abort_physic(modname,abort_message,1) endif if (iaero_venus3.eq.4) then print*,'iaero_venus3= ',iaero_venus3 else abort_message='iaero_venus3 is not 4...' call abort_physic(modname,abort_message,1) endif if (iaero_venusUV.eq.5) then print*,'iaero_venusUV= ',iaero_venusUV else abort_message='iaero_venus1UV is not 5...' call abort_physic(modname,abort_message,1) endif ! cloud model allocate(pbot2(ngrid,naerkind),pbot(ngrid,naerkind)) allocate(ptop(ngrid,naerkind),ptop2(ngrid,naerkind)) allocate(hbot2(ngrid,naerkind),hbot(ngrid,naerkind)) allocate(htop(ngrid,naerkind),htop2(ngrid,naerkind)) allocate(nmode(ngrid,naerkind)) call cloud_haus16(ngrid,pbot2,pbot,hbot2,hbot,ptop,ptop2,htop,htop2,nmode) firstcall=.false. ENDIF ! of IF (firstcall) ! --------------------------------------------------------- ! --------------------------------------------------------- aerosol(:,:,:) = 0. tau_col(:) = 0. do ig=1,ngrid do iaer=1,naerkind ! Opacity calculation ! determine the approximate middle of the mode layer ll=1 DO l=1,nlayer-1 ! high p to low p IF (pplay(ig,l) .gt. (ptop(ig,iaer)+pbot(ig,iaer))/2) ll=l ENDDO ! from there go DOWN for profile N mode_dens = nmode(ig,iaer) ! m-3 DO l=ll,1,-1 hlay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG IF (pplay(ig,l) .le. pbot(ig,iaer)) THEN mode_dens = nmode(ig,iaer) ! m-3 ELSEIF (pplay(ig,l) .gt. pbot(ig,iaer) .and. pplay(ig,l) .le. pbot2(ig,iaer)) THEN mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(hlay/hbot(ig,iaer)) ELSEIF (pplay(ig,l) .gt. pbot2(ig,iaer)) THEN mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(hlay/hbot2(ig,iaer)) ENDIF mode_dens=max(1.e-30,mode_dens) if (iaer.lt.5) then aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* & RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & mode_dens*hlay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l) else ! for UV absorber ! normalized to 0.35 microns (peak of absorption) aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens endif ENDDO ! ... then UP for profile N mode_dens = nmode(ig,iaer) ! m-3 DO l=ll+1,nlayer-1 hlay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG IF (pplay(ig,l) .gt. ptop(ig,iaer)) THEN mode_dens = nmode(ig,iaer) ! m-3 ELSEIF (pplay(ig,l) .le. ptop(ig,iaer) .and. pplay(ig,l) .gt. ptop2(ig,iaer)) THEN mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(hlay/htop(ig,iaer)) ELSEIF (pplay(ig,l) .le. ptop2(ig,iaer)) THEN mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(hlay/htop2(ig,iaer)) ENDIF mode_dens=max(1.e-30,mode_dens) if (iaer.lt.5) then aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* & RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & mode_dens*hlay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l) else ! for UV absorber ! normalized to 0.35 microns (peak of absorption) aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens endif ENDDO if (iaer.eq.5) then ! for UV absorber DO l=1,nlayer tau_col(ig) = tau_col(ig) & + aerosol(ig,l,iaer) ENDDO DO l=1,nlayer aerosol(ig,l,iaer)=aerosol(ig,l,iaer)*0.205/tau_col(ig) ENDDO endif enddo enddo !================================================================== ! if (is_master) then ! ig=10 ! do l=1,nlayer ! write(82,*) l,pplay(ig,l),aerosol(ig,l,1) ! write(83,*) l,pplay(ig,l),aerosol(ig,l,2) ! write(84,*) l,pplay(ig,l),aerosol(ig,l,3) ! write(85,*) l,pplay(ig,l),aerosol(ig,l,4) ! write(86,*) l,pplay(ig,l),aerosol(ig,l,5) ! enddo ! endif ! stop !================================================================== ! -------------------------------------------------------------------------- ! Column integrated visible optical depth in each point (used for diagnostic) tau_col(:)=0.0 do ig=1,ngrid do l=1,nlayer do iaer = 1, naerkind 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 ! -------------------------------------------------------------------------- ! -------------------------------------------------------------------------- subroutine cloud_haus16(ngrid,pbot2,pbot,hbot2,hbot,ptop,ptop2,htop,htop2,nmode) !================================================================== ! Venus clouds (4 modes) ! S. Lebonnois (jan 2016) !================================================================== ! distributions from Haus et al, 2016 ! ------------------------------------ ! ! mode 1 2 2p 3 ! r (microns) 0.30 1.05 1.40 3.65 ! sigma 1.56 1.29 1.23 1.28 ! reff (microns) 0.49 1.23 1.56 4.25 ! nueff 0.21 0.067 0.044 0.063 ! (nueff=exp(ln^2 sigma)-1) ! ! pbot <=> zb ; ptop <=> zb+zc ; hbot <=> Hlo ; htop <=> Hup ! ppbot: N(i-1)=N(i)*(p(i)/p(i-1))**(hlay(i)/hbot) ! N is in m-3 ! ! values in each mode below from Table 1 and text (p.186) ! ! dTau = Qext*[pi*reff**2*exp(-3*ln(1+nueff))]*N*hlay*(-dp)/p ! ! Latitudinal dependence (values from tables 3 and 4): ! mode factor MF12(lat) for modes 1 and 2 ! mode 2' unchanged ! mode factor MF3(lat) for mode 3 ! + for mode 2 only : pbot(lat), ptop(lat), htop(lat), htop2(lat) !================================================================== use radinc_h, only : naerkind USE geometry_mod, ONLY: latitude_deg implicit none INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns real,INTENT(out) :: pbot2(ngrid,naerkind),pbot(ngrid,naerkind) real,INTENT(out) :: ptop(ngrid,naerkind),ptop2(ngrid,naerkind) real,INTENT(out) :: hbot2(ngrid,naerkind),hbot(ngrid,naerkind) real,INTENT(out) :: htop(ngrid,naerkind),htop2(ngrid,naerkind) real,INTENT(out) :: nmode(ngrid,naerkind) ! For latitudinal dependence integer :: i,j,ilat real :: latit,factlat integer,parameter :: nlat_H16 = 16 real :: lat_H16(nlat_H16) real :: MF12(nlat_H16),MF3(nlat_H16) real :: pbotM2(nlat_H16),ptopM2(nlat_H16),htopM2(nlat_H16) !------------------------- !! Equatorial model !------------------------- ! initialization pbot2(:,:) = 1.e8 pbot(:,:) = 1.e8 ptop2(:,:) = 0. ptop(:,:) = 0. hbot2(:,:) = 1. hbot(:,:) = 1. htop2(:,:) = 1. htop(:,:) = 1. nmode(:,:) = 0. ! table of values(lat,mode) ! Mode 1 pbot2(:,1) = 2.0e5 ! Pa hbot2(:,1) = 5.0e3 ! m pbot(:,1) = 1.2e5 hbot(:,1) = 1.0e3 nmode(:,1) = 1.935e8 ! m-3 ptop(:,1) = 1.0e4 htop(:,1) = 3.5e3 ptop2(:,1) = 4.5e2 htop2(:,1) = 2.0e3 ! Mode 2 pbot(:,2) = 1.0e4 hbot(:,2) = 3.0e3 nmode(:,2) = 1.00e8 ! m-3 ptop(:,2) = 8.0e3 htop(:,2) = 3.5e3 ptop2(:,2) = 4.5e2 htop2(:,2) = 2.0e3 ! Mode 2p pbot(:,3) = 1.2e5 hbot(:,3) = 0.1e3 nmode(:,3) = 5.0e7 ! m-3 ptop(:,3) = 2.3e4 htop(:,3) = 1.0e3 ! Mode 3 pbot(:,4) = 1.2e5 hbot(:,4) = 0.5e3 nmode(:,4) = 1.4e7 ! m-3 ptop(:,4) = 3.9e4 htop(:,4) = 1.0e3 ! UV absorber pbot(:,5) = 3.3e4 ! 58 km hbot(:,5) = 1.0e3 nmode(:,5) = 1.00e7 !m-3 ptop(:,5) = 3.7e3 ! 70 km htop(:,5) = 1.0e3 !---------------------------- !! Latitudinal variations !---------------------------- lat_H16(:) = (/0.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90./) MF12(:) = (/0.98,0.98,0.99,1.00,0.98,0.94,0.86,0.81,0.73,0.67,0.64,0.61,0.59,0.47,0.36,0.36/) MF3(:) = (/1.30,1.30,1.26,1.23,1.17,1.13,1.06,1.03,1.04,1.09,1.22,1.51,1.82,2.02,2.09,2.09/) pbotM2(:) = (/1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,9.4e3, & 9.2e3,9.1e3,9.6e3,1.14e4,1.21e4,1.25e4,1.25e4/) ptopM2(:) = (/8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,7.7e3, & 7.5e3,7.5e3,8.0e3,9.2e3,1.0e4,1.06e4,1.06e4/) htopM2(:) = (/3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.4e3, & 3.2e3,2.6e3,2.0e3,2.0e3,0.6e3,0.5e3,0.5e3/) do i=1,ngrid latit = abs(latitude_deg(i)) ilat=2 do j=2,nlat_H16 if (latit.gt.lat_H16(j-1)) ilat=j enddo factlat = (lat_H16(ilat)-latit)/(lat_H16(ilat)-lat_H16(ilat-1)) pbot(i,2) = pbotM2(ilat)*(1.-factlat) + pbotM2(ilat-1)*factlat ptop(i,2) = ptopM2(ilat)*(1.-factlat) + ptopM2(ilat-1)*factlat htop(i,2) = htopM2(ilat)*(1.-factlat) + htopM2(ilat-1)*factlat htop2(i,2) = min(htop2(i,2),htop(i,2)) nmode(i,1) = nmode(i,1)*(MF12(ilat)*(1.-factlat)+MF12(ilat-1)*factlat) nmode(i,2) = nmode(i,2)*(MF12(ilat)*(1.-factlat)+MF12(ilat-1)*factlat) nmode(i,3) = nmode(i,3)*(MF3(ilat) *(1.-factlat)+MF3(ilat-1) *factlat) enddo return end subroutine cloud_haus16