source: trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90 @ 773

!==================================================================
module radii_mod
!==================================================================
!  module to centralize the radii calculations for aerosols
! OK for water but should be extended to other aerosols (CO2,...)
!==================================================================
      
!     water cloud optical properties
      
      real, save ::  rad_h2o
      real, save ::  rad_h2o_ice
      real, save ::  Nmix_h2o
      real, save ::  Nmix_h2o_ice
      real, parameter ::  coef_chaud=0.13
      real, parameter ::  coef_froid=0.09


      contains


!==================================================================
   subroutine su_aer_radii(reffrad,nueffrad)
!==================================================================
!     Purpose
!     -------
!     Compute the effective radii of liquid and icy water particles
!
!     Authors
!     -------
!     Jeremy Leconte (2012)
!
!==================================================================
 ! to use  'getin'
      use ioipsl_getincom 
      use radinc_h, only: naerkind
      use aerosol_mod
      Implicit none

#include "callkeys.h"
#include "dimensions.h"
#include "dimphys.h"
#include "tracer.h"

      real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)
      real, intent(inout) :: nueffrad(ngridmx,nlayermx,naerkind)     !variance     

      logical, save :: firstcall=.true.
      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!

            if(iaer.eq.iaero_co2)then ! CO2 ice
               reffrad(1:ngridmx,1:nlayermx,iaer) = 1.e-4
               nueffrad(1:ngridmx,1:nlayermx,iaer) = 0.1 
            endif

            if(iaer.eq.iaero_h2o)then ! H2O ice
               reffrad(1:ngridmx,1:nlayermx,iaer) = 1.e-5
               nueffrad(1:ngridmx,1:nlayermx,iaer) = 0.1 
            endif

            if(iaer.eq.iaero_dust)then ! dust
               reffrad(1:ngridmx,1:nlayermx,iaer) = 1.e-5
               nueffrad(1:ngridmx,1:nlayermx,iaer) = 0.1 
            endif
 
            if(iaer.eq.iaero_h2so4)then ! H2O ice
               reffrad(1:ngridmx,1:nlayermx,iaer) = 1.e-6
               nueffrad(1:ngridmx,1:nlayermx,iaer) = 0.1 
            endif

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

         enddo


         if (radfixed) then

            write(*,*)"radius of H2O water particles:"
            rad_h2o=13. ! default value
            call getin("rad_h2o",rad_h2o)
            write(*,*)" rad_h2o = ",rad_h2o

            write(*,*)"radius of H2O ice particles:"
            rad_h2o_ice=35. ! default value
            call getin("rad_h2o_ice",rad_h2o_ice)
            write(*,*)" rad_h2o_ice = ",rad_h2o_ice

         else

            write(*,*)"Number mixing ratio of H2O water particles:"
            Nmix_h2o=1.e6 ! default value
            call getin("Nmix_h2o",Nmix_h2o)
            write(*,*)" Nmix_h2o = ",Nmix_h2o

            write(*,*)"Number mixing ratio of H2O ice particles:"
            Nmix_h2o_ice=Nmix_h2o ! default value
            call getin("Nmix_h2o_ice",Nmix_h2o_ice)
            write(*,*)" Nmix_h2o_ice = ",Nmix_h2o_ice
         endif

      print*,'exit su_aer_radii'

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


!==================================================================
   subroutine h2o_reffrad(pq,pt,reffrad,nueffrad)
!==================================================================
!     Purpose
!     -------
!     Compute the effective radii of liquid and icy water particles
!
!     Authors
!     -------
!     Jeremy Leconte (2012)
!
!==================================================================
      use radinc_h, only: naerkind
      use watercommon_h, Only: T_h2O_ice_liq,T_h2O_ice_clouds,rhowater,rhowaterice
      use aerosol_mod, only : iaero_h2o
      Implicit none

#include "callkeys.h"
#include "dimensions.h"
#include "dimphys.h"
#include "tracer.h"
#include "comcstfi.h"

      real, intent(in) :: pq(ngridmx,nlayermx,nqmx) !tracer mixing ratios (kg/kg)
      real, intent(in) :: pt(ngridmx,nlayermx)      !temperature (K)
      real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)
      real, intent(inout) :: nueffrad(ngridmx,nlayermx,naerkind) ! dispersion      

      integer :: ig,l
      real zfice ,zrad,zrad_liq,zrad_ice
      real,external :: CBRT            
      

      if (radfixed) then
         do l=1,nlayermx
            do ig=1,ngridmx
               zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
               zfice = MIN(MAX(zfice,0.0),1.0)
               reffrad(ig,l,iaero_h2o)= rad_h2o * (1.-zfice) + rad_h2o_ice * zfice
               nueffrad(ig,l,iaero_h2o) = coef_chaud * (1.-zfice) + coef_froid * zfice
            enddo
         enddo
      else
         do l=1,nlayermx
            do ig=1,ngridmx
               zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
               zfice = MIN(MAX(zfice,0.0),1.0)
               zrad_liq  = CBRT( 3*pq(ig,l,igcm_h2o_ice)/(4*Nmix_h2o*pi*rhowater) )
               zrad_ice  = CBRT( 3*pq(ig,l,igcm_h2o_ice)/(4*Nmix_h2o_ice*pi*rhowaterice) )
               nueffrad(ig,l,iaero_h2o) = coef_chaud * (1.-zfice) + coef_froid * zfice
               zrad = zrad_liq * (1.-zfice) + zrad_ice * zfice
               reffrad(ig,l,iaero_h2o) = min(max(zrad,1.e-6),100.e-6)
               enddo
            enddo      
      end if

   end subroutine h2o_reffrad
!==================================================================


!==================================================================
   subroutine h2o_cloudrad(pql,reffliq,reffice)
!==================================================================
!     Purpose
!     -------
!     Compute the effective radii of liquid and icy water particles
!
!     Authors
!     -------
!     Jeremy Leconte (2012)
!
!==================================================================
      use watercommon_h, Only: rhowater,rhowaterice
      Implicit none

#include "callkeys.h"
#include "dimensions.h"
#include "dimphys.h"
#include "tracer.h"
#include "comcstfi.h"

      real, intent(in) :: pql(ngridmx,nlayermx) !condensed water mixing ratios (kg/kg)
      real, intent(out) :: reffliq(ngridmx,nlayermx),reffice(ngridmx,nlayermx)     !liquid and ice water particle radii (K)

      real,external :: CBRT            
      

      if (radfixed) then
         reffliq(1:ngridmx,1:nlayermx)= rad_h2o
         reffice(1:ngridmx,1:nlayermx)= rad_h2o_ice
      else
         reffliq(1:ngridmx,1:nlayermx)  = CBRT( 3*pql(1:ngridmx,1:nlayermx)/(4*Nmix_h2o*pi*rhowater) )
         reffliq(1:ngridmx,1:nlayermx)  = min(max(reffliq(1:ngridmx,1:nlayermx),1.e-6),100.e-6)
         reffice(1:ngridmx,1:nlayermx)  = CBRT( 3*pql(1:ngridmx,1:nlayermx)/(4*Nmix_h2o_ice*pi*rhowaterice) )
         reffice(1:ngridmx,1:nlayermx)  = min(max(reffice(1:ngridmx,1:nlayermx),1.e-6),100.e-6)
      end if

   end subroutine h2o_cloudrad
!==================================================================



!==================================================================
   subroutine co2_reffrad(pq,reffrad)
!==================================================================
!     Purpose
!     -------
!     Compute the effective radii of co2 ice particles
!
!     Authors
!     -------
!     Jeremy Leconte (2012)
!
!==================================================================
      use radinc_h, only: naerkind
      use aerosol_mod, only : iaero_co2
      Implicit none

#include "callkeys.h"
#include "dimensions.h"
#include "dimphys.h"
#include "tracer.h"
#include "comcstfi.h"

      real, intent(in) :: pq(ngridmx,nlayermx,nqmx) !tracer mixing ratios (kg/kg)
      real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)

      integer :: ig,l
      real :: zrad   
      real,external :: CBRT            
            
      

      if (radfixed) then
         reffrad(1:ngridmx,1:nlayermx,iaero_co2) = 5.e-5 ! CO2 ice
      else
         do l=1,nlayermx
            do ig=1,ngridmx
               zrad = CBRT( 3*pq(ig,l,igcm_co2_ice)/(4*Nmix_co2*pi*rho_co2) )
               reffrad(ig,l,iaero_co2) = min(max(zrad,1.e-6),100.e-6)
            enddo
         enddo      
      end if

   end subroutine co2_reffrad
!==================================================================



!==================================================================
   subroutine dust_reffrad(reffrad)
!==================================================================
!     Purpose
!     -------
!     Compute the effective radii of dust particles
!
!     Authors
!     -------
!     Jeremy Leconte (2012)
!
!==================================================================
      use radinc_h, only: naerkind
      use aerosol_mod, only : iaero_dust
      Implicit none

#include "callkeys.h"
#include "dimensions.h"
#include "dimphys.h"

      real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)
            
      reffrad(1:ngridmx,1:nlayermx,iaero_dust) = 2.e-6 ! dust

   end subroutine dust_reffrad
!==================================================================


!==================================================================
   subroutine h2so4_reffrad(reffrad)
!==================================================================
!     Purpose
!     -------
!     Compute the effective radii of h2so4 particles
!
!     Authors
!     -------
!     Jeremy Leconte (2012)
!
!==================================================================
      use radinc_h, only: naerkind
      use aerosol_mod, only : iaero_h2so4
      Implicit none

#include "callkeys.h"
#include "dimensions.h"
#include "dimphys.h"

      real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)
                
      reffrad(1:ngridmx,1:nlayermx,iaero_h2so4) = 1.e-6 ! h2so4

   end subroutine h2so4_reffrad
!==================================================================



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