!==================================================================
module radii_mod
!==================================================================
!  module to centralize the radii calculations for aerosols
!==================================================================
      
      use callkeys_mod, only: pres_bottom_tropo,pres_top_tropo, &
                size_tropo,pres_bottom_strato,size_strato

contains


!==================================================================
   subroutine su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
!==================================================================
!     Purpose
!     -------
!     Compute the effective radii of liquid and icy water particles
!
!     Authors
!     -------
!     Jeremy Leconte (2012)
!
!==================================================================
      use ioipsl_getin_p_mod, only: getin_p
      use radinc_h, only: naerkind
      use aerosol_mod, only: iaero_back2lay
      Implicit none

      integer,intent(in) :: ngrid
      integer,intent(in) :: nlayer

      real, intent(out) :: reffrad(ngrid,nlayer,naerkind)      !aerosols radii (K)
      real, intent(out) :: nueffrad(ngrid,nlayer,naerkind)     !variance     

      logical, save :: firstcall=.true.
!$OMP THREADPRIVATE(firstcall)
      integer :: iaer   
      
      print*,'enter su_aer_radii'
          do iaer=1,naerkind
!     these values will change once the microphysics gets to work
!     UNLESS tracer=.false., in which case we should be working with
!     a fixed aerosol layer, and be able to define reffrad in a 
!     .def file. To be improved!


!     WARNING : Titan adapt. (J. Vatant d'Ollone, 2017)
!            - ONLY THE NO AEROSOL CASE FOR NOW SINCE WE COMPUTE THEM ANOTHER WAY !
!            - This routine is just here to keep the code running without unplugging all (yet)
!            - There's only the dummy aerosol case on iaer = 1      
            if(iaer.eq.1)then
               reffrad(1:ngrid,1:nlayer,iaer) = 1.e-4
               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
            endif
            
            if(iaer.eq.iaero_back2lay)then ! Two-layer aerosols
               reffrad(1:ngrid,1:nlayer,iaer) = 2.e-6
               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
            endif

            if(iaer.gt.5)then
               print*,'Error in callcorrk, naerkind is too high (>5).'
               print*,'The code still needs generalisation to arbitrary'
               print*,'aerosol kinds and number.'
               call abort
            endif

         enddo

      print*,'exit su_aer_radii'

   end subroutine su_aer_radii
!==================================================================


!==================================================================
   subroutine back2lay_reffrad(ngrid,reffrad,nlayer,pplev)
!==================================================================
!     Purpose
!     -------
!     Compute the effective radii of particles in a 2-layer model
!
!     Authors
!     -------
!     Sandrine Guerlet (2013)
!
!==================================================================
 
      use aerosol_mod   !! Particle sizes and boundaries of aerosol layers defined there
     Implicit none

      integer,intent(in) :: ngrid

      real, intent(out) :: reffrad(ngrid,nlayer)      ! particle radii (m)
      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
      REAL :: expfactor
      INTEGER l,ig
            
      reffrad(:,:)=1e-6  !!initialization, not important
          DO ig=1,ngrid
            DO l=1,nlayer-1
              IF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
                reffrad(ig,l) = size_tropo
              ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
                expfactor=log(size_strato/size_tropo) / log(pres_bottom_strato/pres_top_tropo)
                reffrad(ig,l)= size_tropo*((pplev(ig,l)/pres_top_tropo)**expfactor)
              ELSEIF (pplev(ig,l) .le. pres_bottom_strato) then
                reffrad(ig,l) = size_strato
              ENDIF
            ENDDO
          ENDDO

   end subroutine back2lay_reffrad
!==================================================================



end module radii_mod
!==================================================================