rrtmgp_sw_cloud_optics.F90

 
 
module rrtmgp_sw_cloud_optics
  use machine,                  only: kind_phys
  use mo_rte_kind,              only: wl
  use mo_cloud_optics,          only: ty_cloud_optics
  use rrtmgp_sw_gas_optics,     only: sw_gas_props
  use radiation_tools,          only: check_error_msg
  use netcdf
#ifdef MPI
  use mpi_f08
#endif
 
  implicit none
  
  type(ty_cloud_optics) :: sw_cloud_props 
  integer :: &
       nrghice_fromfilesw, nbandsw, nsize_liqsw, nsize_icesw, nsizeregsw, &
       ncoeff_extsw, ncoeff_ssa_gsw, nboundsw, npairssw
  real(kind_phys) :: &
       radliq_facsw,          & !< Factor for calculating LUT interpolation indices for liquid   
       radice_facsw
  real(kind_phys), dimension(:,:), allocatable :: &
       lut_extliqsw,          & !< LUT shortwave liquid extinction coefficient  
       lut_ssaliqsw,          & !< LUT shortwave liquid single scattering albedo   
       lut_asyliqsw,          & !< LUT shortwave liquid asymmetry parameter  
       band_limscldsw
  real(kind_phys), dimension(:,:,:), allocatable :: &
       lut_exticesw,          & !< LUT shortwave ice extinction coefficient
       lut_ssaicesw,          & !< LUT shortwave ice single scattering albedo
       lut_asyicesw
  real(kind_phys), dimension(:), allocatable :: &
       pade_sizereg_extliqsw, & !< Particle size regime boundaries for shortwave liquid extinction 
 
       pade_sizereg_ssaliqsw, & 
       pade_sizereg_asyliqsw, & 
       pade_sizereg_exticesw, & 
       pade_sizereg_ssaicesw, & 
       pade_sizereg_asyicesw    
  real(kind_phys), dimension(:,:,:), allocatable :: &
       pade_extliqsw,         & !< PADE coefficients for shortwave liquid extinction
       pade_ssaliqsw,         & !< PADE coefficients for shortwave liquid single scattering albedo
       pade_asyliqsw
  real(kind_phys), dimension(:,:,:,:), allocatable :: &
       pade_exticesw,         & !< PADE coefficients for shortwave ice extinction
       pade_ssaicesw,         & !< PADE coefficients for shortwave ice single scattering albedo
       pade_asyicesw
  real(kind_phys) :: &
       radliq_lwrsw,          & !< Liquid particle size lower bound for LUT interpolation
       radliq_uprsw,          & !< Liquid particle size upper bound for LUT interpolation
       radice_lwrsw,          & !< Ice particle    size upper bound for LUT interpolation
       radice_uprsw
 
  ! Parameters used for rain and snow(+groupel) RRTMGP cloud-optics. *NOTE* Same as in RRTMG
  ! Need to document these magic numbers below.
  real(kind_phys),parameter :: &
       a0r = 3.07e-3,        & !
       a0s = 0.0,            & !
       a1s = 1.5               !  
  real(kind_phys),dimension(:),allocatable :: b0r,b0s,b1s,c0r,c0s
 
contains
  ! ######################################################################################
  ! SUBROUTINE sw_cloud_optics_init
  ! ######################################################################################
  subroutine rrtmgp_sw_cloud_optics_init( rrtmgp_root_dir, rrtmgp_sw_file_clouds,        &
       doGP_cldoptics_PADE, doGP_cldoptics_LUT, nrghice, mpicomm, mpirank, mpiroot,      &
       errmsg, errflg)
 
    ! Inputs
    character(len=128),intent(in) :: &
         rrtmgp_root_dir,    & !< RTE-RRTMGP root directory
         rrtmgp_sw_file_clouds
    logical, intent(in) :: &
         doGP_cldoptics_PADE,& !< Use RRTMGP cloud-optics: PADE approximation?
         doGP_cldoptics_LUT
    integer, intent(inout) :: &
         nrghice
    type(mpi_comm), intent(in) :: &
         mpicomm
    integer, intent(in) :: &
         mpirank,            & !< Current MPI rank
         mpiroot
 
    ! Outputs
    character(len=*), intent(out) :: &
         errmsg
    integer,          intent(out) :: &
         errflg
    
    ! Local variables
    integer :: status,ncid,dimid,varID,mpierr
    character(len=264) :: sw_cloud_props_file
 
    ! Initialize
    errmsg = ''
    errflg = 0
 
    ! Filenames are set in the physics_nml
    sw_cloud_props_file = trim(rrtmgp_root_dir)//trim(rrtmgp_sw_file_clouds)
 
    ! #######################################################################################
    !
    ! Read dimensions for shortwave cloud-optics fields...
    ! (ONLY master processor(0), if MPI enabled)
    !
    ! #######################################################################################
#ifdef MPI
    if (mpirank .eq. mpiroot) then
#endif
       write (*,*) 'Reading RRTMGP shortwave cloud-optics metadata ... '
 
       ! Open file
       status = nf90_open(trim(sw_cloud_props_file), nf90_nowrite, ncid)
 
       ! Read dimensions
       status = nf90_inq_dimid(ncid, 'nband', dimid)
       status = nf90_inquire_dimension(ncid, dimid, len=nbandsw)
       status = nf90_inq_dimid(ncid, 'nrghice', dimid)
       status = nf90_inquire_dimension(ncid, dimid, len=nrghice_fromfilesw)
       status = nf90_inq_dimid(ncid, 'nsize_liq', dimid)
       status = nf90_inquire_dimension(ncid, dimid, len=nsize_liqsw)
       status = nf90_inq_dimid(ncid, 'nsize_ice', dimid)
       status = nf90_inquire_dimension(ncid, dimid, len=nsize_icesw)
       status = nf90_inq_dimid(ncid, 'nsizereg', dimid)
       status = nf90_inquire_dimension(ncid, dimid, len=nsizeregsw)
       status = nf90_inq_dimid(ncid, 'ncoeff_ext', dimid)
       status = nf90_inquire_dimension(ncid, dimid, len=ncoeff_extsw)
       status = nf90_inq_dimid(ncid, 'ncoeff_ssa_g', dimid)
       status = nf90_inquire_dimension(ncid, dimid, len=ncoeff_ssa_gsw)
       status = nf90_inq_dimid(ncid, 'nbound', dimid)
       status = nf90_inquire_dimension(ncid, dimid, len=nboundsw)
       status = nf90_inq_dimid(ncid, 'pair', dimid)
       status = nf90_inquire_dimension(ncid, dimid, len=npairssw) 
#ifdef MPI
    endif ! On master processor
 
    ! Other processors waiting...
    call mpi_barrier(mpicomm, mpierr)
 
    ! #######################################################################################
    !
    ! Broadcast dimensions...
    ! (ALL processors)
    !
    ! #######################################################################################
    call mpi_bcast(nbandsw,            1, mpi_integer, mpiroot, mpicomm, mpierr)
    call mpi_bcast(nsize_liqsw,        1, mpi_integer, mpiroot, mpicomm, mpierr)
    call mpi_bcast(nsize_icesw,        1, mpi_integer, mpiroot, mpicomm, mpierr)
    call mpi_bcast(nsizeregsw,         1, mpi_integer, mpiroot, mpicomm, mpierr)
    call mpi_bcast(ncoeff_extsw,       1, mpi_integer, mpiroot, mpicomm, mpierr)
    call mpi_bcast(ncoeff_ssa_gsw,     1, mpi_integer, mpiroot, mpicomm, mpierr)
    call mpi_bcast(nboundsw,           1, mpi_integer, mpiroot, mpicomm, mpierr)
    call mpi_bcast(npairssw,           1, mpi_integer, mpiroot, mpicomm, mpierr)
#endif
    
    ! Has the number of ice-roughnes categories been provided from the namelist?
    ! If so, override nrghice from cloud-optics file
    if (nrghice .ne. 0) nrghice_fromfilesw = nrghice
#ifdef MPI
    call mpi_bcast(nrghice_fromfilesw, 1, mpi_integer, mpiroot, mpicomm, mpierr)
#endif
 
    ! #######################################################################################
    !
    ! Allocate space for arrays...
    ! (ALL processors)
    !
    ! #######################################################################################
    if (dogp_cldoptics_lut) then
       allocate(lut_extliqsw(nsize_liqsw, nbandsw))
       allocate(lut_ssaliqsw(nsize_liqsw, nbandsw))
       allocate(lut_asyliqsw(nsize_liqsw, nbandsw))
       allocate(lut_exticesw(nsize_icesw, nbandsw, nrghice_fromfilesw))
       allocate(lut_ssaicesw(nsize_icesw, nbandsw, nrghice_fromfilesw))
       allocate(lut_asyicesw(nsize_icesw, nbandsw, nrghice_fromfilesw))
    endif
    if (dogp_cldoptics_pade) then
       allocate(pade_extliqsw(nbandsw, nsizeregsw,  ncoeff_extsw ))
       allocate(pade_ssaliqsw(nbandsw, nsizeregsw,  ncoeff_ssa_gsw))
       allocate(pade_asyliqsw(nbandsw, nsizeregsw,  ncoeff_ssa_gsw))
       allocate(pade_exticesw(nbandsw, nsizeregsw,  ncoeff_extsw,   nrghice_fromfilesw))
       allocate(pade_ssaicesw(nbandsw, nsizeregsw,  ncoeff_ssa_gsw, nrghice_fromfilesw))
       allocate(pade_asyicesw(nbandsw, nsizeregsw,  ncoeff_ssa_gsw, nrghice_fromfilesw))
       allocate(pade_sizereg_extliqsw(nboundsw))
       allocate(pade_sizereg_ssaliqsw(nboundsw))
       allocate(pade_sizereg_asyliqsw(nboundsw))
       allocate(pade_sizereg_exticesw(nboundsw))
       allocate(pade_sizereg_ssaicesw(nboundsw))
       allocate(pade_sizereg_asyicesw(nboundsw))
    endif
    allocate(band_limscldsw(2,nbandsw))
 
    ! #######################################################################################
    !
    ! Read in data ...
    ! (ONLY master processor(0), if MPI enabled) 
    !
    ! #######################################################################################
#ifdef MPI
    if (mpirank .eq. mpiroot) then
#endif 
       if (dogp_cldoptics_lut) then
          write (*,*) 'Reading RRTMGP shortwave cloud data (LUT) ... '
          status = nf90_inq_varid(ncid,'radliq_lwr',varid)
          status = nf90_get_var(ncid,varid,radliq_lwrsw)
          status = nf90_inq_varid(ncid,'radliq_upr',varid)
          status = nf90_get_var(ncid,varid,radliq_uprsw)
          status = nf90_inq_varid(ncid,'radliq_fac',varid)
          status = nf90_get_var(ncid,varid,radliq_facsw)
          status = nf90_inq_varid(ncid,'radice_lwr',varid)
          status = nf90_get_var(ncid,varid,radice_lwrsw)
          status = nf90_inq_varid(ncid,'radice_upr',varid)
          status = nf90_get_var(ncid,varid,radice_uprsw)
          status = nf90_inq_varid(ncid,'radice_fac',varid)
          status = nf90_get_var(ncid,varid,radice_facsw)
          status = nf90_inq_varid(ncid,'lut_extliq',varid)
          status = nf90_get_var(ncid,varid,lut_extliqsw)
          status = nf90_inq_varid(ncid,'lut_ssaliq',varid)
          status = nf90_get_var(ncid,varid,lut_ssaliqsw)
          status = nf90_inq_varid(ncid,'lut_asyliq',varid)
          status = nf90_get_var(ncid,varid,lut_asyliqsw)
          status = nf90_inq_varid(ncid,'lut_extice',varid)
          status = nf90_get_var(ncid,varid,lut_exticesw)
          status = nf90_inq_varid(ncid,'lut_ssaice',varid)
          status = nf90_get_var(ncid,varid,lut_ssaicesw)
          status = nf90_inq_varid(ncid,'lut_asyice',varid)
          status = nf90_get_var(ncid,varid,lut_asyicesw)
          status = nf90_inq_varid(ncid,'bnd_limits_wavenumber',varid)
          status = nf90_get_var(ncid,varid,band_limscldsw)
       endif
       if (dogp_cldoptics_pade) then
          write (*,*) 'Reading RRTMGP shortwave cloud data (PADE) ... '
          status = nf90_inq_varid(ncid,'radliq_lwr',varid)
          status = nf90_get_var(ncid,varid,radliq_lwrsw)
          status = nf90_inq_varid(ncid,'radliq_upr',varid)
          status = nf90_get_var(ncid,varid,radliq_uprsw)
          status = nf90_inq_varid(ncid,'radliq_fac',varid)
          status = nf90_get_var(ncid,varid,radliq_facsw)
          status = nf90_inq_varid(ncid,'radice_lwr',varid)
          status = nf90_get_var(ncid,varid,radice_lwrsw)
          status = nf90_inq_varid(ncid,'radice_upr',varid)
          status = nf90_get_var(ncid,varid,radice_uprsw)
          status = nf90_inq_varid(ncid,'radice_fac',varid)
          status = nf90_get_var(ncid,varid,radice_facsw)
          status = nf90_inq_varid(ncid,'pade_extliq',varid)
          status = nf90_get_var(ncid,varid,pade_extliqsw)
          status = nf90_inq_varid(ncid,'pade_ssaliq',varid)
          status = nf90_get_var(ncid,varid,pade_ssaliqsw)
          status = nf90_inq_varid(ncid,'pade_asyliq',varid)
          status = nf90_get_var(ncid,varid,pade_asyliqsw)
          status = nf90_inq_varid(ncid,'pade_extice',varid)
          status = nf90_get_var(ncid,varid,pade_exticesw)
          status = nf90_inq_varid(ncid,'pade_ssaice',varid)
          status = nf90_get_var(ncid,varid,pade_ssaicesw)
          status = nf90_inq_varid(ncid,'pade_asyice',varid)
          status = nf90_get_var(ncid,varid,pade_asyicesw)
          status = nf90_inq_varid(ncid,'pade_sizreg_extliq',varid)
          status = nf90_get_var(ncid,varid,pade_sizereg_extliqsw)
          status = nf90_inq_varid(ncid,'pade_sizreg_ssaliq',varid)
          status = nf90_get_var(ncid,varid,pade_sizereg_ssaliqsw)
          status = nf90_inq_varid(ncid,'pade_sizreg_asyliq',varid)
          status = nf90_get_var(ncid,varid,pade_sizereg_asyliqsw)
          status = nf90_inq_varid(ncid,'pade_sizreg_extice',varid)
          status = nf90_get_var(ncid,varid,pade_sizereg_exticesw)
          status = nf90_inq_varid(ncid,'pade_sizreg_ssaice',varid)
          status = nf90_get_var(ncid,varid,pade_sizereg_ssaicesw)
          status = nf90_inq_varid(ncid,'pade_sizreg_asyice',varid)
          status = nf90_get_var(ncid,varid,pade_sizereg_asyicesw)
          status = nf90_inq_varid(ncid,'bnd_limits_wavenumber',varid)
          status = nf90_get_var(ncid,varid,band_limscldsw)
       endif
 
       ! Close file
       status = nf90_close(ncid)       
 
#ifdef MPI
    endif ! Master process
 
    ! Other processors waiting...
    call mpi_barrier(mpicomm, mpierr)
 
    ! #######################################################################################
    !
    ! Broadcast data... 
    ! (ALL processors)
    !
    ! #######################################################################################
 
    ! Real scalars
    call mpi_bcast(radliq_facsw, 1, mpi_double_precision, mpiroot, mpicomm, mpierr)
    call mpi_bcast(radice_facsw, 1, mpi_double_precision, mpiroot, mpicomm, mpierr)
    call mpi_bcast(radliq_lwrsw, 1, mpi_double_precision, mpiroot, mpicomm, mpierr)
    call mpi_bcast(radliq_uprsw, 1, mpi_double_precision, mpiroot, mpicomm, mpierr)
    call mpi_bcast(radice_lwrsw, 1, mpi_double_precision, mpiroot, mpicomm, mpierr)
    call mpi_bcast(radice_uprsw, 1, mpi_double_precision, mpiroot, mpicomm, mpierr)
 
    ! Real arrays
    call mpi_bcast(band_limscldsw, size(band_limscldsw),  &
         mpi_double_precision, mpiroot, mpicomm, mpierr)
    if (dogp_cldoptics_lut) then
       call mpi_bcast(lut_extliqsw, size(lut_extliqsw),   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(lut_ssaliqsw, size(lut_ssaliqsw),   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(lut_asyliqsw, size(lut_asyliqsw),   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(lut_exticesw, size(lut_exticesw),   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(lut_ssaicesw, size(lut_ssaicesw),   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(lut_asyicesw, size(lut_asyicesw),   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
    endif
    if (dogp_cldoptics_pade) then
       call mpi_bcast(pade_extliqsw, size(pade_extliqsw),                   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(pade_ssaliqsw, size(pade_ssaliqsw),                   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(pade_asyliqsw, size(pade_asyliqsw),                   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(pade_exticesw, size(pade_exticesw),                   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(pade_ssaicesw, size(pade_ssaicesw),                   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(pade_asyicesw, size(pade_asyicesw),                   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(pade_sizereg_extliqsw, size(pade_sizereg_extliqsw),   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(pade_sizereg_ssaliqsw, size(pade_sizereg_ssaliqsw),   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(pade_sizereg_asyliqsw, size(pade_sizereg_asyliqsw),   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(pade_sizereg_exticesw, size(pade_sizereg_exticesw),   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(pade_sizereg_ssaicesw, size(pade_sizereg_ssaicesw),   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
       call mpi_bcast(pade_sizereg_asyicesw, size(pade_sizereg_asyicesw),   &
            mpi_double_precision, mpiroot, mpicomm, mpierr)
    endif
#endif
 
    ! #######################################################################################
    !   
    ! Initialize RRTMGP DDT's...
    !
    ! #######################################################################################
    if (dogp_cldoptics_lut) then
       call check_error_msg('sw_cloud_optics_init',sw_cloud_props%load(band_limscldsw,      &
            radliq_lwrsw, radliq_uprsw, radliq_facsw, radice_lwrsw, radice_uprsw,           &
            radice_facsw, lut_extliqsw, lut_ssaliqsw, lut_asyliqsw, lut_exticesw,           &
            lut_ssaicesw, lut_asyicesw))
    endif
 
    if (dogp_cldoptics_pade) then
       call check_error_msg('sw_cloud_optics_init', sw_cloud_props%load(band_limscldsw,     &
            pade_extliqsw, pade_ssaliqsw, pade_asyliqsw, pade_exticesw, pade_ssaicesw,      &
            pade_asyicesw, pade_sizereg_extliqsw, pade_sizereg_ssaliqsw,                    &
            pade_sizereg_asyliqsw, pade_sizereg_exticesw, pade_sizereg_ssaicesw,            &
            pade_sizereg_asyicesw))
    endif
 
    call check_error_msg('sw_cloud_optics_init',sw_cloud_props%set_ice_roughness(nrghice_fromfilesw))
 
    ! Initialize coefficients for rain and snow(+groupel) cloud optics
    allocate(b0r(sw_cloud_props%get_nband()),b0s(sw_cloud_props%get_nband()), &
             b1s(sw_cloud_props%get_nband()),c0r(sw_cloud_props%get_nband()), &
             c0s(sw_cloud_props%get_nband()))
    b0r = (/0.496, 0.466,   0.437,   0.416, 0.391, 0.374, 0.352,    &
            0.183, 0.048,   0.012,   0.000, 0.000, 0.000, 0.000/)
    b0s = (/0.460, 0.460,   0.460,   0.460, 0.460, 0.460, 0.460,    &
            0.460, 0.000,   0.000,   0.000, 0.000, 0.000, 0.000/)
    b1s = (/0.000, 0.000,   0.000,   0.000, 0.000, 0.000, 0.000,    &
            0.000, 1.62e-5, 1.62e-5, 0.000, 0.000, 0.000, 0.000/)
    c0r = (/0.980, 0.975,   0.965,   0.960, 0.955, 0.952, 0.950,    &
            0.944, 0.894,   0.884,   0.883, 0.883, 0.883, 0.883/)
    c0s = (/0.970, 0.970,   0.970,   0.970, 0.970, 0.970, 0.970,    &
            0.970, 0.970,   0.970,   0.700, 0.700, 0.700, 0.700/)
 
  end subroutine rrtmgp_sw_cloud_optics_init
end module rrtmgp_sw_cloud_optics