!!!!!  ==========================================================  !!!!!
!!!!!             'module_radiation_driver' descriptions           !!!!!
!!!!!  ==========================================================  !!!!!
!                                                                      !
!   this is the radiation driver module.  it prepares atmospheric      !
!   profiles and invokes main radiation calculations.                  !
!                                                                      !
!   in module 'module_radiation_driver' there are twe externally       !
!   callable subroutine:                                               !
!                                                                      !
!      'radinit'    -- initialization routine                          !
!         input:                                                       !
!           ( si, NLAY, iflip, idate, jdate, ICTM, ISOL, ICO2,         !
!             IAER, IALB, IEMS, ICWP, NP3D, isubcsw, isubclw,          !
!             iovrsw, iovrlw, me )                                     !
!         output:                                                      !
!           ( none )                                                   !
!                                                                      !
!      'grrad'      -- setup and invoke main radiation calls           !
!         input:                                                       !
!          ( prsi,prsl,prslk,tgrs,qgrs,oz,vvl,slmsk,                   !
!            xlon,xlat,tsfc,snowd,sncovr,snoalb,zorl,hprim,            !
!            alvsf,alnsf,alvwf,alnwf,facsf,facwf,fice,tisfc,           !
!            solcon,coszen,coszdg,k1oz,k2oz,facoz,                     !
!            cv,cvt,cvb,iovrsw,iovrlw,fcice,frain,rrime,flgmin,        !
!            icsdsw,icsdlw, np3d,ntcw,ncld,ntoz, NTRAC,NFXR,           !
!            dtlw,dtsw, lsswr,lslwr,lssav,sashal,norad_precip,         !
!            crick_proof, ccnorm,                                      !
!            IX, IM, LM, iflip, me, lprnt, ipt, kdt,                   !
!         output:                                                      !
!            htrsw,topfsw,sfcfsw,sfalb,                                !
!            htrlw,topflw,sfcflw,tsflw,semis,cldcov,cldsa              !
!         input/output:                                                !
!            fluxr                                                     !
!         optional output:                                             !
!            HTRSWB,HTRLWB)                                            !
!                                                                      !
!                                                                      !
!   external modules referenced:                                       !
!       'module machine'                    in 'machine.f'             !
!       'module funcphys'                   in 'funcphys.f'            !
!       'module physcons'                   in 'physcons.f             !
!                                                                      !
!       'module module_radiation_gases'     in 'radiation_gases.f'     !
!       'module module_radiation_aerosols'  in 'radiation_aerosols.f'  !
!       'module module_radiation_surface'   in 'radiation_surface.f'   !
!       'module module_radiation_clouds'    in 'radiation_clouds.f'    !
!                                                                      !
!       'module module_radsw_cntr_para'     in 'radsw_xxxx_param.f'    !
!       'module module_radsw_parameters'    in 'radsw_xxxx_param.f'    !
!       'module module_radsw_main'          in 'radsw_xxxx_main.f'     !
!                                                                      !
!       'module module_radlw_cntr_para'     in 'radlw_xxxx_param.f'    !
!       'module module_radlw_parameters'    in 'radlw_xxxx_param.f'    !
!       'module module_radlw_main'          in 'radlw_xxxx_main.f'     !
!                                                                      !
!    where xxxx may vary according to different scheme selection       !
!                                                                      !
!                                                                      !
!   program history log:                                               !
!     mm-dd-yy    ncep         - created program grrad                 !
!     08-12-03    yu-tai hou   - re-written for modulized radiations   !
!     11-06-03    yu-tai hou   - modified                              !
!     01-18-05    s. moorthi   - NOAH/ICE model changes added          !
!     05-10-05    yu-tai hou   - modified module structure             !
!     12-xx-05    s. moorthi   - sfc lw flux adj by mean temperature   !
!     02-20-06    yu-tai hou   - add time variation for co2 data, and  !
!                                solar const. add sfc emiss change     !
!     03-21-06    s. Moorthi   - added surface temp over ice           !
!     07-28-06    yu-tai hou   - add stratospheric vocanic aerosols    !
!     03-14-07    yu-tai hou   - add generalized spectral band interp  !
!                                for aerosol optical prop. (sw and lw) !
!     04-10-07    yu-tai hou   - spectral band sw/lw heating rates     !
!     05-04-07    yu-tai hou   - make options for clim based and modis !
!                                based (h. wei and c. marshall) albedo !
!     09-05-08    yu-tai hou   - add the initial date and time 'idate' !
!                    and control param 'ICTM' to the passing param list!
!                    to handel different time/date requirements for    !
!                    external data (co2, aeros, solcon, ...)           !
!     10-10-08    yu-tai hou   - add the ICTM=-2 option for combining  !
!                    initial condition data with seasonal cycle from   !
!                    climatology.                                      !
!     03-12-09    yu-tai hou   - use two time stamps to keep tracking  !
!                    dates for init cond and fcst time. remove volcanic!
!                    aerosols data in climate hindcast (ICTM=-2).      !
!     03-16-09    yu-tai hou   - included sub-column clouds approx.    !
!                    control flags isubcsw/isubclw in initialization   !
!                    subroutine. passed auxiliary cloud control arrays !
!                    icsdsw/icsdlw (if isubcsw/isubclw =2, it will be  !
!                    the user provided permutation seeds) to the sw/lw !
!                    radiation calculation programs. also moved cloud  !
!                    overlapping control flags iovrsw/iovrlw from main !
!                    radiation routines to the initialization routines.!
!     04-02-09    yu-tai hou   - modified surface control flag iems to !
!                    have additional function of if the surface-air    !
!                    interface have the same or different temperature  !
!                    for radiation calculations.                       !
!     04-03-09    yu-tai hou   - modified to add lw surface emissivity !
!                    as output variable. changed the sign of sfcnsw to !
!                    be positive value denote solar flux goes into the !
!                    ground (this is needed to reduce sign confusion   !
!                    in other part of model)                           !
!     04-20-09    carlos perez - prepare driver for nmmb.  added option!
!                    of run the gfs's radiation on nmmb                !
!     09-09-09    fanglin yang (thru s.moorthi) added QME5 QME6 to E-20!
!     01-09-10    sarah lu     - added gocart option, revised grrad for!
!                    gocart coupling. calling argument modifed: ldiag3 !
!                    removed; cldcov/fluxr sequence changed; cldcov is !
!                    changed from accumulative to instant field and    !
!                    from input/output to output field                 !
!     01-24-10    sarah lu     - added aod to fluxr, added prslk and   !
!                    oz to setaer input argument (for gocart coupling),!
!                    added tau_gocart to setaer output argument (for,  !
!                    aerosol diag)                                     !
!     07-08-10    s.moorthi - updated the NEMS version for new physics !
!     07-28-10    yu-tai hou   - changed grrad interface to allow all  !
!                    components of sw/lw toa/sfc instantaneous values  !
!                    being passed to the calling program. moved the    !
!                    computaion of sfc net sw flux (sfcnsw) to the     !
!                    calling program. merged carlos' nmmb modification.!
!     07-30-10    s. moorthi - corrected some errors associated with   !
!                    unit changes                                      !
!     12-02-10    s. moorthi/y. hou - removed the use of aerosol flags !
!                    'iaersw' 'iaerlw' from radiations and replaced    !
!                    them by using the runtime variable iaerflg and    !
!                    laswflg defined in module radiation_aerosols.     !
!                    also replaced param nspc in grrad with the use of !
!                    max_num_gridcomp in module radiation_aerosols.    !
!     01-03-11    y. hou     - added sea/land madk 'slmsk' to the      !
!                    argument list of subrotine setaer call for the    !
!                    newly modified horizontal bi-linear interpolation !
!                    in climatological aerosols schem.                 !
!     09-30-11    H.M. Lin   - rename "grrad" to "grrad_nmmb" for using!
!                    in regional to avoid the confliction with global  !
!                    when compile the nems. Also change names of the   !
!                    related subroutines to "astronomy_nmmb" &         !
!                    "solinit_nmmb"                                    !
!     03-15-12    H.M. Lin   - changed grrad interface to allow 'cldsa'!
!                    (fraction of clouds for low, mid, hi, tot, bl)    !
!                     passed to the calling program.                   !
!     09-25-12    y. hou     - added optional extra top layer 'LTP' on !
!                    top of low ceiling models                         !
!                    (** to cover the gases atop, expecially O3)       !
!                                                                      !
!!!!!  ==========================================================  !!!!!
!!!!!                       end descriptions                       !!!!!
!!!!!  ==========================================================  !!!!!



!=============================================!
      module module_radiation_driver_nmmb     !
!.............................................!
!
      use physparam
      use physcons,                 only : con_eps, con_epsm1,          &
     &                                     con_fvirt, PI=>con_pi
      use funcphys,                 only : fpvs, fpkapx

      use MODULE_CONSTANTS,         only : EPSQ

      use module_radiation_astronomy_nmmb,only : sol_init_nmmb,         &
     &                                     sol_update_nmmb, coszmn_nmmb
      use module_radiation_gases_nmmb,only : NF_VGAS, getgases, getozn, &
     &                                     gas_init, gas_update
      use module_radiation_aerosols_nmmb,only : NF_AESW, NF_AELW,setaer,&
     &                                     aer_init, aer_update
!    &,                                    NSPC1                        ! optn for aod output
      use module_radiation_surface_nmmb, only : NF_ALBD,sfc_init,setalb,&
     &                                     setemis
      use module_radiation_clouds_nmmb,  only : NF_CLDS, cld_init,      &
     &                                     progcld1, progcld2, diagcld1 &
     &                                   , progcld8

      use module_radsw_parameters,  only : topfsw_type, sfcfsw_type,    &
     &                                     profsw_type,cmpfsw_type,NBDSW
      use module_radsw_main_nmmb,   only : rswinit,  swrad

      use module_radlw_parameters,  only : topflw_type, sfcflw_type,    &
     &                                     proflw_type, NBDLW
      use module_radlw_main_nmmb,   only : rlwinit,  lwrad

!
      implicit   none
!
      private

!  ---  version tag and last revision date
      character(40), parameter ::                                       &
     &   VTAGRAD='NCEP-Radiation_driver    v5.2  Jan 2013 '
!    &   VTAGRAD='NCEP-Radiation_driver    v5.1  Nov 2012 '
!    &   VTAGRAD='NCEP-Radiation_driver    v5.0  Aug 2012 '

!  ---  constant values
      real (kind=kind_phys) :: QMIN, QME5, QME6
!     parameter (QMIN=1.0e-10, QME5=1.0e-5,  QME6=1.0e-6,  EPSQ=1.0e-12)
      parameter (QMIN=1.0e-10, QME5=1.0e-7,  QME6=1.0e-7)
!     parameter (QMIN=1.0e-10, QME5=1.0e-20, QME6=1.0e-20, EPSQ=1.0e-12)
      real, parameter :: prsmin = 1.0e-6 ! toa pressure minimum value in mb (hpa)

!  ---  control flags set in subr radinit:
      integer :: itsfc  =0             ! flag for lw sfc air/ground interface temp setting

!  ---  data input control variables set in subr radupdate:
      integer :: month0=0,   iyear0=0,   monthd=0
      logical :: loz1st =.true.       ! first-time clim ozone data read flag

!  ---  optional extra top layer on top of low ceiling models (Y. Hou, 2012-09-25)
!     integer, parameter :: LTP = 0   ! do no add an extra top layer
      integer, parameter :: LTP = 1   ! add an extra top layer
      logical :: lextop = (LTP > 0)
!
      integer, parameter :: CHK= CHNK_RRTM      ! CHNK_RRTM is macro defined in configure.nems

!  ---  publicly accessible module programs:

      public radinit_nmmb, radupdate_nmmb, grrad_nmmb, dayparts


! =================
      contains
! =================


!-----------------------------------
      subroutine radinit_nmmb                                           &
!...................................

!  ---  inputs:
     &     ( si, NLAY, me )
!  ---  outputs:
!          ( none )

! =================   subprogram documentation block   ================ !
!                                                                       !
! subprogram:   radinit     initialization of radiation calculations    !
!                                                                       !
!                                                                       !
! program history log:                                                  !
!   08-14-2003   yu-tai hou   created                                   !
!                                                                       !
! usage:        call radinit                                            !
!                                                                       !
! attributes:                                                           !
!   language:  fortran 90                                               !
!   machine:   ibm sp                                                   !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
! input parameters:                                                     !
!   si               : model vertical sigma interface                   !
!   NLAY             : number of model vertical layers                  !
!   iflip            : control flag for direction of vertical index     !
!                     =0: index from toa to surface                     !
!                     =1: index from surface to toa                     !
!   idate(8)         : ncep absolute date and time of initial condition !
!                      (yr, mon, day, t-zone, hr, min, sec, mil-sec)    !
!   jdate(8)         : ncep absolute date and time at fcst time         !
!                      (yr, mon, day, t-zone, hr, min, sec, mil-sec)    !
!   ICTM             :=yyyy#, external data time/date control flag      !
!                     =   -2: same as 0, but superimpose seasonal cycle !
!                             from climatology data set.                !
!                     =   -1: use user provided external data for the   !
!                             forecast time, no extrapolation.          !
!                     =    0: use data at initial cond time, if not     !
!                             available, use latest, no extrapolation.  !
!                     =    1: use data at the forecast time, if not     !
!                             available, use latest and extrapolation.  !
!                     =yyyy0: use yyyy data for the forecast time,      !
!                             no further data extrapolation.            !
!                     =yyyy1: use yyyy data for the fcst. if needed, do !
!                             extrapolation to match the fcst time.     !
!   ISOL             :=0: use a fixed solar constant value              !
!                     =1: use 11-year cycle solar constant table        !
!   ICO2             :=0: use prescribed global mean co2 (old  oper)    !
!                     =1: use observed co2 annual mean value only       !
!                     =2: use obs co2 monthly data with 2-d variation   !
!   IAER             : 3-digit aerosol flag (for volc, lw, sw)          !
!                     =  0: turn all aeros effects off (sw,lw,volc)     !
!                     =  1: use clim tropspheric aerosol for sw only    !
!                     = 10: use clim tropspheric aerosol for lw only    !
!                     = 11: use clim tropspheric aerosol for both sw/lw !
!                     =100: volc aerosol only for both sw and lw        !
!                     =101: volc and clim trops aerosol for sw only     !
!                     =110: volc and clim trops aerosol for lw only     !
!                     =111: volc and clim trops aerosol for both sw/lw  !
!                     = 2: gocart prognostic, without volc forcing      !
!                     =12: gocart prognostic, with volcanic forcing     !
!   IALB             : control flag for surface albedo schemes          !
!                     =0: climatology, based on surface veg types       !
!                     =1: modis retrieval based surface albedo scheme   !
!                     =2: use externally provided albedoes directly.    !
!   IEMS             : ab 2-digit control flag                          !
!                      a =0 set sfc air/ground t same for lw radiation  !
!                        =1 set sfc air/ground t diff for lw radiation  !
!                      b =0 use fixed sfc emissivity=1.0 (black-body)   !
!                        =1 use varying climtology sfc emiss (veg based)!
!                        =2 future development (not yet)                !
!   ICWP             : control flag for cloud generation schemes        !
!                     =-1: use diagnostic cloud scheme (GFDL type       !
!                     =0 : use diagnostic cloud scheme                  !
!                     =1 : use prognostic cloud scheme (default)        !
!   NP3D             :=3: ferrier's microphysics cloud scheme           !
!                     =4: zhao/carr/sundqvist microphysics cloud        !
!                     =5: nmmb ferrier+bmj microphysics scheme          !
!                     =8: Thompson microphysics scheme                  !
!   isubcsw/isubclw  : sub-column cloud approx control flag (sw/lw rad) !
!                     =0: with out sub-column cloud approximation       !
!                     =1: mcica sub-col approx. prescribed random seed  !
!                     =2: mcica sub-col approx. provided random seed    !
!   iovrsw/iovrlw    : control flag for cloud overlap (sw/lw rad)       !
!                     =0: random overlapping clouds                     !
!                     =1: max/ran overlapping clouds                    !
!   me               : print control flag                               !
!                                                                       !
!  outputs: (none)                                                      !
!                                                                       !
!  usage:       call radinit                                            !
!                                                                       !
!  subroutines called:    cldinit, aerinit, rlwinit, rswinit, gasinit   !
!                                                                       !
!  ===================================================================  !
!
      implicit none

!  ---  inputs:
      integer, intent(in) :: NLAY, me

      real (kind=kind_phys), intent(in) :: si(:)

!  ---  outputs: (none, to module variables)

!  ---  locals:
      integer :: RICWP, icldflg_org

!
!===> ...  begin here
!
!  --- ...  set up control variables
      itsfc  = iemsflg / 10             ! sfc air/ground temp control
      loz1st = (ioznflg == 0)           ! first-time clim ozone data read flag

      if (me == 0) then
!       print *,' NEW RADIATION PROGRAM STRUCTURES -- SEP 01 2004'
        print *,' NEW RADIATION PROGRAM STRUCTURES BECAME OPER. ',      &
     &          '  May 01 2007'
        print *, VTAGRAD                !print out version tag
        print *,' - Selected Control Flag settings: ICTMflg=',ictmflg,  &
     &    ' ISOLar =',isolar, ' ICO2flg=',ico2flg,' IAERflg=',iaerflg,  &
     &    ' IALBflg=',ialbflg,' IEMSflg=',iemsflg,' ICLDflg=',icldflg,  &
     &    ' ICMPHYS=',icmphys,' IOZNflg=',ioznflg
        print *,' IVFLIP=',ivflip,' IOVRSW=',iovrsw,' IOVRLW=',iovrlw,  &
     &    ' ISUBCSW=',isubcsw,' ISUBCLW=',isubclw
        print *,' LSASHAL=',lsashal,' LCRICK=',lcrick,' LCNORM=',lcnorm,&
     &    ' LNOPREC=',lnoprec
        print *,' LTP =',LTP,', add extra top layer =',lextop

        if ( ictmflg==0 .or. ictmflg==-2 ) then
          print *,'   Data usage is limited by initial condition!'
          print *,'   No volcanic aerosols'
        endif

        if ( isubclw == 0 ) then
          print *,' - ISUBCLW=',isubclw,' No McICA, use grid ',         &
     &            'averaged cloud in LW radiation'
        elseif ( isubclw == 1 ) then
          print *,' - ISUBCLW=',isubclw,' Use McICA with fixed',        &
     &            'permutation seeds for LW random number generator'
        elseif ( isubclw == 2 ) then
          print *,' - ISUBCLW=',isubclw,' Use McICA with random ',      &
     &            'permutation seeds for LW random number generator'
        else
          print *,' - ERROR!!! ISUBCLW=',isubclw,' is not a ',          &
     &            'valid option '
          stop
        endif

        if ( isubcsw == 0 ) then
          print *,' - ISUBCSW=',isubcsw,' No McICA, use grid ',         &
     &            'averaged cloud in SW radiation'
        elseif ( isubcsw == 1 ) then
          print *,' - ISUBCSW=',isubcsw,' Use McICA with fixed',        &
     &            'permutation seeds for SW random number generator'
        elseif ( isubcsw == 2 ) then
          print *,' - ISUBCSW=',isubcsw,' Use McICA with random ',      &
     &            'permutation seeds for SW random number generator'
        else
          print *,' - ERROR!!! ISUBCSW=',isubcsw,' is not a ',          &
     &            'valid option '
          stop
        endif

        if ( isubcsw /= isubclw ) then
          print *,' - *** Notice *** ISUBCSW /= ISUBCLW !!!',           &
     &            isubcsw, isubclw
        endif
      endif

!  --- ...  call astronomy initialization routine

      call sol_init_nmmb (  me )

!  --- ...  call aerosols initialization routine

      call aer_init ( NLAY, me )

!  --- ...  call co2 and other gases initialization routine

      call gas_init ( me )

!  --- ...  call surface initialization routine

      call sfc_init ( me )

!  --- ...  call cloud initialization routine

      call cld_init ( si, NLAY, me)


!==========================================================
! The following use for NP3D=5 GFDL cloud in radiation
!==========================================================

      if ( icmphys==5 ) then
         ricwp = 0
         icldflg_org = icldflg
         icldflg = ricwp
      endif

!  --- ...  call lw radiation initialization routine

      call rlwinit ( me )

!  --- ...  call sw radiation initialization routine

      call rswinit ( me )
!

!==========================================================
! The following use for NP3D=5 GFDL cloud in radiation
! return icldflg to initial setting
!==========================================================

      if ( icmphys==5 ) then
         icldflg = icldflg_org
      endif

      
      return
!...................................
      end subroutine radinit_nmmb
!-----------------------------------


!-----------------------------------
      subroutine radupdate_nmmb                                              &
!...................................

!  ---  inputs:
     &     ( idate, jdate, deltsw, deltim, lsswr, me,                   &
!  ---  outputs:
     &       slag,sdec,cdec,solcon                                      &
     &     )

! =================   subprogram documentation block   ================ !
!                                                                       !
! subprogram:   radupdate   calls many update subroutines to check and  !
!   update radiation required but time varying data sets and module     !
!   variables.                                                          !
!                                                                       !
! usage:        call radupdate                                          !
!                                                                       !
! attributes:                                                           !
!   language:  fortran 90                                               !
!   machine:   ibm sp                                                   !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
! input parameters:                                                     !
!   idate(8)       : ncep absolute date and time of initial condition   !
!                    (yr, mon, day, t-zone, hr, min, sec, mil-sec)      !
!   jdate(8)       : ncep absolute date and time at fcst time           !
!                    (yr, mon, day, t-zone, hr, min, sec, mil-sec)      !
!   deltsw         : sw radiation calling frequency in seconds          !
!   deltim         : model timestep in seconds                          !
!   lsswr          : logical flags for sw radiation calculations        !
!   me             : print control flag                                 !
!                                                                       !
!  outputs:                                                             !
!   slag           : equation of time in radians                        !
!   sdec, cdec     : sin and cos of the solar declination angle         !
!   solcon         : sun-earth distance adjusted solar constant (w/m2)  !
!                                                                       !
!  module variables:                                                    !
!   isolar   : solar constant cntrl  (see ISOL description)             !
!   ictmfl   : flag for initial condition gh-gas data source            !
!   loz1st   : first-time clim ozone data read flag                     !
!                                                                       !
!  subroutines called: sol_update, aer_update, gas_update               !
!                                                                       !
!  ===================================================================  !
!
      implicit none

!  ---  inputs:
      integer, intent(in) :: idate(:), jdate(:), me
      logical, intent(in) :: lsswr

      real (kind=kind_phys), intent(in) :: deltsw, deltim

!  ---  outputs:
      real (kind=kind_phys), intent(out) :: slag, sdec, cdec, solcon

!  ---  locals:
      integer :: iyear, imon, iday, ihour
      integer :: kyear, kmon, kday, khour

      logical :: lmon_chg       ! month change flag
      logical :: lco2_chg       ! cntrl flag for updating co2 data
      logical :: lsol_chg       ! cntrl flag for updating solar constant
!
!===> ...  begin here
!

!  --- ...  time stamp at fcst time

      iyear = jdate(1)
      imon  = jdate(2)
      iday  = jdate(3)
      ihour = jdate(5)

!  --- ...  set up time stamp used for green house gases (** currently co2 only)

      if ( ictmflg==0 .or. ictmflg==-2 ) then  ! get external data at initial condition time
        kyear = idate(1)
        kmon  = idate(2)
        kday  = idate(3)
        khour = idate(5)
      else                           ! get external data at fcst or specified time
        kyear = iyear
        kmon  = imon
        kday  = iday
        khour = ihour
      endif   ! end if_ictmflg_block

      if ( month0 /= imon ) then
        lmon_chg = .true.
        month0 = imon
      else
        lmon_chg = .false.
      endif

!  --- ...  call astronomy update routine, yearly update, no time interpolation

      if (lsswr) then

        if  ( isolar == 0 .or. isolar == 10 ) then
          lsol_chg = .false.
        elseif ( iyear0/=iyear ) then
          lsol_chg = .true.
        else
          lsol_chg = ( isolar==4 .and. lmon_chg )
        endif
        iyear0 = iyear

        call sol_update_nmmb                                            &
!  ---  input:
     &     ( jdate,kyear,deltsw,deltim,lsol_chg, me,                    &
!  ---  outputs:
     &       slag,sdec,cdec,solcon                                      &
     &     )

      endif  ! end_if_lsswr_block

!  --- ...  call aerosols update routine, monthly update, no time interpolation

      if ( lmon_chg ) then
        call aer_update ( iyear, imon, me )
      endif

!  --- ...  call co2 and other gases update routine

      if ( monthd /= kmon ) then
        monthd = kmon
        lco2_chg = .true.
      else
        lco2_chg = .false.
      endif

      call gas_update ( kyear,kmon,kday,khour,loz1st,lco2_chg, me )

      if ( loz1st ) loz1st = .false.

!  --- ...  call surface update routine (currently not needed)
!     call sfc_update ( iyear, imon, me )

!  --- ...  call clouds update routine (currently not needed)
!     call cld_update ( iyear, imon, me )
!
      return
!...................................
      end subroutine radupdate_nmmb
!-----------------------------------

!-----------------------------------
      subroutine dayparts                                               &
!...................................
!  ---  inputs:
     &     (         XLON,SINLAT,COSLAT,SOLHR,MYPE,dtswav,nrads,        &
     &               veclen ,                                           &
!  ---  outputs:
     &               dp_start,dp_len,dp_day,ndayparts                   &
     &     )

      implicit none
      integer, intent(in) ::  veclen,mype,nrads
      real (kind=kind_phys), intent(in) :: dtswav, solhr
      real (kind=kind_phys), dimension(veclen), intent(in) ::           &
     &       sinlat, coslat, xlon
      integer, intent(out), dimension(veclen) :: dp_start,dp_len
      logical, intent(out), dimension(veclen) :: dp_day
      integer, intent(out) :: ndayparts
! local
      logical isday
      integer i
      real (kind=kind_phys), dimension(veclen) :: coszen_loc, coszdg

      call coszmn_nmmb                                                  &
!  ---  inputs:
     &     ( xlon,sinlat,coslat,solhr,veclen,mype,                      &
     &       dtswav,nrads ,                                             &  ! Extra input
!  ---  outputs:
     &       coszen_loc, coszdg                                         &  !jm
     &      )

      if ( veclen < 1 ) then
        ndayparts = 0
      else
        ndayparts = 1
        dp_day(ndayparts) = (coszen_loc(1) >= 0.0001)
        dp_start(ndayparts) = 1
        dp_len(ndayparts) = 1
        do i = 2, veclen
          isday = (coszen_loc(i) >= 0.0001)
          if ( dp_day(ndayparts) .NEQV. isday ) THEN
              ndayparts=ndayparts+1
              dp_start(ndayparts) = i
              dp_day(ndayparts) = isday
              dp_len(ndayparts) = 1
          else
              dp_len(ndayparts) = dp_len(ndayparts) + 1
          endif
        enddo
      endif
!...................................
      end subroutine dayparts
!-----------------------------------

!-----------------------------------
      subroutine grrad_nmmb                                             &
!...................................

!  ---  inputs:
     &     ( prsi,prsl,prslk,tgrs,qgrs,tracer,vvl,slmsk,                &
     &       xlon,xlat,tsfc,snowd,sncovr,snoalb,zorl,hprim,             &
 !    &       alvsf,alnsf,alvwf,alnwf,facsf,facwf,                      &  ! processed inside
     &       SALBEDO,SM,                                                &  ! input for albedo cal.
     &       fice,tisfc,                                                &
     &       sinlat,coslat,solhr,jdate,solcon,                          &
     &       dtswav,nrads,                                              &  ! extra input
     &       cv,cvt,cvb,fcice,frain,rrime,flgmin,                       &
     &       icsdsw,icsdlw, ntcw,ncld,ntoz, NTRAC,NFXR,                 &
     &       CPATHFAC4LW,                                               &  ! enhance factor of cloud depth for LW
     &       dtlw,dtsw, lsswr,lslwr,lssav,                              &
     &       its,jts, IX, IM, LM, isday, me, lprnt, ipt, kdt,           &
!  ---  additional inputs:                                                 ! GFDL type
     &       tauclouds,cldf,                                            &  ! GFDL type
     &       cld_fraction,                                              &  ! Thompson cloud fraction
!  ---  outputs:
     &       htrsw,topfsw,sfcfsw,sfalb,coszen,coszdg,                   &
     &       htrlw,topflw,sfcflw,tsflw,semis,cldcov,cldsa,              &
!  ---  input/output:
     &       fluxr                                                      &
!! ---  optional outputs:
     &,      HTRSWB,HTRLWB                                              &
     &     )

! =================   subprogram documentation block   ================ !
!                                                                       !
!    this program is the driver of radiation calculation subroutines. * !
!    It sets up profile variables for radiation input, including      * !
!    clouds, surface albedos, atmospheric aerosols, ozone, etc.       * !
!                                                                     * !
!                                                                     * !
!    usage:        call grrad                                         * !
!                                                                     * !
!    subprograms called:                                              * !
!                  setalb, setemis, setaer, getozn, getgases,         * !
!                  progcld1, progcld2, progcld8, diagcds,             * !
!                  swrad, lwrad, fpvs                                 * !
!                                                                     * !
!    attributes:                                                      * !
!      language:   fortran 90                                         * !
!      machine:    ibm-sp, sgi                                        * !
!                                                                     * !
!                                                                     * !
!  ====================  defination of variables  ====================  !
!                                                                       !
!    input variables:                                                   !
!      prsi  (IX,LM+1) : model level pressure in Pa                     !
!      prsl  (IX,LM)   : model layer mean pressure in Pa                !
!      prslk (IX,LM)   : Exner function                                 !
!      tgrs  (IX,LM)   : model layer mean temperature in k              !
!      qgrs  (IX,LM)   : layer specific humidity in gm/gm               !
!      oz  (IX,LM,NTRAC):layer ozone mass mixing ratio                  !
!      vvl   (IX,LM)   : layer mean vertical velocity in Pa/sec         !
!      slmsk (IM)      : sea/land mask array (sea:0,land:1,sea-ice:2)   !
!      xlon,xlat (IM)  : grid longitude/latitude in radians             !
!      tsfc  (IM)      : surface temperature in k                       !
!      snowd (IM)      : snow depth water equivalent in mm              !
!      sncovr(IM)      : snow cover in fraction                         !
!      snoalb(IM)      : maximum snow albedo in fraction                !
!      zorl  (IM)      : surface roughness in cm                        !
!      hprim (IM)      : topographic standard deviation in m            !
!      alvsf (IM)      : mean vis albedo with strong cosz dependency    !
!      alnsf (IM)      : mean nir albedo with strong cosz dependency    !
!      alvwf (IM)      : mean vis albedo with weak cosz dependency      !
!      alnwf (IM)      : mean nir albedo with weak cosz dependency      !
!      facsf (IM)      : fractional coverage with strong cosz dependen  !
!      facwf (IM)      : fractional coverage with weak cosz dependency  !
!      fice  (IM)      : ice fraction over open water grid              !
!      tisfc (IM)      : surface temperature over ice fraction          !
!      solcon          : solar constant (sun-earth distant adjusted)    !
!      coszen(IM)      : mean cos of zenith angle over rad call period  !
!      coszdg(IM)      : daytime mean cosz over rad call period         !
!      k1oz,k2oz,facoz : parameters for climatological ozone            !
!      cv    (IM)      : fraction of convective cloud                   !
!      cvt, cvb (IM)   : convective cloud top/bottom pressure in cb     !
!      iovrsw/iovrlw   : control flag for cloud overlap (sw/lw rad)     !
!                        =0 random overlapping clouds                   !
!                        =1 max/ran overlapping clouds                  !
!      fcice           : fraction of cloud ice  (in ferrier scheme)     !
!      frain           : fraction of rain water (in ferrier scheme)     !
!      rrime           : mass ratio of total to unrimed ice ( >= 1 )    !
!      flgmin          : minimim large ice fraction                     !
!      icsdsw/icsdlw   : auxiliary cloud control arrays passed to main  !
!           (IM)         radiations. if isubcsw/isubclw (input to init) !
!                        are set to 2, the arrays contains provided     !
!                        random seeds for sub-column clouds generators  !
!      np3d            : =3 brad ferrier microphysics scheme            !
!                        =4 zhao/carr/sundqvist microphysics scheme     !
!                        =5 external microphysics scheme provided bulk/ !
!                           grey quantities of cloud fields. (nmmb vars)!
!                        =8 Thompson microphysics scheme                !
!      ntcw            : =0 no cloud condensate calculated              !
!                        >0 array index location for cloud condensate   !
!      ncld            : only used when ntcw .gt. 0                     !
!      ntoz            : =0 climatological ozone profile                !
!                        >0 interactive ozone profile                   !
!      NTRAC           : dimension veriable for array oz                !
!      NFXR            : second dimension of input/output array fluxr   !
!      dtlw, dtsw      : time duration for lw/sw radiation call in sec  !
!      lsswr, lslwr    : logical flags for sw/lw radiation calls        !
!      lssav           : logical flag for store 3-d cloud field         !
!      sashal          : logical flag for Jongil's shallow convection   !
!      norad_precip    : logical flag for not using precip in radiation !
!      crick_proof     : logical flag for eliminating CRICK             !
!      ccnorm          : logical flag for incloud condensate mixing ratio!
!      IX,IM           : horizontal dimention and num of used points    !
!      LM              : vertical layer dimension                       !
!      iflip           : control flag for in/out vertical indexing      !
!                        =0 index from toa to surface                   !
!                        =1 index from surface to toa                   !
!      me              : control flag for parallel process              !
!      lprnt           : control flag for diagnostic print out          !
!      ipt             : index for diagnostic printout point            !
!      kdt             : time-step number                               !
!                                                                       !
!    output variables:                                                  !
!      htrsw (IX,LM)   : total sky sw heating rate in k/sec             !
!      topfsw(IM)      : sw radiation fluxes at toa, components:        !
!                      (check module_radsw_parameters for definition)   !
!       %upfxc           - total sky upward sw flux at toa (w/m**2)     !
!       %dnflx           - total sky downward sw flux at toa (w/m**2)   !
!       %upfx0           - clear sky upward sw flux at toa (w/m**2)     !
!      sfcfsw(IM)      : sw radiation fluxes at sfc, components:        !
!                      (check module_radsw_parameters for definition)   !
!       %upfxc           - total sky upward sw flux at sfc (w/m**2)     !
!       %dnfxc           - total sky downward sw flux at sfc (w/m**2)   !
!       %upfx0           - clear sky upward sw flux at sfc (w/m**2)     !
!       %dnfx0           - clear sky downward sw flux at sfc (w/m**2)   !
!      sfalb (IM)      : mean surface diffused sw albedo                !
!      cldcov(IX,LM)   : 3-d cloud fraction                             !
!      cldsa(IX,5)     : fraction of clouds for low, mid, hi, tot, bl   !
!      htrlw (IX,LM)   : total sky lw heating rate in k/sec             !
!      topflw(IM)      : lw radiation fluxes at top, component:         !
!                        (check module_radlw_paramters for definition)  !
!       %upfxc           - total sky upward lw flux at toa (w/m**2)     !
!       %upfx0           - clear sky upward lw flux at toa (w/m**2)     !
!      sfcflw(IM)      : lw radiation fluxes at sfc, component:         !
!                        (check module_radlw_paramters for definition)  !
!       %upfxc           - total sky upward lw flux at sfc (w/m**2)     !
!       %upfx0           - clear sky upward lw flux at sfc (w/m**2)     !
!       %dnfxc           - total sky downward lw flux at sfc (w/m**2)   !
!       %dnfx0           - clear sky downward lw flux at sfc (w/m**2)   !
!      semis (IM)      : surface lw emissivity in fraction              !
!      tsflw (IM)      : surface air temp during lw calculation in k    !
!                                                                       !
!    input and output variables:                                        !
!      fluxr (IX,NFXR) : to save time accumulated 2-d fields defined as:!
!                 1      - toa total sky upwd lw radiation flux         !
!                 2      - toa total sky upwd sw radiation flux         !
!                 3      - sfc total sky upwd sw radiation flux         !
!                 4      - sfc total sky dnwd sw radiation flux         !
!                 5      - high domain cloud fraction                   !
!                 6      - mid  domain cloud fraction                   !
!                 7      - low  domain cloud fraction                   !
!                 8      - high domain mean cloud top pressure          !
!                 9      - mid  domain mean cloud top pressure          !
!                10      - low  domain mean cloud top pressure          !
!                11      - high domain mean cloud base pressure         !
!                12      - mid  domain mean cloud base pressure         !
!                13      - low  domain mean cloud base pressure         !
!                14      - high domain mean cloud top temperature       !
!                15      - mid  domain mean cloud top temperature       !
!                16      - low  domain mean cloud top temperature       !
!                17      - total cloud fraction                         !
!                18      - boundary layer domain cloud fraction         !
!                19      - sfc total sky dnwd lw radiation flux         !
!                20      - sfc total sky upwd lw radiation flux         !
!                21      - sfc total sky dnwd sw uv-b radiation flux    !
!                22      - sfc clear sky dnwd sw uv-b radiation flux    !
!                23      - toa incoming solar radiation flux            !
!                24      - sfc vis beam dnwd sw radiation flux          !
!                25      - sfc vis diff dnwd sw radiation flux          !
!                26      - sfc nir beam dnwd sw radiation flux          !
!                27      - sfc nir diff dnwd sw radiation flux          !
!                28      - toa clear sky upwd lw radiation flux         !
!                29      - toa clear sky upwd sw radiation flux         !
!                30      - sfc clear sky dnwd lw radiation flux         !
!                31      - sfc clear sky upwd sw radiation flux         !
!                32      - sfc clear sky dnwd sw radiation flux         !
!                33      - sfc clear sky upwd lw radiation flux         !
!optional        34      - aeros opt depth at 550nm (all components)    !
!               ....     - optional for test and future use             !
!                                                                       !
!    optional output variables:                                         !
!      htrswb(IX,LM,NBDSW) : spectral band total sky sw heating rate    !
!      htrlwb(IX,LM,NBDLW) : spectral band total sky lw heating rate    !
!                                                                       !
!                                                                       !
!    definitions of internal variable arrays:                           !
!                                                                       !
!     1. fixed gases:         (defined in 'module_radiation_gases')     !
!          gasvmr(:,:,1)  -  co2 volume mixing ratio                    !
!          gasvmr(:,:,2)  -  n2o volume mixing ratio                    !
!          gasvmr(:,:,3)  -  ch4 volume mixing ratio                    !
!          gasvmr(:,:,4)  -  o2  volume mixing ratio                    !
!          gasvmr(:,:,5)  -  co  volume mixing ratio                    !
!          gasvmr(:,:,6)  -  cf11 volume mixing ratio                   !
!          gasvmr(:,:,7)  -  cf12 volume mixing ratio                   !
!          gasvmr(:,:,8)  -  cf22 volume mixing ratio                   !
!          gasvmr(:,:,9)  -  ccl4 volume mixing ratio                   !
!                                                                       !
!     2. cloud profiles:      (defined in 'module_radiation_clouds')    !
!                ---  for  prognostic cloud  ---                        !
!          clouds(:,:,1)  -  layer total cloud fraction                 !
!          clouds(:,:,2)  -  layer cloud liq water path                 !
!          clouds(:,:,3)  -  mean effective radius for liquid cloud     !
!          clouds(:,:,4)  -  layer cloud ice water path                 !
!          clouds(:,:,5)  -  mean effective radius for ice cloud        !
!          clouds(:,:,6)  -  layer rain drop water path                 !
!          clouds(:,:,7)  -  mean effective radius for rain drop        !
!          clouds(:,:,8)  -  layer snow flake water path                !
!          clouds(:,:,9)  -  mean effective radius for snow flake       !
!                ---  for  diagnostic cloud  ---                        !
!          clouds(:,:,1)  -  layer total cloud fraction                 !
!          clouds(:,:,2)  -  layer cloud optical depth                  !
!          clouds(:,:,3)  -  layer cloud single scattering albedo       !
!          clouds(:,:,4)  -  layer cloud asymmetry factor               !
!                                                                       !
!     3. surface albedo:      (defined in 'module_radiation_surface')   !
!          sfcalb( :,1 )  -  near ir direct beam albedo                 !
!          sfcalb( :,2 )  -  near ir diffused albedo                    !
!          sfcalb( :,3 )  -  uv+vis direct beam albedo                  !
!          sfcalb( :,4 )  -  uv+vis diffused albedo                     !
!                                                                       !
!     4. sw aerosol profiles: (defined in 'module_radiation_aerosols')  !
!          faersw(:,:,:,1)-  sw aerosols optical depth                  !
!          faersw(:,:,:,2)-  sw aerosols single scattering albedo       !
!          faersw(:,:,:,3)-  sw aerosols asymmetry parameter            !
!                                                                       !
!     5. lw aerosol profiles: (defined in 'module_radiation_aerosols')  !
!          faerlw(:,:,:,1)-  lw aerosols optical depth                  !
!          faerlw(:,:,:,2)-  lw aerosols single scattering albedo       !
!          faerlw(:,:,:,3)-  lw aerosols asymmetry parameter            !
!                                                                       !
!     6. sw fluxes at toa:    (defined in 'module_radsw_main')          !
!        (topfsw_type -- derived data type for toa rad fluxes)          !
!          topfsw(:)%upfxc  -  total sky upward flux at toa             !
!          topfsw(:)%dnfxc  -  total sky downward flux at toa           !
!          topfsw(:)%upfx0  -  clear sky upward flux at toa             !
!                                                                       !
!     7. lw fluxes at toa:    (defined in 'module_radlw_main')          !
!        (topflw_type -- derived data type for toa rad fluxes)          !
!          topflw(:)%upfxc  -  total sky upward flux at toa             !
!          topflw(:)%upfx0  -  clear sky upward flux at toa             !
!                                                                       !
!     8. sw fluxes at sfc:    (defined in 'module_radsw_main')          !
!        (sfcfsw_type -- derived data type for sfc rad fluxes)          !
!          sfcfsw(:)%upfxc  -  total sky upward flux at sfc             !
!          sfcfsw(:)%dnfxc  -  total sky downward flux at sfc           !
!          sfcfsw(:)%upfx0  -  clear sky upward flux at sfc             !
!          sfcfsw(:)%dnfx0  -  clear sky downward flux at sfc           !
!                                                                       !
!     9. lw fluxes at sfc:    (defined in 'module_radlw_main')          !
!        (sfcflw_type -- derived data type for sfc rad fluxes)          !
!          sfcflw(:)%upfxc  -  total sky upward flux at sfc             !
!          sfcflw(:)%dnfxc  -  total sky downward flux at sfc           !
!          sfcflw(:)%dnfx0  -  clear sky downward flux at sfc           !
!                                                                       !
!! optional radiation outputs:                                          !
!!   10. sw flux profiles:    (defined in 'module_radsw_main')          !
!!       (profsw_type -- derived data type for rad vertical profiles)   !
!!         fswprf(:,:)%upfxc - total sky upward flux                    !
!!         fswprf(:,:)%dnfxc - total sky downward flux                  !
!!         fswprf(:,:)%upfx0 - clear sky upward flux                    !
!!         fswprf(:,:)%dnfx0 - clear sky downward flux                  !
!!                                                                      !
!!   11. lw flux profiles:    (defined in 'module_radlw_main')          !
!!       (proflw_type -- derived data type for rad vertical profiles)   !
!!         flwprf(:,:)%upfxc - total sky upward flux                    !
!!         flwprf(:,:)%dnfxc - total sky downward flux                  !
!!         flwprf(:,:)%upfx0 - clear sky upward flux                    !
!!         flwprf(:,:)%dnfx0 - clear sky downward flux                  !
!!                                                                      !
!!   12. sw sfc components:   (defined in 'module_radsw_main')          !
!!       (cmpfsw_type -- derived data type for component sfc fluxes)    !
!!         scmpsw(:)%uvbfc  -  total sky downward uv-b flux at sfc      !
!!         scmpsw(:)%uvbf0  -  clear sky downward uv-b flux at sfc      !
!!         scmpsw(:)%nirbm  -  total sky sfc downward nir direct flux   !
!!         scmpsw(:)%nirdf  -  total sky sfc downward nir diffused flux !
!!         scmpsw(:)%visbm  -  total sky sfc downward uv+vis direct flx !
!!         scmpsw(:)%visdf  -  total sky sfc downward uv+vis diff flux  !
!                                                                       !
!  ======================  end of definations  =======================  !
!
      implicit none
 
!  ---  constant parameter
!       The following need to be commented out after the newer version 
!       (Hsin-mu Lin, 2011-04-25)

 !     integer, parameter :: NSPC = 6


!  ---  inputs: (for rank>1 arrays, horizontal dimensioned by IX)
      integer,  intent(in) :: IX,IM, LM, NTRAC, NFXR, me,               &
     &                        ntoz, ntcw, ncld, ipt, kdt,               &
     &                        nrads, cld_fraction
      integer,  intent(in) :: its,jts  
      logical,  intent(in) :: isday  ! true = day, false = night
      integer,  intent(in) :: icsdsw(IX), icsdlw(IX), jdate(8)

      logical,  intent(in) :: lsswr, lslwr, lssav, lprnt

      real (kind=kind_phys), dimension(IX,LM+1), intent(in) :: prsi

      real (kind=kind_phys), dimension(IX,LM), intent(in) :: prsl, &
     &       prslk, tgrs, qgrs, vvl, fcice, frain, rrime
      real (kind=kind_phys), dimension(IX), intent(in) :: flgmin

      real (kind=kind_phys), dimension(IX),      intent(in) ::  slmsk,  &
     &       xlon, xlat, tsfc, snowd, zorl, hprim,                      &
 !    & alvsf, alnsf, alvwf, alnwf, facsf, facwf,                       &
     &       cv, cvt, cvb, fice, tisfc,                                 &
     &       sncovr, snoalb, sinlat, coslat

      real (kind=kind_phys), intent(in) :: solcon, dtlw, dtsw, solhr,   &
     & tracer(IX,LM,NTRAC) &
     & , dtswav, SALBEDO(IX), SM(IX) ! extra input

      real (kind=kind_phys), dimension(IX,LM), intent(in) :: tauclouds, & ! GFDL type
     &       cldf

! --- outputs: (horizontal dimensioned by IX)
      real (kind=kind_phys), dimension(IX,LM),intent(out):: htrsw,htrlw,&
     &       cldcov
      real (kind=kind_phys), dimension(IX,5),intent(out):: cldsa    

      real (kind=kind_phys), dimension(IX),   intent(out):: tsflw,      &
     &       sfalb, semis, coszen, coszdg

      real, intent(in) :: CPATHFAC4LW

      type (topfsw_type), dimension(IX), intent(out) :: topfsw
      type (sfcfsw_type), dimension(IX), intent(out) :: sfcfsw
      type (topflw_type), dimension(IX), intent(out) :: topflw
      type (sfcflw_type), dimension(IX), intent(out) :: sfcflw

!  ---  variables are for both input and output:
      real (kind=kind_phys), intent(inout) :: fluxr(IX,NFXR)

!! ---  optional outputs:
      real (kind=kind_phys), dimension(IX,LM,NBDSW), optional, &
     &                       intent(out) :: htrswb
      real (kind=kind_phys), dimension(IX,LM,NBDLW), optional, &
     &                       intent(out) :: htrlwb

!  ---  local variables: (horizontal dimensioned by CHK)
      type (topfsw_type), dimension(CHK)     :: topfsw_loc
      type (sfcfsw_type), dimension(CHK)     :: sfcfsw_loc
      real (kind=kind_phys), dimension(CHK)  :: coszen_loc    !jm
      integer, dimension(CHK)                :: icsdsw_loc    !jm

      real (kind=kind_phys), dimension(CHK,LM+1+LTP)::plvl, tlvl

      real (kind=kind_phys), dimension(CHK,LM+LTP)  ::          &
     &       plyr, tlyr, qlyr,                                          &
     &       olyr, rhly, qstl, vvel, clw, prslk1, tem2da, tem2db, tvly, &
     &       cldfk1

      real (kind=kind_phys), dimension(CHK,LM+LTP)  ::          &
     &       qc2d, qi2d, qs2d, ni2d

      real (kind=kind_phys), dimension(CHK) :: tsfa, cvt1, cvb1, tem1d,  &
     &       sfcemis, tsfg, tskn

      real (kind=kind_phys), dimension(CHK) ::                           &
     &                       alvsf, alnsf, alvwf, alnwf, facsf, facwf,  &
     &                       CZMEAN

      real (kind=kind_phys), dimension(CHK,LM+LTP,NF_CLDS) :: clouds
      real (kind=kind_phys), dimension(CHK,LM+LTP,NF_CLDS) :: dummys
      real (kind=kind_phys), dimension(CHK,LM+LTP,NF_VGAS) :: gasvmr
      real (kind=kind_phys), dimension(CHK,       NF_ALBD) :: sfcalb
!     real (kind=kind_phys), dimension(CHK,       NSPC1)   :: aerodp      ! optn for aod output
      real (kind=kind_phys), dimension(CHK,LM+LTP,NTRAC)   :: tracer1

      real (kind=kind_phys), dimension(CHK,LM+LTP,NBDSW,NF_AESW):: faersw
      real (kind=kind_phys), dimension(CHK,LM+LTP,NBDLW,NF_AELW):: faerlw

 !     real (kind=kind_phys), dimension(CHK,LM+LTP,NSPC-1) :: tau_gocart

      real (kind=kind_phys), dimension(CHK,LM+LTP) :: htswc
      real (kind=kind_phys), dimension(CHK,LM+LTP) :: htlwc

      real (kind=kind_phys), dimension(CHK,LM+LTP) :: gcice, grain, grime

!! ---  may be used for optional sw/lw outputs:
!!      take out "!!" as needed
!!    real (kind=kind_phys), dimension(CHK,LM+LTP)   :: htsw0
!!    type (profsw_type),    dimension(CHK,LM+1+LTP) :: fswprf
      type (cmpfsw_type),    dimension(CHK)          :: scmpsw
      real (kind=kind_phys), dimension(CHK,LM+LTP,NBDSW) :: htswb

!!    real (kind=kind_phys), dimension(CHK,LM+LTP)   :: htlw0
!!    type (proflw_type),    dimension(CHK,LM+1+LTP) :: flwprf
      real (kind=kind_phys), dimension(CHK,LM+LTP,NBDLW) :: htlwb

!jm      real (kind=kind_phys)::raddt,es(CHK),qs, delt, tem0d, cldsa(CHK,5)&
!jm     & ,SFCALBEDO(CHK), SMX(CHK)
      real (kind=kind_phys)::raddt,es(CHK),qs, delt, tem0d              &
     & ,SFCALBEDO(CHK), SMX(CHK)

      integer j2
      integer :: i, j, k, k1, lv, icec, itop, ibtc, nday, idxday(CHK),   &
     &      mbota(CHK,3), mtopa(CHK,3), LP1, nb, LMK, LMP, kd, lla, llb, &
     &      lya, lyb, kt, kb, np3d

!  ---  for debug test use
!     real (kind=kind_phys) :: temlon, temlat, alon, alat
!     integer :: ipt
!     logical :: lprnt1

      integer ipass,iorder  !jm

!!==========================================================================
!  --- Albedo (ETA era) calculation from 2012 RRTM  (Lin, 20130225)
!
      INTEGER :: IQ,JX(CHK)

!      INTEGER,EXTERNAL :: omp_get_thread_num

 !     REAL (kind=kind_phys) :: ALBD0, ALVD1, ALND1
      REAL (kind=kind_phys) :: ZEN, DZEN, ALB1, ALB2

      REAL (kind=kind_phys), PARAMETER :: TWENTY=20.0,        &
                                          HP537=0.537,        &
                                          ONE=1.,             &
                                          DEGRAD1=180.0/PI,   &
                                          H74E1=74.0,         &
                                          HAF=0.5,            &
                                          HNINETY=90.,        &
                                          FIFTY=50.,          &
                                          QUARTR=0.25,        &
                                          HNINE=9.0,          &
                                          HP1=0.1,            &
                                          H15E1=15.0

      REAL (kind=kind_phys), DIMENSION(IM) :: ALVB,ALNB,ALVD,ALND

      REAL (kind=kind_phys), DIMENSION(20) :: ZA
      REAL (kind=kind_phys), DIMENSION(19) :: DZA
      REAL (kind=kind_phys), DIMENSION(21) :: TRN
      REAL (kind=kind_phys), DIMENSION(21,20) :: ALBD

      DATA TRN/.00,.05,.10,.15,.20,.25,.30,.35,.40,.45,.50,.55,.60,.65,  &
               .70,.75,.80,.85,.90,.95,1.00/

      DATA  ALBD/.061,.062,.072,.087,.115,.163,.235,.318,.395,.472,.542, &
       .604,.655,.693,.719,.732,.730,.681,.581,.453,.425,.061,.062,.070, &
       .083,.108,.145,.198,.263,.336,.415,.487,.547,.595,.631,.656,.670, &
       .652,.602,.494,.398,.370,.061,.061,.068,.079,.098,.130,.174,.228, &
       .290,.357,.424,.498,.556,.588,.603,.592,.556,.488,.393,.342,.325, &
       .061,.061,.065,.073,.086,.110,.150,.192,.248,.306,.360,.407,.444, &
       .469,.480,.474,.444,.386,.333,.301,.290,.061,.061,.065,.070,.082, &
       .101,.131,.168,.208,.252,.295,.331,.358,.375,.385,.377,.356,.320, &
       .288,.266,.255,.061,.061,.063,.068,.077,.092,.114,.143,.176,.210, &
       .242,.272,.288,.296,.300,.291,.273,.252,.237,.266,.220,.061,.061, &
       .062,.066,.072,.084,.103,.127,.151,.176,.198,.219,.236,.245,.250, &
       .246,.235,.222,.211,.205,.200,                                    &
                 .061,.061,.061,.065,.071,.079,.094,.113,.134,.154,.173, &
       .185,.190,.193,.193,.190,.188,.185,.182,.180,.178,.061,.061,.061, &
       .064,.067,.072,.083,.099,.117,.135,.150,.160,.164,.165,.164,.162, &
       .160,.159,.158,.157,.157,.061,.061,.061,.062,.065,.068,.074,.084, &
       .097,.111,.121,.127,.130,.131,.131,.130,.129,.127,.126,.125,.122, &
       .061,.061,.061,.061,.062,.064,.070,.076,.085,.094,.101,.105,.107, &
       .106,.103,.100,.097,.096,.095,.095,.095,.061,.061,.061,.060,.061, &
       .062,.065,.070,.075,.081,.086,.089,.090,.088,.084,.080,.077,.075, &
       .074,.074,.074,.061,.061,.060,.060,.060,.061,.063,.065,.068,.072, &
       .076,.077,.076,.074,.071,.067,.064,.062,.061,.061,.061,.061,.061, &
       .060,.060,.060,.060,.061,.062,.065,.068,.069,.069,.068,.065,.061, &
       .058,.055,.054,.053,.052,.052,                                    &
                 .061,.061,.060,.060,.060,.060,.060,.060,.062,.065,.065, &
       .063,.060,.057,.054,.050,.047,.046,.045,.044,.044,.061,.061,.060, &
       .060,.060,.059,.059,.059,.059,.059,.058,.055,.051,.047,.043,.039, &
       .035,.033,.032,.031,.031,.061,.061,.060,.060,.060,.059,.059,.058, &
       .057,.056,.054,.051,.047,.043,.039,.036,.033,.030,.028,.027,.026, &
       .061,.061,.060,.060,.060,.059,.059,.058,.057,.055,.052,.049,.045, &
       .040,.036,.032,.029,.027,.026,.025,.025,.061,.061,.060,.060,.060, &
       .059,.059,.058,.056,.053,.050,.046,.042,.038,.034,.031,.028,.026, &
       .025,.025,.025,.061,.061,.060,.060,.059,.058,.058,.057,.055,.053, &
       .050,.046,.042,.038,.034,.030,.028,.029,.025,.025,.025/
!
      DATA ZA/90.,88.,86.,84.,82.,80.,78.,76.,74.,70.,66.,62.,58.,54., &
              50.,40.,30.,20.,10.,0.0/
!
      DATA DZA/8*2.0,6*4.0,5*10.0/
!
!  --- end of Albedo (ETA era) calculation from 2012 RRTM
!!==========================================================================

!
!===> ...  begin here
!
      np3d = icmphys

      LP1 = LM + 1               ! num of in/out levels

!===========================================================================
!  --- ...  set local /level/layer indexes corresponding to in/out variables
!           ( for optional extra top layer)

      LMK = LM + LTP             ! num of local layer
      LMP = LMK + 1              ! num of local level

      if ( lextop ) then
        if ( ivflip == 1 ) then   ! vertical from sfc upward
          kd = 0                   ! index diff between in/out and local
          kt = 1                   ! index diff between lyr and upper bound
          kb = 0                   ! index diff between lyr and lower bound
          lla = LMK                ! local index at the 2nd level from top
          llb = LMP                ! local index at toa level
          lya = LM                 ! local index for the 2nd layer from top
          lyb = LP1                ! local index for the top layer
        else                     ! vertical from toa downward
          kd = 1                   ! index diff between in/out and local
          kt = 0                   ! index diff between lyr and upper bound
          kb = 1                   ! index diff between lyr and lower bound
          lla = 2                  ! local index at the 2nd level from top
          llb = 1                  ! local index at toa level
          lya = 2                  ! local index for the 2nd layer from top
          lyb = 1                  ! local index for the top layer
        endif                    ! end if_ivflip_block
      else
        kd = 0
        if ( ivflip == 1 ) then  ! vertical from sfc upward
          kt = 1                   ! index diff between lyr and upper bound
          kb = 0                   ! index diff between lyr and lower bound
        else                     ! vertical from toa downward
          kt = 0                   ! index diff between lyr and upper bound
          kb = 1                   ! index diff between lyr and lower bound
        endif                    ! end if_ivflip_block
      endif   ! end if_lextop_block
!
!  --- End of index setting for optional extra top layer
!===========================================================================

      raddt = min(dtsw, dtlw)

!  --- ...  for debug test
!     alon = 120.0
!     alat = 29.5
!     ipt = 0
!     do i = 1, IM
!       temlon = xlon(i) * 57.29578
!       if (temlon < 0.0) temlon = temlon + 360.0
!       temlat = xlat(i) * 57.29578
!       lprnt1 = abs(temlon-alon) < 1.1 .and. abs(temlat-alat) < 1.1
!       if ( lprnt1 ) then
!         ipt = i
!         exit
!       endif
!     enddo

!     print *,' in grrad : raddt=',raddt

!  --- ...  setup surface ground temp and ground/air skin temp if required

      if ( itsfc == 0 ) then            ! use same sfc skin-air/ground temp
        do i = 1, IM
          tskn(i) = tsfc(i)
          tsfg(i) = tsfc(i)
        enddo
      else                              ! use diff sfc skin-air/ground temp
        do i = 1, IM
!!        tskn(i) = ta  (i)               ! not yet
!!        tsfg(i) = tg  (i)               ! not yet
          tskn(i) = tsfc(i)
          tsfg(i) = tsfc(i)
        enddo
      endif

!  --- ...  prepare atmospheric profiles for radiation input
!           convert pressure unit from Pa to mb

!     if (im > ipt) then
!       write(0,*)' prsi=',prsi(ipt,1:10)
!       write(0,*)' prsi=',prsl(ipt,1:10)
!       write(0,*)' tgrs=',tgrs(ipt,1:10)
!     endif

!=============================================================
! ** NOTE: For NMMB, the conversion is "from cb to mb"
!=============================================================

      do k = 1, LM
        k1 = k + kd      ! kd: index diff between in/out and local

!DEC$ SIMD
        do i = 1, IM
          plvl(i,k1) = 10.0 * prsi(i,k)   ! cb to mb
          plyr(i,k1) = 10.0 * prsl(i,k)   ! cb to mb
!         plvl(i,k1) = 0.01 * prsi(i,k)   ! Pa to mb
!         plyr(i,k1) = 0.01 * prsl(i,k)   ! Pa to mb
          tlyr(i,k1) = tgrs(i,k)
          prslk1(i,k1) = prslk(i,k)
        ENDDO

!  --- ...  compute relative humidity

          ! Note for qs (**Saturation specific humidity): 
          !   ** prsl unit in NMMB is kPa
          !   it can be done in es (**Saturation vapor pressure) for 0.001*fpvs
          !   or multiple every prsl by 1000
          !   (Hsin-mu Lin, 20110927)

! note this loop won't vectorize because of the call to fpvs

        DO i = 1, IM
          es(i) = min( prsl(i,k), 0.001 * fpvs( tgrs(i,k) ) ) ! fpvs in pa
        ENDDO

!DEC$ SIMD
        DO i = 1, IM
          qs = max( QMIN,con_eps*es(i)/(prsl(i,k) + con_epsm1*es(i)) )
          rhly(i,k1) = max( 0.0, min( 1.0, max(QMIN, qgrs(i,k))/qs ) )
          qstl(i,k1) = qs
        enddo

        do j = 1, NTRAC
          do i = 1, IM
             tracer1(i,k1,j) = tracer(i,k,j)
          enddo
        enddo
      enddo

      do i = 1, IM
        plvl(i,LP1+kd) = 10.0 * prsi(i,LP1)  ! cb to mb
!       plvl(i,LP1+kd) = 0.01 * prsi(i,LP1)  ! Pa to mb
      enddo

!======================================================================
!  --- ...  values for extra top layer

      if ( lextop ) then                 ! values for extra top layer
!DEC$ SIMD
        do i = 1, IM
          plvl(i,llb) = prsmin
          if ( plvl(i,lla)<=prsmin ) plvl(i,lla) = 2.0*prsmin
          plyr(i,lyb) = 0.5 * plvl(i,lla)
          tlyr(i,lyb) = tlyr(i,lya)
          rhly(i,lyb) = rhly(i,lya)
          qstl(i,lyb) = qstl(i,lya)
        enddo

! loop won't vectorize because of call to fpkapx
        do i = 1, IM
          prslk1(i,lyb) = fpkapx(plyr(i,lyb)*100.0) ! fpkapx in pa
        enddo

!  ---  note: may need to take care the top layer mount

         do i = 1, IM
           tracer1(i,lyb,1) = tracer1(i,lya,1)    ! 1st tracer element is ozone
         enddo

         do j = 2, NTRAC     ! This loop variable should now be correct
           do i = 1, IM
             tracer1(i,lyb,j) = 0.       !  enforce no cloud above domain top
           enddo
         enddo

      endif
!
!======================================================================

!======================================================================
!  --- ...  extra variables needed for using Thompson cloud fraction in
!           progcld2 & progcld8 for rrtm. For non Thompson cloud fraction
!           the value will be "0" as initialization in module_RA_RRTM
!
      do k = 1, LM
        k1 = k + kd

        do i = 1, IM
           cldfk1(i,k1)= cldf(i,k)
        enddo
      enddo

      !========================================
      !  --- ...  values for extra top layer

      if ( lextop ) then
        do i = 1, IM
           cldfk1(i,lyb) = cldf(i,lya)
        enddo
      endif

!
!======================================================================

!======================================================================
!  --- ...  extra variables needed for ferrier's microphysics

      if (np3d == 3) then
        do k = 1, LM
          k1 = k + kd

          do i = 1, IM
            gcice(i,k1)= fcice(i,k)
            grain(i,k1)= frain(i,k)
            grime(i,k1)= rrime(i,k)
          enddo
        enddo

        !========================================
        !  --- ...  values for extra top layer

        if ( lextop ) then
          do i = 1, IM
            gcice(i,lyb) = 1.0   ! fcice(i,lya)
            grain(i,lyb) = 0.0   ! frain(i,lya)
            grime(i,lyb) = 1.0   ! rrime(i,lya)
          enddo
        endif
      endif   ! if_np3d

!
!======================================================================

!  --- ...  get layer ozone mass mixing ratio

      if (ntoz > 0) then            ! interactive ozone generation

        do k = 1, LMK
          do i = 1, IM
            olyr(i,k) = max( QMIN, tracer1(i,k,ntoz) )
          enddo
        enddo

      else                          ! climatological ozone

!     print *,' in grrad : calling getozn'

        call getozn                                                     &
!  ---  inputs:
     &     ( prslk1,xlat,                                               &
     &       IM, LMK,                                                   &
!  ---  outputs:
     &       olyr                                                       &
     &     )

      endif                            ! end_if_ntoz

!  --- ...  compute cosin of zenith angle

      call coszmn_nmmb                                                  &
!  ---  inputs:
     &     ( xlon,sinlat,coslat,solhr,IM, me,                           &
     &       dtswav,nrads ,                                             &  ! Extra input
!  ---  outputs:
     &       coszen_loc, coszdg                                         &  !jm
     &      )

!  --- ...  set up non-prognostic gas volume mixing ratioes

      call getgases                                                     &
!  ---  inputs:
     &    ( plvl, xlon, xlat,                                           &
     &      IM, LMK,                                                    &
!  ---  outputs:
     &      gasvmr                                                      &
     &     )

!  --- ...  get temperature at layer interface, and layer moisture

      do k = 2, LMK
        do i = 1, IM
          tem2da(i,k) = log( plyr(i,k) )
          tem2db(i,k) = log( plvl(i,k) )
        enddo
      enddo

      if (ivflip == 0) then               ! input data from toa to sfc

        do i = 1, IM
          tem1d (i)   = QME6
          tem2da(i,1) = log( plyr(i,1) )

        !---------------------------------------------
        ! take care for the model top (20130410, Lin)
        !---------------------------------------------

          if ( plvl(i,1) == 0.0 ) then
             tem2db(i,1) = -14.0
          else
             tem2db(i,1) = log( plvl(i,1) )
          endif

         ! tem2db(i,1) = 1.0

          tsfa  (i)   = tlyr(i,LMK)                  ! sfc layer air temp
          tlvl(i,1)   = tlyr(i,1)
          tlvl(i,LMP) = tskn(i)
        enddo

        do k = 1, LM
          k1 = k + kd      ! kd: index diff between in/out and local

          do i = 1, IM
            qlyr(i,k1) = max( tem1d(i), qgrs(i,k) )
            tem1d(i)   = min( QME5, qlyr(i,k1) )
            tvly(i,k1) = tgrs(i,k) * (1.0 + con_fvirt*qlyr(i,k1))          ! virtual temp in K
          enddo
        enddo

        !============================================
        !  --- ...  values for extra top layer

        if ( lextop ) then
          do i = 1, IM
            qlyr(i,lyb) = qlyr(i,lya)
            tvly(i,lyb) = tvly(i,lya)
          enddo
        endif

        !============================================

        do k = 2, LMK
          do i = 1, IM
            tlvl(i,k) = tlyr(i,k) + (tlyr(i,k-1) - tlyr(i,k))           &
     &                * (tem2db(i,k)   - tem2da(i,k))                   &
     &                / (tem2da(i,k-1) - tem2da(i,k))
          enddo
        enddo

      else                               ! input data from sfc to toa

        do i = 1, IM
          tem1d (i)   = QME6
          tem2da(i,1) = log( plyr(i,1) )
          tem2db(i,1) = log( plvl(i,1) )
          tsfa  (i)   = tlyr(i,1)                    ! sfc layer air temp
          tlvl(i,1)   = tskn(i)
          tlvl(i,LMP) = tlyr(i,LMK)
        enddo

        do k = LM, 1, -1
          do i = 1, IM
            qlyr(i,k) = max( tem1d(i), qgrs(i,k) )
            tem1d(i)  = min( QME5, qlyr(i,k) )
            tvly(i,k) = tgrs(i,k) * (1.0 + con_fvirt*qlyr(i,k))          ! virtual temp in K
          enddo
        enddo

        !============================================
        !  --- ...  values for extra top layer

        if ( lextop ) then
          do i = 1, IM
            qlyr(i,lyb) = qlyr(i,lya)
            tvly(i,lyb) = tvly(i,lya)
          enddo
        endif

        !============================================

        do k = 1, LMK-1
          do i = 1, IM
            tlvl(i,k+1) = tlyr(i,k) + (tlyr(i,k+1) - tlyr(i,k))         &
     &                  * (tem2db(i,k+1) - tem2da(i,k))                 &
     &                  / (tem2da(i,k+1) - tem2da(i,k))
          enddo
        enddo

      endif                              ! end_if_ivflip

!  ---  check for daytime points
!jm with modification in RA_RRTM, the grrad_nmm routine will 
!jm only be called for all day or all night points, the number
!jm of which is passed in as IM.  
      nday = im

!  --- ...  setup aerosols property profile for radiation

!check  print *,' in grrad : calling setaer '

      call setaer                                                     &
!  ---  inputs:
     &   ( plvl,plyr,prslk1,tvly,rhly,slmsk,tracer1,xlon,xlat,        &
     &     IM,LMK,LMP, lsswr,lslwr,                                   &
!  ---   extra 2 vars for old code that has been replaced by 'tvly'
!        can be eliminated after aggrements being reached
     &     tlyr,qlyr,                                                 &
!  ---  outputs:
     &     faersw,faerlw                                              &
 !   &,    tau_gocart                                                 &
     &   )

!  --- ...  obtain cloud information for radiation calculations

      if (ntcw > 0) then                   ! prognostic cloud scheme

        do k = 1, LMK
          do i = 1, IM
            clw(i,k) = 0.0
          enddo

          do j = 1, ncld
            lv = ntcw + j - 1
            do i = 1, IM
               clw(i,k) = clw(i,k) + tracer1(i,k,lv)    ! cloud condensate amount
            enddo
          enddo
        enddo

        do k = 1, LMK
          do i = 1, IM
            if ( clw(i,k) < EPSQ ) clw(i,k) = 0.0
          enddo
        enddo

        if (np3d == 4) then              ! zhao/moorthi's prognostic cloud scheme

          call progcld1                                                 &
!  ---  inputs:
     &     ( plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw,                    &
     &       xlat,xlon,slmsk,                                           &
     &       IM, LMK, LMP,                                              &
!  ---  outputs:
     &       clouds,cldsa,mtopa,mbota                                   &
     &      )

        elseif (np3d == 3) then          ! ferrier's microphysics

!     print *,' in grrad : calling progcld2'
          call progcld2                                                 &
!  ---  inputs:
     &     ( plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw,                    &
     &       xlat,xlon,slmsk, gcice,grain,grime,flgmin,                 &
     &       IM, LMK, LMP,                                              &
     &       cldfk1, cld_fraction,                                      &
!  ---  outputs:
     &       clouds,cldsa,mtopa,mbota                                   &
     &      )

        elseif (np3d == 5) then         ! nmmb ferrier+bmj  (GFDL type diagnostic)
!
! ======================================================================
!  The original GFS RRTM will use "call progcld3".
! ** After merging the nmmb & GFS RRTM, "np3d=5" is temporarily used
!    for the GFDL diagnostic cloud type which is in the original meso
!    nmmb                                    (Hsin-mu Lin 12/28/2011)
! ======================================================================
!
!         call progcld3                                                 &
!  ---  inputs:
!    &     ( plyr,plvl,tlyr,qlyr,qstl,rhly,clw,                         &
!    &       xlat,xlon,slmsk, fcice,frain,rrime,flgmin,                 &
!    &       IM, LMK, LMP, iflip, iovrsw,                               &
!  ---  outputs:
!    &       clouds,cldsa,mtopa,mbota                                   &
!    &      )

 !     print *,' in grrad : !!! need to develop progcld3 for nmmb',      &
 !    &        ' temporarily using progcld2 !!!'

          do k = 1, LM                                          ! GFDL type
            k1 = k + kd      ! kd: index diff between in/out and local

            do i = 1, IM                                        ! GFDL type
              clouds(i,k1,1) = cldf(i,k)                        ! GFDL type
              clouds(i,k1,2) = tauclouds(i,k)                   ! GFDL type
              clouds(i,k1,3) = 0.99                             ! GFDL type
              clouds(i,k1,4) = 0.84                             ! GFDL type
            enddo                                               ! GFDL type
          enddo                                                 ! GFDL type

          !============================================
          !  --- ...  values for extra top layer

          if ( lextop ) then
            do i = 1, IM
              clouds(i,lyb,1) = clouds(i,lya,1)
              clouds(i,lyb,2) = clouds(i,lya,2)
              clouds(i,lyb,3) = clouds(i,lya,3)
              clouds(i,lyb,4) = clouds(i,lya,4)
            enddo
          endif

          !============================================

!..This is far from optimal!  The specific array element numbers are
!.. potentially flexible for cloud water, cloud ice, snow mixing ratios
!.. and cloud ice number concentration.  Therefore, numerical values of
!.. 4, 5, 6, 9 are hard-wired when they should be made flexible.

        elseif (np3d == 8) then            ! Thompson microphysics

!         write(*,*) ' DEBUG-GT, double-check radiation flags ',iswcliq , iswcice, ilwcliq, ilwcice
!         write(*,*) ' DEBUG-GT, double-check NTRAC number ', ntrac
!         write(*,*) ' DEBUG-GT, double-check LMK (LM+LTP) ', LMK

          do k = 1, LMK
            do i = 1, IM
              qc2d(i,k) = tracer1(i,k,4)
              if (qc2d(i,k) .LE. EPSQ) qc2d(i,k)=0.0
              qi2d(i,k) = tracer1(i,k,5)
              if (qi2d(i,k) .LE. EPSQ) qi2d(i,k)=0.0
              qs2d(i,k) = tracer1(i,k,6)
              if (qs2d(i,k) .LE. EPSQ) qs2d(i,k)=0.0
              ni2d(i,k) = MAX(1.E-8, tracer1(i,k,9))
            enddo
          enddo


          call progcld8 (plyr,plvl, tlyr, qlyr, qc2d, qi2d, qs2d, ni2d, &
     &                   xlat, IM, LMK, LMP,                            &
     &                   cldfk1, cld_fraction,                          &
     &                   clouds, cldsa, mtopa, mbota)

        endif                              ! end if_np3d
!
      else                                 ! diagnostic cloud scheme

        do i = 1, IM
          cvt1(i) = 10.0 * cvt(i)
          cvb1(i) = 10.0 * cvb(i)
!         cvt1(i) = 0.01 * cvt(i)
!         cvb1(i) = 0.01 * cvb(i)
        enddo

        do k = 1, LM
          k1 = k + kd      ! kd: index diff between in/out and local

          do i = 1, IM
            vvel(i,k1) = 10.0 * vvl (i,k)
!           vvel(i,k1) = 0.01 * vvl (i,k)
          enddo
        enddo

        !============================================
        !  --- ...  values for extra top layer

        if ( lextop ) then
          do i = 1, IM
            vvel(i,lyb) = vvel(i,lya)
          enddo
        endif

        !============================================

!  ---  compute diagnostic cloud related quantities

        call diagcld1                                                   &
!  ---  inputs:
     &     ( plyr,plvl,tlyr,rhly,vvel,cv,cvt1,cvb1,                     &
     &       xlat,xlon,slmsk,                                           &
     &       IM, LMK, LMP,                                              &
!  ---  outputs:
     &       clouds,cldsa,mtopa,mbota                                   &
     &      )


      endif                                ! end_if_ntcw

!==========================================================================
!  Albedo (ETA era) calculation from 2012 RRTM  (Lin, 20130225)
!==========================================================================
      IF (ialbflg == 2) THEN

         do i = 1, IM

            SMX(i) = SM(i)
            SFCALBEDO(i) = SALBEDO(i)

!..... THE FOLLOWING CODE GETS ALBEDO FROM PAYNE,1972 TABLES IF
!         1) OPEN SEA POINT (SLMSK=1);  2) KALB=0

            IQ=INT(TWENTY*HP537+ONE)
            CZMEAN(i)=coszen_loc(i)   !jm

            IF(CZMEAN(i).GT.0.0 .AND. SMX(i).GT.0.5) THEN
               ZEN=DEGRAD1*ACOS(MAX(CZMEAN(i),0.0))

               IF(ZEN.GE.H74E1) JX(I)=INT(HAF*(HNINETY-ZEN)+ONE)
               IF(ZEN.LT.H74E1 .AND. ZEN.GE.FIFTY) &
                 JX(I)=INT(QUARTR*(H74E1-ZEN)+HNINE)
               IF(ZEN .LT. FIFTY) JX(I)=INT(HP1*(FIFTY-ZEN)+H15E1)

               DZEN=-(ZEN-ZA(JX(I)))/DZA(JX(I))

               ALB1=ALBD(IQ,JX(I))+DZEN*(ALBD(IQ,JX(I)+1)&
                     -ALBD(IQ,JX(I)))
               ALB2=ALBD(IQ+1,JX(I))+DZEN*(ALBD(IQ+1,JX(I)+1)&
                     -ALBD(IQ+1,JX(I)))

               SFCALBEDO(I)=ALB1+TWENTY*(ALB2-ALB1)*(HP537-TRN(IQ)) ! BSF
            ENDIF

!.....     VISIBLE AND NEAR IR DIFFUSE ALBEDO
            ALVD(I) = SFCALBEDO(I) ! BSF
            ALND(I) = SFCALBEDO(I) ! BSF

!.....     VISIBLE AND NEAR IR DIRECT BEAM ALBEDO
            ALVB(I) = SFCALBEDO(I) ! BSF
            ALNB(I) = SFCALBEDO(I) ! BSF
!
!--- Remove diurnal variation of land surface albedos (Ferrier, 6/28/05)
!--- Turn back on to mimic NAM 8/17/05
!
!.....     VISIBLE AND NEAR IR DIRECT BEAM ALBEDO, IF NOT OCEAN NOR SNOW
!        ..FUNCTION OF COSINE SOLAR ZENITH ANGLE..
!
!=== The following line commented out by the "fixed" are used instead of the
!    original version (09/2012)
!
!fixed            IF (SMX .LT. 0.5) THEN
!fixed             IF (SFCALBEDO .LE. 0.5) THEN
!fixed             ALBD0=-18.0 * (0.5 - ACOS(CZMEAN(I))/PI)
!fixed             ALBD0=EXP (ALBD0)
!fixed             ALVD1=(ALVD(I) - 0.054313) / 0.945687
!fixed             ALND1=(ALND(I) - 0.054313) / 0.945687
!fixed             ALVB(I)=ALVD1 + (1.0 - ALVD1) * ALBD0
!fixed             ALNB(I)=ALND1 + (1.0 - ALND1) * ALBD0
!fixed! !-- Put in an upper limit on beam albedos
!fixed             ALVB(I)=MIN(0.5,ALVB(I))
!fixed             ALNB(I)=MIN(0.5,ALNB(I))
!fixed            ENDIF
!fixed           ENDIF

!!! WE INTRODUCE HERE DIRECT AND DIFFUSE ALBEDO... FOR THIS OPTION, THERE IS A CHANGE IN GRRAD.f

            alvsf(i) = ALVB(I) !For this option ALVSF1 is direct visible albedo
            alnsf(i) = ALNB(I) !For this option ALNSF1 is direct nir albedo
            alvwf(i) = ALVD(I) !For this option ALVWF1 is diffuse visible albedo
            alnwf(i) = ALND(I) !For this option ALNWF1 is diffuse nir albedo
            facsf(i) = 0.      !not used with this option
            facwf(i) = 0.      !not used for this option

         ENDDO
      ENDIF    !---- end of IALB=2  GFDL TYPE RADIATION

!
! --- End of Albedo (ETA era) calculation from 2012 RRTM ------------------
!==========================================================================



!  --- ...  start radiation calculations 
!           remember to set heating rate unit to k/sec!

!#ifdef ORIG_LWSW_ORDERING
      iorder = 0
!#else
!      iorder = mod(omp_get_thread_num(),2)   ! every other thread
!      iorder = mod(omp_get_thread_num()/2,2)    ! alternate by groups of 4
!#endif
      DO ipass=0,1
      IF( ipass .eq. iorder ) then

      if (lsswr) then

!  ---  setup surface albedo for sw radiation, incl xw (nov04) sea-ice

        call setalb                                                     &
!  ---  inputs:
     &     ( slmsk,snowd,sncovr,snoalb,zorl,coszen_loc,tsfg,tsfa,hprim, &  !jm
     &       alvsf,alnsf,alvwf,alnwf,facsf,facwf,fice,tisfc,            &
     &       IM,                                                        &
!  ---  outputs:
     &       sfcalb                                                     &
     &     )

!  --- lu [+4L]: derive SFALB from vis- and nir- diffuse surface albedo

        do i = 1, IM
          sfalb(i) = max(0.01, 0.5 * (sfcalb(i,2) + sfcalb(i,4)))
        enddo

!        if (nday > 0) then
        if (isday) then

!     print *,' in grrad : calling swrad'


!jm pad out vector
#if 1
        if ( nday < CHK ) then
          coszen_loc(nday+1:CHK) = coszen_loc(nday)
          do k = 1,lmp
            plvl   (nday+1:CHK,k)  = plvl  (nday,k)
            tlvl   (nday+1:CHK,k)  = tlvl  (nday,k)
          enddo
          do k = 1,lmk
            plyr   (nday+1:CHK,k)  = plyr  (nday,k)
            tlyr   (nday+1:CHK,k)  = tlyr  (nday,k)
            qlyr   (nday+1:CHK,k)  = qlyr  (nday,k)
            olyr   (nday+1:CHK,k)  = olyr  (nday,k)
          enddo
          do j = 1,4
            sfcalb (nday+1:CHK,j)  = sfcalb  (nday,j)
          enddo
          do j = 1,9
            do k = 1,lmk
              gasvmr    (nday+1:CHK,k,j)  = gasvmr  (nday,k,j)
              clouds    (nday+1:CHK,k,j)  = clouds  (nday,k,j)
            enddo
          enddo
          do j = 1,3
            do j2 = 1,nbdsw
              do k = 1,lmk
                faersw (nday+1:CHK,k,j2,j)  = faersw (nday,k,j2,j)
              enddo
            enddo
          enddo
        endif
#endif

          if ( present(htrswb) ) then

            call swrad                                                  &
!  ---  inputs:
     &     ( plyr,plvl,tlyr,tlvl,qlyr,olyr,gasvmr,                      &
     &       clouds,icsdsw_loc,faersw,sfcalb,                           &
     &       coszen_loc,solcon, nday,idxday,                            &
     &       IM, LMK, LMP, lprnt,                                       &
!  ---  outputs:
     &       htswc,topfsw_loc,sfcfsw_loc                                &
!! ---  optional:
!!   &,      HSW0=htsw0,FLXPRF=fswprf                                   &
     &,      HSWB=htswb,FDNCMP=scmpsw                                   &
     &,      ITS=its,JTS=jts                                            &
     &     )

            do k = 1, LM
              k1 = k + kd      ! kd: index diff between in/out and local

              do j = 1, NBDSW
                do i = 1, IM
                  htrswb(i,k,j) = htswb(i,k1,j)
                enddo
              enddo
            enddo

          else

            call swrad                                                  &
!  ---  inputs:
     &     ( plyr,plvl,tlyr,tlvl,qlyr,olyr,gasvmr,                      &
     &       clouds,icsdsw_loc,faersw,sfcalb,                           &
     &       coszen_loc,solcon, nday,idxday,                            &
     &       IM, LMK, LMP, lprnt,                                       &
!  ---  outputs:
     &       htswc,topfsw_loc,sfcfsw_loc                                &
!! ---  optional:
!!   &,      HSW0=htsw0,FLXPRF=fswprf,HSWB=htswb                        &
     &,      FDNCMP=scmpsw                                              &
     &,      ITS=its,JTS=jts                                            &
     &     )

          endif


          do k = 1, LM
            k1 = k + kd      ! kd: index diff between in/out and local

            do i = 1, IM
              htrsw(i,k) = htswc(i,k1)
            enddo
          enddo

        else                   ! if_isday_block

          do k = 1, LM
            do i = 1, IM
              htrsw(i,k) = 0.0
            enddo
          enddo

          sfcfsw_loc = sfcfsw_type( 0.0, 0.0, 0.0, 0.0 )
          topfsw_loc = topfsw_type( 0.0, 0.0, 0.0 )
          scmpsw = cmpfsw_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 )

!! ---  optional:
!!        fswprf= profsw_type( 0.0, 0.0, 0.0, 0.0 )

          if ( present(htrswb) ) then
            do j = 1, NBDSW
              do k = 1, LM
                do i = 1, IM
                  htrswb(i,k,j) = 0.0
                enddo
              enddo
            enddo
          endif

        endif                  ! end_if_isday

        topfsw(1:IM) = topfsw_loc(1:IM)
        sfcfsw(1:IM) = sfcfsw_loc(1:IM)

      endif                                ! end_if_lsswr

      ELSE  ! ipass == iorder

      if (lslwr) then

!  ---  setup surface emissivity for lw radiation

        call setemis                                                    &
!  ---  inputs:
     &     ( xlon,xlat,slmsk,snowd,sncovr,zorl,tsfg,tsfa,hprim,         &
     &       IM,                                                        &
!  ---  outputs:
     &       sfcemis                                                    &
     &     )

!     print *,' in grrad : calling lwrad'

!==========================================================================
!  Special for lwrad to enhence the emissivity
!  it is similar to *CPATHFAC4LW to odcld in radlw (Hsin-Mu Lin, 20140520)
!==========================================================================

        if (np3d == 3) then          ! ferrier's microphysics
          do k = 1, LMK
            do i = 1, IM
              dummys(i,k,2) = clouds(i,k,2)
              dummys(i,k,4) = clouds(i,k,4)
              clouds(i,k,2) = dummys(i,k,2)*CPATHFAC4LW  ! cloud liquid path 
              clouds(i,k,4) = dummys(i,k,4)*CPATHFAC4LW  ! cloud ice path
            enddo
          enddo
        endif

!==========================================================================

        if ( present(htrlwb) ) then
          call lwrad                                                    &
!  ---  inputs:
     &     ( plyr,plvl,tlyr,tlvl,qlyr,olyr,gasvmr,                      &
     &       clouds,icsdlw,faerlw,sfcemis,tsfg,                         &
     &       IM, LMK, LMP, lprnt,                                       &
!  ---  outputs:
     &       htlwc,topflw,sfcflw                                        &
!! ---  optional:
!!   &,      HLW0=htlw0,FLXPRF=flwprf                                   &
     &,      HLWB=htlwb                                                 &
     &,      ITS=its,JTS=jts                                            &
     &     )

          do k = 1, LM
            k1 = k + kd      ! kd: index diff between in/out and local

            do j = 1, NBDLW
              do i = 1, IM
                htrlwb(i,k,j) = htlwb(i,k1,j)
              enddo
            enddo
          enddo

        else

!          do i = 1, IM
!      if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k sfcemis",3i7,e25.15)')its+i-1,jts,k,sfcemis(i)
!      if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k tsfg",3i7,e25.15)')its+i-1,jts,k,tsfg(i)
!      if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k icsdlw",4i7)')its+i-1,jts,k,icsdlw(i)
!        do k = 1, LM
!      if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k plyr",3i7,e25.15)')its+i-1,jts,k,plyr(i,k)
!      if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k plvl",3i7,e25.15)')its+i-1,jts,k,plvl(i,k)
!      if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k tlyr",3i7,e25.15)')its+i-1,jts,k,tlyr(i,k)
!      if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k tlvl",3i7,e25.15)')its+i-1,jts,k,tlvl(i,k)
!      if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k qlyr",3i7,e25.15)')its+i-1,jts,k,qlyr(i,k)
!      if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k olyr",3i7,e25.15)')its+i-1,jts,k,olyr(i,k)
!        enddo
!          enddo

          call lwrad                                                    &
!  ---  inputs:
     &     ( plyr,plvl,tlyr,tlvl,qlyr,olyr,gasvmr,                      &
     &       clouds,icsdlw,faerlw,sfcemis,tsfg,                         &
     &       IM, LMK, LMP, lprnt,                                       &
!  ---  outputs 
     &       htlwc,topflw,sfcflw                                        &
!! ---  optional:
!!   &,      HLW0=htlw0,FLXPRF=flwprf,HLWB=htlwb                        &
     &,      ITS=its,JTS=jts                                            &
     &     )

        endif

        do i = 1, IM
          semis (i) = sfcemis(i)
!  ---  save surface air temp for diurnal adjustment at model t-steps
          tsflw (i) = tsfa(i)
        enddo

        do k = 1, LM
          k1 = k + kd      ! kd: index diff between in/out and local

          do i = 1, IM
            htrlw(i,k) = htlwc(i,k1)
!      if (k.eq.1.and.its+i-1.eq.140.and.jts.eq.47)write(0,'("two IM,its,jts",2i7,e25.15)')its+i-1,jts,htrlw(i,1)
          enddo
        enddo

!==========================================================================
!   return *1.5 clouds to original value  (Hsin-Mu Lin, 20140520)
!==========================================================================

        if (np3d == 3) then          ! ferrier's microphysics
          do k = 1, LMK
            do i = 1, IM
              clouds(i,k,2) = dummys(i,k,2)  ! cloud liquid path
              clouds(i,k,4) = dummys(i,k,4)  ! cloud ice path
            enddo
          enddo
        endif

!==========================================================================

      endif                                ! end_if_lslwr

      ENDIF  ! iorder
      ENDDO  ! jm ipass

!  --- ...  collect the fluxr data for wrtsfc

      if (lssav) then

!  ---  save lw toa and sfc fluxes

        if (lslwr) then
          do i = 1, IM
!  ---  lw total-sky fluxes
            fluxr(i,1 ) = fluxr(i,1 ) + dtlw * topflw(i)%upfxc   ! total sky top lw up
            fluxr(i,19) = fluxr(i,19) + dtlw * sfcflw(i)%dnfxc   ! total sky sfc lw dn
            fluxr(i,20) = fluxr(i,20) + dtlw * sfcflw(i)%upfxc   ! total sky sfc lw up
!  ---  lw clear-sky fluxes
            fluxr(i,28) = fluxr(i,28) + dtlw * topflw(i)%upfx0   ! clear sky top lw up
            fluxr(i,30) = fluxr(i,30) + dtlw * sfcflw(i)%dnfx0   ! clear sky sfc lw dn
            fluxr(i,33) = fluxr(i,33) + dtlw * sfcflw(i)%upfx0   ! clear sky sfc lw up
          enddo
        endif

!  ---  save sw toa and sfc fluxes with proper diurnal sw wgt. coszen=mean cosz over daylight
!       part of sw calling interval, while coszdg= mean cosz over entire interval

        if (lsswr) then
          do i = 1, IM
            if (coszen_loc(i) > 0.) then
!  ---  sw total-sky fluxes
              tem0d = dtsw * coszdg(i) / coszen_loc(i)
              fluxr(i,2 ) = fluxr(i,2)  + topfsw(i)%upfxc * tem0d  ! total sky top sw up
              fluxr(i,3 ) = fluxr(i,3)  + sfcfsw(i)%upfxc * tem0d  ! total sky sfc sw up
              fluxr(i,4 ) = fluxr(i,4)  + sfcfsw(i)%dnfxc * tem0d  ! total sky sfc sw dn
!  ---  sw uv-b fluxes
              fluxr(i,21) = fluxr(i,21) + scmpsw(i)%uvbfc * tem0d  ! total sky uv-b sw dn
              fluxr(i,22) = fluxr(i,22) + scmpsw(i)%uvbf0 * tem0d  ! clear sky uv-b sw dn
!  ---  sw toa incoming fluxes
              fluxr(i,23) = fluxr(i,23) + topfsw(i)%dnfxc * tem0d  ! top sw dn
!  ---  sw sfc flux components
              fluxr(i,24) = fluxr(i,24) + scmpsw(i)%visbm * tem0d  ! uv/vis beam sw dn
              fluxr(i,25) = fluxr(i,25) + scmpsw(i)%visdf * tem0d  ! uv/vis diff sw dn
              fluxr(i,26) = fluxr(i,26) + scmpsw(i)%nirbm * tem0d  ! nir beam sw dn        !?
              fluxr(i,27) = fluxr(i,27) + scmpsw(i)%nirdf * tem0d  ! nir diff sw dn        !?
!  ---  sw clear-sky fluxes
              fluxr(i,29) = fluxr(i,29) + topfsw(i)%upfx0 * tem0d  ! clear sky top sw up
              fluxr(i,31) = fluxr(i,31) + sfcfsw(i)%upfx0 * tem0d  ! clear sky sfc sw up
              fluxr(i,32) = fluxr(i,32) + sfcfsw(i)%dnfx0 * tem0d  ! clear sky sfc sw dn
            endif
          enddo
        endif

!  ---  save total cloud and bl cloud

        if (lsswr .or. lslwr) then
          do i = 1, IM
            fluxr(i,17) = fluxr(i,17) + raddt * cldsa(i,4)             !?
            fluxr(i,18) = fluxr(i,18) + raddt * cldsa(i,5)             !?
          enddo

!  ---  save cld frac,toplyr,botlyr and top temp, note that the order
!       of h,m,l cloud is reversed for the fluxr output.
!  ---  save interface pressure (Pa) of top/bot

          do j = 1, 3
            do i = 1, IM
              tem0d = raddt * cldsa(i,j)
              itop  = mtopa(i,j) - kd
              ibtc  = mbota(i,j) - kd

            !==========================================
            !=== Zavisa mentioned on 01/15/2013
            !    for our current run, there are no problem for np3d=3
            !    without this modification, while np3d=5 need this

              itop=min(max(itop,1),lm) !zj quick fix
              ibtc=min(max(ibtc,1),lm) !zj quick fix
 
            !==========================================

              fluxr(i, 8-j) = fluxr(i, 8-j) + tem0d
              fluxr(i,11-j) = fluxr(i,11-j) + prsi(i,itop+kt) * tem0d
              fluxr(i,14-j) = fluxr(i,14-j) + prsi(i,ibtc+kb) * tem0d
              fluxr(i,17-j) = fluxr(i,17-j) + tgrs(i,itop) * tem0d
            enddo
          enddo
        endif

        do k = 1, LM
          k1 = k + kd       ! kd: index diff between in/out and local

          do i = 1, IM
            cldcov(i,k) = clouds(i,k1,1)
          enddo
        enddo

!  ---  save optional vertically integrated aerosol optical depth at
!       wavelenth of 550nm aerodp(:,1), and other optional aod for
!       individual species aerodp(:,2:NSPC1)

!       if ( laswflg ) then
!         if ( NFXR > 33 ) then
!           do i = 1, IM
!             fluxr(i,34) = fluxr(i,34) + dtsw*aerodp(i,1)  ! total aod at 550nm (all species)
!           enddo

!           if ( lspcodp ) then
!             do j = 2, NSPC1
!               k = 33 + j

!               do i = 1, IM
!                 fluxr(i,k) = fluxr(i,k) + dtsw*aerodp(i,j) ! aod at 550nm for indiv species
!               enddo
!             enddo
!           endif     ! end_if_lspcodp
!         else
!           print *,'  !Error! Need to increase array fluxr size NFXR ',&
!    &              ' to be able to output aerosol optical depth'
!           stop
!         endif     ! end_if_nfxr
!       endif       ! end_if_laswflg

      endif                                ! end_if_lssav

      do i = 1, IM
        coszen(i) = coszen_loc(i)   ! jm
      enddo

!
      return
!...................................
      end subroutine grrad_nmmb
!-----------------------------------


!
!.............................................!
      end module module_radiation_driver_nmmb !
!=============================================!