diff --git a/CHANGELOG.md b/CHANGELOG.md index daa19ef4..3947fee3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- Addition of SST and FRACI forecast (ocean) boundary conditions generation capability in `pre/prepare_ocnExtData` - Added EASE grid option for remapping of land restarts in remap_restarts.py package (facilitates use of package in GEOSldas setup script) - Added support for SLES15, NAS site and log for remap_lake_landice_saltwater in remap_restarts.py package - Added "land_only" option for remapping of restarts @@ -56,6 +57,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [2.0.2] - 2023-06-28 + ### Fixed - Updates to fix dates for comparison plots. Updates for plot levels. diff --git a/pre/prepare_ocnExtData/CMakeLists.txt b/pre/prepare_ocnExtData/CMakeLists.txt index 19b827c0..e67d9938 100644 --- a/pre/prepare_ocnExtData/CMakeLists.txt +++ b/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) diff --git a/pre/prepare_ocnExtData/daily_clim_SST_FRACI_eight.F90 b/pre/prepare_ocnExtData/daily_clim_SST_FRACI_eight.F90 new file mode 100644 index 00000000..22ca2762 --- /dev/null +++ b/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 diff --git a/pre/prepare_ocnExtData/extract_day_SST_FRACI_eight.F90 b/pre/prepare_ocnExtData/extract_day_SST_FRACI_eight.F90 new file mode 100644 index 00000000..93e90ca5 --- /dev/null +++ b/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 diff --git a/pre/prepare_ocnExtData/forecast_SST_FRACI_eight.F90 b/pre/prepare_ocnExtData/forecast_SST_FRACI_eight.F90 new file mode 100644 index 00000000..b754f4fa --- /dev/null +++ b/pre/prepare_ocnExtData/forecast_SST_FRACI_eight.F90 @@ -0,0 +1,231 @@ +PROGRAM forecast_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) + real :: sst0(nlat, nlon), ice0(nlat, nlon) + real :: sst_daily_clim(nlat, nlon), ice_daily_clim(nlat, nlon) + real :: anom_sst0(nlat, nlon), anom_ice0(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) :: format_str = "(I8)" + character (len=8) :: date_str, today_str, tom_str + character (len = 200) :: sst_file, ice_file + + integer :: i, iarg, argc, iargc + character (len=255) :: argv + character (len=255) :: method + integer :: numArgs + integer :: year_s, month_s, day_s, tom + integer :: tom_year, tom_mon, tom_day + integer :: today_year, today_mon, today_day + + integer :: fcst_nDays +!---- + + print *, " " + print *, "-----------------------------------------------------------" + print *, " Generating SST and FRACI for forecasts " + print *, "-----------------------------------------------------------" + print *, " " + + argc = iargc() + if (argc < 3) then + print *, "Need minmum 3 inputs, see below examples. Try again!" + print *, "gen_forecast_bcs.x -year 2010 -month 1 -day 1" + print *, " " + print *, "** 3 additional optional inputs are allowed ** " + print *, "gen_forecast_bcs.x -year 2010 -month 1 -day 1 -fcst_nDays NN" + print *, " " + print *, "gen_forecast_bcs.x -year 2010 -month 1 -day 1 -input_data_path xx" + print *, " " + print *, "gen_forecast_bcs.x -year 2010 -month 1 -day 1 -fcst_nDays NN -input_data_path xx" + print *, " " + print *, "gen_forecast_bcs.x -year 2010 -month 1 -day 1 -fcst_nDays NN -input_data_path xx -method yy" + print *, " " + STOP + end if + + ! Following variables have defaults that can be overridden by user inputs. + fcst_nDays = 15 ! default length of each forecast: 15 days. + ! GMAO OPS data path + data_path= '/discover/nobackup/projects/gmao/share/dao_ops/fvInput/g5gcm/bcs/realtime/OSTIA_REYNOLDS/2880x1440/' + method = 'persistence' + + !-- 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("-fcst_nDays") + iarg = iarg + 1 + call getarg ( iarg, argv ) + read(argv, '(I14)') fcst_nDays + case ("-input_data_path") + iarg = iarg + 1 + call getarg ( iarg, argv ) + read(argv, *) data_path + data_path = trim(data_path) + case ("-method") + iarg = iarg + 1 + call getarg ( iarg, argv ) + read(argv, *) method + method = trim(method) + case default + end select + end do + print *, "Inputs: " + print *, year_s, month_s, day_s, fcst_nDays, data_path, method + print *, " " + + ! --- Read forecast start date data --- + ! -- + !-- 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) + + ! read start date sst + CALL read_bin_SST_ICE( sst_file, extract_date, date_out, nlon_out, nlat_out, lon, lat, sst0) + ! read start date fraci + CALL read_bin_SST_ICE( ice_file, extract_date, date_out, nlon_out, nlat_out, lon, lat, ice0) + ! -- + + print *, " " + print *, "For: ", date_out + print *, "Read: ", sst_file, ice_file + print *, " " + + ! day 1: [year_s, month_s, day_s --> year_s+1, month_s+1, day_s+1] + ! day 2: [year_s+1, month_s+1, day_s+1 --> year_s+2, month_s+2, day_s+2] + ! .. + print *, "Write out: ", trim(method), " BCs data..." + print *, " " + + today_year=year_s; today_mon=month_s; today_day=day_s + do i = 1, fcst_nDays + + tom = INCYMD( today_year*10000+today_mon*100+today_day, 1) ! tomorrow = today + (1 day) + tom_year = int(tom/10000) + tom_mon = int((tom-tom_year*10000)/100) + tom_day = tom - (tom_year*10000+tom_mon*100) + + write (today_str, format_str) today_year*10000+today_mon*100+today_day + write (tom_str, format_str) tom_year*10000 +tom_mon*100 +tom_day + print *, trim(today_str), "---->", trim(tom_str) + + if (method == "persistence") then + sst = sst0; ice = ice0 + else if (method == "persist_anomaly") then + call get_daily_clim( today_str, nlat_out, nlon_out, sst_daily_clim, ice_daily_clim) + if (i > 1) then + sst = sst_daily_clim + anom_sst0 + ice = ice_daily_clim + anom_ice0 + else + sst = sst0; ice = ice0 ! forecast day-1: just write out the same field(s). + anom_sst0 = sst0 - sst_daily_clim + anom_ice0 = ice0 - ice_daily_clim + end if + else + print *, "Unknown method: ", method, "Exiting. Fix and try again." + STOP + end if + + CALL write_bin( 'bcs_2880x1440', today_year, today_mon, today_day, tom_year, tom_mon, tom_day, today_str, & + nlat, nlon, transpose(sst), transpose(ice)) + + CALL write_netcdf( 'sst_ice', today_year, today_mon, today_day, today_str, & + nlat, nlon, transpose(sst), transpose(ice)) + + ! done with today; increment it to tomorrow + today_year=tom_year; today_mon=tom_mon; today_day=tom_day + end do + + print *, " " + print *, "Done." + print *, " " + +!--------------------------------------------------------------------------- + CONTAINS + SUBROUTINE get_daily_clim(today_str, nlat_out, nlon_out, sst_daily_clim, ice_daily_clim) + + integer :: date_out, nlon_out, nlat_out + real :: sst_daily_clim(nlat, nlon), ice_daily_clim(nlat, nlon) + character (len = 200) :: clim_data_path, sst_clim_pref, ice_clim_pref + character (len=8) :: today_str + + clim_data_path = "/discover/nobackup/projects/gmao/advda/sakella/future_sst_fraci/data/binFiles/" + sst_clim_pref = "daily_clim_mean_sst_0001" !clim_year = 1 ! from daily_clim_SST_FRACI_eight.F90 + ice_clim_pref = "daily_clim_mean_ice_0001" + + sst_file = trim(clim_data_path)//"/"//trim(sst_clim_pref)//trim(today_str(5:8))//".bin" + ice_file = trim(clim_data_path)//"/"//trim(ice_clim_pref)//trim(today_str(5:8))//".bin" + print *, "Reading daily clim from: " + print *, sst_file + print *, ice_file + + CALL read_daily_clim(sst_file, sst_daily_clim, nlat, nlon) ! read sst anomaly + CALL read_daily_clim(ice_file, ice_daily_clim, nlat, nlon) ! read fraci anomaly + print *, nlat, nlon + END SUBROUTINE get_daily_clim + + SUBROUTINE read_daily_clim(fName, bcs_field, nlat, nlon) + integer, intent(in) :: nlat, nlon + character (len=*), intent(in) :: fName + real, intent(out) :: bcs_field(nlat, nlon) + + real year1,month1,day1,hour1,min1,sec1, year2,month2,day2,hour2,min2,sec2 + real dum1,dum2, field(nlon,nlat) + integer rc + + open (10,file=fName,form='unformatted',access='sequential', STATUS = 'old') ! these files only 1 record (header+data) + rc = 0 + do while (rc.eq.0) + ! read header + read (10,iostat=rc) year1,month1,day1,hour1,min1,sec1,& + year2,month2,day2,hour2,min2,sec2,dum1,dum2 + ! read data + if( rc.eq.0 ) read (10,iostat=rc) field + bcs_field = TRANSPOSE(field) + close(10) + end do + END SUBROUTINE read_daily_clim +!--------------------------------------------------------------------------- +END PROGRAM forecast_SST_FRACI_eight diff --git a/pre/prepare_ocnExtData/gen_wts_regrid_ostia_geos.py b/pre/prepare_ocnExtData/gen_wts_regrid_ostia_geos.py new file mode 100755 index 00000000..4e4d79e9 --- /dev/null +++ b/pre/prepare_ocnExtData/gen_wts_regrid_ostia_geos.py @@ -0,0 +1,22 @@ +#!/usr/bin/env python + +import xarray as xr +import xesmf as xe + +ostia_path = "/discover/nobackup/sakella/projects/bcs_pert/data/data_OSTIA" +geos_path = "/discover/nobackup/sakella/projects/bcs_pert/tests" + +ds_ostia = xr.open_dataset(ostia_path + "/" + "20230101-UKMO-L4HRfnd-GLOB-v01-fv02-OSTIA.nc").squeeze() +ds_geos = xr.open_dataset(geos_path + "/" + "sst_fraci_20230101.nc4").squeeze() + +print("Generating weights...") +regridder = xe.Regridder(ds_ostia, ds_geos, "conservative", periodic=True) +print("Done.") + +regridder.filename = 'wts_cons_ostia_to_geos_eight.nc' +regridder.to_netcdf() +print("Saved weights for reuse to: ", regridder.filename) + +# use generated weights +# regridder = xe.Regridder(ds_ostia, ds_geos, 'conservative', weights='wts_cons_ostia_to_geos_eight.nc') +#ds_ostia_on_geos_grid = regridder(ds_ostia, keep_attrs=True) diff --git a/pre/prepare_ocnExtData/hermes_m_tick.f90 b/pre/prepare_ocnExtData/hermes_m_tick.f90 new file mode 100644 index 00000000..dd8b741b --- /dev/null +++ b/pre/prepare_ocnExtData/hermes_m_tick.f90 @@ -0,0 +1,146 @@ +! -- +! Copied from and renamed +! GMAO_Shared/main/GMAO_hermes/m_tick.f90 +! -- + +module m_tick +private +public tick +public INCYMD +contains + subroutine tick (nymd,nhms,ndt) +!*********************************************************************** +! Purpose +! Tick the Date (nymd) and Time (nhms) by NDT (seconds) +! +!*********************************************************************** +!* GODDARD LABORATORY FOR ATMOSPHERES * +!*********************************************************************** + + IF(NDT.NE.0) THEN + NSEC = NSECF(NHMS) + NDT + + IF (NSEC.GT.86400) THEN + DO WHILE (NSEC.GT.86400) + NSEC = NSEC - 86400 + NYMD = INCYMD (NYMD,1) + ENDDO + ENDIF + + IF (NSEC.EQ.86400) THEN + NSEC = 0 + NYMD = INCYMD (NYMD,1) + ENDIF + + IF (NSEC.LT.00000) THEN + DO WHILE (NSEC.LT.0) + NSEC = 86400 + NSEC + NYMD = INCYMD (NYMD,-1) + ENDDO + ENDIF + + NHMS = NHMSF (NSEC) + ENDIF + + RETURN + END + + FUNCTION INCYMD (NYMD,M) +!*********************************************************************** +! PURPOSE +! INCYMD: NYMD CHANGED BY ONE DAY +! MODYMD: NYMD CONVERTED TO JULIAN DATE +! DESCRIPTION OF PARAMETERS +! NYMD CURRENT DATE IN YYMMDD FORMAT +! M +/- 1 (DAY ADJUSTMENT) +! +!*********************************************************************** +!* GODDARD LABORATORY FOR ATMOSPHERES * +!*********************************************************************** + + INTEGER NDPM(12) + DATA NDPM /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/ + LOGICAL LEAP + DATA NY00 / 1900 / + LEAP(NY) = MOD(NY,4).EQ.0 .AND. (NY.NE.0 .OR. MOD(NY00,400).EQ.0) + +!*********************************************************************** +! + NY = NYMD / 10000 + NM = MOD(NYMD,10000) / 100 + ND = MOD(NYMD,100) + M + + IF (ND.EQ.0) THEN + NM = NM - 1 + IF (NM.EQ.0) THEN + NM = 12 + NY = NY - 1 + ENDIF + ND = NDPM(NM) + IF (NM.EQ.2 .AND. LEAP(NY)) ND = 29 + ENDIF + + IF (ND.EQ.29 .AND. NM.EQ.2 .AND. LEAP(NY)) GO TO 20 + + IF (ND.GT.NDPM(NM)) THEN + ND = 1 + NM = NM + 1 + IF (NM.GT.12) THEN + NM = 1 + NY = NY + 1 + ENDIF + ENDIF + + 20 CONTINUE + INCYMD = NY*10000 + NM*100 + ND + RETURN + +!*********************************************************************** +! E N T R Y M O D Y M D +!*********************************************************************** + + ENTRY MODYMD (NYMD) + NY = NYMD / 10000 + NM = MOD(NYMD,10000) / 100 + ND = MOD(NYMD,100) + + 40 CONTINUE + IF (NM.LE.1) GO TO 60 + NM = NM - 1 + ND = ND + NDPM(NM) + IF (NM.EQ.2 .AND. LEAP(NY)) ND = ND + 1 + GO TO 40 + + 60 CONTINUE + MODYMD = ND + RETURN + END + + function nhmsf (nsec) +!*********************************************************************** +! Purpose +! Converts Total Seconds to NHMS format +! +!*********************************************************************** +!* GODDARD LABORATORY FOR ATMOSPHERES * +!*********************************************************************** + implicit none + integer nhmsf, nsec + nhmsf = nsec/3600*10000 + mod(nsec,3600)/60*100 + mod(nsec,60) + return + end + + function nsecf (nhms) +!*********************************************************************** +! Purpose +! Converts NHMS format to Total Seconds +! +!*********************************************************************** +!* GODDARD LABORATORY FOR ATMOSPHERES * +!*********************************************************************** + implicit none + integer nhms, nsecf + nsecf = nhms/10000*3600 + mod(nhms,10000)/100*60 + mod(nhms,100) + return + end +end module m_tick diff --git a/pre/prepare_ocnExtData/plot_diff_ostia_native_regridded_binned.py b/pre/prepare_ocnExtData/plot_diff_ostia_native_regridded_binned.py new file mode 100755 index 00000000..568f9b4c --- /dev/null +++ b/pre/prepare_ocnExtData/plot_diff_ostia_native_regridded_binned.py @@ -0,0 +1,48 @@ +#!/usr/bin/env python + +import xarray as xr +import xesmf as xe + +import matplotlib.pyplot as plt + +import cartopy.crs as ccrs +import cartopy.feature as cfeature +import sys + +ostia_path = "/discover/nobackup/sakella/projects/bcs_pert/data/data_OSTIA" +geos_path = "/discover/nobackup/sakella/projects/bcs_pert/tests" + +date_str = sys.argv[1]#'20230101' +ds_ostia = xr.open_dataset(ostia_path + "/" + date_str + "-UKMO-L4HRfnd-GLOB-v01-fv02-OSTIA.nc").squeeze() +ds_geos = xr.open_dataset(geos_path + "/" + "sst_fraci_" + date_str + ".nc4").squeeze() + +regridder = xe.Regridder(ds_ostia, ds_geos, 'conservative', weights='wts_cons_ostia_to_geos_eight.nc') +ds_ostia_on_geos_grid = regridder(ds_ostia, keep_attrs=True) +# + +fig = plt.figure(figsize=(16,4)) + +ax1 = plt.subplot(121, projection=ccrs.PlateCarree( central_longitude=0)) +ax1.coastlines() +ax1.add_feature(cfeature.LAND, facecolor='white', zorder=100) +ax1.add_feature(cfeature.LAKES, edgecolor='gray', zorder=50) +im1 = ax1.pcolormesh(ds_geos['lon'], ds_geos['lat'], ds_ostia_on_geos_grid.analysed_sst-ds_geos.SST, transform=ccrs.PlateCarree(), cmap=plt.cm.bwr, vmin=-0.25, vmax=0.25) +gl = ax1.gridlines(crs=ccrs.PlateCarree(), linewidth=1, color='0.85', alpha=0.5, linestyle='-', draw_labels=True) +gl.top_labels = False +cbar=fig.colorbar(im1, extend='both', shrink=0.5, ax=ax1, pad=0.1) +cbar.set_label(r'diff SST') +ax1.set_title(r'%s'%(date_str)) + +ax2 = plt.subplot(122, projection=ccrs.PlateCarree()) +ax2.coastlines() +ax2.add_feature(cfeature.LAND, facecolor='white', zorder=100) +ax2.add_feature(cfeature.LAKES, edgecolor='gray', zorder=50) +im2 = ax2.pcolormesh(ds_geos['lon'], ds_geos['lat'], ds_ostia_on_geos_grid.sea_ice_fraction-ds_geos.FRACI, transform=ccrs.PlateCarree(), cmap=plt.cm.bwr, vmin=-0.1, vmax=0.1) +gl = ax2.gridlines(crs=ccrs.PlateCarree(), linewidth=1, color='0.85', alpha=0.5, linestyle='-', draw_labels=True) +gl.top_labels = False +cbar=fig.colorbar(im2, extend='both', shrink=0.5, ax=ax2, pad=0.1) +cbar.set_label(r'diff FRACI') +ax2.set_title(r'%s'%(date_str)) + +fName_fig = 'consr_binned_sst_ice_diff_' + date_str + '.png' +plt.savefig(fName_fig, dpi=100) diff --git a/pre/prepare_ocnExtData/proc_SST_FRACI_eight.F90 b/pre/prepare_ocnExtData/proc_SST_FRACI_eight.F90 index 2de2d495..06f85944 100644 --- a/pre/prepare_ocnExtData/proc_SST_FRACI_eight.F90 +++ b/pre/prepare_ocnExtData/proc_SST_FRACI_eight.F90 @@ -255,14 +255,14 @@ PROGRAM proc_SST_FRACI_eight READ( tomrw_Day, 99) tomrw_iDay if (bin_output == 'yes') then - call write_bin(today_iYear, today_iMon, today_iDay, & + call write_bin('bcs_2880x1440', today_iYear, today_iMon, today_iDay, & tomrw_iYear, tomrw_iMon, tomrw_iDay, & today, & NLAT_out, NLON_out, & ostia_SST_eigth, ostia_ICE_eigth) endif - call write_netcdf(today_iYear, today_iMon, today_iDay, & + call write_netcdf('sst_ice', today_iYear, today_iMon, today_iDay, & today, & NLAT_out, NLON_out, & ostia_SST_eigth, ostia_ICE_eigth) diff --git a/pre/prepare_ocnExtData/sst_ice_helpers.F90 b/pre/prepare_ocnExtData/sst_ice_helpers.F90 index 4304d275..e11bb4bf 100644 --- a/pre/prepare_ocnExtData/sst_ice_helpers.F90 +++ b/pre/prepare_ocnExtData/sst_ice_helpers.F90 @@ -17,6 +17,7 @@ module sst_ice_helpers public read_input_quart public write_bin public write_netcdf +public read_bin_SST_ICE public check contains @@ -974,7 +975,7 @@ END SUBROUTINE check ! ! Write out a binary file with `HEADER` and data, see below for its format. ! - SUBROUTINE write_bin( today_year, today_mon, today_day, & + SUBROUTINE write_bin( fileName_pref, today_year, today_mon, today_day, & tomrw_year, tomrw_mon, tomrw_day, & today, nlat, nlon, sst, ice) !--------------------------------------------------------------------------- @@ -983,6 +984,8 @@ SUBROUTINE write_bin( today_year, today_mon, today_day, & INTEGER, INTENT(IN) :: today_year, today_mon, today_day, & tomrw_year, tomrw_mon, tomrw_day, & nlat, nlon + + CHARACTER (LEN=*),INTENT(IN) :: fileName_pref CHARACTER (LEN=*),INTENT(IN) :: today REAL, INTENT(IN) :: sst(nlon, nlat), & @@ -1001,8 +1004,8 @@ SUBROUTINE write_bin( today_year, today_mon, today_day, & !--------------------------------------------------------------------------- ! Write out SST and ICE concentration data ! - fileName_sst = 'Ostia_sst_' // today //'.bin' - fileName_ice = 'Ostia_ice_' // today //'.bin' + fileName_sst = fileName_pref // '_sst_' // today //'.bin' + fileName_ice = fileName_pref // '_ice_' // today //'.bin' OPEN (UNIT = 991, FILE = fileName_sst, FORM = 'unformatted', STATUS = 'new') OPEN (UNIT = 992, FILE = fileName_ice, FORM = 'unformatted', STATUS = 'new') @@ -1014,8 +1017,8 @@ SUBROUTINE write_bin( today_year, today_mon, today_day, & CLOSE(991) CLOSE(992) - PRINT*, "Finished writing:" - PRINT*, fileName_sst, fileName_ice + !PRINT*, "Finished writing:" + !PRINT*, fileName_sst, fileName_ice !--------------------------------------------------------------------------- END SUBROUTINE write_bin @@ -1023,7 +1026,7 @@ END SUBROUTINE write_bin ! ! Write out a netcdf file, see below for its format. ! - SUBROUTINE write_netcdf( today_year, today_mon, today_day, & + SUBROUTINE write_netcdf( fileName_pref, today_year, today_mon, today_day, & today, nlat, nlon, sst, ice) !--------------------------------------------------------------------------- IMPLICIT NONE @@ -1035,7 +1038,8 @@ SUBROUTINE write_netcdf( today_year, today_mon, today_day, & REAL, INTENT(IN) :: sst(nlon, nlat), & ice(nlon, nlat) - CHARACTER (LEN = 40) :: fileName + CHARACTER (LEN=*),INTENT(IN) :: fileName_pref + CHARACTER (LEN=60) :: fileName integer :: ncid integer, parameter :: ndims = 3 @@ -1079,7 +1083,7 @@ SUBROUTINE write_netcdf( today_year, today_mon, today_day, & file_hist = "File created on " // trim(date_creation) // trim(time_creation) // ". In yyyymmddHHMMSS.millisec" write(DTSTR,10) today_year, today_mon, today_day - 10 FORMAT(I4,'-', I2, '-', I0.2) + 10 FORMAT(I0.4,'-', I0.2, '-', I0.2) TIME_UNITS = "days since " // trim(DTSTR) //" 12:00:00" dlat = 180./real(nlat); dlon = 360./real(nlon) @@ -1096,7 +1100,7 @@ SUBROUTINE write_netcdf( today_year, today_mon, today_day, & sst_out(:,:, 1) = sst; fraci_out(:,:, 1) = ice - fileName = 'sst_fraci_' // today //'.nc4' + fileName = fileName_pref // '_' // today //'.nc' call check( nf90_create(fileName, nf90_clobber, ncid)) ! Create the file. @@ -1146,9 +1150,108 @@ SUBROUTINE write_netcdf( today_year, today_mon, today_day, & call check( nf90_close(ncid) ) ! Close the file. !--------------------------------------------------------------------------- - PRINT*, "Finished writing ... ", fileName + !PRINT*, "Finished writing ... ", fileName END SUBROUTINE write_netcdf !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! +! Input: +! - File name containing (binary formatted) SST or ICE at (1440, 2880) resolution. +! - Date (format): YYYYMMDD (4 digit year, 2 digit month, 2 digit day). +! +! Output: +! - Date (= above input date). +! - Array of lat (1440), lon (2880), field (1440, 2880) for user input date; field = SST or ICE. +! +! ..................................................................... + SUBROUTINE read_bin_SST_ICE( fileName, nymd_in, nymd_out, NLON, NLAT, LON, LAT, bcs_field) + + IMPLICIT NONE + + CHARACTER (LEN = *), INTENT(IN) :: fileName + INTEGER, INTENT(IN) :: nymd_in + + INTEGER, INTENT(OUT) :: nymd_out + INTEGER, INTENT(OUT) :: NLON + INTEGER, INTENT(OUT) :: NLAT + + REAL, INTENT(OUT) :: LON(2880) + REAL, INTENT(OUT) :: LAT(1440) + REAL, INTENT(OUT) :: bcs_field(1440, 2880) + +! LOCAL VARS + real year1,month1,day1,hour1,min1,sec1 + real year2,month2,day2,hour2,min2,sec2 + real dum1,dum2 + real dLat,dLon + real field(2880,1440) + + integer nymd1,nhms1 + integer nymd2,nhms2 + + integer rc, ix +! .................................................................... + + !print *, 'Input date: ', nymd_in + + open (10,file=fileName,form='unformatted',access='sequential', STATUS = 'old') +! .................................................................... + rc = 0 + do while (rc.eq.0) + ! read header + read (10,iostat=rc) year1,month1,day1,hour1,min1,sec1, & + year2,month2,day2,hour2,min2,sec2,dum1,dum2 + if( rc.eq.0 ) then + read (10,iostat=rc) field + + NLON = nint(dum1) + NLAT = nint(dum2) + IF( (NLON /= 2880) .or. (NLAT /= 1440)) THEN + PRINT *, 'ERROR in LAT/LON dimension in file: ', fileName, 'NLON=', NLON, 'NLAT=', NLAT + EXIT + END IF + dLat = 180./NLAT + dLon = 360./NLON + + ! start date + nymd1 = nint( year1*10000 ) + nint (month1*100) + nint( day1 ) + nhms1 = nint( hour1*10000 ) + nint ( min1*100) + nint( sec1 ) + + ! end date + nymd2 = nint( year2*10000 ) + nint (month2*100) + nint( day2 ) + nhms2 = nint( hour2*10000 ) + nint ( min2*100) + nint( sec2 ) + end if + + !print *, nymd1, nymd2 + + if (nymd1 .eq. nymd_in) then ! Found the date for which data is requested. + nymd_out = nymd1 + +! print *, 'Yay! Found data for: ', nymd_out +!! print *, bcs_field + + LAT(1) = -90. + dLat/2. + DO ix = 2, NLAT + LAT(ix) = LAT(ix-1) + dLat + END DO + + LON(1) = -180. + dLon/2. + DO ix = 2, NLON + LON(ix) = LON(ix-1) + dLon + END DO + + bcs_field = TRANSPOSE(field) + close(10) + RETURN + else +! print *, "Date in sequential file: ", nymd1, "Continue looking: ", nymd_in + end if + end do +! .................................................................... + close(10) +!--------------------------------------------------------------------------- + + END SUBROUTINE read_bin_SST_ICE +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! end module sst_ice_helpers