text:
!>
!! This module is the interface between nwp_nh_interface to the radiation schemes
!! (RRTM or Ritter-Geleyn).
!!
!! @author Thorsten Reinhardt, AGeoBw, Offenbach
!!
!! @par Revision History
!! Initial release by Thorsten Reinhardt, AGeoBw, Offenbach (2011-01-13)
!!
!! @par Copyright and License
!!
!! This code is subject to the DWD and MPI-M-Software-License-Agreement in
!! its most recent form.
!! Please see the file LICENSE in the root of the source tree for this code.
!! Where software is supplied by third parties, it is indicated in the
!! headers of the routines.
!!

!----------------------------
#include "omp_definitions.inc"
!----------------------------
#if defined __xlC__
@PROCESS SPILL(1058)
#endif
MODULE mo_nwp_rrtm_interface

  USE mo_atm_phy_nwp_config,   ONLY: atm_phy_nwp_config, iprog_aero, icpl_aero_conv
  USE mo_nwp_tuning_config,    ONLY: tune_dust_abs
  USE mo_grid_config,          ONLY: l_limited_area
  USE mo_exception,            ONLY: message,  finish, message_text
  USE mo_ext_data_types,       ONLY: t_external_data
  USE mo_parallel_config,      ONLY: nproma, p_test_run, test_parallel_radiation
  USE mo_run_config,           ONLY: msg_level, iqv, iqc, iqi
  USE mo_impl_constants,       ONLY: min_rlcell_int, io3_ape, nexlevs_rrg_vnest, &
                                     iss, iorg, ibc, iso4, idu, MAX_CHAR_LENGTH
  USE mo_impl_constants_grf,   ONLY: grf_bdywidth_c, grf_ovlparea_start_c
  USE mo_kind,                 ONLY: wp
  USE mo_loopindices,          ONLY: get_indices_c
  USE mo_nwp_lnd_types,        ONLY: t_lnd_prog
  USE mo_model_domain,         ONLY: t_patch, p_patch_local_parent
  USE mo_phys_nest_utilities,  ONLY: upscale_rad_input, downscale_rad_output
  USE mo_nonhydro_types,       ONLY: t_nh_diag
  USE mo_nwp_phy_types,        ONLY: t_nwp_phy_diag
  USE mo_o3_util,              ONLY: calc_o3_clim, calc_o3_gems
  USE mo_radiation,            ONLY: radiation, radiation_nwp
  USE mo_radiation_config,     ONLY: irad_o3, irad_aero
  USE mo_radiation_rg_par,     ONLY: aerdis
  USE mo_aerosol_util,         ONLY: tune_dust
  USE mo_lrtm_par,             ONLY: nbndlw
  USE mo_sync,                 ONLY: global_max, global_min

  USE mo_rrtm_data_interface,  ONLY: t_rrtm_data, recv_rrtm_input, send_rrtm_output
  USE mo_timer,                ONLY: timer_start, timer_stop, timers_level,    &
    &                                timer_radiaton_recv, timer_radiaton_comp, &
    &                                timer_radiaton_send, timer_preradiaton
  USE mtime,                     ONLY: datetime, newDatetime, deallocateDatetime
  USE mo_bcs_time_interpolation, ONLY: t_time_interpolation_weights,         &
    &                                  calculate_time_interpolation_weights

#ifdef COUP_OAS_ICON
  USE oas_icon_define,         ONLY: oas_rcv_field_icon
#endif

  IMPLICIT NONE

  PRIVATE

  !> module name string
  CHARACTER(LEN=*), PARAMETER :: modname = 'mo_nwp_rrtm_interface'

  PUBLIC :: nwp_ozon_aerosol
  PUBLIC :: nwp_rrtm_radiation
  PUBLIC :: nwp_rrtm_radiation_reduced
  PUBLIC :: nwp_rrtm_radiation_repartition

  REAL(wp), PARAMETER::  &
    & zaeops = 0.05_wp,   &
    & zaeopl = 0.2_wp,    &
    & zaeopu = 0.1_wp,    &
    & zaeopd = 1.9_wp,    &
    & ztrpt  = 30.0_wp,   &
    & ztrbga = 0.03_wp  / (101325.0_wp - 19330.0_wp), &
    & zvobga = 0.007_wp /  19330.0_wp , &
    & zstbga = 0.045_wp /  19330.0_wp!, &
!      & zaeadk(1:3) = (/0.3876E-03_wp,0.6693E-02_wp,0.8563E-03_wp/)

CONTAINS

  !---------------------------------------------------------------------------------------
  !>
  !! @par Revision History
  !! Initial release by Thorsten Reinhardt, AGeoBw, Offenbach (2011-01-13)
  !!
  SUBROUTINE nwp_ozon_aerosol ( p_sim_time, mtime_datetime, pt_patch, ext_data, &
    & pt_diag,prm_diag,zaeq1,zaeq2,zaeq3,zaeq4,zaeq5,zduo3 )

!    CHARACTER(len=MAX_CHAR_LENGTH), PARAMETER::  &
!      &  routine = 'mo_nwp_rad_interface:'

    REAL(wp),INTENT(in)         :: p_sim_time

    TYPE(datetime), POINTER, INTENT(in)    :: mtime_datetime
    TYPE(t_patch), TARGET,   INTENT(in)    :: pt_patch     !<grid/patch info.
    TYPE(t_external_data),   INTENT(inout) :: ext_data
    TYPE(t_nh_diag), TARGET, INTENT(in)    :: pt_diag     !<the diagnostic variables
    TYPE(t_nwp_phy_diag),    INTENT(inout) :: prm_diag

    REAL(wp), INTENT(out) :: &
      & zaeq1(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq2(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq3(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq4(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq5(nproma,pt_patch%nlev,pt_patch%nblks_c)

    ! for Ritter-Geleyn radiation:
    REAL(wp), OPTIONAL, INTENT(out) :: zduo3(nproma,pt_patch%nlev,pt_patch%nblks_c)

    ! for ozone:
    REAL(wp):: &
      & zptop32(nproma,pt_patch%nblks_c), &
      & zo3_hm (nproma,pt_patch%nblks_c), &
      & zo3_top (nproma,pt_patch%nblks_c), &
      & zpbot32(nproma,pt_patch%nblks_c), &
      & zo3_bot (nproma,pt_patch%nblks_c)
    ! for aerosols:
    REAL(wp):: &
      & zsign(nproma,pt_patch%nlevp1), &
      & zvdaes(nproma,pt_patch%nlevp1), &
      & zvdael(nproma,pt_patch%nlevp1), &
      & zvdaeu(nproma,pt_patch%nlevp1), &
      & zvdaed(nproma,pt_patch%nlevp1), &
      & zaetr_top(nproma,pt_patch%nblks_c), zaetr_bot, zaetr,       &
      & zaeqdo   (nproma,pt_patch%nblks_c), zaeqdn,                 &
      & zaequo   (nproma,pt_patch%nblks_c), zaequn,                 &
      & zaeqlo   (nproma,pt_patch%nblks_c), zaeqln,                 &
      & zaeqso   (nproma,pt_patch%nblks_c), zaeqsn, zw

    ! Local scalars:
    INTEGER:: jc,jk,jb
    INTEGER:: jg                !domain id
    INTEGER:: nlev, nlevp1      !< number of full and half levels

    INTEGER:: rl_start, rl_end
    INTEGER:: i_startblk, i_endblk    !> blocks
    INTEGER:: i_startidx, i_endidx    !< slices
    INTEGER:: i_nchdom                !< domain index

    INTEGER:: imo1,imo2 !for Tegen aerosol time interpolation

    REAL(wp) :: wfac, ncn_bg

    TYPE(datetime), POINTER :: current_time_hours
    TYPE(t_time_interpolation_weights) :: current_time_interpolation_weights

    i_nchdom  = MAX(1,pt_patch%n_childdom)
    jg        = pt_patch%id

    ! number of vertical levels
    nlev   = pt_patch%nlev
    nlevp1 = pt_patch%nlevp1

    IF (timers_level > 3) CALL timer_start(timer_preradiaton)

    !-------------------------------------------------------------------------
    !> Radiation setup
    !-------------------------------------------------------------------------

    !O3
    SELECT CASE (irad_o3)
    CASE(io3_ape)
      ! APE ozone: do nothing since everything is already
      ! set in mo_nwp_phy_init
    CASE (6)
      CALL calc_o3_clim(                             &
        & kbdim      = nproma,                       & ! in
        & jg         = jg,                           &
        & p_inc_rad  = atm_phy_nwp_config(jg)%dt_rad,& ! in
        & z_sim_time = p_sim_time,                   & ! in
        & pt_patch   = pt_patch,                     & ! in
        & zvio3      = prm_diag%vio3,                & !inout
        & zhmo3      = prm_diag%hmo3  )                !inout
    CASE (7,9,79)
      CALL calc_o3_gems(pt_patch,mtime_datetime,pt_diag,prm_diag,ext_data)
    CASE(10)
      !CALL message('mo_nwp_rg_interface:irad_o3=10', &
      !  &          'Ozone used for radiation is calculated by ART')
    END SELECT

    IF ( irad_aero == 6  .OR. irad_aero == 9) THEN
      current_time_hours => newDatetime(mtime_datetime)
      current_time_hours%time%minute = 0
      current_time_hours%time%second = 0
      current_time_hours%time%ms = 0
      current_time_interpolation_weights = calculate_time_interpolation_weights(current_time_hours)
      imo1 = current_time_interpolation_weights%month1
      imo2 = current_time_interpolation_weights%month2
      zw = current_time_interpolation_weights%weight2
      CALL deallocateDatetime(current_time_hours)
    ENDIF

    rl_start = 1
    rl_end   = min_rlcell_int

    i_startblk = pt_patch%cells%start_blk(rl_start,1)
    i_endblk   = pt_patch%cells%end_blk(rl_end,i_nchdom)

!$OMP PARALLEL
!$OMP DO PRIVATE(jb,jc,jk,i_endidx, &
!$OMP       zsign,zvdaes, zvdael, zvdaeu, zvdaed, zaeqsn, zaeqln, &
!$OMP zaequn,zaeqdn,zaetr_bot,zaetr,wfac,ncn_bg )  ICON_OMP_DEFAULT_SCHEDULE
    DO jb = i_startblk, i_endblk

      CALL get_indices_c(pt_patch, jb, i_startblk, i_endblk, &
        &                        i_startidx, i_endidx, rl_start, rl_end)

      IF ( irad_aero == 5 ) THEN ! Tanre aerosols

        DO jk = 2, nlevp1
          DO jc = 1,i_endidx
            zsign(jc,jk) = pt_diag%pres_ifc(jc,jk,jb) / 101325._wp
          ENDDO
        ENDDO

        ! The routine aerdis is called to recieve some parameters for the vertical
        ! distribution of background aerosol.
        CALL aerdis ( &
          & kbdim  = nproma,      & !in
          & jcs    = 1,           & !in
          & jce    = i_endidx,    & !in
          & klevp1 = nlevp1,      & !in
          & petah  = zsign(1,1),  & !in
          & pvdaes = zvdaes(1,1), & !out
          & pvdael = zvdael(1,1), & !out
          & pvdaeu = zvdaeu(1,1), & !out
          & pvdaed = zvdaed(1,1) )  !out

        ! top level
        DO jc = 1,i_endidx
          zaeqso   (jc,jb) = zaeops*prm_diag%aersea(jc,jb)*zvdaes(jc,1)
          zaeqlo   (jc,jb) = zaeopl*prm_diag%aerlan(jc,jb)*zvdael(jc,1)
          zaequo   (jc,jb) = zaeopu*prm_diag%aerurb(jc,jb)*zvdaeu(jc,1)
          zaeqdo   (jc,jb) = zaeopd*prm_diag%aerdes(jc,jb)*zvdaed(jc,1)
          zaetr_top(jc,jb) = 1.0_wp
        ENDDO

        ! loop over layers
        DO jk = 1,nlev
          DO jc = 1,i_endidx
            zaeqsn         = zaeops*prm_diag%aersea(jc,jb)*zvdaes(jc,jk+1)
            zaeqln         = zaeopl*prm_diag%aerlan(jc,jb)*zvdael(jc,jk+1)
            zaequn         = zaeopu*prm_diag%aerurb(jc,jb)*zvdaeu(jc,jk+1)
            zaeqdn         = zaeopd*prm_diag%aerdes(jc,jb)*zvdaed(jc,jk+1)
            zaetr_bot      = zaetr_top(jc,jb) &
              & * ( MIN (1.0_wp, pt_diag%temp_ifc(jc,jk,jb)/pt_diag%temp_ifc(jc,jk+1,jb)) )**ztrpt

            zaetr          = SQRT(zaetr_bot*zaetr_top(jc,jb))
            zaeq1(jc,jk,jb)= (1._wp-zaetr) &
              & * (ztrbga* pt_diag%dpres_mc(jc,jk,jb)+zaeqln-zaeqlo(jc,jb)+zaeqdn-zaeqdo(jc,jb))
            zaeq2(jc,jk,jb)   = (1._wp-zaetr) * ( zaeqsn-zaeqso(jc,jb) )
            zaeq3(jc,jk,jb)   = (1._wp-zaetr) * ( zaequn-zaequo(jc,jb) )
            zaeq4(jc,jk,jb)   =     zaetr  *   zvobga*pt_diag%dpres_mc(jc,jk,jb)
            zaeq5(jc,jk,jb)   =     zaetr  *   zstbga*pt_diag%dpres_mc(jc,jk,jb)

            zaetr_top(jc,jb) = zaetr_bot
            zaeqso(jc,jb)    = zaeqsn
            zaeqlo(jc,jb)    = zaeqln
            zaequo(jc,jb)    = zaequn
            zaeqdo(jc,jb)    = zaeqdn
          ENDDO
        ENDDO

      ELSE IF ((irad_aero == 6) .OR. (irad_aero == 9)) THEN ! Tegen aerosol climatology

        IF (iprog_aero == 0) THEN ! purely climatological aerosol
!DIR$ IVDEP
          DO jc = 1,i_endidx
            prm_diag%aerosol(jc,iss,jb) = ext_data%atm_td%aer_ss(jc,jb,imo1) + &
              & ( ext_data%atm_td%aer_ss(jc,jb,imo2)   - ext_data%atm_td%aer_ss(jc,jb,imo1)   ) * zw
            prm_diag%aerosol(jc,iorg,jb) = ext_data%atm_td%aer_org(jc,jb,imo1) + &
              & ( ext_data%atm_td%aer_org(jc,jb,imo2)  - ext_data%atm_td%aer_org(jc,jb,imo1)  ) * zw
            prm_diag%aerosol(jc,ibc,jb) = ext_data%atm_td%aer_bc(jc,jb,imo1) + &
              & ( ext_data%atm_td%aer_bc(jc,jb,imo2)   - ext_data%atm_td%aer_bc(jc,jb,imo1)   ) * zw
            prm_diag%aerosol(jc,iso4,jb) = ext_data%atm_td%aer_so4(jc,jb,imo1) + &
              & ( ext_data%atm_td%aer_so4(jc,jb,imo2)  - ext_data%atm_td%aer_so4(jc,jb,imo1)  ) * zw
            prm_diag%aerosol(jc,idu,jb) = ext_data%atm_td%aer_dust(jc,jb,imo1) + &
              & ( ext_data%atm_td%aer_dust(jc,jb,imo2) - ext_data%atm_td%aer_dust(jc,jb,imo1) ) * zw
          ENDDO
        ELSE ! simple prognostic aerosol; fill extra variables for climatology needed for relaxation equation
!DIR$ IVDEP
          DO jc = 1,i_endidx
            prm_diag%aercl_ss(jc,jb) = ext_data%atm_td%aer_ss(jc,jb,imo1) + &
              & ( ext_data%atm_td%aer_ss(jc,jb,imo2)   - ext_data%atm_td%aer_ss(jc,jb,imo1)   ) * zw
            prm_diag%aercl_or(jc,jb) = ext_data%atm_td%aer_org(jc,jb,imo1) + &
              & ( ext_data%atm_td%aer_org(jc,jb,imo2)  - ext_data%atm_td%aer_org(jc,jb,imo1)  ) * zw
            prm_diag%aercl_bc(jc,jb) = ext_data%atm_td%aer_bc(jc,jb,imo1) + &
              & ( ext_data%atm_td%aer_bc(jc,jb,imo2)   - ext_data%atm_td%aer_bc(jc,jb,imo1)   ) * zw
            prm_diag%aercl_su(jc,jb) = ext_data%atm_td%aer_so4(jc,jb,imo1) + &
              & ( ext_data%atm_td%aer_so4(jc,jb,imo2)  - ext_data%atm_td%aer_so4(jc,jb,imo1)  ) * zw
            prm_diag%aercl_du(jc,jb) = ext_data%atm_td%aer_dust(jc,jb,imo1) + &
              & ( ext_data%atm_td%aer_dust(jc,jb,imo2) - ext_data%atm_td%aer_dust(jc,jb,imo1) ) * zw
          ENDDO
        ENDIF

        DO jk = 2, nlevp1
          DO jc = 1,i_endidx
            zsign(jc,jk) = pt_diag%pres_ifc(jc,jk,jb) / 101325._wp
          ENDDO
        ENDDO

        ! The routine aerdis is called to recieve some parameters for the vertical
        ! distribution of background aerosol.
        CALL aerdis ( &
          & kbdim  = nproma,      & !in
          & jcs    = 1,           & !in
          & jce    = i_endidx,    & !in
          & klevp1 = nlevp1,      & !in
          & petah  = zsign(1,1),  & !in
          & pvdaes = zvdaes(1,1), & !out
          & pvdael = zvdael(1,1), & !out
          & pvdaeu = zvdaeu(1,1), & !out
          & pvdaed = zvdaed(1,1) )  !out

        ! top level
        DO jc = 1,i_endidx
          zaeqso   (jc,jb) = prm_diag%aerosol(jc,iss,jb)                                  *zvdaes(jc,1)
          zaeqlo   (jc,jb) = ( prm_diag%aerosol(jc,iorg,jb)+prm_diag%aerosol(jc,iso4,jb) )*zvdael(jc,1)
          zaequo   (jc,jb) =  prm_diag%aerosol(jc,ibc,jb)                                 *zvdaeu(jc,1)
          zaeqdo   (jc,jb) =  prm_diag%aerosol(jc,idu,jb)                                 *zvdaed(jc,1)
          zaetr_top(jc,jb) = 1.0_wp
        ENDDO

        ! loop over layers
        DO jk = 1,nlev
          DO jc = 1,i_endidx
            zaeqsn         =  prm_diag%aerosol(jc,iss,jb)                                 * zvdaes(jc,jk+1)
            zaeqln         =  (prm_diag%aerosol(jc,iorg,jb)+prm_diag%aerosol(jc,iso4,jb)) * zvdael(jc,jk+1)
            zaequn         =  prm_diag%aerosol(jc,ibc,jb)                                 * zvdaeu(jc,jk+1)
            zaeqdn         =  prm_diag%aerosol(jc,idu,jb)                                 * zvdaed(jc,jk+1)
            zaetr_bot      = zaetr_top(jc,jb) &
              & * ( MIN (1.0_wp, pt_diag%temp_ifc(jc,jk,jb)/pt_diag%temp_ifc(jc,jk+1,jb)) )**ztrpt

            zaetr          = SQRT(zaetr_bot*zaetr_top(jc,jb))
            zaeq1(jc,jk,jb)= (1.0_wp-zaetr)*( ztrbga*pt_diag%dpres_mc(jc,jk,jb) &
              &            + zaeqln - zaeqlo(jc,jb) )
            zaeq2(jc,jk,jb)   = (1._wp-zaetr) * ( zaeqsn-zaeqso(jc,jb) )
            zaeq3(jc,jk,jb)   = (1.0_wp-zaetr)*(zaeqdn-zaeqdo(jc,jb))
            zaeq4(jc,jk,jb)   = (1.0_wp-zaetr)*(zaequn-zaequo(jc,jb))
            zaeq5(jc,jk,jb)   =     zaetr  *   zstbga*pt_diag%dpres_mc(jc,jk,jb)

            zaetr_top(jc,jb) = zaetr_bot
            zaeqso(jc,jb)    = zaeqsn
            zaeqlo(jc,jb)    = zaeqln
            zaequo(jc,jb)    = zaequn
            zaeqdo(jc,jb)    = zaeqdn

          ENDDO
        ENDDO

      ELSE !no aerosols

        zaeq1(1:i_endidx,:,jb) = 0.0_wp
        zaeq2(1:i_endidx,:,jb) = 0.0_wp
        zaeq3(1:i_endidx,:,jb) = 0.0_wp
        zaeq4(1:i_endidx,:,jb) = 0.0_wp
        zaeq5(1:i_endidx,:,jb) = 0.0_wp

      ENDIF ! irad_aero

      ! Compute cloud number concentration depending on aerosol climatology if
      ! aerosol-microphysics or aerosol-convection coupling is turned on
      IF (atm_phy_nwp_config(jg)%icpl_aero_gscp == 1 .OR. icpl_aero_conv == 1) THEN

        DO jk = 1,nlev
!DIR$ IVDEP
          DO jc = 1, i_endidx
            wfac = MAX(1._wp,(MIN(8._wp,0.8_wp*pt_diag%pres_sfc(jc,jb)/pt_diag%pres(jc,jk,jb))))**2
            ncn_bg = MIN(prm_diag%cloud_num(jc,jb),50.e6_wp)
            prm_diag%acdnc(jc,jk,jb) = (ncn_bg+(prm_diag%cloud_num(jc,jb)-ncn_bg)*(EXP(1._wp-wfac)))
          END DO
        END DO

      ENDIF

      IF ( irad_o3 == 6 ) THEN ! Old GME ozone climatology

        ! 3-dimensional O3
        ! top level
        ! Loop starts with 1 instead of i_startidx because the start index is missing in RRTM
        DO jc = 1,i_endidx
          zptop32  (jc,jb) = (SQRT(pt_diag%pres_ifc(jc,1,jb)))**3
          zo3_hm   (jc,jb) = (SQRT(prm_diag%hmo3(jc,jb)))**3
          zo3_top  (jc,jb) = prm_diag%vio3(jc,jb)*zptop32(jc,jb)/(zptop32(jc,jb)+zo3_hm(jc,jb))
        ENDDO
        ! loop over layers
        DO jk = 1,nlev
!DIR$ IVDEP
          DO jc = 1,i_endidx
            zpbot32  (jc,jb) = (SQRT(pt_diag%pres_ifc(jc,jk+1,jb)))**3
            zo3_bot  (jc,jb) = prm_diag%vio3(jc,jb)* zpbot32(jc,jb)    &
              /( zpbot32(jc,jb) + zo3_hm(jc,jb))
            !O3 content
            ext_data%atm%o3(jc,jk,jb) = (zo3_bot(jc,jb) - zo3_top(jc,jb))/pt_diag%dpres_mc(jc,jk,jb)
            ! store previous bottom values in arrays for top of next layer
            zo3_top (jc,jb) = zo3_bot (jc,jb)
          ENDDO
        ENDDO

      ENDIF

      ! Needed for RG radiation only
      IF (PRESENT(zduo3)) THEN
        DO jk = 1,nlev
          DO jc = 1,i_endidx
            zduo3(jc,jk,jb) = ext_data%atm%o3(jc,jk,jb)*pt_diag%dpres_mc(jc,jk,jb)
          ENDDO
        ENDDO
      ENDIF
    ENDDO !jb
!$OMP END DO NOWAIT
!$OMP END PARALLEL

    IF (timers_level > 3) CALL timer_stop(timer_preradiaton)

  END SUBROUTINE nwp_ozon_aerosol
  !---------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------
  !>
  !! @par Revision History
  !! Initial release by Thorsten Reinhardt, AGeoBw, Offenbach (2011-01-13)
  !!
  SUBROUTINE nwp_rrtm_radiation ( current_date, pt_patch, ext_data,                      &
    &  zaeq1, zaeq2, zaeq3, zaeq4, zaeq5, pt_diag, prm_diag, lnd_prog, irad )

    CHARACTER(len=MAX_CHAR_LENGTH), PARAMETER::  &
      &  routine = modname//'::nwp_rrtm_radiation'

    TYPE(datetime), POINTER, INTENT(in) :: current_date

    TYPE(t_patch),        TARGET,INTENT(in) :: pt_patch     !<grid/patch info.
    TYPE(t_external_data),INTENT(in):: ext_data

    REAL(wp), INTENT(in) :: &
      & zaeq1(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq2(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq3(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq4(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq5(nproma,pt_patch%nlev,pt_patch%nblks_c)

    TYPE(t_nh_diag), TARGET, INTENT(in)  :: pt_diag     !<the diagnostic variables
    TYPE(t_nwp_phy_diag),       INTENT(inout):: prm_diag
    TYPE(t_lnd_prog),           INTENT(inout):: lnd_prog

    INTEGER, INTENT(IN) :: irad ! To distinguish between RRTM (1) and PSRAD (3)

    REAL(wp):: aclcov(nproma,pt_patch%nblks_c), dust_tunefac(nproma,nbndlw)

    ! Local scalars:
    INTEGER:: jb
    INTEGER:: jg                !domain id
    INTEGER:: nlev, nlevp1      !< number of full and half levels

    INTEGER:: rl_start, rl_end
    INTEGER:: i_startblk, i_endblk    !> blocks
    INTEGER:: i_startidx, i_endidx    !< slices
    INTEGER:: i_nchdom                !< domain index
    REAL(wp) :: t_sf(nproma,pt_patch%nblks_c)

    i_nchdom  = MAX(1,pt_patch%n_childdom)
    jg        = pt_patch%id

    ! number of vertical levels
    nlev   = pt_patch%nlev
    nlevp1 = pt_patch%nlevp1

    !-------------------------------------------------------------------------
    !> Radiation
    !-------------------------------------------------------------------------

    rl_start = grf_bdywidth_c+1
    rl_end   = min_rlcell_int

    i_startblk = pt_patch%cells%start_blk(rl_start,1)
    i_endblk   = pt_patch%cells%end_blk(rl_end,i_nchdom)

    IF (msg_level >= 12) &
      &           CALL message(routine, 'RRTM radiation on full grid')

    IF (p_test_run) THEN
      CALL get_indices_c(pt_patch, i_startblk, i_startblk, i_endblk, &
           &                         i_startidx, i_endidx, rl_start, rl_end)
      aclcov(1:i_startidx-1,i_startblk) = 0
      prm_diag%lwflxall(1:i_startidx-1,:,i_startblk) = 0
      prm_diag%trsolall(1:i_startidx-1,:,i_startblk) = 0
      prm_diag%lwflx_up_sfc_rs(1:i_startidx-1,i_startblk) = 0
      prm_diag%trsol_up_toa(1:i_startidx-1,i_startblk) = 0
      prm_diag%trsol_up_sfc(1:i_startidx-1,i_startblk) = 0
      prm_diag%trsol_par_sfc(1:i_startidx-1,i_startblk) = 0
      prm_diag%trsol_dn_sfc_diff(1:i_startidx-1,i_startblk) = 0
      prm_diag%trsolclr_sfc(1:i_startidx-1,i_startblk) = 0
    END IF

#ifdef COUP_OAS_ICON
      t_sf(:,:) = oas_rcv_field_icon(:,:,8)
#else
      t_sf = lnd_prog%t_g(:,:)
#endif

!$OMP PARALLEL PRIVATE(jb,i_startidx,i_endidx,dust_tunefac)
!$OMP DO ICON_OMP_GUIDED_SCHEDULE
    DO jb = i_startblk, i_endblk

      CALL get_indices_c(pt_patch, jb, i_startblk, i_endblk, &
        &                         i_startidx, i_endidx, rl_start, rl_end)

      ! It may happen that an MPI patch contains only nest boundary points
      ! In this case, no action is needed
      IF (i_startidx > i_endidx) CYCLE

      prm_diag%tsfctrad(1:i_endidx,jb) = t_sf(1:i_endidx,jb)

      IF (tune_dust_abs > 0._wp) THEN
!DIR$ NOINLINE
        CALL tune_dust(pt_patch%cells%center(:,jb)%lat,pt_patch%cells%center(:,jb)%lon,i_endidx,dust_tunefac)
      ELSE
        dust_tunefac(:,:) = 1._wp
      ENDIF

      CALL radiation_nwp(               &
                              !
                              ! input
                              ! -----
                              !
        & current_date                     ,&!< in current date
                              ! indices and dimensions
        & jg         =jg                   ,&!< in domain index
        & jb         =jb                   ,&!< in block index
        & irad       =irad                 ,&!< in option for radiation scheme (RRTM/PSRAD)
        & jce        =i_endidx             ,&!< in  end   index for loop over block
        & kbdim      =nproma               ,&!< in  dimension of block over cells
        & klev       =nlev                 ,&!< in  number of full levels = number of layers
        & klevp1     =nlevp1               ,&!< in  number of half levels = number of layer ifcs
                              !
        & ktype      =prm_diag%ktype(:,jb) ,&!< in     type of convection
                              !
                              ! surface: albedo + temperature
        & zland      =ext_data%atm%fr_land_smt(:,jb)   ,&!< in     land fraction
        & zglac      =ext_data%atm%fr_glac_smt(:,jb)   ,&!< in     land glacier fraction
                              !
        & cos_mu0    =prm_diag%cosmu0  (:,jb) ,&!< in  cos of zenith angle mu0
        & alb_vis_dir=prm_diag%albvisdir(:,jb) ,&!< in surface albedo for visible range, direct
        & alb_nir_dir=prm_diag%albnirdir(:,jb) ,&!< in surface albedo for near IR range, direct
        & alb_vis_dif=prm_diag%albvisdif(:,jb),&!< in surface albedo for visible range, diffuse
        & alb_nir_dif=prm_diag%albnirdif(:,jb),&!< in surface albedo for near IR range, diffuse
        & emis_rad   =ext_data%atm%emis_rad(:,jb),&!< in longwave surface emissivity
        & tk_sfc     =prm_diag%tsfctrad(:,jb) ,&!< in surface temperature
                              !
                              ! atmosphere: pressure, tracer mixing ratios and temperature
        & pp_hl      =pt_diag%pres_ifc  (:,:,jb)     ,&!< in  pres at half levels at t-dt [Pa]
        & pp_fl      =pt_diag%pres      (:,:,jb)     ,&!< in  pres at full levels at t-dt [Pa]
        & tk_fl      =pt_diag%temp      (:,:,jb)     ,&!< in  temperature at full level at t-dt
        & qm_vap     =prm_diag%tot_cld  (:,:,jb,iqv) ,&!< in  water vapor mass mix ratio at t-dt
        & qm_liq     =prm_diag%tot_cld  (:,:,jb,iqc) ,&!< in cloud water mass mix ratio at t-dt
        & qm_ice     =prm_diag%tot_cld  (:,:,jb,iqi) ,&!< in cloud ice mass mixing ratio at t-dt
        & qm_o3      =ext_data%atm%o3   (:,:,jb)     ,&!< in o3 mass mixing ratio at t-dt
        & cdnc       =prm_diag%acdnc    (:,:,jb)     ,&!< in  cloud droplet numb conc. [1/m**3]
        & cld_frc    =prm_diag%clc      (:,:,jb)     ,&!< in  cloud fraction [m2/m2]
        & zaeq1      = zaeq1(:,:,jb)                 ,&!< in aerosol continental
        & zaeq2      = zaeq2(:,:,jb)                 ,&!< in aerosol maritime
        & zaeq3      = zaeq3(:,:,jb)                 ,&!< in aerosol urban
        & zaeq4      = zaeq4(:,:,jb)                 ,&!< in aerosol volcano ashes
        & zaeq5      = zaeq5(:,:,jb)                 ,&!< in aerosol stratospheric background
        & dust_tunefac = dust_tunefac (:,:)          ,&!< in LW tuning factor for dust aerosol
        & dt_rad     = atm_phy_nwp_config(jg)%dt_rad ,&
                              ! output
                              ! ------
                              !
        & cld_cvr    =  aclcov             (:,jb),  &     !< out cloud cover in a column [m2/m2]
        & flx_lw_net =  prm_diag%lwflxall(:,:,jb), &      !< out terrestrial flux, all sky, net down
        & trsol_net  =  prm_diag%trsolall(:,:,jb), &      !< out solar transmissivity, all sky, net down
        & flx_uplw_sfc = prm_diag%lwflx_up_sfc_rs(:,jb),& !< out longwave upward flux at surface
        & trsol_up_toa = prm_diag%trsol_up_toa(:,jb), &   !< out upward solar transmissivity at TOA
        & trsol_up_sfc = prm_diag%trsol_up_sfc(:,jb), &   !< out upward solar transmissivity at surface
        & trsol_par_sfc = prm_diag%trsol_par_sfc(:,jb), & !< out downward transmissivity for photosynthetically active rad. at surface
        & trsol_dn_sfc_diffus = prm_diag%trsol_dn_sfc_diff(:,jb), &  !< out downward diffuse solar transmissivity at surface
        & trsol_clr_sfc = prm_diag%trsolclr_sfc(:,jb)  )  !< out clear-sky net transmissvity at surface

      ENDDO ! blocks

!$OMP END DO NOWAIT
!$OMP END PARALLEL

  END SUBROUTINE nwp_rrtm_radiation
  !---------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------
  !>
  !!
  !! @par Revision History
  !! Initial release by Thorsten Reinhardt, AGeoBw, Offenbach (2011-01-13)
  !!
  SUBROUTINE nwp_rrtm_radiation_reduced ( current_date, pt_patch, pt_par_patch, ext_data, &
    &                                     zaeq1,zaeq2,zaeq3,zaeq4,zaeq5,    &
    &                                     pt_diag,prm_diag,lnd_prog,irad    )

    CHARACTER(len=MAX_CHAR_LENGTH), PARAMETER::  &
      &  routine = modname//'::nwp_rrtm_radiation_reduced'

    TYPE(datetime), POINTER, INTENT(in) :: current_date

    TYPE(t_patch),        TARGET,INTENT(in) :: pt_patch     !<grid/patch info.
    TYPE(t_patch),        TARGET,INTENT(in) :: pt_par_patch !<grid/patch info (parent grid)
    TYPE(t_external_data),INTENT(in):: ext_data
    REAL(wp),             INTENT(in) ::               &
      & zaeq1(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq2(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq3(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq4(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq5(nproma,pt_patch%nlev,pt_patch%nblks_c)

    TYPE(t_nh_diag), TARGET,    INTENT(inout):: pt_diag     !<the diagnostic variables
    TYPE(t_nwp_phy_diag),       INTENT(inout):: prm_diag
    TYPE(t_lnd_prog),           INTENT(inout):: lnd_prog

    INTEGER, INTENT(IN) :: irad ! To distinguish between RRTM (1) and PSRAD (3)

    REAL(wp):: aclcov(nproma,pt_patch%nblks_c), dust_tunefac(nproma,nbndlw)
    ! For radiation on reduced grid
    ! These fields need to be allocatable because they have different dimensions for
    ! the global grid and nested grids, and for runs with/without MPI parallelization
    ! Input fields
    REAL(wp), ALLOCATABLE, TARGET:: zrg_fr_land  (:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_fr_glac  (:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_emis_rad (:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_cosmu0   (:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_albvisdir(:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_albnirdir(:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_albvisdif(:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_albnirdif(:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_albdif   (:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_tsfc     (:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_rtype    (:,:) ! type of convection (integer)
    INTEGER,  ALLOCATABLE, TARGET:: zrg_ktype    (:,:) ! type of convection (real)

    REAL(wp), ALLOCATABLE, TARGET:: zrg_pres_ifc (:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zlp_pres_ifc (:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_pres     (:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_temp     (:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_o3       (:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_acdnc    (:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_tot_cld  (:,:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zlp_tot_cld  (:,:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_clc      (:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_aeq1(:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_aeq2(:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_aeq3(:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_aeq4(:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_aeq5(:,:,:)
    ! Output fields
    REAL(wp), ALLOCATABLE, TARGET:: zrg_aclcov   (:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_lwflxall (:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_trsolall (:,:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_lwflx_up_sfc   (:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_trsol_up_toa   (:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_trsol_up_sfc   (:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_trsol_par_sfc  (:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_trsol_dn_sfc_diff(:,:)
    REAL(wp), ALLOCATABLE, TARGET:: zrg_trsol_clr_sfc   (:,:)

    ! Pointer to parent patach or local parent patch for reduced grid
    TYPE(t_patch), POINTER       :: ptr_pp

    ! Variables for debug output
    REAL(wp) :: max_albvisdir, min_albvisdir, max_albvisdif, min_albvisdif, &
                max_albdif, min_albdif, max_tsfc, min_tsfc, max_psfc, min_psfc

    REAL(wp), DIMENSION(:), ALLOCATABLE :: max_pres_ifc, max_pres, max_temp, max_acdnc, &
        max_qv, max_qc, max_qi, max_cc, min_pres_ifc, min_pres, min_temp, min_acdnc, &
        min_qv, min_qc, min_qi, min_cc

    REAL(wp), DIMENSION(pt_patch%nlevp1) :: max_lwflx, min_lwflx, max_swtrans, min_swtrans

    ! Local scalars:
    INTEGER:: jk,jb
    INTEGER:: jg                      !domain id
    INTEGER:: nlev, nlevp1, nlev_rg   !< number of full and half levels
    INTEGER:: nblks_par_c, nblks_lp_c !nblks for reduced grid

    INTEGER:: rl_start, rl_end
    INTEGER:: i_startblk, i_endblk    !> blocks
    INTEGER:: i_startidx, i_endidx    !< slices
    INTEGER:: i_nchdom                !< domain index
    INTEGER:: i_chidx

    i_nchdom  = MAX(1,pt_patch%n_childdom)
    jg        = pt_patch%id

    ! number of vertical levels
    nlev   = pt_patch%nlev
    nlevp1 = pt_patch%nlevp1

    !-------------------------------------------------------------------------
    !> Radiation
    !-------------------------------------------------------------------------

      ! section for computing radiation on reduced grid

      IF (msg_level >= 12) &
        &       CALL message(routine, 'RRTM radiation on reduced grid')

      i_chidx     =  pt_patch%parent_child_index

      IF (jg == 1 .AND. .NOT. l_limited_area) THEN
        ptr_pp => pt_par_patch
        nblks_par_c = pt_par_patch%nblks_c
        nblks_lp_c  =  p_patch_local_parent(jg)%nblks_c
      ELSE ! Nested domain with MPI parallelization
        ptr_pp      => p_patch_local_parent(jg)
        nblks_par_c =  ptr_pp%nblks_c
        nblks_lp_c  =  ptr_pp%nblks_c
      ENDIF

      ! Add extra layer for atmosphere above model top if requested
      IF (atm_phy_nwp_config(jg)%latm_above_top) THEN
        IF (jg == 1 .OR. pt_patch%nshift == 0) THEN
          nlev_rg = nlev + 1
        ELSE ! add a specified number levels up to the top of the parent domain in case of vertical nesting
          nlev_rg = MIN(nlev+nexlevs_rrg_vnest, pt_par_patch%nlev)
        ENDIF
      ELSE
        nlev_rg = nlev
      ENDIF

      ALLOCATE (zrg_cosmu0   (nproma,nblks_par_c),     &
        zrg_fr_land  (nproma,nblks_par_c),             &
        zrg_fr_glac  (nproma,nblks_par_c),             &
        zrg_emis_rad (nproma,nblks_par_c),             &
        zrg_albvisdir(nproma,nblks_par_c),             &
        zrg_albnirdir(nproma,nblks_par_c),             &
        zrg_albvisdif(nproma,nblks_par_c),             &
        zrg_albnirdif(nproma,nblks_par_c),             &
        zrg_albdif   (nproma,nblks_par_c),             &
        zrg_tsfc     (nproma,nblks_par_c),             &
        zrg_rtype    (nproma,nblks_par_c),             &
        zrg_ktype    (nproma,nblks_par_c),             &
        zrg_pres_ifc (nproma,nlev_rg+1,nblks_par_c),   &
        zlp_pres_ifc (nproma,nlev_rg+1,nblks_lp_c ),   &
        zrg_pres     (nproma,nlev_rg  ,nblks_par_c),   &
        zrg_temp     (nproma,nlev_rg  ,nblks_par_c),   &
        zrg_o3       (nproma,nlev_rg  ,nblks_par_c),   &
        zrg_aeq1     (nproma,nlev_rg  ,nblks_par_c),   &
        zrg_aeq2     (nproma,nlev_rg  ,nblks_par_c),   &
        zrg_aeq3     (nproma,nlev_rg  ,nblks_par_c),   &
        zrg_aeq4     (nproma,nlev_rg  ,nblks_par_c),   &
        zrg_aeq5     (nproma,nlev_rg  ,nblks_par_c),   &
        zrg_acdnc    (nproma,nlev_rg  ,nblks_par_c),   &
        zrg_tot_cld  (nproma,nlev_rg  ,nblks_par_c,3), &
        zlp_tot_cld  (nproma,nlev_rg  ,nblks_lp_c,3),  &
        zrg_clc      (nproma,nlev_rg  ,nblks_par_c),   &
        zrg_aclcov   (nproma,          nblks_par_c),   &
        zrg_lwflx_up_sfc   (nproma,    nblks_par_c),   &
        zrg_trsol_up_toa   (nproma,    nblks_par_c),   &
        zrg_trsol_up_sfc   (nproma,    nblks_par_c),   &
        zrg_trsol_par_sfc  (nproma,    nblks_par_c),   &
        zrg_trsol_dn_sfc_diff(nproma,  nblks_par_c),   &
        zrg_trsol_clr_sfc    (nproma,  nblks_par_c),   &
        zrg_lwflxall (nproma,nlev_rg+1,nblks_par_c),   &
        zrg_trsolall (nproma,nlev_rg+1,nblks_par_c)    )

      rl_start = 1 ! SR radiation is not set up to handle boundaries of nested domains
      rl_end   = min_rlcell_int

      i_startblk = pt_patch%cells%start_blk(rl_start,1)
      i_endblk   = pt_patch%cells%end_blk(rl_end,i_nchdom)

      IF (p_test_run) THEN
        CALL get_indices_c(pt_patch, i_startblk, i_startblk, i_endblk, &
             &                         i_startidx, i_endidx, rl_start, rl_end)
        zrg_lwflxall(:,:,:) = 0._wp
        zrg_trsolall(:,:,:) = 0._wp
        zrg_lwflx_up_sfc(1:i_startidx-1,i_startblk) = 0
        zrg_trsol_up_toa(1:i_startidx-1,i_startblk) = 0
        zrg_trsol_up_sfc(1:i_startidx-1,i_startblk) = 0
        zrg_trsol_par_sfc(1:i_startidx-1,i_startblk) = 0
        zrg_trsol_dn_sfc_diff(1:i_startidx-1,i_startblk) = 0
        zrg_trsol_clr_sfc(1:i_startidx-1,i_startblk) = 0
      ENDIF

      DO jb = i_startblk, i_endblk

        CALL get_indices_c(pt_patch, jb, i_startblk, i_endblk, &
          &                       i_startidx, i_endidx, rl_start, rl_end)

        ! Loop starts with 1 instead of i_startidx because the start index is missing in RRTM
        prm_diag%tsfctrad(1:i_endidx,jb) = lnd_prog%t_g(1:i_endidx,jb)

      ENDDO ! blocks

      ! For limited-area radiation grids, tsfc needs to be filled with air temp along the nest boundary
      ! because the surface scheme is not active on the boundary points
      IF (jg == 1 .AND. l_limited_area) THEN
        rl_start = 1
        rl_end   = grf_bdywidth_c
        i_startblk = pt_patch%cells%start_blk(rl_start,1)
        i_endblk   = pt_patch%cells%end_blk(rl_end,i_nchdom)

        DO jb = i_startblk, i_endblk

          CALL get_indices_c(pt_patch, jb, i_startblk, i_endblk, &
            &                       i_startidx, i_endidx, rl_start, rl_end)

          prm_diag%tsfctrad(i_startidx:i_endidx,jb) = pt_diag%temp(i_startidx:i_endidx,nlev,jb)

        ENDDO ! blocks
      ENDIF

      CALL upscale_rad_input(pt_patch%id, pt_par_patch%id,              &
        & nlev_rg, ext_data%atm%fr_land_smt, ext_data%atm%fr_glac_smt,  &
        & ext_data%atm%emis_rad, prm_diag%cosmu0,                       &
        & prm_diag%albvisdir, prm_diag%albnirdir, prm_diag%albvisdif,   &
        & prm_diag%albnirdif, prm_diag%albdif, prm_diag%tsfctrad,       &
        & prm_diag%ktype, pt_diag%pres_ifc, pt_diag%pres,               &
        & pt_diag%temp,prm_diag%acdnc, prm_diag%tot_cld, prm_diag%clc,  &
        & ext_data%atm%o3, zaeq1, zaeq2, zaeq3, zaeq4, zaeq5,           &
        & zrg_fr_land, zrg_fr_glac, zrg_emis_rad,                       &
        & zrg_cosmu0, zrg_albvisdir, zrg_albnirdir, zrg_albvisdif,      &
        & zrg_albnirdif, zrg_albdif, zrg_tsfc, zrg_rtype, zrg_pres_ifc, &
        & zrg_pres, zrg_temp, zrg_acdnc, zrg_tot_cld, zrg_clc, zrg_o3,  &
        & zrg_aeq1, zrg_aeq2, zrg_aeq3, zrg_aeq4, zrg_aeq5,             &
        & zlp_pres_ifc, zlp_tot_cld, prm_diag%buffer_rrg)

      rl_start = grf_ovlparea_start_c
      rl_end   = min_rlcell_int

      i_startblk = ptr_pp%cells%start_blk(rl_start,i_chidx)
      i_endblk   = ptr_pp%cells%end_blk(rl_end,i_chidx)

      ! Debug output of radiation input fields
      IF (msg_level >= 16) THEN

        ALLOCATE(max_pres_ifc(nlev_rg), max_pres(nlev_rg), max_temp(nlev_rg), max_acdnc(nlev_rg), &
                 max_qv(nlev_rg), max_qc(nlev_rg), max_qi(nlev_rg), max_cc(nlev_rg),              &
                 min_pres_ifc(nlev_rg), min_pres(nlev_rg), min_temp(nlev_rg), min_acdnc(nlev_rg), &
                 min_qv(nlev_rg), min_qc(nlev_rg), min_qi(nlev_rg), min_cc(nlev_rg)               )

        max_albvisdir = 0._wp
        min_albvisdir = 1.e10_wp
        max_albvisdif = 0._wp
        min_albvisdif = 1.e10_wp
        max_albdif = 0._wp
        min_albdif = 1.e10_wp
        max_tsfc = 0._wp
        min_tsfc = 1.e10_wp
        max_psfc = 0._wp
        min_psfc = 1.e10_wp
        max_pres_ifc = 0._wp
        max_pres    = 0._wp
        max_temp     = 0._wp
        max_acdnc    = 0._wp
        max_qv = 0._wp
        max_qc = 0._wp
        max_qi = 0._wp
        max_cc  = 0._wp
        min_pres_ifc = 1.e10_wp
        min_pres    = 1.e10_wp
        min_temp     = 1.e10_wp
        min_acdnc    = 1.e10_wp
        min_qv = 1.e10_wp
        min_qc = 1.e10_wp
        min_qi = 1.e10_wp
        min_cc  = 1.e10_wp

        DO jb = i_startblk, i_endblk

         CALL get_indices_c(ptr_pp, jb, i_startblk, i_endblk, &
                            i_startidx, i_endidx, rl_start, rl_end)

         max_albvisdir = MAX(max_albvisdir,MAXVAL(zrg_albvisdir(i_startidx:i_endidx,jb)))
         min_albvisdir = MIN(min_albvisdir,MINVAL(zrg_albvisdir(i_startidx:i_endidx,jb)))
         max_albvisdif = MAX(max_albvisdif,MAXVAL(zrg_albvisdif(i_startidx:i_endidx,jb)))
         min_albvisdif = MIN(min_albvisdif,MINVAL(zrg_albvisdif(i_startidx:i_endidx,jb)))
         max_albdif    = MAX(max_albdif,   MAXVAL(zrg_albdif(i_startidx:i_endidx,jb)))
         min_albdif    = MIN(min_albdif,   MINVAL(zrg_albdif(i_startidx:i_endidx,jb)))
         max_tsfc = MAX(max_tsfc,MAXVAL(zrg_tsfc(i_startidx:i_endidx,jb)))
         min_tsfc = MIN(min_tsfc,MINVAL(zrg_tsfc(i_startidx:i_endidx,jb)))
         max_psfc = MAX(max_psfc,MAXVAL(zrg_pres_ifc(i_startidx:i_endidx,nlev_rg+1,jb)))
         min_psfc = MIN(min_psfc,MINVAL(zrg_pres_ifc(i_startidx:i_endidx,nlev_rg+1,jb)))
         DO jk = 1, nlev_rg
          max_pres_ifc(jk) = MAX(max_pres_ifc(jk),MAXVAL(zrg_pres_ifc(i_startidx:i_endidx,jk,jb)))
          max_pres(jk)    = MAX(max_pres(jk),MAXVAL(zrg_pres     (i_startidx:i_endidx,jk,jb)))
          max_temp(jk)     = MAX(max_temp(jk),MAXVAL(zrg_temp     (i_startidx:i_endidx,jk,jb)))
          max_acdnc(jk)    = MAX(max_acdnc(jk),MAXVAL(zrg_acdnc    (i_startidx:i_endidx,jk,jb)))
          max_qv(jk) = MAX(max_qv(jk),MAXVAL(zrg_tot_cld(i_startidx:i_endidx,jk,jb,iqv)))
          max_qc(jk) = MAX(max_qc(jk),MAXVAL(zrg_tot_cld(i_startidx:i_endidx,jk,jb,iqc)))
          max_qi(jk) = MAX(max_qi(jk),MAXVAL(zrg_tot_cld(i_startidx:i_endidx,jk,jb,iqi)))
          max_cc(jk)  = MAX(max_cc(jk),MAXVAL(zrg_clc(i_startidx:i_endidx,jk,jb)))
          min_pres_ifc(jk) = MIN(min_pres_ifc(jk),MINVAL(zrg_pres_ifc(i_startidx:i_endidx,jk,jb)))
          min_pres(jk)    = MIN(min_pres(jk),MINVAL(zrg_pres     (i_startidx:i_endidx,jk,jb)))
          min_temp(jk)     = MIN(min_temp(jk),MINVAL(zrg_temp     (i_startidx:i_endidx,jk,jb)))
          min_acdnc(jk)    = MIN(min_acdnc(jk),MINVAL(zrg_acdnc    (i_startidx:i_endidx,jk,jb)))
          min_qv(jk) = MIN(min_qv(jk),MINVAL(zrg_tot_cld(i_startidx:i_endidx,jk,jb,iqv)))
          min_qc(jk) = MIN(min_qc(jk),MINVAL(zrg_tot_cld(i_startidx:i_endidx,jk,jb,iqc)))
          min_qi(jk) = MIN(min_qi(jk),MINVAL(zrg_tot_cld(i_startidx:i_endidx,jk,jb,iqi)))
          min_cc(jk)  = MIN(min_cc(jk),MINVAL(zrg_clc(i_startidx:i_endidx,jk,jb)))
         ENDDO
        ENDDO ! blocks

        max_albvisdir = global_max(max_albvisdir)
        min_albvisdir = global_min(min_albvisdir)
        max_albvisdif = global_max(max_albvisdif)
        min_albvisdif = global_min(min_albvisdif)
        max_tsfc = global_max(max_tsfc)
        min_tsfc = global_min(min_tsfc)
        max_psfc = global_max(max_psfc)
        min_psfc = global_min(min_psfc)
        max_pres_ifc = global_max(max_pres_ifc)
        max_pres    = global_max(max_pres)
        max_temp     = global_max(max_temp)
        max_acdnc    = global_max(max_acdnc)
        max_qv = global_max(max_qv)
        max_qc = global_max(max_qc)
        max_qi = global_max(max_qi)
        max_cc  = global_max(max_cc)
        min_pres_ifc = global_min(min_pres_ifc)
        min_pres    = global_min(min_pres)
        min_temp     = global_min(min_temp)
        min_acdnc    = global_min(min_acdnc)
        min_qv = global_min(min_qv)
        min_qc = global_min(min_qc)
        min_qi = global_min(min_qi)
        min_cc  = global_min(min_cc)

        WRITE(message_text,'(a,4f12.8)') 'max/min alb = ', max_albvisdir, min_albvisdir, &
          max_albvisdif, min_albvisdif
        CALL message(routine, TRIM(message_text))

        WRITE(message_text,'(a,2f10.3,2f10.2)') 'max/min sfc temp/pres = ', max_tsfc, min_tsfc, &
          max_psfc, min_psfc
        CALL message(routine, TRIM(message_text))

        WRITE(message_text,'(a)') 'max/min pres_ifc, pres, temp, acdnc'
        CALL message(routine, TRIM(message_text))

        DO jk = 1, nlev_rg
          WRITE(message_text,'(i4,4f10.2,2f10.3,2f12.1)') jk,max_pres_ifc(jk), min_pres_ifc(jk), &
            max_pres(jk), min_pres(jk), max_temp(jk), min_temp(jk), max_acdnc(jk), min_acdnc(jk)
          CALL message(routine, TRIM(message_text))
        ENDDO

        WRITE(message_text,'(a)') 'max/min QV, QC, QI, CC'
        CALL message(routine, TRIM(message_text))

        DO jk = 1, nlev_rg
          WRITE(message_text,'(i4,8e13.5)') jk,max_qv(jk), min_qv(jk), max_qc(jk), min_qc(jk), &
             max_qi(jk), min_qi(jk), max_cc(jk), min_cc(jk)
          CALL message(routine, TRIM(message_text))
        ENDDO

        DEALLOCATE(max_pres_ifc, max_pres, max_temp, max_acdnc, max_qv, max_qc, max_qi, max_cc, &
                   min_pres_ifc, min_pres, min_temp, min_acdnc, min_qv, min_qc, min_qi, min_cc)

      ENDIF ! msg_level >= 16

!$OMP PARALLEL
!$OMP DO PRIVATE(jb,jk,i_startidx,i_endidx,dust_tunefac) ICON_OMP_GUIDED_SCHEDULE
      DO jb = i_startblk, i_endblk

        CALL get_indices_c(ptr_pp, jb, i_startblk, i_endblk, &
          &                         i_startidx, i_endidx, rl_start, rl_end)

        ! It may happen that an MPI patch contains only nest boundary points
        ! In this case, no action is needed
        IF (i_startidx > i_endidx) CYCLE

        ! Unfortunately, the coding of SR radiation is not compatible with the presence
        ! of nested domains. Therefore, the normally unused elements of the first block
        ! need to be filled with dummy values
        IF ( (jg > 1 .OR. l_limited_area) .AND. jb == i_startblk) THEN
          zrg_fr_land   (1:i_startidx-1,jb) = zrg_fr_land   (i_startidx,jb)
          zrg_fr_glac   (1:i_startidx-1,jb) = zrg_fr_glac   (i_startidx,jb)
          zrg_emis_rad  (1:i_startidx-1,jb) = zrg_emis_rad  (i_startidx,jb)
          zrg_cosmu0    (1:i_startidx-1,jb) = zrg_cosmu0    (i_startidx,jb)
          zrg_albvisdir (1:i_startidx-1,jb) = zrg_albvisdir (i_startidx,jb)
          zrg_albnirdir (1:i_startidx-1,jb) = zrg_albnirdir (i_startidx,jb)
          zrg_albvisdif (1:i_startidx-1,jb) = zrg_albvisdif (i_startidx,jb)
          zrg_albnirdif (1:i_startidx-1,jb) = zrg_albnirdif (i_startidx,jb)
          zrg_albdif    (1:i_startidx-1,jb) = zrg_albdif    (i_startidx,jb)
          zrg_tsfc      (1:i_startidx-1,jb) = zrg_tsfc      (i_startidx,jb)
          zrg_rtype     (1:i_startidx-1,jb) = zrg_rtype     (i_startidx,jb)
          zrg_pres_ifc (1:i_startidx-1,nlev_rg+1,jb) = zrg_pres_ifc (i_startidx,nlev_rg+1,jb)
          DO jk = 1, nlev_rg
            zrg_pres_ifc (1:i_startidx-1,jk,jb) = zrg_pres_ifc (i_startidx,jk,jb)
            zrg_pres     (1:i_startidx-1,jk,jb) = zrg_pres     (i_startidx,jk,jb)
            zrg_temp     (1:i_startidx-1,jk,jb) = zrg_temp     (i_startidx,jk,jb)
            zrg_o3       (1:i_startidx-1,jk,jb) = zrg_o3       (i_startidx,jk,jb)
            zrg_aeq1     (1:i_startidx-1,jk,jb) = zrg_aeq1     (i_startidx,jk,jb)
            zrg_aeq2     (1:i_startidx-1,jk,jb) = zrg_aeq2     (i_startidx,jk,jb)
            zrg_aeq3     (1:i_startidx-1,jk,jb) = zrg_aeq3     (i_startidx,jk,jb)
            zrg_aeq4     (1:i_startidx-1,jk,jb) = zrg_aeq4     (i_startidx,jk,jb)
            zrg_aeq5     (1:i_startidx-1,jk,jb) = zrg_aeq5     (i_startidx,jk,jb)
            zrg_acdnc    (1:i_startidx-1,jk,jb) = zrg_acdnc    (i_startidx,jk,jb)
            zrg_tot_cld  (1:i_startidx-1,jk,jb,iqv) = zrg_tot_cld(i_startidx,jk,jb,iqv)
            zrg_tot_cld  (1:i_startidx-1,jk,jb,iqc) = zrg_tot_cld(i_startidx,jk,jb,iqc)
            zrg_tot_cld  (1:i_startidx-1,jk,jb,iqi) = zrg_tot_cld(i_startidx,jk,jb,iqi)
            zrg_clc      (1:i_startidx-1,jk,jb) = zrg_clc(i_startidx,jk,jb)
          ENDDO
        ENDIF

        IF (tune_dust_abs > 0._wp) THEN
!DIR$ NOINLINE
          CALL tune_dust(ptr_pp%cells%center(:,jb)%lat,ptr_pp%cells%center(:,jb)%lon,i_endidx,dust_tunefac)
        ELSE
          dust_tunefac(:,:) = 1._wp
        ENDIF

        ! Type of convection is required as INTEGER field
        zrg_ktype(1:i_endidx,jb) = NINT(zrg_rtype(1:i_endidx,jb))

        CALL radiation_nwp(               &
                                !
                                ! input
                                ! -----
                                !
          & current_date                     ,&!< in current date
                                ! indices and dimensions
          & jg          =jg                  ,&!< in domain index
          & jb          =jb                  ,&!< in block index
          & irad        =irad                ,&!< in option for radiation scheme (RRTM/PSRAD)
          & jce         =i_endidx            ,&!< in  end   index for loop over block
          & kbdim       =nproma              ,&!< in  dimension of block over cells
          & klev        =nlev_rg             ,&!< in  number of full levels = number of layers
          & klevp1      =nlev_rg+1           ,&!< in  number of half levels = number of layer ifcs
                                !
          & ktype       =zrg_ktype(:,jb)     ,&!< in type of convection
                                !
          & zland       =zrg_fr_land (:,jb)  ,&!< in land mask,     1. over land
          & zglac       =zrg_fr_glac (:,jb)  ,&!< in glacier mask,  1. over land ice
                                !
          & cos_mu0     =zrg_cosmu0  (:,jb)  ,&!< in    cos of zenith angle mu0
          & alb_vis_dir=zrg_albvisdir(:,jb)  ,&!< in    surface albedo for visible range, direct
          & alb_nir_dir=zrg_albnirdir(:,jb)  ,&!< in    surface albedo for near IR range, direct
          & alb_vis_dif=zrg_albvisdif(:,jb)  ,&!< in    surface albedo for visible range, diffuse
          & alb_nir_dif=zrg_albnirdif(:,jb)  ,&!< in    surface albedo for near IR range, diffuse
          & emis_rad   =zrg_emis_rad(:,jb)   ,&!< in longwave surface emissivity
          & tk_sfc     =zrg_tsfc     (:,jb)       ,&!< in    surface temperature
                                !
                                ! atmosphere: pressure, tracer mixing ratios and temperature
          & pp_hl      =zrg_pres_ifc(:,:,jb)    ,&!< in    pressure at half levels at t-dt [Pa]
          & pp_fl      =zrg_pres    (:,:,jb)    ,&!< in    pressure at full levels at t-dt [Pa]
          & tk_fl      =zrg_temp    (:,:,jb)    ,&!< in    temperature at full level at t-dt
          & qm_vap     =zrg_tot_cld (:,:,jb,iqv),&!< in    water vapor mass mixing ratio at t-dt
          & qm_liq     =zrg_tot_cld (:,:,jb,iqc),&!< in    cloud water mass mixing ratio at t-dt
          & qm_ice     =zrg_tot_cld (:,:,jb,iqi),&!< in    cloud ice mass mixing ratio at t-dt
          & qm_o3      = zrg_o3     (:,:,jb)    ,&!< in    O3
          & cdnc       =zrg_acdnc   (:,:,jb)    ,&!< in    cloud droplet numb. conc. [1/m**3]
          & cld_frc    =zrg_clc    (:,:,jb)     ,&!< in    cld_frac = cloud fraction [m2/m2]
          & zaeq1      = zrg_aeq1(:,:,jb)       ,&!< in aerosol continental
          & zaeq2      = zrg_aeq2(:,:,jb)       ,&!< in aerosol maritime
          & zaeq3      = zrg_aeq3(:,:,jb)       ,&!< in aerosol urban
          & zaeq4      = zrg_aeq4(:,:,jb)       ,&!< in aerosol volcano ashes
          & zaeq5      = zrg_aeq5(:,:,jb)       ,&!< in aerosol stratospheric background
          & dust_tunefac = dust_tunefac (:,:)   ,&!< in LW tuning factor for dust aerosol
          & dt_rad     = atm_phy_nwp_config(jg)%dt_rad ,&
                                !
                                ! output
                                ! ------
                                !
          & cld_cvr    =  zrg_aclcov    (:,jb), &      !< out   cloud cover in a column [m2/m2]
          & flx_lw_net =  zrg_lwflxall(:,:,jb), &      !< out terrestrial flux, all sky, net down
          & trsol_net  =  zrg_trsolall(:,:,jb), &      !< out solar transmissivity, all sky, net down
          & flx_uplw_sfc = zrg_lwflx_up_sfc(:,jb), &   !< out longwave upward flux at surface
          & trsol_up_toa = zrg_trsol_up_toa(:,jb), &   !< out upward solar transmissivity at TOA
          & trsol_up_sfc = zrg_trsol_up_sfc(:,jb), &   !< out upward solar transmissivity at surface
          & trsol_par_sfc = zrg_trsol_par_sfc(:,jb), & !< downward transmissivity for photosynthetically active rad. at surface
          & trsol_dn_sfc_diffus = zrg_trsol_dn_sfc_diff(:,jb), &  !< out downward diffuse solar transmissivity at surface
          & trsol_clr_sfc = zrg_trsol_clr_sfc(:,jb)   )  !< out clear-sky net transmissvity at surface (used with reduced grid only)

      ENDDO ! blocks

!$OMP END DO NOWAIT
!$OMP END PARALLEL

      CALL downscale_rad_output(pt_patch%id, pt_par_patch%id,                                     &
        &  nlev_rg, zrg_aclcov, zrg_lwflxall, zrg_trsolall, zrg_trsol_clr_sfc, zrg_lwflx_up_sfc,  &
        &  zrg_trsol_up_toa, zrg_trsol_up_sfc, zrg_trsol_par_sfc, zrg_trsol_dn_sfc_diff,          &
        &  zrg_tsfc, zrg_albdif, zrg_emis_rad, zrg_cosmu0, zrg_tot_cld, zlp_tot_cld, zrg_pres_ifc,&
        &  zlp_pres_ifc, prm_diag%tsfctrad, prm_diag%albdif, aclcov, prm_diag%lwflxall,           &
        &  prm_diag%trsolall, prm_diag%lwflx_up_sfc_rs, prm_diag%trsol_up_toa,                    &
        &  prm_diag%trsol_up_sfc, prm_diag%trsol_par_sfc, prm_diag%trsol_dn_sfc_diff,             &
        &  prm_diag%trsolclr_sfc )

      ! Debug output of radiation output fields
      IF (msg_level >= 16) THEN
        max_lwflx = 0._wp
        min_lwflx = 1.e10_wp
        max_swtrans = 0._wp
        min_swtrans = 1.e10_wp

        rl_start = grf_bdywidth_c + 1
        rl_end   = min_rlcell_int

        i_startblk = pt_patch%cells%start_blk(rl_start,1)
        i_endblk   = pt_patch%cells%end_blk(rl_end,i_nchdom)

        DO jb = i_startblk, i_endblk
         CALL get_indices_c(pt_patch, jb, i_startblk, i_endblk, &
                            i_startidx, i_endidx, rl_start, rl_end)

         DO jk = 1, nlevp1
          max_lwflx(jk)   = MAX(max_lwflx(jk),  MAXVAL(prm_diag%lwflxall(i_startidx:i_endidx,jk,jb)))
          max_swtrans(jk) = MAX(max_swtrans(jk),MAXVAL(prm_diag%trsolall(i_startidx:i_endidx,jk,jb)))
          min_lwflx(jk)   = MIN(min_lwflx(jk),  MINVAL(prm_diag%lwflxall(i_startidx:i_endidx,jk,jb)))
          min_swtrans(jk) = MIN(min_swtrans(jk),MINVAL(prm_diag%trsolall(i_startidx:i_endidx,jk,jb)))
         ENDDO
        ENDDO ! blocks

        max_lwflx = global_max(max_lwflx)
        min_lwflx = global_min(min_lwflx)
        max_swtrans = global_max(max_swtrans)
        min_swtrans = global_min(min_swtrans)

        WRITE(message_text,'(a)') 'max/min LW flux, SW transmissivity'
        CALL message(routine, TRIM(message_text))

        DO jk = 1, nlevp1
          WRITE(message_text,'(i4,2f10.3,2f10.7)') jk,max_lwflx(jk), min_lwflx(jk), &
            max_swtrans(jk), min_swtrans(jk)
          CALL message(routine, TRIM(message_text))
        ENDDO

      ENDIF ! msg_level >= 16

      DEALLOCATE (zrg_cosmu0, zrg_albvisdir, zrg_albnirdir, zrg_albvisdif, zrg_albnirdif, &
        zrg_albdif, zrg_tsfc, zrg_pres_ifc, zrg_pres, zrg_temp, zrg_o3, zrg_ktype,        &
        zrg_aeq1,zrg_aeq2,zrg_aeq3,zrg_aeq4,zrg_aeq5, zrg_acdnc, zrg_tot_cld, zrg_clc,    &
        zrg_aclcov, zrg_lwflxall, zrg_trsolall, zrg_lwflx_up_sfc, zrg_trsol_up_toa,       &
        zrg_trsol_up_sfc, zrg_trsol_par_sfc, zrg_trsol_dn_sfc_diff, zrg_trsol_clr_sfc,    &
        zrg_fr_land, zrg_fr_glac, zrg_emis_rad, zlp_pres_ifc, zlp_tot_cld)

  END SUBROUTINE nwp_rrtm_radiation_reduced
  !---------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------
  !>
  !! @par Revision History
  !! Initial release by Thorsten Reinhardt, AGeoBw, Offenbach (2011-01-13)
  !!
  SUBROUTINE nwp_rrtm_radiation_repartition (current_date, pt_patch, ext_data, &
    &  zaeq1, zaeq2, zaeq3, zaeq4, zaeq5, pt_diag, prm_diag, lnd_prog )

    CHARACTER(len=MAX_CHAR_LENGTH), PARAMETER::  &
      &  routine = modname//'::nwp_rrtm_radiation_repartition'

    TYPE(datetime), POINTER, INTENT(in) :: current_date

    TYPE(t_patch),        TARGET,INTENT(in) :: pt_patch     !<grid/patch info.
    TYPE(t_external_data),INTENT(in)        :: ext_data

    REAL(wp), INTENT(in) :: &
      & zaeq1(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq2(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq3(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq4(nproma,pt_patch%nlev,pt_patch%nblks_c), &
      & zaeq5(nproma,pt_patch%nlev,pt_patch%nblks_c)

    TYPE(t_nh_diag), TARGET, INTENT(in)   :: pt_diag     !<the diagnostic variables
    TYPE(t_nwp_phy_diag),    INTENT(inout):: prm_diag
    TYPE(t_lnd_prog),        INTENT(inout):: lnd_prog

    ! Local scalars:
    INTEGER:: jc,jb
    INTEGER:: jg                !domain id
    INTEGER:: nlev, nlevp1      !< number of full and half levels

    INTEGER:: rl_start, rl_end
    INTEGER:: i_startblk, i_endblk    !> blocks
    INTEGER:: i_startidx, i_endidx    !< slices
    INTEGER:: i_nchdom                !< domain index

    TYPE(t_rrtm_data), POINTER :: rrtm_data
    INTEGER:: return_status
    REAL(wp), POINTER :: &
        & test_aclcov  (:,    :),  &
        & test_lwflxclr(:, :, :),  &
        & test_trsolclr(:, :, :),  &
        & test_lwflxall(:, :, :),  &
        & test_trsolall(:, :, :)
    REAL(wp) :: check_diff

    CHARACTER(*), PARAMETER :: method_name = "nwp_rrtm_radiation_repartition"

    IF (timers_level > 3) CALL timer_start(timer_preradiaton)
    !-------------------------------------------------------------------------
    IF (msg_level >= 12) &
      &  CALL message(routine, 'RRTM radiation on redistributed grid')

    i_nchdom  = MAX(1,pt_patch%n_childdom)
    jg        = pt_patch%id

    ! number of vertical levels
    nlev   = pt_patch%nlev
    nlevp1 = pt_patch%nlevp1

    !-------------------------------------------------------------------------
    !> Radiation
    !-------------------------------------------------------------------------

    rl_start = grf_bdywidth_c+1
    rl_end   = min_rlcell_int

    i_startblk = pt_patch%cells%start_blk(rl_start,1)
    i_endblk   = pt_patch%cells%end_blk(rl_end,i_nchdom)

    IF (test_parallel_radiation) THEN

      ! allocate temp arrays for comparing the results
      ! from the direct radiation call and the redistributed
      ! radiation call
      ALLOCATE( &
        & test_aclcov  (nproma,         pt_patch%nblks_c),  &
        & test_lwflxclr(nproma, nlevp1, pt_patch%nblks_c),  &
        & test_trsolclr(nproma, nlevp1, pt_patch%nblks_c),  &
        & test_lwflxall(nproma, nlevp1, pt_patch%nblks_c),  &
        & test_trsolall(nproma, nlevp1, pt_patch%nblks_c),  &
        & STAT=return_status)

!$OMP PARALLEL
!$OMP DO PRIVATE(jb,jc,i_startidx,i_endidx) ICON_OMP_GUIDED_SCHEDULE
      DO jb = i_startblk, i_endblk

        CALL get_indices_c(pt_patch, jb, i_startblk, i_endblk, &
          &                         i_startidx, i_endidx, rl_start, rl_end)

        ! It may happen that an MPI patch contains only nest boundary points
        ! In this case, no action is needed
        IF (i_startidx > i_endidx) CYCLE

        prm_diag%tsfctrad(1:i_endidx,jb) = lnd_prog%t_g(1:i_endidx,jb)

        CALL radiation(               &
                                !
                                ! input
                                ! -----
                                !
          & current_date                     ,&!< in current date
                                ! indices and dimensions
          & jg         =jg                   ,&!< in domain index
          & jb         =jb                   ,&!< in block index
          & jce        =i_endidx             ,&!< in  end   index for loop over block
          & kbdim      =nproma               ,&!< in  dimension of block over cells
          & klev       =nlev                 ,&!< in  number of full levels = number of layers
          & klevp1     =nlevp1               ,&!< in  number of half levels = number of layer ifcs
                                !
          & ktype      =prm_diag%ktype(:,jb) ,&!< in     type of convection
                                !
                                ! surface: albedo + temperature
          & zland      =ext_data%atm%fr_land_smt(:,jb)   ,&!< in     land fraction
          & zglac      =ext_data%atm%fr_glac_smt(:,jb)   ,&!< in     land glacier fraction
                                !
          & cos_mu0    =prm_diag%cosmu0  (:,jb) ,&!< in  cos of zenith angle mu0
          & alb_vis_dir=prm_diag%albvisdir(:,jb) ,&!< in surface albedo for visible range, direct
          & alb_nir_dir=prm_diag%albnirdir(:,jb) ,&!< in surface albedo for near IR range, direct
          & alb_vis_dif=prm_diag%albvisdif(:,jb) ,&!< in surface albedo for visible range, diffuse
          & alb_nir_dif=prm_diag%albnirdif(:,jb) ,&!< in surface albedo for near IR range, diffuse
          & emis_rad=ext_data%atm%emis_rad(:,jb) ,&!< in longwave surface emissivity
          & tk_sfc     =prm_diag%tsfctrad(:,jb)  ,&!< in surface temperature
                                !
                                ! atmosphere: pressure, tracer mixing ratios and temperature
          & pp_hl      =pt_diag%pres_ifc  (:,:,jb)     ,&!< in  pres at half levels at t-dt [Pa]
          & pp_fl      =pt_diag%pres      (:,:,jb)     ,&!< in  pres at full levels at t-dt [Pa]
          & tk_fl      =pt_diag%temp      (:,:,jb)     ,&!< in  temperature at full level at t-dt
          & qm_vap     =prm_diag%tot_cld  (:,:,jb,iqv) ,&!< in  water vapor mass mix ratio at t-dt
          & qm_liq     =prm_diag%tot_cld  (:,:,jb,iqc) ,&!< in cloud water mass mix ratio at t-dt
          & qm_ice     =prm_diag%tot_cld  (:,:,jb,iqi) ,&!< in cloud ice mass mixing ratio at t-dt
          & qm_o3      =ext_data%atm%o3   (:,:,jb)     ,&!< in o3 mass mixing ratio at t-dt
          & cdnc       =prm_diag%acdnc    (:,:,jb)     ,&!< in  cloud droplet numb conc. [1/m**3]
          & cld_frc    =prm_diag%clc      (:,:,jb)     ,&!< in  cloud fraction [m2/m2]
          & zaeq1      = zaeq1(:,:,jb)                 ,&!< in aerosol continental
          & zaeq2      = zaeq2(:,:,jb)                 ,&!< in aerosol maritime
          & zaeq3      = zaeq3(:,:,jb)                 ,&!< in aerosol urban
          & zaeq4      = zaeq4(:,:,jb)                 ,&!< in aerosol volcano ashes
          & zaeq5      = zaeq5(:,:,jb)                 ,&!< in aerosol stratospheric background
          & dt_rad     = atm_phy_nwp_config(jg)%dt_rad ,&
                                ! output
                                ! ------
                                !
          & cld_cvr        = test_aclcov             (:,jb),&!< out cloud cover in a column [m2/m2]
          & flx_lw_net_clr = test_lwflxclr(:,:,jb),&!< out terrestrial flux, clear sky, net down
          & trm_sw_net_clr = test_trsolclr(:,:,jb),&!< out sol. transmissivity, clear sky, net down
          & flx_lw_net     = test_lwflxall(:,:,jb),&!< out terrestrial flux, all sky, net down
          & trm_sw_net     = test_trsolall(:,:,jb),&!< out solar transmissivity, all sky, net down
          & opt_halo_cosmu0 = .FALSE. )

      ENDDO ! blocks
!$OMP END DO NOWAIT
!$OMP END PARALLEL

    ENDIF !test_parallel_radiation

    IF (timers_level > 3) THEN
      CALL timer_stop(timer_preradiaton)
      CALL timer_start(timer_radiaton_recv)
    ENDIF

    ! this maybe used, so we fill it while not an output of the radiation
    prm_diag%tsfctrad(:,:) = lnd_prog%t_g(:,:)
    CALL recv_rrtm_input( &
          & ktype      =prm_diag%ktype(:,:)             ,&!< in     type of convection
          & zland      =ext_data%atm%fr_land_smt(:,:)   ,&!< in     land fraction
          & zglac      =ext_data%atm%fr_glac_smt(:,:)   ,&!< in     land glacier fraction
                                !
          & cos_mu0    =prm_diag%cosmu0  (:,:) ,&!< in  cos of zenith angle mu0
          ! these are actualy computed from albvisdif
!          & alb_vis_dir=albvisdir        (:,:) ,&!< in surface albedo for visible range, direct
!          & alb_nir_dir=albnirdir        (:,:) ,&!< in surface albedo for near IR range, direct
          & alb_vis_dif=prm_diag%albvisdif(:,:),&!< in surface albedo for visible range, diffuse
!          & alb_nir_dif=prm_diag%albnirdif(:,:) ,&!< in surface albedo for near IR range, diffuse
          & emis_rad=ext_data%atm%emis_rad(:,:),&!< in longwave surface emissivity
          & tk_sfc     =prm_diag%tsfctrad(:,:) ,&!< in surface temperature
                                !
                                ! atmosphere: pressure, tracer mixing ratios and temperature
          & pp_hl      =pt_diag%pres_ifc  (:,:,:)     ,&!< in  pres at half levels at t-dt [Pa]
          & pp_fl      =pt_diag%pres      (:,:,:)     ,&!< in  pres at full levels at t-dt [Pa]
          & tk_fl      =pt_diag%temp      (:,:,:)     ,&!< in  temperature at full level at t-dt
          & qm_vap     =prm_diag%tot_cld  (:,:,:,iqv) ,&!< in  water vapor mass mix ratio at t-dt
          & qm_liq     =prm_diag%tot_cld  (:,:,:,iqc) ,&!< in cloud water mass mix ratio at t-dt
          & qm_ice     =prm_diag%tot_cld  (:,:,:,iqi) ,&!< in cloud ice mass mixing ratio at t-dt
          & qm_o3      =ext_data%atm%o3   (:,:,:)     ,&!< in o3 mass mixing ratio at t-dt
          & cdnc       =prm_diag%acdnc    (:,:,:)     ,&!< in  cloud droplet numb conc. [1/m**3]
          & cld_frc    =prm_diag%clc      (:,:,:)     ,&!< in  cloud fraction [m2/m2]
          & zaeq1      = zaeq1(:,:,:)                 ,&!< in aerosol continental
          & zaeq2      = zaeq2(:,:,:)                 ,&!< in aerosol maritime
          & zaeq3      = zaeq3(:,:,:)                 ,&!< in aerosol urban
          & zaeq4      = zaeq4(:,:,:)                 ,&!< in aerosol volcano ashes
          & zaeq5      = zaeq5(:,:,:)                 ,&!< in aerosol stratospheric background
          & patch      = pt_patch                     ,&!< in
          & rrtm_data  = rrtm_data)                     !< out, pointer to rrtm input values

    IF (timers_level > 3) THEN
      CALL timer_stop(timer_radiaton_recv)
      CALL timer_start(timer_radiaton_comp)
    ENDIF

!$OMP PARALLEL
!$OMP DO PRIVATE(jb,jc,i_startidx,i_endidx) ICON_OMP_GUIDED_SCHEDULE
    DO jb = 1, rrtm_data%no_of_blocks

      i_endidx = MERGE(rrtm_data%block_size, rrtm_data%end_index, &
           jb /= rrtm_data%no_of_blocks)

      !Calculate direct albedo from diffuse albedo and solar zenith angle
      !formula as in Ritter-Geleyn's fesft
      DO jc = 1, i_endidx
        rrtm_data%albedo_vis_dir(jc,jb) =  ( 1.0_wp                                                           &
          &  + 0.5_wp * (rrtm_data%cosmu0(jc,jb) * (1.0_wp/rrtm_data%albedo_vis_dif(jc,jb) - 1.0_wp))) &
          & / (1.0_wp + (rrtm_data%cosmu0(jc,jb) * (1.0_wp/rrtm_data%albedo_vis_dif(jc,jb) - 1.0_wp)))**2
      ENDDO
!      IF (i_startidx > 1) albvisdir(1:i_startidx-1,jb) = albvisdir(i_startidx,jb)

      ! no distiction between vis and nir albedo
      rrtm_data%albedo_nir_dir(1:i_endidx,jb) = rrtm_data%albedo_vis_dir(1:i_endidx,jb)
      rrtm_data%albedo_nir_dif(1:i_endidx,jb) = rrtm_data%albedo_vis_dif(1:i_endidx,jb)

      CALL radiation(               &
                              !
                              ! input
                              ! -----
                              !
        & current_date                       ,&!< in current date
                              ! indices and dimensions
        & jg         =jg                     ,&!< in domain index
        & jb         =jb                     ,&!< in block index
        & jce        = i_endidx              ,&!< in  end   index for loop over block
        & kbdim      = rrtm_data%block_size  ,&!< in  dimension of block over cells
        & klev       = rrtm_data%full_levels ,&!< in  number of full levels =  number of layers
        & klevp1     = rrtm_data%half_levels ,&!< in  number of half levels =  number of layer ifcs
                              !
        & ktype      = rrtm_data%convection_type(:,jb) ,&!< in     type of convection
                              !
                              ! surface: albedo + temperature
        & zland      = rrtm_data%fr_land_smt(:,jb)     ,&!< in     land fraction
        & zglac      = rrtm_data%fr_glac_smt(:,jb)     ,&!< in     land glacier fraction
                              !
        & cos_mu0    = rrtm_data%cosmu0         (:,jb) ,&!< in  cos of zenith angle mu0
        & alb_vis_dir= rrtm_data%albedo_vis_dir (:,jb) ,&!< in surface albedo for visible range, direct
        & alb_nir_dir= rrtm_data%albedo_nir_dir (:,jb) ,&!< in surface albedo for near IR range, direct
        & alb_vis_dif= rrtm_data%albedo_vis_dif (:,jb) ,&!< in surface albedo for visible range, diffuse
        & alb_nir_dif= rrtm_data%albedo_nir_dif (:,jb) ,&!< in surface albedo for near IR range, diffuse
        & emis_rad   = rrtm_data%emis_rad       (:,jb) ,&!< in longwave surface emissivity
        & tk_sfc     = rrtm_data%tsfctrad       (:,jb) ,&!< in surface temperature
                              !
                              ! atmosphere: pressure, tracer mixing ratios and temperature
        & pp_hl      = rrtm_data%pres_ifc  (:,:,jb)    ,&!< in  pres at half levels at t-dt [Pa]
        & pp_fl      = rrtm_data%pres      (:,:,jb)    ,&!< in  pres at full levels at t-dt [Pa]
        & tk_fl      = rrtm_data%temp      (:,:,jb)    ,&!< in  temperature at full level at t-dt
        & qm_vap     = rrtm_data%qm_vapor  (:,:,jb)    ,&!< in  water vapor mass mix ratio at t-dt
        & qm_liq     = rrtm_data%qm_liquid (:,:,jb)    ,&!< in cloud water mass mix ratio at t-dt
        & qm_ice     = rrtm_data%qm_ice    (:,:,jb)    ,&!< in cloud ice mass mixing ratio at t-dt
        & qm_o3      = rrtm_data%qm_o3     (:,:,jb)    ,&!< in o3 mass mixing ratio at t-dt
        & cdnc       = rrtm_data%acdnc     (:,:,jb)    ,&!< in  cloud droplet numb conc. [1/m**3]
        & cld_frc    = rrtm_data%cld_frc   (:,:,jb)    ,&!< in  cloud fraction [m2/m2]
        & zaeq1      = rrtm_data%zaeq1     (:,:,jb)    ,&!< in aerosol continental
        & zaeq2      = rrtm_data%zaeq2     (:,:,jb)    ,&!< in aerosol maritime
        & zaeq3      = rrtm_data%zaeq3     (:,:,jb)    ,&!< in aerosol urban
        & zaeq4      = rrtm_data%zaeq4     (:,:,jb)    ,&!< in aerosol volcano ashes
        & zaeq5      = rrtm_data%zaeq5     (:,:,jb)    ,&!< in aerosol stratospheric background
        & dt_rad     = atm_phy_nwp_config(jg)%dt_rad   ,&
                              ! output
                              ! ------
                              !
        & cld_cvr        = rrtm_data%aclcov  (:,  jb),&!< out cloud cover in a column [m2/m2]
        & flx_lw_net_clr = rrtm_data%lwflxclr(:,:,jb),&!< out terrestrial flux, clear sky, net down
        & trm_sw_net_clr = rrtm_data%trsolclr(:,:,jb),&!< out sol. transmissivity, clear sky, net down
        & flx_lw_net     = rrtm_data%lwflxall(:,:,jb),&!< out terrestrial flux, all sky, net down
        & trm_sw_net     = rrtm_data%trsolall(:,:,jb),&!< out solar transmissivity, all sky, net down
        & opt_halo_cosmu0 = .FALSE. )
    ENDDO
!$OMP END DO NOWAIT
!$OMP END PARALLEL

    IF (timers_level > 3) THEN
      CALL timer_stop(timer_radiaton_comp)
      CALL timer_start(timer_radiaton_send)
    ENDIF

    ! aclcov is also output but not used
    CALL send_rrtm_output(        &
      & rrtm_data               , &
      & prm_diag%lwflxall(:,:,:), &!< out terrestrial flux, all sky, net down
      & prm_diag%trsolall(:,:,:), &!< out sol. transmissivity, all sky, net down
      & prm_diag%lwflxall(:,:,:), &!< out terrestrial flux, all sky, net down
      & prm_diag%trsolall(:,:,:))  !< out solar transmissivity, all sky, net down

    IF (timers_level > 3) &
      & CALL timer_stop(timer_radiaton_send)

    IF (test_parallel_radiation) THEN
      DO jb = i_startblk, i_endblk

        CALL get_indices_c(pt_patch, jb, i_startblk, i_endblk, &
          &                         i_startidx, i_endidx, rl_start, rl_end)

        DO jc = i_startidx, i_endidx

          check_diff = MAXVAL(ABS(prm_diag%lwflxall(jc,:,jb) - test_lwflxall(jc,:,jb)))
          IF (check_diff > 0.0_wp) THEN
            write(0,*) " jc,jb=", jc,jb
            write(0,*) " prm_diag%lwflxall=", prm_diag%lwflxall(jc,:,jb)
            write(0,*) " test_lwflxall=", test_lwflxall(jc,:,jb)
            CALL finish(method_name,"lwflxall differs")
          ENDIF

          check_diff = MAXVAL(ABS(prm_diag%trsolall(jc,:,jb) - test_trsolall(jc,:,jb)))
          IF (check_diff > 0.0_wp) THEN
            write(0,*) " jc,jb=", jc,jb
            write(0,*) " prm_diag%trsolall=", prm_diag%trsolall(jc,:,jb)
            write(0,*) " test_trsolall=", test_trsolall(jc,:,jb)
            CALL finish(method_name,"trsolall differs")
          ENDIF

       ENDDO
     ENDDO

     DEALLOCATE(test_aclcov, test_lwflxclr, test_trsolclr, test_lwflxall, test_trsolall)

   ENDIF ! test_parallel_radiation

  END SUBROUTINE nwp_rrtm_radiation_repartition
  !---------------------------------------------------------------------------------------

END MODULE mo_nwp_rrtm_interface

max_stars_repo_name: niklaswr/TSMP

max_issues_repo_name: niklaswr/TSMP

max_forks_repo_name: niklaswr/TSMP


text: Sagehen Creek Field Station is a UC Nature Reserves UC Nature Reserve administered by UC Berkeley. It is located about 20 miles north of Lake Tahoe and serves as a hub for Sierra Nevada Mountains Sierra Nevada research and education. During the summer, UC Davis uses the facilities for programs in entomology, botany and Nature and Culture.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      MODULE m_inwint
      CONTAINS
      SUBROUTINE inwint(
     >                  e,fl,ki,fkap,cs,cis,s,z,h,dd,rn,rnot,msh,vr,
     X                  a,b,
     <                  kj)
c
c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
c-----x  perform the  inward integration.   this routine is a     x----
c-----x  derivative of start2/diff2.   much unnecessary indexing  x----
c-----x  has been removed among other things.    dale koelling    x----
c
c     adams' procedure. more stable for GGA. Feb.20,1998. T.A.
c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
      use m_juDFT
      IMPLICIT NONE

c     .. Arguments ..
      INTEGER,INTENT (IN)  :: msh,ki
      INTEGER,INTENT (OUT) :: kj
      REAL,   INTENT (IN)  :: e,fl,fkap,cs,cis,s,z,h,dd,rn,rnot
      REAL,   INTENT (IN)  :: vr(msh)
c     ..
      REAL,   INTENT (INOUT) :: a(msh),b(msh)
c     ..
c     .. Local Scalars ..
      REAL atk,btk,d,da1,da2,da3,da4,db1,db2,db3,db4,eabsv,eps,f1,f2,f3,
     +     f4,f5,fg1,fg2,fg3,g,g1,g2,g3,g4,g5,h3,q11,q22,r,rp12,rp21
      INTEGER I,k,l,n
      REAL :: da(5),db(5)
c     ..
c     .. intrinsic functions ..
      INTRINSIC abs,mod,sqrt
c     ..
c     .. equivalences ..
cta+
      EQUIVALENCE (da1,da(1)), (da2,da(2)), (da3,da(3)), (da4,da(4)),
     &  (da5,da(5))
      EQUIVALENCE (db1,db(1)), (db2,db(2)), (db3,db(3)), (db4,db(4)),
     &  (db5,db(5))

      REAL c1,c5,c9,c19, t9,t37,t55,t59, da5,db5
cta-
c     ..
c     .. data statements ..
c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
c-----x this version modified to accept positive energies.        x----
c-----x                              d.d.koelling   7/7/77      x----x
c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
c---* eps    value used to determine the practical infinity
      DATA f1/1.045833333e0/,f2/2.691666666e0/,f3/1.1e0/
      DATA f4/.4416666666e0/,f5/.07916666666e0/
      DATA g1/.9666666666e0/,g2/4.133333333e0/,g3/.8e0/
      DATA g4/.1333333333e0/,g5/.03333333333e0/
      DATA fg1/.9333333333e0/,fg2/4.266666666e0/,fg3/1.6e0/
      DATA eps/75.0e0/
c     ..
      h3 = -h/3.0e0
      n = msh
cta+
      t55=55.e0/8.e0
      t59=59.e0/8.e0
      t37=37.e0/8.e0
      t9=9.e0/8.e0

      c19=19.e0/8.e0
      c9=9.e0/8.e0
      c5=5.e0/8.e0
      c1=1.e0/8.e0
cta-

c     e: energy
c     ki: the grid point corresponding to classical terning point
cc        at which potential exceeded the energy.
cc        ki=249 for Fe and l=0.
c     kj: undefined at the beginning. finally the farthest grid
cc        point which can be considered as infinity for the
cc        level to be obtained.

      d = 1.0e0/dd
      g = sqrt(fkap**2- (cis*z)**2)
      q11 = s*fkap - g
      q22 = -g - s*fkap
      kj = ki + 10
      r = rnot* (dd**kj)

   10 kj = kj + 1
      IF (kj>msh)  CALL juDFT_error("inward integration failed",calledby
     +     ="inwint",hint
     +     ="This often indicates that the potential is unphysical")

c     eps is a big energy. thus if the next statement is fulfilled,
cc    the grid is considered to be infinitely far.

      IF (r* (vr(kj)-e*r).gt.eps) GOTO 20
      r = r*dd
      IF (kj.lt.n) go to 10
      b(kj) = mod(fl,2.0e0)
      a(kj) = 1.0e0 - b(kj)
      go to 30
   20 a(kj) = 1.0e0
      eabsv = abs(e)
      b(kj) = s*sqrt(eabsv/ (2*cs*cs+e))
   30 continue
c  30 do 40 l = 1,4

c     kj=338 for Fe l=0 and e=-333.
c     a(kj)=1.0

      do 40 l = 1,4
         k = kj - l
         a(k) = a(kj)
         b(k) = b(kj)
   40 continue

c     a(k)=1.0 for k=kj-4,..,kj.

      do 70 i = 1,6
         k = kj + 1
         r = rn*d** (n-k)
         do 60 l = 1,5
            k = k - 1
            r = r*d
            rp21 = cis* (e*r-vr(k))
            rp12 = -rp21 - 2.0e0*cs*r
            da(l) = h3* (q11*a(k)+rp12*b(k))
            db(l) = h3* (q22*b(k)+rp21*a(k))
   60    continue
         l = kj - 1
         a(l) = a(kj) + f1*da1 + f2*da2 + f4*da4 - (f3*da3+f5*da(5))
         b(l) = b(kj) + f1*db1 + f2*db2 + f4*db4 - (f3*db3+f5*db(5))
         l = l - 1
         a(l) = a(kj) + g1*da1 + g2*da2 + g3*da3 + g4*da4 - g5*da(5)
         b(l) = b(kj) + g1*db1 + g2*db2 + g3*db3 + g4*db4 - g5*db(5)
         l = l - 1
         a(l) = a(kj) + 1.0125e0*da1 + 3.825e0*da2 + 2.7e0*da3 +
     +          1.575e0*da4 - .1125e0*da(5)
         b(l) = b(kj) + 1.0125e0*db1 + 3.825e0*db2 + 2.7e0*db3 +
     +          1.575e0*db4 - .1125e0*db(5)
         l = l - 1
         a(l) = a(kj) + fg1*da1 + fg2*da2 + fg3*da3 + fg2*da4 +
     +          fg1*da(5)
         b(l) = b(kj) + fg1*db1 + fg2*db2 + fg3*db3 + fg2*db4 +
     +          fg1*db(5)
   70 continue

c     a(kj)=1.0, a(kj-1),..,a(kj-4) have been obtained.

c     adams' modification here on.

      k = kj - 3
      r = rn*d** (n-k)
   80 k = k - 1
      r = r*d
      rp21 = cis* (e*r-vr(k))
      rp12 = -rp21 - 2.0e0*cs*r

cta+
c     adams' predictor.
      atk=a(k+1) + t55*da4-t59*da3+t37*da2-t9*da1
      btk=b(k+1) + t55*db4-t59*db3+t37*db2-t9*db1

      da1 = da2
      da2 = da3
      da3 = da4

      db1 = db2
      db2 = db3
      db3 = db4

cta-
      da4=h3*(q11*atk+rp12*btk)
      db4=h3*(q22*btk+rp21*atk)

cta+
c     adams interpolation formula 25.5.5 (Handbook of math.functions).
      a(k)=a(k+1) + c9*da4+c19*da3-c5*da2+c1*da1
      b(k)=b(k+1) + c9*db4+c19*db3-c5*db2+c1*db1
cta-

      da4 = h3* (q11*a(k)+rp12*b(k))
      db4 = h3* (q22*b(k)+rp21*a(k))

      if(k.gt.ki) go to 80

      return
      END SUBROUTINE
      end

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
!***************************************************************************************************
! eos_bahcall.f90 10/18/17
! EoS Replacement based on Bahcall ().
! This file contains routines which calculate EoS quantites needed to calculate screening
! corrections for reaction rates.
!***************************************************************************************************

Module xnet_eos
  Implicit None

Contains

  Subroutine eos_initialize
    !-----------------------------------------------------------------------------------------------
    ! This routine initializes the EoS.
    !-----------------------------------------------------------------------------------------------
    Implicit None
    Return
  End Subroutine eos_initialize

  Subroutine eos_interface(rho,t9,y,ye,cv,efermkt,defermktdt9)
    !-----------------------------------------------------------------------------------------------
    ! This routine updates the Equation of State for changes in temperature and density.
    !
    ! Bahcall introduced a fit for the plasma parameters needed to calculate the screening factors
    ! for nuclear reaction based on the easier to calculate efermkt (the ratio of the electron
    ! Fermi Energy to the themal energy) instead of the more general expresions based on etae (the
    ! ratio of the electron chemical potential to the thermal energy), which generally require a
    ! complex Equation of State to calculate.   For basic cases, XNet uses Bahcall's approach
    ! to simplify the
    !
    ! Note: the factors efermkt and defermktdt9 (the derivative of efermkt wrt to temperature) in the
    ! subroutine interface are stored in XNet's conditions module as etae and detaedt9.  For more
    ! general EoS, etae and detaedt9 are their proper selves.
    !-----------------------------------------------------------------------------------------------
    Use nuclear_data, Only: zz53
    Use xnet_constants, Only: asig, avn, bok, e2, pi, pi2, third, emass, clt, ele_en, hbar
    Use xnet_controls, Only: idiag, iheat, iscrn, lun_diag
    Use xnet_types, Only: dp
    Use xnet_abundances, Only: y_moment
    Implicit None

    ! Input variables
    Real(dp), Intent(in) :: t9, rho, y(:)

    ! Output variables
    Real(dp), Intent(out) :: ye, cv, efermkt, defermktdt9

    ! Local variables
    Real(dp), Parameter :: a1 = 0.898004, b1 = 0.96786, c1 = 0.220703, d1 = -0.86097
    Real(dp), Parameter :: a2 = -0.86602540378, b2 = 0.29561, c2 = 1.9885
    Real(dp) :: ytot, bkt, abar, zbar, z2bar, zibar, z53bar
    Real(dp) :: etae_mb, ae, gam, gam14, gam32, gamc2, rel_ef
    Real(dp) :: eion, erad, ecoul
    Real(dp) :: deiondt9, deraddt9, decouldt9

    ! Calculate Ye and other needed moments of the abundance distribution
    Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar)

    If ( iscrn > 0 .or. iheat > 0 ) Then

      bkt = bok*t9
      etae_mb = log(rho*avn*ye * (2.0*pi*hbar**2 / (emass*bkt))**1.5)
      rel_ef = hbar * (3.0*pi2*rho*avn*ye)**third / (emass*clt)
      efermkt = ele_en * (sqrt(1.0 + rel_ef**2) - 1.0) / bkt
      defermktdt9 = -efermkt/t9

      ! Calculate energy derivatives (assume etot = eion + erad + ecoul) (ignore degenerate electron energy)
      eion = 1.5*bkt*ytot
      deiondt9 = 1.5*ytot*bok

      erad = asig*bkt*bkt*bkt*bkt / (rho*avn)
      deraddt9 = 4.0*erad / t9

      ae = (3.0 / (4.0*pi*avn*rho*ye))**third ! electron-sphere radius
      z53bar = sum(zz53 * y) / ytot
      gam = z53bar*e2 / (ae*bkt) ! ionic Coulomb coupling parameter
      If ( gam >= 1.0 ) Then
        gam14 = gam**0.25
        ecoul = (a1*gam + b1*gam14 + c1/gam14 + d1) * bkt * ytot
        decouldt9 = ecoul/t9 - bok * ytot * (a1*gam + 0.25 * (b1*gam14 - c1/gam14))
      Else
        gam32 = gam*sqrt(gam)
        gamc2 = gam**c2
        ecoul = (a2*gam32 + b2*gamc2) * bkt * ytot
        decouldt9 = ecoul/t9 - bok * ytot * (1.5*a2*gam32 + b2*c2*gamc2)
      EndIf
      cv = deiondt9 + deraddt9 + decouldt9
    Else
      efermkt = 0.0
      defermktdt9 = 0.0
      cv = 0.0
    EndIf

    If ( idiag >= 3 ) Write(lun_diag,"(a,6es23.15)") 'EOS',t9,rho,ye,cv,efermkt,defermktdt9

    Return
  End Subroutine eos_interface

  Subroutine eos_screen(t9,rho,y,efermkt,defermktdt9,ztilde,zinter,lambda0,gammae,dztildedt9)
    !-----------------------------------------------------------------------------------------------
    ! This routine Bahcall's approach to calculate the factors needed for screening from the input
    ! temperature, density and composition.
    !
    ! Note: the efermkt and defermktdt9 in the subroutine interface are stored in XNet's conditions
    ! module as etae and detaedt9.  For more general EoS, etae and detaedt9 are their proper selves.
    !-----------------------------------------------------------------------------------------------
    Use xnet_constants, Only: avn, bok, clt, e2, ele_en, emass, hbar, pi, pi2, third, two3rd, &
      & thbim2, twm2bi
    Use xnet_controls, Only: idiag, iheat, lun_diag
    Use xnet_types, Only: dp
    Use xnet_abundances, Only: y_moment
    Implicit None

    ! Input variables
    Real(dp), Intent(in) :: t9, rho, y(:), efermkt, defermktdt9

    ! Output variables
    Real(dp), Intent(out) :: ztilde, zinter, lambda0, gammae, dztildedt9

    ! Local variables
    Real(dp) :: ye, ytot, bkt, abar, zbar, z2bar, zibar
    Real(dp) :: sratio, ae, dsratiodefermkt

    ! Calculate Ye and other needed moments of the abundance distribution
    Call y_moment(y,ye,ytot,abar,zbar,z2bar,zibar)

    ! Calculate ratio f'/f for electrons (Salpeter, Eq. 24)
    Call salpeter_ratio(efermkt,sratio,dsratiodefermkt)
    ztilde = sqrt(z2bar + zbar*sratio)
    If ( iheat > 0 ) Then
      dztildedt9 = 0.5*zbar/ztilde * dsratiodefermkt*defermktdt9
    Else
      dztildedt9 = 0.0
    EndIf

    ! Calculate plasma quantities
    bkt = bok*t9
    lambda0 = sqrt(4.0*pi*rho*avn*ytot) * (e2/bkt)**1.5 ! DGC, Eq. 3
    ae = (3.0 / (4.0*pi*avn*rho*ye))**third ! electron-sphere radius
    gammae = e2 / (ae*bkt) ! electron Coulomb coupling parameter
    zinter = zibar / (ztilde**thbim2 * zbar**twm2bi)
    If ( idiag >= 2 ) Write(lun_diag,"(a14,9es23.15)") 'EOS Screen', &
      & t9,rho,ye,z2bar,zbar,sratio,ztilde,ztilde*lambda0,gammae

    Return
  End Subroutine eos_screen

  Subroutine salpeter_ratio(efmkt,ratio,dratiodefmkt)
    !-----------------------------------------------------------------------------------------------
    ! This routine calculates the Salpeter (1954) ratio f'/f(eta) needed for electron screening, using
    ! a fit to Figure 24 in that paper. efmkt is the ratio of electron chemical potential to kT.
    !-----------------------------------------------------------------------------------------------
    Use xnet_constants, Only: ln_10
    Use xnet_controls, Only: iheat
    Use xnet_types, Only: dp
    Implicit None

    ! Input variables
    Real(dp), Intent(in) :: efmkt

    ! Output variables
    Real(dp), Intent(out) :: ratio, dratiodefmkt

    ! Local variables
    Real(dp), Parameter :: a0 = 1.578347, c0 = 0.75793, c1 = -0.54621
    Real(dp), Parameter ::  c2 = -0.30964,  c3 = 0.12535,  c4 = 0.1203,  c5 = -0.012857,  c6 = -0.014768
    Real(dp), Parameter :: dc2 = -0.61928, dc3 = 0.37605, dc4 = 0.4812, dc5 = -0.064285, dc6 = -0.088608
    Real(dp) :: lefmkt

    lefmkt = log10(efmkt)
    If ( lefmkt <= -a0 ) Then ! Bahcall uses limit of -2, but this gives ratio slightly above 1.
      ratio = 1.0
    ElseIf ( lefmkt >= 1.5 ) Then
      ratio = 0.0
    Else
      ratio = c0 + (c1*lefmkt) + (c2*lefmkt**2) + (c3*lefmkt**3) + (c4*lefmkt**4) + (c5*lefmkt**5) + (c6*lefmkt**6)
    Endif

    If ( iheat > 0 .and. lefmkt > a0 .and. lefmkt < 1.5 ) Then
      dratiodefmkt = c1 + dc2*lefmkt + dc3*lefmkt**2 + dc4*lefmkt**3 + dc5*lefmkt**4 + dc6*lefmkt**5
      dratiodefmkt = dratiodefmkt / (ln_10*efmkt)
    Else
      dratiodefmkt = 0.0
    EndIf

    Return
  End Subroutine salpeter_ratio

End Module xnet_eos

max_stars_repo_name: vishaltiwari/XNet

max_issues_repo_name: vishaltiwari/XNet

max_forks_repo_name: vishaltiwari/XNet


text:
c *******************************************************************
      subroutine tc4nosedr(winsl1,winss1,winsl2,winss2,
     +   side,sl1,ss1,nl1,ns1,sl2,ss2,nl2,ns2,
     +   vrdi_unit,winnl,winns,bottom,right,zoom1,zoom2,
     +   lbox,pair,graphicsplane)
c move correlation window & get cursor location.
c side=0 means first pass for new tiepoint
c side=1 left side
c side=2 right side

      integer zoom1,zoom2
      integer four,five,zero,one,xcurs,ycurs
      integer oldxcurs,oldycurs,value
      integer two,three,gstatus,pair
      integer xdclocation,xdipolyline,xdxswitch,xdcset
      integer winnl,winns,winsl1,winss1,winsl2,winss2,side
      integer sl1,ss1,sl2,ss2,bottom,right
      integer   xbox(5),ybox(5),vrdi_unit
      byte      i255, i0
      data one/1/,four/4/,i255/255/,five/5/,zero/0/,two/2/
      data three/3/, i0/0/
      integer graphicsplane

      pair=0
      side=0
      nr=right/6
      xcurs=right/4
      ycurs=(nl1/zoom1)/2
      oldxcurs=xcurs
      oldycurs=ycurs
      gstatus=xdcset(vrdi_unit,one,xcurs,ycurs)
      if(gstatus.ne.1)then
          call prnt(4,1,gstatus,'xdcset error #.')
          call prnt(4,1,xcurs,'xcurs value=.')
          call prnt(4,1,ycurs,'ycurs value=.')
c         call abend
      endif
      winsl1=ycurs*zoom1+sl1-1
      winss1=xcurs*zoom1+ss1-1
      call get_square_coord(xbox,ybox,xcurs,ycurs,
     +         				winns/zoom1,winnl/zoom1)
      gstatus=xdipolyline(vrdi_unit,graphicsplane,i255,five,xbox,ybox)
      if(gstatus.ne.1)then
          call prnt(4,1,gstatus,'xdipolyline error #.')
c         call abend
      endif

      winsl2=(ycurs*zoom2)+sl2-1
      winss2=(xcurs*zoom2)+ss2-1
      winsl1=(ycurs*zoom1)+sl1-1
      winss1=(xcurs*zoom1)+ss1-1

c main loop
ccc 30    call vwait(100)

30	continue
      value=0

c     Check for QUIT condition (left button)
      gstatus=xdxswitch(vrdi_unit,one,one,value)
      if(gstatus.ne.1)then
          call prnt(4,1,gstatus,'xdxswitch error #.')
c         call abend
      endif
      if(value.eq.1)then
         if(side.eq.0)then
            gstatus=xdipolyline(vrdi_unit,graphicsplane,i0,
     +        				five,xbox,ybox)
            if(gstatus.ne.1)then
                call prnt(4,1,gstatus,'xdipolyline error #.')
c               call abend
            endif
         endif
         return
      endif

c     Check for LEFT side condition (left button)
      gstatus=xdxswitch(vrdi_unit,one,two,value)
      if(gstatus.ne.1)then
         call prnt(4,1,gstatus,'xdxswitch error #.')
c         call abend
      endif
      if(value.eq.1)then
         value=0
         side=1
         xcurs=(winss1-ss1+1)/zoom1
         ycurs=(winsl1-sl1+1)/zoom1
         gstatus=xdcset(vrdi_unit,one,xcurs,ycurs)
         if(gstatus.ne.1)then
             call prnt(4,1,gstatus,'xdcset error #.')
c            call abend
         endif
         call get_square_coord(xbox,ybox,xcurs,ycurs,
     +         				winns/zoom1,winnl/zoom1)
        goto 30
      endif

c     Check for RIGHT side condition (third button)
      gstatus=xdxswitch(vrdi_unit,one,three,value)
      if(gstatus.ne.1)then
          call prnt(4,1,gstatus,'xdxswitch error #.')
c         call abend
      endif
      if(value.eq.1)then
         value=0
         side=2
         xcurs=(winss2-ss2+1)/zoom2 + (right/2)
         ycurs=(winsl2-sl2+1)/zoom2
         gstatus=xdcset(vrdi_unit,one,xcurs,ycurs)
         if(gstatus.ne.1)then
             call prnt(4,1,gstatus,'xdcset error #.')
c            call abend
         endif
         call get_square_coord(xbox,ybox,xcurs,ycurs,
     +         				winns/zoom2,winnl/zoom2)
         goto 30
      endif

c read new cursor position.
      gstatus=xdclocation(vrdi_unit,one,xcurs,ycurs)
      if(gstatus.ne.1)then
          call prnt(4,1,gstatus,'xdlocation error #.')
c         call abend
      endif

c restrict correlation area to be within the display area.
      if(side.eq.2)then
         if(xcurs.lt.right/2) xcurs=right/2
         if(xcurs+(winns/zoom2)-1.gt.right)
     +        xcurs=right-(winns/zoom2)+1
      else
         if(xcurs.lt.1) xcurs=1
         if(xcurs+(winns/zoom1)-1.gt.right/2-1)
     +        xcurs=right/2-1-winns/zoom1+1
      endif
      if(ycurs+(winnl/zoom1)-1.gt.bottom)
     +        ycurs=bottom-winnl/zoom1+1
      if(ycurs.lt.1)ycurs=1

      if(xcurs.eq.oldxcurs.and.ycurs.eq.oldycurs) goto 30
      oldxcurs=xcurs
      oldycurs=ycurs
c erase the old box location
      gstatus=xdipolyline(vrdi_unit,graphicsplane,i0,five,xbox,ybox)
      if(gstatus.ne.1)then
          call prnt(4,1,gstatus,'xdipolyline error #.')
c         call abend
      endif
c update picture coordinates
      if(side.eq.2)then
         winsl2=ycurs*zoom2+sl2-1
         winss2=(xcurs-right/2)*zoom2+ss2-1
         call get_square_coord(xbox,ybox,xcurs,ycurs,
     +        				winns/zoom2,winnl/zoom2)
      else
         winsl1=ycurs*zoom1+sl1-1
         winss1=xcurs*zoom1+ss1-1
         call get_square_coord(xbox,ybox,xcurs,ycurs,
     +         				winns/zoom1,winnl/zoom1)
      endif
c draw in box
      gstatus=xdipolyline(vrdi_unit,graphicsplane,i255,five,xbox,ybox)
      if(gstatus.ne.1)then
          call prnt(4,1,gstatus,'xdipolyline error #.')
c         call abend
      endif

c Check if cursor is in one of the control areas below left picture.
      if((ycurs.le.lbox).or.(xcurs.ge.3*nr)
     +    .or.(ycurs.gt.lbox+100) ) goto 30
      if(xcurs.lt.nr)then
         pair=-1
         return
      else if(xcurs.lt.2*nr)then
         pair=1
         return
      else
         side=0
         return
      endif

      goto 30
      end

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
*     MB02FD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2020 NICONET e.V.
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO, ONE
      PARAMETER        ( ZERO = 0.0D0, ONE = 1.0D0 )
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          ITMAX, KMAX, NMAX
      PARAMETER        ( ITMAX = 10, KMAX = 20, NMAX = 20 )
      INTEGER          LDR, LDT, LDWORK
      PARAMETER        ( LDR = NMAX*KMAX, LDT = KMAX,
     $                   LDWORK = ( NMAX + 1 )*KMAX )
*     .. Local Scalars ..
      INTEGER          I, INFO, IT, J, K, LEN, M, N, P, PIT, POS, POSR,
     $                 S1, SCIT
      CHARACTER        TYPET
      DOUBLE PRECISION NNRM
*     .. Local Arrays .. (Dimensioned for TYPET = 'R'.)
      INTEGER          S(ITMAX)
      DOUBLE PRECISION DWORK(LDWORK), R(LDR, NMAX*KMAX),
     $                 T(LDT, NMAX*KMAX), V(NMAX*KMAX), W(NMAX*KMAX),
     $                 Z(NMAX*KMAX)
*     .. External Functions ..
      LOGICAL          LSAME
      DOUBLE PRECISION DNRM2
      EXTERNAL         DNRM2, LSAME
*     .. External Subroutines ..
      EXTERNAL         DAXPY, DCOPY, DGEMV, DLASET, DSCAL, DTRMV, MB02FD
*
*     .. Executable Statements ..
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) N, K, IT
      TYPET = 'R'
      M = N*K
      IF( N.LE.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99993 ) N
      ELSE IF( K.LE.0 .OR. K.GT.KMAX ) THEN
         WRITE ( NOUT, FMT = 99992 ) K
      ELSE IF( IT.LE.0 .OR. IT.GT.ITMAX ) THEN
         WRITE ( NOUT, FMT = 99991 ) IT
      ELSE
         READ ( NIN, FMT = * ) ( S(I), I = 1, IT )
         READ ( NIN, FMT = * ) ( ( T(I,J), J = 1,M ), I = 1,K )
         P   = 0
         POS = 1
         WRITE ( NOUT, FMT = 99997 )
         DO 90  SCIT = 1, IT
            CALL MB02FD( TYPET, K, N, P, S(SCIT), T(1,POS), LDT,
     $                   R(POS,POS), LDR, DWORK, LDWORK, INFO )
            IF ( INFO.NE.0 ) THEN
               WRITE ( NOUT, FMT = 99998 ) INFO
               STOP
            END IF
            S1 = S(SCIT) + P
            IF ( S1.EQ.0 ) THEN
*              Estimate the 2-norm of the Toeplitz matrix with 5 power
*              iterations.
               LEN = N*K
               CALL DLASET( 'All', LEN, 1, ONE, ONE, V, 1 )
               DO 30  PIT = 1, 5
                  DO 10  I = 1, N
                     CALL DGEMV( 'NoTranspose', K, LEN-(I-1)*K, ONE, T,
     $                           LDT, V((I-1)*K+1), 1, ZERO,
     $                           W((I-1)*K+1), 1 )
   10             CONTINUE
                  DO 20 I = 1, N-1
                     CALL DGEMV( 'Transpose', K, (N-I)*K, ONE,
     $                           T(1,K+1), LDT, V((I-1)*K+1), 1,
     $                           ONE, W(I*K+1), 1 )
   20             CONTINUE
                  CALL DCOPY( LEN, W, 1, V, 1 )
                  NNRM = DNRM2( LEN, V, 1 )
                  CALL DSCAL( LEN, ONE/NNRM, V, 1 )
   30          CONTINUE
            ELSE
*              Estimate the 2-norm of the Schur complement with 5 power
*              iterations.
               LEN = ( N - S1 )*K
               CALL DLASET( 'All', LEN, 1, ONE, ONE, V, 1 )
               DO 80  PIT = 1, 5
                  POSR = ( S1 - 1 )*K + 1
                  DO 40  I = 1, N - S1
                     CALL DGEMV( 'NoTranspose', K, LEN-(I-1)*K, ONE,
     $                           T(1,POSR+K), LDT, V((I-1)*K+1), 1,
     $                           ZERO, W((I-1)*K+1), 1 )
   40             CONTINUE
                  DO 50  I = 1, N - S1
                     CALL DTRMV( 'Upper', 'NoTranspose', 'NonUnit', K,
     $                           R(POSR,POSR), LDR, V((I-1)*K+1), 1 )
                     CALL DGEMV( 'NoTranspose', K, LEN-I*K, ONE,
     $                           R(POSR,POSR+K), LDR, V(I*K+1), 1, ONE,
     $                           V((I-1)*K+1), 1 )
   50             CONTINUE
                  CALL DLASET( 'All', LEN, 1, ZERO, ZERO, Z, 1 )
                  DO 60  I = 1, N - S1
                     CALL DGEMV( 'Transpose', K, LEN-I*K, ONE,
     $                           R(POSR,POSR+K), LDR, V((I-1)*K+1), 1,
     $                           ONE, Z(I*K+1), 1 )
                     CALL DTRMV( 'Upper', 'Transpose', 'NonUnit', K,
     $                           R(POSR,POSR), LDR, V((I-1)*K+1), 1 )
                     CALL DAXPY( K, ONE, V((I-1)*K+1), 1, Z((I-1)*K+1),
     $                           1 )
   60             CONTINUE
                  CALL DLASET( 'All', LEN, 1, ZERO, ZERO, V, 1 )
                  DO 70  I = 1, N - S1
                     CALL DGEMV( 'Transpose', K, LEN-(I-1)*K, ONE,
     $                           T(1,POSR+K), LDT, W((I-1)*K+1), 1,
     $                           ONE, V((I-1)*K+1), 1 )
   70             CONTINUE
                  CALL DAXPY( LEN, -ONE, Z, 1, V, 1 )
                  NNRM = DNRM2( LEN, V, 1 )
                  CALL DSCAL( LEN, -ONE/NNRM, V, 1 )
   80          CONTINUE
               POS = ( S1 - 1 )*K + 1
               P   = S1
            END IF
            WRITE ( NOUT, FMT = 99995 ) P*K, NNRM
   90    CONTINUE
         WRITE ( NOUT, FMT = 99996 )
         DO 100  I = 1, P*K
            WRITE ( NOUT, FMT = 99994 ) ( R(I,J), J = 1, M )
  100    CONTINUE
      END IF
      STOP
*
99999 FORMAT (' MB02FD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB02FD = ',I2)
99997 FORMAT ('   Incomplete Cholesky factorization ',
     $         //'   rows    norm(Schur complement)',/)
99996 FORMAT (/' The upper ICC factor of the block Toeplitz matrix is '
     $       )
99995 FORMAT (I4,5X,F8.4)
99994 FORMAT (20(1X,F8.4))
99993 FORMAT (/' N is out of range.',/' N = ',I5)
99992 FORMAT (/' K is out of range.',/' K = ',I5)
99991 FORMAT (/' IT is out of range.',/' IT = ',I5)
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: bnavigator/SLICOT-Reference

max_forks_repo_name: bnavigator/SLICOT-Reference