Misc fixes identified by Claude (batch 2)

src/aerosol/radiative_aerosol.F90

@@ -271,7 +271,7 @@
 
    ! bin name
    if (present(bin_name)) then
-      bin_name = bins%names(m_idx)
+      bin_name = bins%names(mm)
    end if
 
    ! name of interstitial number mixing ratio
@@ -459,7 +459,7 @@
             if (present(mode_idx)) mode_idx = i
             if (present(spec_idx)) spec_idx = ispec
             found_spectype = .true.
             exit
          end if
       end do
 
@@ -574,8 +574,8 @@
 
 !================================================================================================
 
 integer function rad_aer_get_idx(list_idx, aer_name)
    use string_utils,   only: to_str
    use cam_abortutils, only: endrun
    use radiative_aerosol_definitions, only: N_DIAG, aerlist_t, bulk_aerosol_list
 
@@ -1328,8 +1328,8 @@
 !------------------------------------------------------------------------
 ! Return the physprop ID for a mode in the modal aerosol list
 !------------------------------------------------------------------------
 integer function rad_aer_mode_physprop_id(list_idx, mode_idx)
    use string_utils,   only: to_str
    use cam_abortutils, only: endrun
    use radiative_aerosol_definitions, only: N_DIAG, modelist_t, modal_aerosol_list
 
@@ -1357,8 +1357,8 @@
 !------------------------------------------------------------------------
 ! Return the physprop ID for an aerosol in the bulk aerosol list
 !------------------------------------------------------------------------
 integer function rad_aer_bulk_physprop_id(list_idx, aer_idx)
    use string_utils,   only: to_str
    use cam_abortutils, only: endrun
    use radiative_aerosol_definitions, only: N_DIAG, aerlist_t, bulk_aerosol_list
 
@@ -1386,8 +1386,8 @@
 !------------------------------------------------------------------------
 ! Return the physprop ID for a bin in the sectional (CARMA) aerosol list
 !------------------------------------------------------------------------
 integer function rad_aer_bin_physprop_id(list_idx, bin_idx)
    use string_utils,   only: to_str
    use cam_abortutils, only: endrun
    use radiative_aerosol_definitions, only: N_DIAG, binlist_t, sectional_aerosol_list
 

src/control/cam_comp.F90

@@ -174,7 +174,8 @@ subroutine cam_init(caseid, ctitle, model_doi_url,                         &
            ref_tod, stop_ymd, stop_tod, curr_ymd, curr_tod,                   &
            perpetual_run, perpetual_ymd, initial_run_in)
 
-      dtime_phys = 0.0_r8
+      ! dtime_phys is seeded in phys_init, after the registry allocation
+      ! resets it; mark it here so the IC read does not try to fill it.
       call mark_as_initialized('timestep_for_physics')
 
       is_first_timestep = .true.

src/cpl/nuopc/atm_comp_nuopc.F90

@@ -417,6 +417,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc)
        call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO)
     end if
 
+    call shr_log_getLogUnit(shrlogunit)
     call shr_log_setLogUnit(iulog)
 
     !----------------------------------------------------------------------------
@@ -738,7 +739,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc)
           allocate(ownedElemCoords(spatialDim*numOwnedElements), stat=ierr, errmsg=tempc1)
           call check_allocate(ierr, subname, 'ownedElemCoords(spatialDim*numOwnedElements)', &
                               file=__FILE__, line=__LINE__, errmsg=tempc1)
-          call ESMF_MeshGet(model_mesh, ownedElemCoords=ownedElemCoords)
+          call ESMF_MeshGet(model_mesh, ownedElemCoords=ownedElemCoords, rc=rc)
           if (ChkErr(rc,__LINE__,u_FILE_u)) return
           do n = 1,lsize
              lonMesh(n) = ownedElemCoords(2*n-1)
@@ -916,6 +917,7 @@ subroutine DataInitialize(gcomp, rc)
           call ESMF_LogWrite('CAM - Initialize-Data-Dependency Returning to mediator without doing tphysbc', &
                ESMF_LOGMSG_INFO, rc=rc)
           if (ChkErr(rc,__LINE__,u_FILE_u)) return
+          call shr_log_setLogUnit (shrlogunit)
           return
        end if
 
@@ -1103,7 +1105,7 @@ subroutine ModelAdvance(gcomp, rc)
     ! Determine current time
     !--------------------------------
 
-    call ESMF_ClockGet(clock, currTime=currTime)
+    call ESMF_ClockGet(clock, currTime=currTime, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
     call ESMF_ClockGetNextTime(clock, nextTime=nextTime, rc=rc)
@@ -1262,7 +1264,7 @@ subroutine ModelAdvance(gcomp, rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
     call shr_cal_ymd2date(yr_sync, mon_sync, day_sync, ymd_sync)
 
-    if ((ymd /= ymd_sync) .and. (tod /= tod_sync))then
+    if ((ymd /= ymd_sync) .or. (tod /= tod_sync))then
        write(iulog,*)' cam ymd=',ymd     ,'  cam tod= ',tod
        write(iulog,*)'sync ymd=',ymd_sync,' sync tod= ',tod_sync
        call shr_sys_abort(subname//': CAM clock is not in sync with master Sync Clock')
@@ -1440,7 +1442,7 @@ subroutine ModelFinalize(gcomp, rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
     ! Determine if time to write restart
-    call ESMF_ClockGet(clock, currTime=currTime)
+    call ESMF_ClockGet(clock, currTime=currTime, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
     call ESMF_ClockGetNextTime(clock, nextTime=nextTime, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
@@ -2002,7 +2004,7 @@ subroutine cam_write_srfrest(gcomp, yr_spec, mon_spec, day_spec, sec_spec, rc)
           else if (lrank == 2) then
 
              ! There is an output variable for each element of the undistributed dimension
-             call ESMF_FieldGet(lfield, ungriddedUBound=ungriddedUBound, rc=rc)
+             call ESMF_FieldGet(lfield, ungriddedUBound=ungriddedUBound, gridToFieldMap=gridToFieldMap, rc=rc)
              if (chkerr(rc,__LINE__,u_FILE_u)) return
              call ESMF_FieldGet(lfield, farrayPtr=fldptr2d, rc=rc)
              if (ChkErr(rc,__LINE__,u_FILE_u)) return

src/cpl/nuopc/atm_stream_ndep.F90

@@ -80,11 +80,11 @@
     stream_ndep_year_align    = 1 ! align stream_ndep_year_first with this model year
 
     ! For now variable list in stream data file is hard-wired
     stream_varlist_ndep = (/'NDEP_NHx_month', 'NDEP_NOy_month'/)
 
     ! Read ndep_stream namelist
     if (masterproc) then
        open(newunit=nu_nml, file=trim(nlfile), status='old', iostat=nml_error, iomsg=nml_errmsg)
        if (nml_error /= 0) then
           call endrun(subName//': ERROR opening '//trim(nlfile)//' : '//&
                       errMsg(sourcefile, __LINE__)//' : '//trim(nml_errmsg))
@@ -168,7 +168,7 @@
          model_clock         = model_clock,                         &
          model_mesh          = model_mesh,                          &
          stream_meshfile     = trim(stream_ndep_mesh_filename),     &
          stream_filenames    = (/trim(stream_ndep_data_filename)/), &
          stream_yearFirst    = stream_ndep_year_first,              &
          stream_yearLast     = stream_ndep_year_last,               &
          stream_yearAlign    = stream_ndep_year_align,              &
@@ -211,6 +211,7 @@
     character(len=CS) :: ndepunits! ndep units
     !-----------------------------------------------------------------------
 
+    ndepunits = ''
     call cam_pio_openfile( File, trim(stream_fldFileName_ndep), PIO_NOWRITE)
     call pio_seterrorhandling(File, PIO_BCAST_ERROR, err_handling)
     ierr = pio_inq_varid(File, stream_varlist_ndep(1), vardesc)
@@ -219,6 +220,10 @@
             trim(stream_fldFileName_ndep)//errMsg(sourcefile, __LINE__))
     else
        ierr = PIO_get_att(File, vardesc, "units", ndepunits)
+       if (ierr /= PIO_NOERR) then
+          call endrun(' ERROR reading units of variable: '//trim(stream_varlist_ndep(1))//' in file: '// &
+               trim(stream_fldFileName_ndep)//errMsg(sourcefile, __LINE__))
+       end if
     end if
     call pio_seterrorhandling(File, err_handling)
     call cam_pio_closefile(File)

src/data/ref_pres.F90

@@ -15,8 +15,8 @@
    use cam_abortutils,  only: endrun
 
    implicit none
    public
    save
 
    ! Reference pressures (Pa)
    !> \section arg_table_ref_pres  Argument Table
@@ -198,7 +198,7 @@
       ! do_molec_diff
       call mark_as_initialized("do_molecular_diffusion")
       ! nbot_molec
-      call mark_as_initialized("index_of_pressure_at_bottom_of_molecular_diffusion")
+      call mark_as_initialized("vertical_layer_index_at_bottom_of_molecular_diffusion")
       ! gravity_wave_taper_bot_press
       call mark_as_initialized("largest_model_top_pressure_that_allows_tapering_gravity_wave_drag_at_model_top")
       ! nbot_gravity_wave_top_taper

src/dynamics/none/dyn_grid.F90

@@ -383,15 +383,28 @@ subroutine model_grid_init()
             call endrun(errormsg)
          end if
          if ((num_lats > 1) .and. (dimlens(1) == num_lats)) then
-            allocate(local_areas(num_lats), stat=iret, errmsg=errormsg)
-            call check_allocate(iret, subname, 'local_areas(num_lats)', &
+            ! Read the global area / weight variable and keep only this
+            ! task's latitudes, mirroring the latitude read above, so that
+            ! local_areas is indexed by the task-local latitude index.
+            start(1) = ((col_start - 1) / num_lons) + 1
+            kount(1) = ((col_end - 1) / num_lons) + 1 - start(1) + 1
+            allocate(local_areas(kount(1)), stat=iret, errmsg=errormsg)
+            call check_allocate(iret, subname, 'local_areas(kount)', &
                                 file=__FILE__, line=__LINE__, errmsg=errormsg)
-            start(1) = 1
-            kount(1) = num_lats
-            iret = pio_get_var(fh_ini, vardesc, start, kount, local_areas)
+            allocate(temp_arr(num_lats), stat=iret, errmsg=errormsg)
+            call check_allocate(iret, subname, 'temp_arr(num_lats)', &
+                                file=__FILE__, line=__LINE__, errmsg=errormsg)
+            iret = pio_get_var(fh_ini, vardesc, [1], [num_lats], temp_arr)
             call cam_pio_handle_error(iret,                                   &
                  subname//': Unable to read '//trim(var_name))
+            lindex = start(1) + kount(1) - 1
+            local_areas(1:kount(1)) = temp_arr(start(1):lindex)
+            deallocate(temp_arr)
          else if (dimlens(1) == num_global_columns) then
+            if (grid_is_latlon) then
+               ! iodesc is only created for the unstructured read above
+               call endrun(subname//': per-column grid areas on a lat/lon grid are not supported')
+            end if
             allocate(local_areas(num_local_columns), stat=iret, errmsg=errormsg)
             call check_allocate(iret, subname, 'local_areas(num_local_columns)', &
                                 file=__FILE__, line=__LINE__, errmsg=errormsg)
@@ -401,6 +414,7 @@ subroutine model_grid_init()
             write(errormsg, '(a,3(a,i0))') subname,                           &
                  'Unsupported number of grid areas, ', dimlens(1),            &
                  ', num_lats = ', num_lats, ', num_cols = ', num_global_columns
+            call endrun(errormsg)
          end if
       else
          call endrun(subname//'Unable to find grid areas')
@@ -473,7 +487,7 @@ subroutine set_dyn_col_values()
       do lindex = 1, num_local_columns
          if (grid_is_latlon) then
             lat_index = ((global_col_offset + lindex - 1) / num_lons) + 1 - lat1
-            lon_index = MOD(lindex - 1, num_lons) + 1
+            lon_index = MOD(global_col_offset + lindex - 1, num_lons) + 1
          else
             lat_index = lindex
             lon_index = lindex
@@ -615,7 +629,7 @@ subroutine find_dimension(file, dim_names, found_name, dim_len)
       character(len=128)               :: errmsg
       character(len=*), parameter      :: subname = 'find_dimension'
 
-      num_dims = len(dim_names)
+      num_dims = size(dim_names)
       found_name = ''
       do index = 1, num_dims
          ierr = PIO_Inq_DimID(file, trim(dim_names(index)), dimid)
@@ -631,7 +645,7 @@ subroutine find_dimension(file, dim_names, found_name, dim_len)
       if (len_trim(found_name) == 0) then
          write(errmsg, '(2a,20("  ",a))') subname,                            &
               ": Did not find any of these dimensions on initial data file:", &
-              (trim(dim_names(index)), index=1, len(dim_names))
+              (trim(dim_names(index)), index=1, size(dim_names))
          call endrun(errmsg)
       end if
    end subroutine find_dimension

src/dynamics/se/dp_coupling.F90

@@ -70,7 +70,6 @@ subroutine d_p_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_out)
 
    real(r8),  allocatable :: ps_tmp(:,:)         ! temp array to hold ps
    real(r8),  allocatable :: dp3d_tmp(:,:,:)     ! temp array to hold dp3d
-   real(r8),  allocatable :: dp3d_tmp_tmp(:,:)
    real(r8),  allocatable :: phis_tmp(:,:)       ! temp array to hold phis
    real(r8),  allocatable :: T_tmp(:,:,:)        ! temp array to hold T
    real(r8),  allocatable :: uv_tmp(:,:,:,:)     ! temp array to hold u and v
@@ -130,10 +129,6 @@ subroutine d_p_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_out)
    call check_allocate(ierr, subname, 'dp3d_tmp(nphys_pts,pver,nelemd)', &
                        file=__FILE__, line=__LINE__, errmsg=errmsg)
 
-   allocate(dp3d_tmp_tmp(nphys_pts,pver), stat=ierr, errmsg=errmsg)
-   call check_allocate(ierr, subname, 'dp3d_tmp_tmp(nphys_pts,pver)', &
-                       file=__FILE__, line=__LINE__, errmsg=errmsg)
-
    allocate(phis_tmp(nphys_pts,nelemd), stat=ierr, errmsg=errmsg)
    call check_allocate(ierr, subname, 'phis_tmp(nphys_pts,nelemd)', &
                        file=__FILE__, line=__LINE__, errmsg=errmsg)
@@ -676,7 +671,7 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
    call shr_vmath_log(phys_state%pintdry(1:pcols,1), &
                       phys_state%lnpintdry(1:pcols,1), pcols)
 
-   !$omp parallel do num_threads(horz_num_threads) private (k, i)
+   ! pintdry recurrence is loop-carried in k, do not parallelize:
    do k = 1, nlev
       do i = 1, pcols
          ! Calculate dry pressure variables for rest of column:
@@ -724,7 +719,8 @@ subroutine derived_phys_dry(cam_runtime_opts, phys_state, phys_tend)
       phys_state%pint(i,1) = phys_state%pintdry(i,1)
    end do
 
-   !$omp parallel do num_threads(horz_num_threads) private (k, i)
+   ! pint recurrence is loop-carried in k and ps accumulates over k,
+   ! do not parallelize:
    do k = 1, nlev
       do i=1, pcols
          ! Calculate wet (total) pressure variables for rest of column:

src/dynamics/se/stepon.F90

@@ -357,7 +357,7 @@ subroutine diag_dynvar_ic(elem, fvm)
    real(r8), allocatable :: fld_fvm(:,:,:,:,:), fld_gll(:,:,:,:,:)
    real(r8), allocatable :: fld_2d(:,:)
    logical,  allocatable :: llimiter(:)
-   real(r8)              :: qtmp(np,np,nlev), dp_ref(np,np,nlev), ps_ref(np,np)
+   real(r8)              :: dp_ref(np,np,nlev), ps_ref(np,np)
    real(r8), allocatable :: factor_array(:,:,:)
 
    character(len=*), parameter :: subname = 'diag_dynvar_ic'
@@ -375,8 +375,6 @@ subroutine diag_dynvar_ic(elem, fvm)
      tfname = trim(cnst_diag_name_gll(m_cnst))//'_gll'
      if (is_history_field_active(tfname)) then
        do ie = 1, nelemd
-         qtmp(:,:,:) =  elem(ie)%state%Qdp(:,:,:,m_cnst,tl_qdp)/&
-              elem(ie)%state%dp3d(:,:,:,tl_f)
          do j = 1, np
            do i = 1, np
              ftmp(i+(j-1)*np,:,1) = elem(ie)%state%Qdp(i,j,:,m_cnst,tl_qdp)/&

src/dynamics/tests/initial_conditions/ic_baroclinic.F90

@@ -173,7 +173,7 @@ subroutine bc_wav_set_ic(vcoord,latvals, lonvals, zint, U, V, T, PS, PHIS, &
       verbose_use = .true.
     end if
 
-    if(masterproc .and. verbose .and. present(PS)) then
+    if(masterproc .and. verbose_use .and. present(PS)) then
       write(iulog,*) subname, ': Model top (in km) is at z= ',ztop/1000.0_r8
     end if
 
@@ -673,7 +673,7 @@ END FUNCTION evaluate_exponential
   !-----------------------------------------------------------------------
   !    Stream function perturbation function
   !-----------------------------------------------------------------------
-  real(r8) FUNCTION evaluate_streamfunction(z,lon_local,lat_local)
+  real(r8) FUNCTION evaluate_streamfunction(lon_local,lat_local,z)
 
     real(r8), INTENT(IN)  :: lon_local
     real(r8), INTENT(IN)  :: lat_local

src/dynamics/tests/initial_conditions/ic_us_standard_atm.F90

@@ -134,6 +134,7 @@ subroutine us_std_atm_set_ic(latvals, lonvals, zint, U, V, T, PS, PHIS_IN, &
    end if
 
    if (present(T)) then
+      nlev = size(T, 2)
       allocate(pmid(nlev), stat=iret)
       call check_allocate(iret, subname, 'pmid(nlev)', &
                           file=__FILE__, line=__LINE__)
@@ -229,11 +230,11 @@ subroutine us_std_atm_set_ic(latvals, lonvals, zint, U, V, T, PS, PHIS_IN, &
       const_props => cam_model_const_properties()
 
       do m = 1, ncnst
-         if (m_cnst(m) == m_cnst_ix_q) then
+         if (m == m_cnst_ix_q) then
             ! No water vapor in profile
             do k = 1, nlev
                where(mask_use)
-                  Q(:,k,m_cnst(m)) = 0.0_r8
+                  Q(:,k,m) = 0.0_r8
                end where
             end do
             if(masterproc .and. verbose_use) then

src/dynamics/utils/dyn_thermo.F90

@@ -1398,14 +1398,14 @@ subroutine get_exner(tracer,mixing_ratio,active_species_idx,&
             deallocate(poverp0_phys)
          end if
 
+         !Set output variables back to dynamics kind:
+         exner = real(exner_phys, kind_dyn)
+
          !Deallocate variables:
          deallocate(tracer_phys)
          deallocate(dp_dry_phys)
          deallocate(exner_phys)
 
-         !Set output variables back to dynamics kind:
-         exner = real(exner_phys, kind_dyn)
-
       end if !kind check
 
    end subroutine get_exner

src/history/cam_hist_file.F90

@@ -697,7 +697,7 @@ subroutine config_set_beg_time(this, day, sec)
       class(hist_file_t), intent(inout) :: this
       integer, intent(in) :: day
       integer, intent(in) :: sec
-      integer, parameter :: seconds_per_day = 86400._r8
+      real(r8), parameter :: seconds_per_day = 86400._r8
 
       this%beg_time = day + (sec/seconds_per_day)
 
@@ -1559,6 +1559,10 @@ subroutine config_write_field(this, field, split_file_index, restart, &
          allocate(field_data(end_dims(1) - beg_dims(1) + 1, field_shape(2)), stat=ierr, errmsg=errmsg)
          call check_allocate(ierr, subname, 'field_data', file=__FILE__, line=__LINE__-1, errmsg=errmsg)
       end if
+      ! Points not touched by hist_field_norm_value (fields with zero
+      ! accumulated samples) must be written as fill, not as uninitialized
+      ! memory; this matches the _FillValue attribute set on the variable.
+      field_data = field%fill_value()
       ! Shape of array
       dimind = field%dimensions()
 
@@ -1857,6 +1861,12 @@ subroutine read_namelist_entry(unitn, hfile_config, hist_inst_fields,      &
            masterprocid, mpicom, ierr)
       call MPI_Bcast(hist_interp_type, flen, MPI_CHARACTER,                   &
            masterprocid, mpicom, ierr)
+      call MPI_Bcast(hist_interp_out, 1, MPI_LOGICAL, masterprocid,           &
+           mpicom, ierr)
+      call MPI_Bcast(hist_interp_nlat, 1, MPI_INTEGER, masterprocid,          &
+           mpicom, ierr)
+      call MPI_Bcast(hist_interp_nlon, 1, MPI_INTEGER, masterprocid,          &
+           mpicom, ierr)
       call MPI_Bcast(hist_filename_spec, CL, MPI_CHARACTER,                   &
            masterprocid, mpicom, ierr)
       ! Configure the history file

src/history/cam_history.F90

@@ -129,7 +129,7 @@ subroutine history_write_files()
                   write_history = .true.
                end if
             case('second')
-               nstep_freq = out_frq_mult / dtime
+               nstep_freq = nint(real(out_frq_mult, r8) / dtime)
                if (mod(nstep, nstep_freq) == 0) then
                   write_history = .true.
                end if
@@ -582,6 +582,10 @@ subroutine history_add_field_nd(diagnostic_name, standard_name, dimnames, avgfla
          num_levels = 1
       end if
 
+      ! cam_grid_get_array_bounds only fills mapped dimensions; default the
+      ! rest to a degenerate extent of 1 so no undefined values are stored.
+      dimbounds(:,1) = 1
+      dimbounds(:,2) = 1
       call cam_grid_get_array_bounds(grid_decomp, dimbounds)
 
       call cam_grid_dimensions(grid_decomp, grid_dims, rank)

src/physics/utils/gravity_wave_drag_ridge_read.F90

@@ -132,7 +132,9 @@ subroutine gravity_wave_drag_ridge_read_file()
 
     ! Local variables
     type(file_desc_t), pointer    :: fh_topo
-    type(file_desc_t), pointer    :: fh_rdggm
+    ! fh_rdggm is a locally-owned file handle (unlike fh_topo, which is
+    ! borrowed from cam_initfiles), so it is not a pointer.
+    type(file_desc_t)             :: fh_rdggm
     integer                       :: errflg
     character(len=512)            :: errmsg
     character(len=*), parameter   :: subname = 'gravity_wave_drag_ridge_read_file'
@@ -144,7 +146,6 @@ subroutine gravity_wave_drag_ridge_read_file()
     errflg = 0
 
     nullify(fh_topo)
-    nullify(fh_rdggm)
 
     has_gbxar_from_topo = .false.
     call mark_as_initialized('number_of_ridges_in_ridge_gravity_wave_drag')
@@ -301,6 +302,9 @@ subroutine gravity_wave_drag_ridge_read_file()
         ! Convert from m2 to km2
         rdg_gbxarg = rdg_gbxarg * (rearth/1000._kind_phys) * (rearth/1000._kind_phys)
       else
+        ! rdg_gbxar was already read from the topo file and converted to
+        ! km2 above; copy it so the gamma grid box area is not left zero.
+        rdg_gbxarg(:) = rdg_gbxar(:)
         if(masterproc) then
           write(iulog,*) trim(subname) // ': Using GBXAR from topo file, skipping gamma file GBXAR'
         end if
@@ -347,8 +351,6 @@ subroutine gravity_wave_drag_ridge_read_file()
       end if
 
       call cam_pio_closefile(fh_rdggm)
-      deallocate(fh_rdggm)
-      nullify(fh_rdggm)
 
       ! Mark variables as initialized so they are not read from ic file.
       call mark_as_initialized('grid_box_area_for_gamma_ridge_gravity_wave_drag')

src/physics/utils/orbital_data.F90

@@ -85,8 +85,10 @@ subroutine orbital_data_advance(calendar_day, latitudes, longitudes, use_rad_uni
 
     ! Compute the solar zenith angle [radians]
     do i = 1, size(latitudes)
-      solar_zenith_angle(i) = acos(shr_orb_cosz(calendar_day, latitudes(i), &
-                                                longitudes(i), solar_declination))
+      ! Clamp to [-1,1] so roundoff in the cosine cannot make acos return NaN
+      solar_zenith_angle(i) = acos(max(-1.0_kind_phys, min(1.0_kind_phys, &
+                                   shr_orb_cosz(calendar_day, latitudes(i), &
+                                                longitudes(i), solar_declination))))
     end do
 
     ! Compute the cosine of solar zenith angle for radiation [radians]