rrtmgp_sw_cloud_optics.F90

 
 
module rrtmgp_sw_cloud_optics
  use mo_rte_kind,              only: wl, wp
  use mo_cloud_optics_rrtmgp,   only: ty_cloud_optics => ty_cloud_optics_rrtmgp
  use rrtmgp_sw_gas_optics,     only: sw_gas_props
  use radiation_tools,          only: check_error_msg
  use netcdf
  use mpi_f08
 
  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(wp), 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(wp), dimension(:,:,:), allocatable :: &
       lut_exticesw,          & !< LUT shortwave ice extinction coefficient
       lut_ssaicesw,          & !< LUT shortwave ice single scattering albedo
       lut_asyicesw
  real(wp) :: &
       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(wp),parameter :: &
       a0r = 3.07e-3,        & !
       a0s = 0.0,            & !
       a1s = 1.5               !  
  real(wp),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,        &
       nrghice, mpicomm, mpirank, mpiroot, is_initialized, errmsg, errflg)
 
    ! Inputs
    character(len=128),intent(in) :: &
         rrtmgp_root_dir,    & !< RTE-RRTMGP root directory
         rrtmgp_sw_file_clouds
    integer, intent(inout) :: &
         nrghice
    type(mpi_comm), intent(in) :: &
         mpicomm
    integer, intent(in) :: &
         mpirank,            & !< Current MPI rank
         mpiroot
 
    ! Outputs
    logical,          intent(inout) :: &
         is_initialized
    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
 
    if (is_initialized) return
 
    ! 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)
    !
    ! #######################################################################################
    if (mpirank .eq. mpiroot) then
       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) 
    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)
    
    ! 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
    call mpi_bcast(nrghice_fromfilesw, 1, mpi_integer, mpiroot, mpicomm, mpierr)
 
    ! #######################################################################################
    !
    ! Allocate space for arrays...
    ! (ALL processors)
    !
    ! #######################################################################################
    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))
    allocate(band_limscldsw(2,nbandsw))
 
    ! #######################################################################################
    !
    ! Read in data ...
    ! (ONLY master processor(0), if MPI enabled) 
    !
    ! #######################################################################################
    if (mpirank .eq. mpiroot) 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,'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,'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)
 
       ! Close file
       status = nf90_close(ncid)       
 
    endif ! Master process
 
    ! Other processors waiting...
    call mpi_barrier(mpicomm, mpierr)
 
    ! #######################################################################################
    !
    ! Broadcast data... 
    ! (ALL processors)
    !
    ! #######################################################################################
 
    ! Real scalars
#ifdef RTE_USE_SP
    call mpi_bcast(radliq_lwrsw, 1, mpi_real,             mpiroot, mpicomm, mpierr)
    call mpi_bcast(radliq_uprsw, 1, mpi_real,             mpiroot, mpicomm, mpierr)
    call mpi_bcast(radice_lwrsw, 1, mpi_real,             mpiroot, mpicomm, mpierr)
    call mpi_bcast(radice_uprsw, 1, mpi_real,             mpiroot, mpicomm, mpierr)
#else
    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)
#endif
 
    ! Real arrays
#ifdef RTE_USE_SP
    call mpi_bcast(band_limscldsw, size(band_limscldsw),  mpi_real,             mpiroot, mpicomm, mpierr)
    call mpi_bcast(lut_extliqsw,   size(lut_extliqsw),    mpi_real,             mpiroot, mpicomm, mpierr)
    call mpi_bcast(lut_ssaliqsw,   size(lut_ssaliqsw),    mpi_real,             mpiroot, mpicomm, mpierr)
    call mpi_bcast(lut_asyliqsw,   size(lut_asyliqsw),    mpi_real,             mpiroot, mpicomm, mpierr)
    call mpi_bcast(lut_exticesw,   size(lut_exticesw),    mpi_real,             mpiroot, mpicomm, mpierr)
    call mpi_bcast(lut_ssaicesw,   size(lut_ssaicesw),    mpi_real,             mpiroot, mpicomm, mpierr)
    call mpi_bcast(lut_asyicesw,   size(lut_asyicesw),    mpi_real,             mpiroot, mpicomm, mpierr)
#else
    call mpi_bcast(band_limscldsw, size(band_limscldsw),  mpi_double_precision, mpiroot, mpicomm, mpierr)
    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
 
    ! #######################################################################################
    !   
    ! Initialize RRTMGP DDT's...
    !
    ! #######################################################################################
    call check_error_msg('sw_cloud_optics_init',sw_cloud_props%load(band_limscldsw,         &
            radliq_lwrsw, radliq_uprsw, radice_lwrsw, radice_uprsw,                         &
            lut_extliqsw, lut_ssaliqsw, lut_asyliqsw,                                       &
            lut_exticesw, lut_ssaicesw, lut_asyicesw))
 
    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/)
 
    is_initialized = .true.
 
  end subroutine rrtmgp_sw_cloud_optics_init
end module rrtmgp_sw_cloud_optics