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

module aeropacity_mod

implicit none

contains

      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,zzlev,pt,pq, &
         dtau_aer,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,micro_indx
       use comcstfi_mod, only: g, pi, mugaz, avocado
       use geometry_mod, only: latitude
       use callkeys_mod, only: kastprof, callmufi
       use mp2m_diagnostics
       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,nlayer,naerkind) \ 3d extinction coefficients
!     QREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths
!
!     Output
!     ------
!     dtau_aer           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) :: zzlev(ngrid,nlayer)   ! Altitude at the layer boundaries.
      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
      REAL,INTENT(IN) :: pt(ngrid,nlayer) ! mid-layer temperature (K)
      REAL,INTENT(OUT) :: dtau_aer(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)

      ! Microphysical tracers
      real sig
      real m0as(ngrid,nlayer)
      real m0af(ngrid,nlayer)

      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
         if (callmufi) then
            ! Convert intensive microphysical tracers to extensive [m-2]
            m0as(:,:) = pq(:,:,micro_indx(1)) * (pplev(:,1:nlayer) - pplev(:,2:nlayer+1)) / g
            m0af(:,:) = pq(:,:,micro_indx(3)) * (pplev(:,1:nlayer) - pplev(:,2:nlayer+1)) / g
            
            ! Spherical aerosols
            sig = 0.2
            dtau_aer(:,:,1) = m0as(:,:) * QREFvis3d(:,:,1) * pi * mp2m_rc_sph(:,:)**2 * exp(2*sig**2)
            ! Fractal aerosols
            sig = 0.35
            dtau_aer(:,:,2) = m0af(:,:) * QREFvis3d(:,:,2) * pi * mp2m_rc_fra(:,:)**2 * exp(2*sig**2)
            
            ! write(*,*) 'dtau_as  :', MINVAL(dtau_aer(:,:,1)), '-', MAXVAL(dtau_aer(:,:,1))
            ! write(*,*) 'dtau_af  :', MINVAL(dtau_aer(:,:,2)), '-', MAXVAL(dtau_aer(:,:,2))
         
         else
            do iaer = 1, naerkind
               ! 1. Initialization
               dtau_aer(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
                     ! Eq. 2.37 - Madeleine's PhD (2011).
                     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)
                     dtau_aer(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)
            enddo ! end iaer = 1, naerkind
         endif ! end callmufi
      endif ! 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) + dtau_aer(ig,l,iaer)
         end do
      end do
   end do

   do ig=1,ngrid
      do l=1,nlayer
         do iaer = 1, naerkind
            if(dtau_aer(ig,l,iaer).gt.1.e3)then
               print*,'WARNING: dtau_aer=',dtau_aer(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*,'dtau_aer=',dtau_aer(ig,:,:)
         print*,'QREFvis3d=',QREFvis3d(ig,:,:)
         print*,'reffrad=',reffrad(ig,:,:)
      endif
   end do
  end subroutine aeropacity

end module aeropacity_mod