Bug fix for Consine Solar Zenith Angle calculation

This commit implements this bug fix the way CESM
fixed it in their latest tag(cesm1_4_beta03).
This commit replaces the flag sol_insolation_bug_fix with
use_rad_dt_cosz. This flag is .false.
by default as turning this true is an answer changing
change (not climate changing). See Issue #139 for more
details.

This commit replaces sol_insolation_bug_fix logical flag
in the cam_inparm block of atm_in namelist

Other modifications:
1. Modified a comment for ndrop.F90 bug fix
2. Removed zenith.F90 from cam/src/physics/cam and added
   orbit.F90 with module "orbit" and subroutine "zenith" in
   cam/src/utils/.Basically moved "zenith" subroutine in
   orbit module.
3. Nag compiler complained about #error CPP directive.Removed
   "#" from "#error" to make it "error" to fix that. This
   fix retains the compile time failure functionality.
4. Shortened lines longer than 132 characters in spelt_mod.F90
5. Made emin and emax intent-inout in viscosity_mod.F90 as they
   are used without being initialized to a value.

Fixes #139
[NML]
AG-189

models/atm/cam/bld/build-namelist

@@ -779,9 +779,9 @@ else {
 add_default($nl, 'ncdata', %atts);
 
 #BSINGH - (when .true.)calculate solar zenith angle averaged over a time step.
-#In current default method solar zenith angle is held constant over time.
+#In default model, solar zenith angle is held constant over time.
 #This flag is kept to test b4b with the default code(this must be REMOVED when we decide to make it default)
-add_default($nl,'sol_insolation_bug_fix','val'=>'.false.');
+add_default($nl,'use_rad_dt_cosz','val'=>'.false.');
 
 
 #

models/atm/cam/bld/namelist_files/namelist_defaults_cam.xml

@@ -1095,8 +1095,8 @@
 <fix_g1_err_ndrop>             .false.  </fix_g1_err_ndrop>
 
 <!-- BSINGH - Flag (when .TRUE.)calculates solar zenith angle averaged over a time step.
-     In current default method solar zenith angle is held constant over time -->
-<sol_insolation_bug_fix>             .false.  </sol_insolation_bug_fix>
+     In default model solar zenith angle is held constant over time -->
+<use_rad_dt_cosz>             .false.  </use_rad_dt_cosz>
 
 
 

models/atm/cam/bld/namelist_files/namelist_definition.xml

@@ -3137,15 +3137,15 @@ Return fluxes per band in addition to the total fluxes.
 Default: FALSE
 </entry>
 
-<entry id="sol_insolation_bug_fix" type="logical"   category="radiation"
+<entry id="use_rad_dt_cosz" type="logical"  category="radiation"
        group="cam_inparm" valid_values="" >
+If true, then use the radiation timestep for all cosz calculations.
 NOTE: This logical flag is *temporary*. It should be REMOVED after we decide
-to call subroutine avg_slr_insol by default from shr_orb_mod.F90.
-When .TRUE. the model calculate solar zenith angle averaged over a time step.
-In current default method solar zenith angle is held constant over time.
-Default: .FALSE.
+to make this functionality default
+Default: FALSE
 </entry>
 
+
 <entry id="mode_defs" type="char*256(60)"  category="radiation"
        group="rad_cnst_nl" valid_values="" >
 Definitions for the aerosol modes that may be used in the rad_climate and

models/atm/cam/src/chemistry/mozart/mo_gas_phase_chemdr.F90

@@ -215,6 +215,7 @@ subroutine gas_phase_chemdr(lchnk, ncol, imozart, q, &
     use infnan,            only : nan, assignment(=)
     use rate_diags,        only : rate_diags_calc
     use mo_mass_xforms,    only : mmr2vmr, vmr2mmr, h2o_to_vmr, mmr2vmri
+    use orbit,             only : zenith
 !
 ! LINOZ
 !

models/atm/cam/src/chemistry/mozart/mo_waccm_hrates.F90

@@ -117,6 +117,7 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf )
       use short_lived_species,only: get_short_lived_species
       use physics_buffer,    only : physics_buffer_desc
       use phys_control,      only : waccmx_is
+      use orbit,             only : zenith
 
 !-----------------------------------------------------------------------
 !        ... dummy arguments

models/atm/cam/src/control/runtime_opts.F90

@@ -277,7 +277,7 @@ module runtime_opts
 logical :: spectralflux  ! calculate fluxes (up and down) per band. Default: FALSE
 
 !BSINGH - Flag to add solar insolation bug fix in /models/csm_share/shr/shr_orb_mod.F90
-logical :: sol_insolation_bug_fix  ! **temporary** flag for fixing solar insolation bug 
+logical :: use_rad_dt_cosz  ! if true, uses the radiation dt for all cosz calculations
 
 ! SCM Options
 logical  :: single_column
@@ -380,7 +380,6 @@ subroutine read_namelist(single_column_in, scmlon_in, scmlat_in, nlfilename_in )
 #if ( defined OFFLINE_DYN )
    use metdata,             only: metdata_readnl
 #endif
-   use shr_orb_mod,         only: shr_orb_mod_init
 
 !---------------------------Arguments-----------------------------------
 
@@ -472,7 +471,7 @@ subroutine read_namelist(single_column_in, scmlon_in, scmlat_in, nlfilename_in )
   namelist /cam_inparm/ print_energy_errors
 
   ! radiative heating calculation options
-  namelist /cam_inparm/ iradsw, iradlw, iradae, irad_always, spectralflux, sol_insolation_bug_fix
+  namelist /cam_inparm/ iradsw, iradlw, iradae, irad_always, spectralflux, use_rad_dt_cosz
 
   ! scam
   namelist /cam_inparm/ iopfile,scm_iop_srf_prop,scm_relaxation, &
@@ -513,7 +512,8 @@ subroutine read_namelist(single_column_in, scmlon_in, scmlat_in, nlfilename_in )
       iradlw_out      = iradlw,     &
       iradae_out      = iradae,     &
       irad_always_out = irad_always, &
-      spectralflux_out = spectralflux )
+      spectralflux_out = spectralflux,&
+      use_rad_dt_cosz_out = use_rad_dt_cosz )
 
    if (present(single_column_in)) then
       call scam_default_opts(scmlat_out=scmlat,scmlon_out=scmlon, &
@@ -696,21 +696,13 @@ subroutine read_namelist(single_column_in, scmlon_in, scmlat_in, nlfilename_in )
    call check_energy_setopts( &
       print_energy_errors_in = print_energy_errors )
 
-   !BSINGH(04/24/2015): call shr_orb_mod_init to set sol_insolation_bug_fix.
-   !This is a temporary variable which should be removed after we decide
-   !to make this fix defualt in the code
-   !This call also initializes radiation time step using 'dtime' and iradsw'
-   !Inserted here to avoid reading these namelist variables again in 
-   !/models/csm_share/shr/avg_slr_insol.F90
-   call shr_orb_mod_init(sol_insolation_bug_fix, dtime, iradsw)
-
    call radiation_setopts( dtime, nhtfrq(1), &
       iradsw_in      = iradsw,     &
       iradlw_in      = iradlw,     &
       iradae_in      = iradae,     &
       irad_always_in = irad_always, &
-      spectralflux_in = spectralflux )
-
+      spectralflux_in = spectralflux,&
+      use_rad_dt_cosz_in = use_rad_dt_cosz )
 ! 
 ! Set runtime options for single column mode
 !
@@ -969,7 +961,7 @@ subroutine distnl
    call mpibcast (spectralflux,1, mpilog, 0, mpicom)
 
    !BSINGH
-   call mpibcast (sol_insolation_bug_fix,1, mpilog, 0, mpicom)
+    call mpibcast (use_rad_dt_cosz,1,mpilog,0,mpicom)
 
 end subroutine distnl
 #endif

models/atm/cam/src/dynamics/se/share/cuda_mod.F90

@@ -24,8 +24,10 @@ module cuda_mod
 #if USE_CUDA_FORTRAN
 
 ! NP > 4 is not supported due to shared memory constraints
+!BSINGH -  NAG compiler doesn't support "#error" CPP directive.
+!Therefore, getting rid of "#" preceding "error" so that the code still fails at compile time
 #if NP > 4
-#error CUDA Fortran build only supported with NP <= 4
+error CUDA Fortran build only supported with NP <= 4
 #endif
 
 #define PAD 1

models/atm/cam/src/dynamics/se/share/spelt_mod.F90

@@ -1160,7 +1160,7 @@ subroutine spelt_runair(elem,spelt,hybrid,deriv,tstep,tl,nets,nete)
              endif
              if (mod(j,2)==1) then            ! works only for nip=3!!!
                fluxval(i,j,2) =  dt6 * sga*(contrauv(i,j,2,1)* slval_air(i,j,1) + & 
-                        4.0D0 * contrauv(i,j,2,2)*slval_air(i,j,2) + contrauv(i,j,2,3)*slval_air(i,j,3) )                                
+                        4.0D0 * contrauv(i,j,2,2)*slval_air(i,j,2) + contrauv(i,j,2,3)*slval_air(i,j,3) )
              endif
           end do
         end do 
@@ -1248,7 +1248,7 @@ subroutine spelt_runair(elem,spelt,hybrid,deriv,tstep,tl,nets,nete)
                                            flux(2) - &
                                            flux(3) - &
                                            flux(4) ) / (spelt(ie)%area_sphere(icell,jcell))
-              spelt(ie)%c(i+1,j+1,k,itr,tl%np1) = spelt(ie)%c(i+1,j+1,k,itr,tl%np1)/spelt(ie)%c(i+1,j+1,k,1,tl%np1)                             
+              spelt(ie)%c(i+1,j+1,k,itr,tl%np1) = spelt(ie)%c(i+1,j+1,k,itr,tl%np1)/spelt(ie)%c(i+1,j+1,k,1,tl%np1)
             end do
           end do
         end do      
@@ -1374,7 +1374,7 @@ subroutine spelt_runair_old(elem,spelt,hybrid,deriv,tstep,tl,nets,nete)
             endif
             if (mod(j,2)==1) then            ! works only for nip=3!!!
                fluxval(i,j,2) =  dt6 * sga*(contrauv(i,j,2,1)* slvalone(i,j,1) + & 
-                        4.0D0 * contrauv(i,j,2,2)*slvalone(i,j,2) + contrauv(i,j,2,3)*slvalone(i,j,3) )                                
+                        4.0D0 * contrauv(i,j,2,2)*slvalone(i,j,2) + contrauv(i,j,2,3)*slvalone(i,j,3) )
             endif
           end do
         end do
@@ -1515,7 +1515,8 @@ subroutine spelt_runair_old(elem,spelt,hybrid,deriv,tstep,tl,nets,nete)
 !                  write(*,*) flux
 !                endif                             
 !                 spelt(ie)%c(i+1,j+1,k,itr,tl%np1)=spelt(ie)%c(i+1,j+1,k,itr,tl%np1)/spelt(ie)%c(i+1,j+1,k,1,tl%np1)
-!                 if (isnan(spelt(ie)%c(i+1,j+1,k,itr,tl%np1))) write(*,*) 'ERROR2',spelt(ie)%c(i+1,j+1,k,itr,tl%np1),spelt(ie)%c(i+1,j+1,k,1,tl%np1)
+!                if (isnan(spelt(ie)%c(i+1,j+1,k,itr,tl%np1))) write(*,*) &
+!                     'ERROR2',spelt(ie)%c(i+1,j+1,k,itr,tl%np1),spelt(ie)%c(i+1,j+1,k,1,tl%np1)
 !                 
             end do
           end do

models/atm/cam/src/dynamics/se/share/viscosity_mod.F90

@@ -387,7 +387,7 @@ subroutine biharmonic_wk_scalar_minmax(elem,qtens,deriv,edgeq,hybrid,nets,nete,e
 real (kind=real_kind), dimension(np,np,nlev,qsize,nets:nete) :: qtens
 type (EdgeBuffer_t)  , intent(inout) :: edgeq
 type (derivative_t)  , intent(in) :: deriv
-real (kind=real_kind), intent(out), dimension(nlev,qsize,nets:nete) :: emin,emax
+real (kind=real_kind), intent(inout), dimension(nlev,qsize,nets:nete) :: emin,emax !BSINGH: emin and emax should be intent-inout
 
 ! local
 integer :: k,kptr,i,j,ie,ic,q

models/atm/cam/src/physics/cam/radiation.F90

@@ -63,6 +63,7 @@ module radiation
                            ! run continuously from the start of an
                            ! initial or restart run
 logical :: spectralflux  = .false. ! calculate fluxes (up and down) per band.
+logical :: use_rad_dt_cosz  = .false. ! if true, uses the radiation dt for all cosz calculations
 
 !Physics buffer indices
 integer :: qrs_idx      = 0
@@ -78,6 +79,8 @@ module radiation
 integer :: nmxrgn_idx = -1
 integer :: pmxrgn_idx = -1
 
+real(r8) :: dt_avg=0  ! time step to use for the shr_orb_cosz calculation, if use_rad_dt_cosz set to true
+
 !===============================================================================
 contains
 !===============================================================================
@@ -98,7 +101,7 @@ end subroutine radiation_register
 
 !================================================================================================
 
-subroutine radiation_defaultopts(iradsw_out, iradlw_out, iradae_out, irad_always_out, spectralflux_out)
+subroutine radiation_defaultopts(iradsw_out, iradlw_out, iradae_out, irad_always_out, spectralflux_out, use_rad_dt_cosz_out)
 !----------------------------------------------------------------------- 
 ! Purpose: Return default runtime options
 !-----------------------------------------------------------------------
@@ -108,20 +111,22 @@ subroutine radiation_defaultopts(iradsw_out, iradlw_out, iradae_out, irad_always
    integer, intent(out), optional :: iradae_out
    integer, intent(out), optional :: irad_always_out
    logical, intent(out), optional :: spectralflux_out
+   logical, intent(out), optional :: use_rad_dt_cosz_out
    !-----------------------------------------------------------------------
 
    if ( present(iradsw_out) )      iradsw_out = iradsw
    if ( present(iradlw_out) )      iradlw_out = iradlw
    if ( present(iradae_out) )      iradae_out = iradae
    if ( present(irad_always_out) ) irad_always_out = irad_always
    if ( present(spectralflux_out) ) spectralflux_out = spectralflux
+   if ( present(use_rad_dt_cosz_out) ) use_rad_dt_cosz_out = use_rad_dt_cosz
 
 end subroutine radiation_defaultopts
 
 !================================================================================================
 
 subroutine radiation_setopts(dtime, nhtfrq, iradsw_in, iradlw_in, iradae_in, &
-   irad_always_in, spectralflux_in)
+   irad_always_in, spectralflux_in, use_rad_dt_cosz_in)
 !----------------------------------------------------------------------- 
 ! Purpose: Set runtime options
 ! *** NOTE *** This routine needs information about dtime (init by dycore) 
@@ -137,7 +142,7 @@ subroutine radiation_setopts(dtime, nhtfrq, iradsw_in, iradlw_in, iradae_in, &
    integer, intent(in), optional :: iradae_in
    integer, intent(in), optional :: irad_always_in
    logical, intent(in), optional :: spectralflux_in
-
+   logical, intent(in), optional :: use_rad_dt_cosz_in
    ! Local
    integer :: ntspdy   ! no. timesteps per day
    integer :: nhtfrq1  ! local copy of input arg nhtfrq
@@ -148,6 +153,7 @@ subroutine radiation_setopts(dtime, nhtfrq, iradsw_in, iradlw_in, iradae_in, &
    if ( present(iradae_in) )      iradae = iradae_in
    if ( present(irad_always_in) ) irad_always = irad_always_in
    if ( present(spectralflux_in) ) spectralflux = spectralflux_in
+   if ( present(use_rad_dt_cosz_in) ) use_rad_dt_cosz = use_rad_dt_cosz_in
 
    ! Convert iradsw, iradlw and irad_always from hours to timesteps if necessary
    if (iradsw      < 0) iradsw      = nint((-iradsw     *3600._r8)/dtime)
@@ -346,6 +352,12 @@ subroutine radiation_init()
     call radlw_init(gravit, stebol)
     call radae_init(gravit, epsilo, stebol, pstd, mwdry, mwco2, mwo3)
 
+    ! Set the radiation timestep for cosz calculations if requested using the adjusted iradsw value from radiation
+    if (use_rad_dt_cosz)  then
+       dtime  = get_step_size()
+       dt_avg = iradsw*dtime
+    end if
+
     ! Get physics buffer indices
     cld_idx    = pbuf_get_index('CLD')
     concld_idx = pbuf_get_index('CONCLD')
@@ -548,6 +560,7 @@ subroutine radiation_tend(state,ptend, pbuf, &
     use interpolate_data, only: vertinterp
     use radiation_data,   only: output_rad_data
     use cloud_cover_diags,only: cloud_cover_diags_out
+    use orbit,            only: zenith
 
 
     ! Arguments
@@ -703,7 +716,7 @@ subroutine radiation_tend(state,ptend, pbuf, &
     ! Cosine solar zenith angle for current time step
     call get_rlat_all_p(lchnk, ncol, clat)
     call get_rlon_all_p(lchnk, ncol, clon)
-    call zenith (calday, clat, clon, coszrs, ncol)
+    call zenith (calday, clat, clon, coszrs, ncol, dt_avg)
 
     call output_rad_data(  pbuf, state, cam_in, landm, coszrs )