Standalone rad

individual-components/am4.f90

@@ -693,7 +693,7 @@ subroutine override_data_3d(name, buffer, time_level, domain)
   real, dimension(:, :, :), allocatable :: override_buffer
   integer, dimension(3) :: count_
   type(FmsNetcdfDomainFile_t) :: dataset
-  integer :: i
+  integer :: i, zmin, zmax
   integer, dimension(3) :: start
 
   do i = 1, size(override_variables)
@@ -711,6 +711,11 @@ subroutine override_data_3d(name, buffer, time_level, domain)
   if (override_z_lower .lt. 0 .or. override_z_upper .lt. 0) then
     call error_mesg("override_data_3d", "you must set the overrize z limits.", fatal)
   endif
+  zmin = max(override_z_lower, 1)
+  zmax = min(override_z_upper, size(buffer, 3))
+  if (zmin .gt. zmax) then
+    call error_mesg("override_data_3d", "override z lower limit must be less than the upper limit.", fatal)
+  endif
   if (.not. open_file(dataset, override_path, "read", domain)) then
     call error_mesg("override_data_3d", "cannot find "//trim(override_path)//".", fatal)
   endif
@@ -719,7 +724,7 @@ subroutine override_data_3d(name, buffer, time_level, domain)
   allocate(override_buffer(size(buffer, 1), size(buffer, 2), size(buffer, 3)))
   call read_data(dataset, name, override_buffer, unlim_dim_level=time_level)
   call close_file(dataset)
-  buffer(:, :, override_z_lower:override_z_upper) = override_buffer(:, :, override_z_lower:override_z_upper)
+  buffer(:, :, zmin:zmax) = override_buffer(:, :, zmin:zmax)
   deallocate(override_buffer)
 end subroutine override_data_3d
 

individual-components/radiation.f90

@@ -27,11 +27,12 @@ program main
 use solar_constant, only: SolarConstant
 use solar_spectrum, only: SolarSpectrum
 use time_interp_external2_mod, only: time_interp_external_init
-use time_manager_mod, only: get_date, julian, print_time, set_calendar_type, time_manager_init, &
-                            time_type, operator(+), operator(-)
+use time_manager_mod, only: get_date, julian, noleap, print_time, set_calendar_type, time_manager_init, &
+                            time_type, get_time, set_time, operator(+), operator(-), operator(/), set_date
 use tracer_manager_mod, only: get_number_tracers, get_tracer_index, &
                               tracer_manager_end, tracer_manager_init
 use utilities, only: catch_error, integrate
+use, intrinsic :: ieee_arithmetic, only: ieee_is_nan
 implicit none
 
 
@@ -76,10 +77,16 @@ program main
 real(kind=wp), dimension(:, :), allocatable :: swabs_integral
 integer :: t
 type(time_type) :: time
+type(time_type) :: forcing_time
 type(time_type) :: time_next
 type(time_type) :: timestep
 real :: top_level_pressure
-
+integer :: invalid_timestep
+integer, allocatable :: block_status(:)
+integer :: retry
+integer :: lookback, orig_t
+logical :: found_good_time
+integer :: max_lookback
 !MPP timers.
 integer :: aerosol_optics_clock
 integer :: cloud_optics_clock
@@ -100,16 +107,20 @@ program main
 type(block_control_type) :: column_blocking
 
 !Runtime options.
+integer :: forcing_year = -1 !Year of the forcing for solar, ghg, and aerosol
 integer :: near_infrared_cutoff = 14600 !Wavenumber [cm-1] that distinguishes the visible from the near-infrared.
 integer :: nxblocks = 1
 integer :: nyblocks = 1
+logical :: remove_last_timestep = .false.
 character(len=256) :: solar_constant_path = ""
 character(len=256) :: solar_spectrum_path = ""
-namelist /standalone_radiation_nml/ near_infrared_cutoff, &
+namelist /standalone_radiation_nml/ forcing_year, & 
+                                    near_infrared_cutoff, &
                                     nxblocks, &
                                     nyblocks, &
                                     solar_constant_path, &
-                                    solar_spectrum_path
+                                    solar_spectrum_path, &
+                                    remove_last_timestep
 
 !Start up mpi.
 call fms_init()
@@ -141,31 +152,22 @@ program main
 !Set the model time.
 if (trim(atm(1)%calendar) .eq. "julian") then
   call set_calendar_type(julian)
+elseif (trim(atm(1)%calendar) .eq. "noleap") then
+  call set_calendar_type(noleap)
 else
-  call error_mesg("main", "only julian calendar supported.", fatal)
+  call error_mesg("main", "only julian and noleap calendar supported.", fatal)
 endif
+time = get_cal_time(atm(1)%time(1), atm(1)%time_units, atm(1)%calendar)
+time = normalize_time(time)
 if (atm(1)%num_times .gt. 1) then
   dt = atm(1)%time(2) - atm(1)%time(1)
-  do t = 3, atm(1)%num_times
-    !Check to see if the input dataset has a constant timestep.
-    if (abs((atm(1)%time(t) - atm(1)%time(t - 1)) - dt) .gt. 1.e-10) then
-      call error_mesg("main", "timestep is not constant in the input datasets.", fatal)
-    endif
-  enddo
 else
-  !If the input dataset only has one time level, we can't determine the
-  !model radiatoin timestep that was used, so set it to zero.
   dt = 0.
 endif
 
-!Model diagnostics are output at the end of a timestep, yet calculated
-!using the time at the beginning of a timestep.  To mimic that, subtract
-!one timestep off of the first time found in the input dataset.
-time = get_cal_time(max(atm(1)%time(1) - dt, 0.), atm(1)%time_units, atm(1)%calendar)
-timestep = get_cal_time(atm(1)%time(1), atm(1)%time_units, atm(1)%calendar)
-timestep = timestep - time
-write(*, *) atm(1)%time, dt
-
+if (mpp_pe() .eq. mpp_root_pe()) then
+  call print_time(time, "Starting radiation timestep: ")
+endif
 
 !Read in the solar data.
 call solar_flux_constant%create("solar_flux", trim(solar_constant_path))
@@ -268,7 +270,6 @@ program main
   call aerosol_species_diags(i)%create(time, axes, radiation_context%aerosol_optics%family(i)%name, &
                                        radiation_context%aerosol_optics%family(i)%name)
 enddo
-call diag_manager_set_time_end(get_cal_time(atm(1)%time(atm(1)%num_times), atm(1)%time_units, atm(1)%calendar))
 
 !Calculate variables need for dealing with the land spectral decomposition.
 num_bands = size(shortwave_band_limits, 2)
@@ -300,33 +301,94 @@ program main
 endif
 deallocate(shortwave_band_limits)
 
+! Allocate block_status once for the run to avoid repeated allocation/race conditions
+allocate(block_status(num_blocks))
+
 !Main loop.
-do t = 1, atm(1)%num_times
+!The last timestep of offline input files might be invalid
+!Remove the last timestep if needed.
+if (remove_last_timestep) then
+        invalid_timestep = 1   
+else
+        invalid_timestep = 0
+endif
+do t = 1, atm(1)%num_times-invalid_timestep
   !Calculate the current time.
-  time_next = time + timestep
+  time = get_cal_time(atm(1)%time(t), atm(1)%time_units, atm(1)%calendar)
+  time = normalize_time(time)
+  time_next = get_cal_time(atm(1)%time(t) + dt, atm(1)%time_units, atm(1)%calendar)
+  time_next = normalize_time(time_next)
   call print_time(time, "Running timestep: ")
 
   !Read in the atmospheric properies.
   call read_time_slice(atm, t, column_blocking)
 
+  if (forcing_year .ge. 0) then
+      call get_date(time, date(1), date(2), date(3), date(4), date(5), date(6))
+      forcing_time = set_date(forcing_year, date(2), date(3), date(4), date(5), date(6))
+  else
+      forcing_time = time
+  endif
+
   !Time interpolation.
-  call solar_flux_constant%update(time)
-  call radiation_context%update(time)
-
-!$omp parallel do private(block_) default(shared)
-  do block_ = 1, num_blocks
-    call radiation_scheme(radiation_context, atm(block_), column_blocking, num_layers, block_, &
-                          aerosol_optics_clock, cloud_optics_clock, flux_solver_clock, &
-                          gas_optics_clock, radiation_driver_clock, h2o, o3, last_infrared_band, &
-                          aerosol_species_diags, cloud_diags, flux_diags, time, time_next, solar_flux_constant, &
-                          surface_albedo_weight, all_, clean, clean_clear, clear, &
-                          olr_integral, swabs_integral)
+  call solar_flux_constant%update(forcing_time)
+  call radiation_context%update(forcing_time)
+
+  ! Try the current timestep and, if it fails, search backwards to
+  ! the most recent earlier timestep that yields a successful radiation run.
+  orig_t = t
+  found_good_time = .false.
+  ! Limit lookback to at most 5 earlier timesteps to avoid excessive searching.
+  max_lookback = min(5, orig_t - 1)
+  do lookback = 0, max_lookback
+    if (lookback == 0) then
+      call read_time_slice(atm, orig_t, column_blocking)
+    else
+      call read_time_slice(atm, orig_t - lookback, column_blocking)
+    endif
+
+    block_status(:) = 0
+    !$omp parallel do private(block_) default(shared)
+    do block_ = 1, num_blocks
+      call radiation_scheme(radiation_context, atm(block_), column_blocking, num_layers, block_, &
+                            aerosol_optics_clock, cloud_optics_clock, flux_solver_clock, &
+                            gas_optics_clock, radiation_driver_clock, h2o, o3, last_infrared_band, &
+                            aerosol_species_diags, cloud_diags, flux_diags, time, time_next, solar_flux_constant, &
+                            surface_albedo_weight, all_, clean, clean_clear, clear, &
+                            olr_integral, swabs_integral, block_status(block_))
+    enddo
+    !$omp end parallel do
+
+    if (.not. any(block_status /= 0)) then
+      if (lookback > 0 .and. mpp_pe() .eq. mpp_root_pe()) then
+        write(logfile_handle, *) 'Warning: radiation failed at t=', orig_t, ' recovered using t=', orig_t - lookback
+      endif
+      found_good_time = .true.
+      exit
+    else
+      if (orig_t - lookback <= 1) then
+        ! No earlier time to try.
+        exit
+      endif
+      ! Otherwise, continue to next earlier timestep and try again.
+    endif
   enddo
 
+  if (.not. found_good_time) then
+    if (mpp_pe() .eq. mpp_root_pe()) then
+      write(logfile_handle, *) 'Error: radiation failed for timestep', orig_t, 'and all earlier times; aborting.'
+    endif
+    call error_mesg('main', 'radiation failed for all earlier timesteps; aborting.', fatal)
+  endif
+
   !Write out diagnostics.
-  call diag_send_complete(timestep)
+  call diag_manager_set_time_end(time)
+  call diag_send_complete(time)
   time = time_next
 enddo
+  if (allocated(block_status)) then
+    deallocate(block_status)
+  endif
 
 !Clean up.
 deallocate(ak, bk)
@@ -352,13 +414,32 @@ program main
 
 contains
 
+function normalize_time(time)
+  type(time_type), intent(in) :: time
+  type(time_type) :: normalize_time
+  integer :: seconds, days, ticks, residual
+  integer, parameter :: time_step = 1800
+
+  call get_time(time, seconds, days, ticks)
+  residual = mod(seconds, time_step)
+  if (residual .ge. time_step / 2) then
+    seconds = seconds + time_step - residual
+    if (seconds .ge. 86400) then
+      seconds = 0
+      days = days + 1
+    endif
+  else
+    seconds = seconds - residual
+  endif
+  normalize_time = set_time(seconds, days, ticks)
+end function normalize_time
 
 subroutine radiation_scheme(radiation_context, atm, column_blocking, num_layers, block_, &
                             aerosol_optics_clock, cloud_optics_clock, flux_solver_clock, &
                             gas_optics_clock, radiation_driver_clock, h2o, o3, last_infrared_band, &
                             aerosol_species_diags, cloud_diags, flux_diags, time, time_next, solar_flux_constant, &
                             surface_albedo_weight, all_, clean, clean_clear, clear, &
-                            olr_integral, swabs_integral)
+                            olr_integral, swabs_integral, status)
 
   type(RadiationContext), intent(inout) :: radiation_context
   type(Atmosphere_t), intent(in), target :: atm
@@ -386,6 +467,7 @@ subroutine radiation_scheme(radiation_context, atm, column_blocking, num_layers,
   integer, intent(in) :: clear
   real(kind=wp), dimension(:, :), intent(inout) :: olr_integral
   real(kind=wp), dimension(:, :), intent(inout) :: swabs_integral
+  integer, intent(out) :: status
 
   integer :: band, column, i, n, num_bands, num_columns, num_lat, num_levels, num_lon, s
   real, dimension(:, :, :), allocatable :: aerosol_relative_humidity
@@ -447,6 +529,33 @@ subroutine radiation_scheme(radiation_context, atm, column_blocking, num_layers,
   num_columns = num_lon*num_lat
   num_levels = num_layers + 1
 
+  ! Initialize status and perform quick NaN checks on essential fields.
+  status = 0
+  if (any(ieee_is_nan(atm%ppmv))) then
+    status = 1
+    return
+  endif
+  if (any(ieee_is_nan(atm%layer_temperature))) then
+    status = 1
+    return
+  endif
+  if (any(ieee_is_nan(atm%level_temperature))) then
+    status = 1
+    return
+  endif
+  if (any(ieee_is_nan(atm%level_pressure))) then
+    status = 1
+    return
+  endif
+  if (any(ieee_is_nan(atm%layer_pressure))) then
+    status = 1
+    return
+  endif
+  if (any(ieee_is_nan(atm%surface_temperature))) then
+    status = 1
+    return
+  endif
+
   !Get gpoint limits
   call radiation_context%longwave_gas_optics%gpoint_limits(longwave_gpoint_limits)
   call radiation_context%shortwave_gas_optics%gpoint_limits(shortwave_gpoint_limits)
@@ -486,12 +595,18 @@ subroutine radiation_scheme(radiation_context, atm, column_blocking, num_layers,
   allocate(surface_emissivity(num_bands, num_columns))
   surface_emissivity(:, :) = 1.
 
-  !Start the radiation timer.
-  call mpp_clock_begin(radiation_driver_clock)
-
   !Update gas concentrations.
   water_vapor(1:num_columns, 1:num_layers) => atm%ppmv(1:num_lon, 1:num_lat, 1:num_layers, h2o)
   ozone(1:num_columns, 1:num_layers) => atm%ppmv(1:num_lon, 1:num_lat, 1:num_layers, o3)
+  ! If ozone VMRs are out of the valid range [0,1], signal failure so
+  ! the driver can retry the timestep using the previous input.
+  if (any(ozone < 0.0_wp .or. ozone > 1.0_wp)) then
+    if (mpp_pe() .eq. mpp_root_pe()) then
+      write(logfile_handle, *) 'Warning: ozone values out of [0,1] found in block', block_, '; requesting retry with previous timestep'
+    endif
+    status = 1
+    return
+  endif
   call radiation_context%update_concentrations(water_vapor, ozone, block_)
 
   !Calculate gas optics.
@@ -501,6 +616,8 @@ subroutine radiation_scheme(radiation_context, atm, column_blocking, num_layers,
   level_temperature(1:num_columns, 1:num_levels) => atm%level_temperature(1:num_lon, 1:num_lat, 1:num_levels)
   surface_temperature(1:num_columns) => atm%surface_temperature(1:num_lon, 1:num_lat)
   zenith(1:num_columns) => atm%solar_zenith_angle(1:num_lon, 1:num_lat)
+  !Start the radiation timer.
+  call mpp_clock_begin(radiation_driver_clock)
   call mpp_clock_begin(gas_optics_clock)
   call radiation_context%calculate_gas_optics(layer_pressure, level_pressure, &
                                               layer_temperature, level_temperature, &