Addition of new features to generate ocean boundary conditions for atmosphere-only **forecasts**

CHANGELOG.md

@@ -52,6 +52,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [2.0.2] - 2023-06-28
 
+### Added
+
+- Addition of SST and FRACI forecast (ocean) boundary conditions generation capability in `pre/prepare_ocnExtData`
+
 ### Fixed
 
 - Updates to fix dates for comparison plots. Updates for plot levels.

pre/prepare_ocnExtData/CMakeLists.txt

@@ -2,6 +2,7 @@ esma_set_this()
 
 set ( srcs
   sst_ice_helpers.F90
+  hermes_m_tick.f90
   )
 
 esma_add_library (${this} SRCS ${srcs} DEPENDENCIES MAPL NetCDF::NetCDF_Fortran)
@@ -19,6 +20,24 @@ ecbuild_add_executable (
   LIBS ${this})
 set_target_properties(sst_sic_eight.x PROPERTIES Fortran_MODULE_DIRECTORY ${include_${this}})
 
+ecbuild_add_executable (
+  TARGET cal_daily_clim.x
+  SOURCES daily_clim_SST_FRACI_eight.F90
+  LIBS ${this})
+set_target_properties(cal_daily_clim.x PROPERTIES Fortran_MODULE_DIRECTORY ${include_${this}})
+
+ecbuild_add_executable (
+  TARGET extract_daily_sst_ice.x
+  SOURCES extract_day_SST_FRACI_eight.F90
+  LIBS ${this})
+set_target_properties(extract_daily_sst_ice.x PROPERTIES Fortran_MODULE_DIRECTORY ${include_${this}})
+
+ecbuild_add_executable (
+  TARGET gen_forecast_bcs.x
+  SOURCES forecast_SST_FRACI_eight.F90
+  LIBS ${this})
+set_target_properties(gen_forecast_bcs.x PROPERTIES Fortran_MODULE_DIRECTORY ${include_${this}})
+
 if (USE_F2PY)
    find_package(F2PY2)
    if (F2PY2_FOUND)

pre/prepare_ocnExtData/daily_clim_SST_FRACI_eight.F90

@@ -0,0 +1,171 @@
+PROGRAM daily_clim_SST_FRACI_eight
+!---------------------------------------------------------------------------
+  USE m_tick, only: INCYMD
+  USE sst_ice_helpers, only: read_bin_SST_ICE, &
+                             write_bin,        &
+                             write_netcdf
+
+  IMPLICIT NONE
+
+  INTEGER                 :: month, day
+  INTEGER                 :: clim_year, year_start, year_end
+  CHARACTER (LEN = 200)   :: data_path
+  CHARACTER (LEN = 200)   :: sst_file_pref, ice_file_pref
+  CHARACTER (LEN = 200)   :: sst_file_suff, ice_file_suff
+
+  INTEGER                 :: nymd_in, nYears
+  INTEGER                 :: nymd_out, NLON, NLAT
+  REAL                    :: LON(2880), LAT(1440)
+
+  INTEGER                 :: year
+  INTEGER                 :: clim_tom, clim_tom_year, clim_tom_mon, clim_tom_day
+  CHARACTER (LEN = 200)   :: year_str, sst_file, ice_file
+  CHARACTER (LEN = 8)     :: clim_date_str
+  REAL                    :: mean_sst(1440, 2880), mean_ice(1440, 2880)
+  REAL                    :: sdev_sst(1440, 2880), sdev_ice(1440, 2880)
+  REAL                    ::      sst(1440, 2880),      ice(1440, 2880)
+  REAL                    ::     dsst(1440, 2880),     dice(1440, 2880)
+
+  LOGICAL                 :: verbose
+
+  character (len=4) :: arg
+  integer :: numArgs
+  numArgs = iargc()
+
+  verbose = .false.
+
+! Inputs
+  if (numArgs < 4)  then
+    print *, "Need 4 inputs, in following format. Try again!"
+    print *, "start_year end_year month day. Last 2 inputs are for the climatology of MM/DD."
+    STOP
+  else ! read user inputs
+    call getarg(1, arg)
+    read(arg, '(I4)') year_start !2007
+    call getarg(2, arg)
+    read(arg, '(I4)') year_end   !2023
+    call getarg(3, arg)
+    read(arg, '(I2)') month
+    call getarg(4, arg)
+    read(arg, '(I2)') day
+  end if
+
+! - Set inputs
+  clim_year  = 1 ! Set, can't be "0" since GEOS/ESMF does not like that!
+! GMAO OPS has following fixed
+  data_path= '/discover/nobackup/projects/gmao/share/dao_ops/fvInput/g5gcm/bcs/realtime/OSTIA_REYNOLDS/2880x1440/'
+  sst_file_pref = 'dataoceanfile_OSTIA_REYNOLDS_SST.2880x1440.'
+  ice_file_pref = 'dataoceanfile_OSTIA_REYNOLDS_ICE.2880x1440.'
+  sst_file_suff = '.data'
+  ice_file_suff = '.data'
+
+! Check input range
+  if (year_start < 2007) then
+    print *, "Invalid starting year: ", year_start, "Must be 2007 or later, try again."
+    STOP
+  endif
+  if ( (month < 1) .or. (month > 12)) then
+    print *, "Invalid month: ", month, "Try again."
+    STOP
+  endif
+  if ( (day < 1) .or. (day > 31)) then
+    print *, "Invalid day: ", day, "Try again."
+    STOP
+  endif
+  if(verbose) print *, "User inputs: ", year_start, year_end, month, day
+!--
+
+  print *, " "
+! mean and std dev over year_start: year_end
+! 1. mean
+  nYears = 0
+  do year = year_start, year_end
+    write (year_str,'(I4.0)') year
+    sst_file = trim(data_path)//"/"//trim(sst_file_pref)//trim(year_str)//trim(sst_file_suff)
+    ice_file = trim(data_path)//"/"//trim(ice_file_pref)//trim(year_str)//trim(ice_file_suff)
+
+    nymd_in = year*10000 + month*100 + day
+    if(verbose) then
+      print *, year, year_str, nymd_in
+      print *, sst_file, ice_file
+      print *, nymd_in
+    endif
+
+    CALL read_bin_SST_ICE( sst_file, nymd_in, nymd_out, NLON, NLAT, LON, LAT, sst) ! Read GMAO OPS dataset
+    CALL read_bin_SST_ICE( ice_file, nymd_in, nymd_out, NLON, NLAT, LON, LAT, ice)
+
+    if ( (nymd_out-nymd_in) .eq. 0) then ! Found data for matching dates
+      print *, "Read SST and FRACI for: ", nymd_out
+      if (verbose) print *, "Found data in file, number of years =", nYears
+      nYears = nYears + 1
+
+      if (nYears == 1) then
+        mean_sst = sst
+        mean_ice = ice
+      else
+        mean_sst = mean_sst + sst 
+        mean_ice = mean_ice + ice
+      endif
+    endif
+  end do
+
+  print *, " "
+  print *, "Number of years= ", nYears
+
+  mean_sst = mean_sst/REAL(nYears) ! nYears should be _promoted_ to real, being safe!
+  mean_ice = mean_ice/REAL(nYears)
+
+! 2. sdev
+  nYears = 0
+  do year = year_start, year_end
+
+    write (year_str,'(I4.0)') year
+    sst_file = trim(data_path)//"/"//trim(sst_file_pref)//trim(year_str)//trim(sst_file_suff)
+    ice_file = trim(data_path)//"/"//trim(ice_file_pref)//trim(year_str)//trim(ice_file_suff)
+
+    nymd_in = year*10000 + month*100 + day
+    CALL read_bin_SST_ICE( sst_file, nymd_in, nymd_out, NLON, NLAT, LON, LAT, sst)
+    CALL read_bin_SST_ICE( ice_file, nymd_in, nymd_out, NLON, NLAT, LON, LAT, ice)
+
+    if ( (nymd_out-nymd_in) .eq. 0) then
+      nYears = nYears + 1
+
+      if (nYears == 1) then
+        dsst = (sst - mean_sst)**2.
+        dice = (ice - mean_ice)**2.
+      else
+        dsst = dsst + (sst - mean_sst)**2.
+        dice = dice + (ice - mean_ice)**2.
+      endif
+    endif
+  end do
+
+  print *, "Number of years= ", nYears
+! based on a sample not population
+  sdev_sst = sqrt( dsst/REAL(nYears-1))
+  sdev_ice = sqrt( dice/REAL(nYears-1))
+! --
+
+  print *, " "
+  clim_tom = INCYMD( clim_year*10000+month*100+day, 1)
+  clim_tom_year = int(clim_tom/10000)
+  clim_tom_mon  = int((clim_tom-clim_tom_year*10000)/100)
+  clim_tom_day = clim_tom - (clim_tom_year*10000+clim_tom_mon*100)
+  if(verbose) print *, "Tomorrow day:", clim_tom_year, clim_tom_mon, clim_tom_day
+
+  write (clim_date_str,'(I0.8)') clim_year*10000+month*100+day
+  if(verbose) print *, clim_date_str
+
+! Only write mean in binary format- good enough.
+  print *, "Write output ..."
+  CALL write_bin(    'daily_clim_mean', clim_year, month, day, clim_tom_year, clim_tom_mon, clim_tom_day, clim_date_str, 1440, 2880, transpose(mean_sst), transpose(mean_ice))
+
+! mean
+  CALL write_netcdf( 'daily_clim_mean_sst_fraci', clim_year, month, day, clim_date_str, 1440, 2880, transpose(mean_sst), transpose(mean_ice))
+! std dev
+  CALL write_netcdf( 'daily_clim_sdev_sst_fraci', clim_year, month, day, clim_date_str, 1440, 2880, transpose(sdev_sst), transpose(sdev_ice))
+  print *, "Done."
+  print *, " "
+
+!---------------------------------------------------------------------------
+END PROGRAM daily_clim_SST_FRACI_eight

pre/prepare_ocnExtData/extract_day_SST_FRACI_eight.F90

@@ -0,0 +1,130 @@
+PROGRAM extract_day_SST_FRACI_eight
+!---------------------------------------------------------------------------
+  USE m_tick, only: INCYMD
+  USE sst_ice_helpers, only: read_bin_SST_ICE, &
+                             write_bin,        &
+                             write_netcdf
+
+  IMPLICIT NONE
+
+  integer, parameter :: nlon = 2880
+  integer, parameter :: nlat = 1440
+  real               :: lon(nlon), lat(nlat)
+  real               :: sst(nlat, nlon), ice(nlat, nlon)
+
+  integer               :: extract_date
+  integer               :: date_out, nlon_out, nlat_out
+
+  character (len = 200) :: data_path
+  character (len = 200) :: sst_file_pref, ice_file_pref
+  character (len = 200) :: sst_file_suff, ice_file_suff
+
+  character (len=8)     :: date_str, format_str = "(I8)"
+  character (len = 200) :: sst_file, ice_file
+
+  integer :: i, iarg, argc, iargc
+  character (len=255) :: argv
+  integer :: numArgs
+  integer :: year_s, month_s, day_s
+  integer :: tom, tom_year, tom_mon, tom_day
+
+  logical :: save_bin
+!----
+
+  argc = iargc()
+  if (argc < 3)  then
+    print *, "Need minmum 3 inputs, see below examples. Try again!"
+    print *, "extract_daily_sst_ice.x -year 2010 -month 1 -day 1"
+    print *, " "
+    print *, "** 2 additional optional inputs are allowed ** "
+    print *, "extract_daily_sst_ice.x -year 2010 -month 1 -day 1 -save_bin"
+    print *, " "
+    print *, "extract_daily_sst_ice.x -year 2010 -month 1 -day 1 -input_data_path xx"
+    print *, " "
+    print *, "extract_daily_sst_ice.x -year 2010 -month 1 -day 1 -save_bin -input_data_path xx"
+    print *, " "
+    STOP
+  end if
+
+  ! Following 2 variables have defaults that can be overridden by user inputs.
+  save_bin = .false.
+  ! GMAO OPS data path
+  data_path= '/discover/nobackup/projects/gmao/share/dao_ops/fvInput/g5gcm/bcs/realtime/OSTIA_REYNOLDS/2880x1440/'
+
+  !-- read user inputs
+  iarg = 0
+  do i = 1, 99
+    iarg = iarg + 1
+    if ( iarg > argc ) exit
+    call getarg(iarg, argv)
+
+    select case (argv)
+      case ("-year")
+        iarg = iarg + 1
+        call getarg ( iarg, argv )
+        read(argv, '(I14)') year_s
+      case ("-month")
+        iarg = iarg + 1
+        call getarg ( iarg, argv )
+        read(argv, '(I14)') month_s
+      case ("-day")
+        iarg = iarg + 1
+        call getarg ( iarg, argv )
+        read(argv, '(I14)') day_s
+      case ("-save_bin")
+        save_bin = .true.
+      case ("-input_data_path")
+        iarg = iarg + 1
+        call getarg ( iarg, argv )
+        read(argv, *) data_path
+        data_path = trim(data_path)
+      case default
+    end select
+  end do 
+  print *, year_s, month_s, day_s, save_bin, data_path
+
+  !-- form a string of input date
+  extract_date = year_s*10000+month_s*100+day_s
+  write (date_str, format_str) extract_date
+
+  ! GMAO OPS has following fixed
+  sst_file_pref = 'dataoceanfile_OSTIA_REYNOLDS_SST.2880x1440.'
+  ice_file_pref = 'dataoceanfile_OSTIA_REYNOLDS_ICE.2880x1440.'
+  sst_file_suff = '.data'
+  ice_file_suff = '.data'
+
+  sst_file = trim(data_path)//"/"//trim(sst_file_pref)//trim(date_str(1:4))//trim(sst_file_suff)
+  ice_file = trim(data_path)//"/"//trim(ice_file_pref)//trim(date_str(1:4))//trim(ice_file_suff)
+
+  print *, " "
+  print *, "Read: ", sst_file, ice_file
+  print *, " "
+
+  ! read sst
+  CALL read_bin_SST_ICE( sst_file, extract_date, date_out, nlon_out, nlat_out, lon, lat, sst)
+  ! read sea ice
+  CALL read_bin_SST_ICE( ice_file, extract_date, date_out, nlon_out, nlat_out, lon, lat, ice)
+
+  print *, date_out
+  print *, " "
+  print *, "Write out the data..."
+
+  ! binary write
+  if (save_bin) then
+    tom = INCYMD( year_s*10000+month_s*100+day_s, 1) ! next day
+
+    tom_year = int(tom/10000)
+    tom_mon  = int((tom-tom_year*10000)/100)
+    tom_day  = tom - (tom_year*10000+tom_mon*100)
+    !print *, "Tomorrow day:", tom_year, tom_mon, tom_day
+    CALL write_bin( 'bcs_2880x1440', year_s, month_s, day_s, tom_year, tom_mon, tom_day, date_str,  nlat, nlon, transpose(sst), transpose(ice))
+  end if
+
+  ! netcdf
+  CALL write_netcdf( 'sst_ice', year_s, month_s, day_s, date_str, nlat, nlon, transpose(sst), transpose(ice))
+
+  print *, "Done."
+  print *, " "
+
+!---------------------------------------------------------------------------
+END PROGRAM extract_day_SST_FRACI_eight