source: trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90 @ 3572

module aeropacity_mod

implicit none

contains

      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pt, pq, &
         aerosol,reffrad,nueffrad, QREFvis3d,QREFir3d,tau_col)

       use radinc_h, only : L_TAUMAX,naerkind
       use aerosol_mod, only: iaero_haze, i_haze
       USE tracer_h, only: noms,rho_n2,rho_ice,rho_q,mmol
       use comcstfi_mod, only: g, pi, mugaz, avocado
       use geometry_mod, only: latitude
       use callkeys_mod, only: kastprof
       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)
!     Generic n-layer aerosol - J. Vatant d'Ollone (2020)
!     Radiative Generic Condensable Species - Lucas Teinturier (2022)
!
!     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,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
      INTEGER,INTENT(IN) :: nq     ! number of tracers
      REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
      REAL,INTENT(IN) :: pt(ngrid,nlayer) ! mid-layer temperature (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(IN) :: QREFir3d(ngrid,nlayer,naerkind)
      REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth

      real aerosol0, obs_tau_col_aurora, pm
      real pcloud_deck, cloud_slope

      real dp_strato(ngrid)
      real dp_tropo(ngrid)
      real dp_layer(ngrid)

      INTEGER l,ig,iq,iaer,ia

      LOGICAL,SAVE :: firstcall=.true.
!$OMP THREADPRIVATE(firstcall)
      REAL CBRT
      EXTERNAL CBRT

      INTEGER,SAVE :: i_n2ice=0      ! n2 ice
!$OMP THREADPRIVATE(i_n2ice)
      CHARACTER(LEN=20) :: tracername ! to temporarily store text

      real CLFtot
      integer igen_ice,igen_gas ! to store the index of generic tracer
      logical dummy_bool ! dummy boolean just in case we need one
      !  for venus clouds
      real      :: p_bot,p_top,h_bot,h_top,mode_dens,h_lay

      ! identify tracers
      IF (firstcall) THEN
        ia =0
        write(*,*) "Tracers found in aeropacity:"
        do iq=1,nq
          tracername=noms(iq)
          if (tracername.eq."n2_ice") then
            i_n2ice=iq
          write(*,*) "i_n2ice=",i_n2ice

          endif
        enddo

        firstcall=.false.
      ENDIF ! of IF (firstcall)


!     ---------------------------------------------------------
!==================================================================
!    Haze aerosols
!==================================================================

      if (iaero_haze.ne.0) then
        iaer=iaero_haze
!       1. Initialization
         aerosol(1:ngrid,1:nlayer,iaer)=0.0
!       2. Opacity calculation
         DO ig=1, ngrid
               DO l=1,nlayer-1 ! to stop the rad tran bug
                  ! if fractal, radius doit etre equivalent sphere radius
                  aerosol0 =                         &
                       (  0.75 * QREFvis3d(ig,l,iaer) /        &
                       ( rho_q(i_haze) * reffrad(ig,l,iaer) )  ) *   &
                       ( pq(ig,l,i_haze) + 1.E-10 ) *         &
                       ( pplev(ig,l) - pplev(ig,l+1) ) / g
                  aerosol0           = max(aerosol0,1.e-10)
                  aerosol0           = min(aerosol0,L_TAUMAX)
                  aerosol(ig,l,iaer) = aerosol0
               ENDDO
         ENDDO
         !QREF est le meme dans toute la colonne par def si size uniforme
         !print*, 'TB17: QREFvis3d=',QREFvis3d(1,:,1)
         !print*, 'TB17: rho_q=',rho_q(i_haze)
         !print*, 'TB17: reffrad=',reffrad(1,:,1)
         !print*, 'TB17: pq=',pq(1,:,i_haze)
         !print*, 'TB17: deltap=',pplev(1,1) - pplev(1,nlayer)
   end if ! if haze aerosols

! --------------------------------------------------------------------------
! 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
  end subroutine aeropacity

end module aeropacity_mod