From 9c8266d37ce8fcd3087a29136894bf4789f3bd76 Mon Sep 17 00:00:00 2001 From: alex-huth Date: Wed, 10 Jan 2024 13:08:10 -0500 Subject: [PATCH 1/3] fixed noted units of calving --- src/icebergs.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/icebergs.F90 b/src/icebergs.F90 index ebff51c..c228aef 100644 --- a/src/icebergs.F90 +++ b/src/icebergs.F90 @@ -5076,7 +5076,7 @@ subroutine icebergs_run(bergs, time, calving, uo, vo, ui, vi, tauxa, tauya, ssh, ! Arguments type(icebergs), pointer :: bergs !< Container for all types and memory type(time_type), intent(in) :: time !< Model time - real, dimension(:,:), intent(inout) :: calving !< Calving (kg/s). This field is updated with melt by bergs. + real, dimension(:,:), intent(inout) :: calving !< Calving (kg m-2 s-1). This field is updated with melt by bergs. real, dimension(:,:), intent(inout) :: calving_hflx !< Calving heat flux (W/m2) real, dimension(:,:), intent(in) :: uo !< Ocean zonal velocity (m/s) real, dimension(:,:), intent(in) :: vo !< Ocean meridional velocity (m/s) From 89cf3e63fde1ddbc613eecbff528178be4a89a47 Mon Sep 17 00:00:00 2001 From: alex-huth Date: Fri, 2 Feb 2024 11:22:09 -0500 Subject: [PATCH 2/3] Use masked-out SST to determine units of SST, preventing an error where some PEs erroneously trigger (or not) a K to C SST conversion due to unrealistic SST values under ice shelves (which occur on zero-thickness ocean surface layers and are unused otherwise) --- src/icebergs.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/icebergs.F90 b/src/icebergs.F90 index c228aef..7a896fe 100644 --- a/src/icebergs.F90 +++ b/src/icebergs.F90 @@ -5337,7 +5337,7 @@ subroutine icebergs_run(bergs, time, calving, uo, vo, ui, vi, tauxa, tauya, ssh, endif call mpp_update_domains(grd%ssh, grd%domain) - max_SST = maxval(sst(:,:)) + max_SST = maxval(sst(:,:) * grd%msk(grd%isc:grd%iec,grd%jsc:grd%jec)) if (max_SST > 120.0) then ! The input sst is in degrees Kelvin, otherwise the water would be boiling. grd%sst(grd%isc:grd%iec,grd%jsc:grd%jec) = sst(:,:)-273.15 ! Note convert from Kelvin to Celsius else ! The input sst is already in degrees Celsius. From 9d77ef15f2ab647ddce47521ddf0e228e1bfe423 Mon Sep 17 00:00:00 2001 From: alex-huth Date: Wed, 7 Feb 2024 13:50:17 -0500 Subject: [PATCH 3/3] Several miscellaneous updates: Modified the icebergs driver and driver_data modules to be FMS2-compatible; updated the default solo-mode make commands and build environment for gaea c5; updated the namelist file for the a68a test case so that the data directory will properly be assigned on any system by default; fixed an allocation error for the iceberg bond trajectory 'broken' (which saves whether a bond is broken or not); fixed a boolean error in icebergs_fms2io.F90 that was preventing iceberg bonds from being written. All solo-mode tests passed. --- build/a68_dims__genmod.f90 | 12 -- build/env | 28 +++-- build/handle_err__genmod.f90 | 10 -- build/mkmkf | 8 +- build/mkmkf_fms2 | 4 +- driver/driver_data.F90 | 187 ++-------------------------- driver/driver_data_fms.F90 | 185 ++++++++++++++++++++++++++++ driver/driver_data_fms2.F90 | 229 +++++++++++++++++++++++++++++++++++ driver/icebergs_driver.F90 | 10 +- src/icebergs_fms2io.F90 | 2 +- src/icebergs_framework.F90 | 13 +- tests/a68_test/long_run.nml | 6 +- 12 files changed, 462 insertions(+), 232 deletions(-) delete mode 100644 build/a68_dims__genmod.f90 delete mode 100644 build/handle_err__genmod.f90 create mode 100644 driver/driver_data_fms.F90 create mode 100644 driver/driver_data_fms2.F90 diff --git a/build/a68_dims__genmod.f90 b/build/a68_dims__genmod.f90 deleted file mode 100644 index 3e781da..0000000 --- a/build/a68_dims__genmod.f90 +++ /dev/null @@ -1,12 +0,0 @@ - !COMPILER-GENERATED INTERFACE MODULE: Fri Oct 1 13:21:29 2021 - ! This source file is for reference only and may not completely - ! represent the generated interface used by the compiler. - MODULE A68_DIMS__genmod - INTERFACE - SUBROUTINE A68_DIMS(DATA_DIR,NI,NJ) - CHARACTER(LEN=1000) :: DATA_DIR - INTEGER(KIND=4), INTENT(OUT) :: NI - INTEGER(KIND=4), INTENT(OUT) :: NJ - END SUBROUTINE A68_DIMS - END INTERFACE - END MODULE A68_DIMS__genmod diff --git a/build/env b/build/env index 8dd47c1..c751803 100644 --- a/build/env +++ b/build/env @@ -1,16 +1,18 @@ -#source /opt/modules/default/etc/modules.sh -module use -a /ncrc/home2/fms/local/modulefiles -setenv MODULEPATH /usw/eslogin/modulefiles-c4:/sw/eslogin-c4/modulefiles:/opt/cray/pe/ari/modulefiles:/opt/cray/ari/modulefiles:/opt/cray/pe/modulefiles:/opt/cray/modulefiles:/opt/modulefiles:/sw/common/modulefiles -#module unload PrgEnv-pgi PrgEnv-intel PrgEnv-gnu darshan; module load PrgEnv-gnu; module unload netcdf gcc; module load gcc/7.3.0 cray-hdf5 cray-netcdf - -module unload PrgEnv-pgi -module unload PrgEnv-pathscale -module unload PrgEnv-intel -module unload PrgEnv-gnu -module unload PrgEnv-cray +#module unload PrgEnv-pgi +#module unload PrgEnv-pathscale +#module unload PrgEnv-intel +#module unload PrgEnv-gnu +#module unload PrgEnv-cray module load PrgEnv-intel -module swap intel intel/18.0.6.288 -module unload netcdf +module load intel-classic/2022.0.2 +#module swap intel intel/18.0.6.288 +#module load intel +#module unload netcdf +module load cray-hdf5 module load cray-netcdf -module load cray-hdf5 \ No newline at end of file + +module switch cray-libsci/22.10.1.2 +#module load cray-hdf5/1.12.1.3 +#module load cray-netcdf/4.8.1.3 +#module load cray-hdf5 diff --git a/build/handle_err__genmod.f90 b/build/handle_err__genmod.f90 deleted file mode 100644 index 0dd846d..0000000 --- a/build/handle_err__genmod.f90 +++ /dev/null @@ -1,10 +0,0 @@ - !COMPILER-GENERATED INTERFACE MODULE: Fri Oct 1 13:21:29 2021 - ! This source file is for reference only and may not completely - ! represent the generated interface used by the compiler. - MODULE HANDLE_ERR__genmod - INTERFACE - SUBROUTINE HANDLE_ERR(STATUS) - INTEGER(KIND=4), INTENT(IN) :: STATUS - END SUBROUTINE HANDLE_ERR - END INTERFACE - END MODULE HANDLE_ERR__genmod diff --git a/build/mkmkf b/build/mkmkf index 6aa5125..595590d 100644 --- a/build/mkmkf +++ b/build/mkmkf @@ -1,3 +1,7 @@ +# rm -f path_names +# ../../mkmf/bin/list_paths ../../FMS/{mpp,diag_manager,time_manager,include,memutils,constants,platform,fms,random_numbers,mosaic,exchange} ../src/ ../driver/ +# ../../mkmf/bin/mkmf -t ../../mkmf/templates/ncrc-intel.mk -c "-Duse_libMPI -Duse_netCDF" -p bergs.x path_names + rm -f path_names -../../mkmf/bin/list_paths ../../FMS/{mpp,diag_manager,time_manager,include,memutils,constants,platform,fms,random_numbers,mosaic,exchange} ../src/ ../driver/ -../../mkmf/bin/mkmf -t ../../mkmf/templates/ncrc-intel.mk -c "-Duse_libMPI -Duse_netCDF" -p bergs.x path_names \ No newline at end of file +../../mkmf/bin/list_paths -l ../src/ ../driver/ +../../mkmf/bin/mkmf -t ../../mkmf/templates/ncrc-intel.mk -o "-I../../../build/fms" -p bergs.x -l '-L../../../build/fms -lfms' -c "-Duse_libMPI -Duse_netCDF -I../../../build/fms" path_names \ No newline at end of file diff --git a/build/mkmkf_fms2 b/build/mkmkf_fms2 index 84932e2..58877cd 100644 --- a/build/mkmkf_fms2 +++ b/build/mkmkf_fms2 @@ -1,3 +1,3 @@ rm -f path_names -../../mkmf/bin/list_paths ../src/ ../driver/ -../../mkmf/bin/mkmf -t ../../mkmf/templates/ncrc-intel.mk -c "-Duse_libMPI -Duse_netCDF -Duse_fms2_io -I../../../build/fms2" -p bergs.x -l '-L../../../build/fms2 -lfms' path_names \ No newline at end of file +../../mkmf/bin/list_paths -l ../src/ ../driver/ +../../mkmf/bin/mkmf -t ../../mkmf/templates/ncrc-intel.mk -o "-I../../../build/fms2" -p bergs.x -l '-L../../../build/fms2 -lfms' -c "-Duse_libMPI -Duse_netCDF -DUSE_FMS2_IO -I../../../build/fms2" path_names \ No newline at end of file diff --git a/driver/driver_data.F90 b/driver/driver_data.F90 index b60b015..12fa0ab 100644 --- a/driver/driver_data.F90 +++ b/driver/driver_data.F90 @@ -1,184 +1,11 @@ +!> Handles reading/writing of restart files and trajectory-based diagnostic files module driver_data +! This file is part of NOAA-GFDL/icebergs. See LICENSE.md for the license. -use fms_mod, only : read_data -use fms_mod, only : file_exist -use fms_mod, only : field_exist -use fms_mod, only : error_mesg -use fms_mod, only : FATAL -use mpp_domains_mod, only : domain2D -use mpp_mod, only : mpp_pe -use mpp_mod, only : mpp_root_pe -use mpp_mod, only: mpp_sync -use constants_mod, only: pi -use netcdf -!include 'netcdf.inc' - -implicit none - -contains - - !> Determine dimensions of A68 background grid - subroutine a68_dims(data_dir,ni,nj) - ! Arguments - character(len=*) :: data_dir - integer, intent(out):: ni,nj - ! Local variables - character(len=1000):: infile - integer(kind=4) :: ncid - integer :: status - - infile=trim(trim(data_dir)//'a68_experiment_ll_p125_grid.nc') - - !Open netCDF file - call handle_err(nf90_open(infile, nf90_nowrite, ncid)) - !Inquire about the dimensions - call handle_err(nf90_inquire_dimension(ncid,1,len=nj)) - call handle_err(nf90_inquire_dimension(ncid,2,len=ni)) - !Close netCDF file - call handle_err(nf90_close(ncid)) - end subroutine a68_dims - - !> Handle for errors - subroutine handle_err(status) - integer, intent (in) :: status - - if(status /= nf90_noerr) then - print *, trim(nf90_strerror(status)) - stop "Stopped" - end if - end subroutine handle_err - - !>A68 grid setup - subroutine a68_prep(data_dir,mpp_domain,lon,lat,dx,dy,area,REarth) - ! Arguments - character(len=*) :: data_dir - type(domain2D) :: mpp_domain - real :: lon(:,:) !< Longitude or position in x - real :: lat(:,:) !< Latitude or position in y - real :: dx(:,:) !< Length of northern edge of cell (m) - real :: dy(:,:) !< Length of eatern edge of cell (m) - real :: area(:,:) !< Area of cell (m2) - real :: REarth - ! Local variables - character(len=1000) :: filename - real, parameter :: gres=0.125 - - !Grid specs: lon,lat,dx,dy,area - filename=trim(trim(data_dir)//'a68_experiment_ll_p125_grid.nc') - call read_a68_data(mpp_domain,filename,'longitude',lon) - lon=lon+360 - call read_a68_data(mpp_domain,filename,'latitude',lat) - - call haversine_dist_and_area(REarth,gres,lon,lat,dx,dy,area) - end subroutine a68_prep - - !> Haversine distance and area - subroutine haversine_dist_and_area(REarth,gres,lon1,lat1,dx,dy,area) - ! Arguments - real, intent(in) :: REarth !Earth radius (m) - real, intent(in) :: gres !grid resolution (degrees) - real, intent(in) :: lon1(:,:) !grid longitude - real, intent(in) :: lat1(:,:) !grid latitude - real, intent(out) :: dx(:,:) !grid dx (m) - real, intent(out) :: dy(:,:) !grid dy (m) - real, intent(out) :: area(:,:) !grid area (m^2) - ! Local variables - real, parameter :: pi_180=pi/180. - real, dimension(size(lon1,1),size(lon1,2)) :: lon2,lat2,p1,p2,dp,dm,a,c - - !calculate dx - lon2=lon1-gres - lat2=lat1 - p1 = lat1 * pi_180; p2 = lat2 * pi_180; dp = (lat2-lat1) * pi_180; dm = (lon2-lon1) * pi_180 - a = sin(0.5*dp)**2. + cos(p1) * cos(p2) * sin(0.5*dm)**2.; c = 2. * atan2(sqrt(a), sqrt(1-a)) - dx = REarth * c - - !calculate dy - lat2=lat1-gres - lon2=lon1 - p1 = lat1 * pi_180; p2 = lat2 * pi_180; dp = (lat2-lat1) * pi_180; dm = (lon2-lon1) * pi_180 - a = sin(0.5*dp)**2. + cos(p1) * cos(p2) * sin(0.5*dm)**2.; c = 2. * atan2(sqrt(a), sqrt(1-a)) - dy = REarth * c - - !calculate area - area=pi_180*(REarth**2.)*abs(sin(lat1*pi_180)-sin(lat2*pi_180))*abs(gres) - end subroutine haversine_dist_and_area - - - !> Read in static (not time-varying) data for A68 experiment from file - subroutine read_a68_data(domain,filename,var,field) - ! Arguments - type(domain2D) :: domain !< Parallel decomposition - character(len=*), intent(in) :: filename, var - real :: field(:,:) - - call mpp_sync() - if (mpp_pe().eq.mpp_root_pe()) write(*,'(2a)') 'icebergs_driver: Reading ',var - if (file_exist(filename)) then - if (field_exist(filename, var)) then - call read_data(filename, var, field, domain) - else - if (mpp_pe().eq.mpp_root_pe()) then - call error_mesg('icebergs_driver:','Variable missing in file!', FATAL) - endif - endif - else - if (mpp_pe().eq.mpp_root_pe()) then - call error_mesg('icebergs_driver:','File for variable not found!', FATAL) - endif - endif - call mpp_sync() - end subroutine read_a68_data - - !> Calls read_a68_data for wind velocity, ocean velocity, and ssh - subroutine a68_prep_3d(data_dir,mpp_domain,tauxa_hr,tauya_hr,uo_hr,vo_hr,ssh_hr) - ! Arguments - character(len=*) :: data_dir - - type(domain2D) :: mpp_domain - real :: tauxa_hr(:,:,:) !< Zonal jra55 wind velocity (m/s, hourly) - real :: tauya_hr(:,:,:) !< Meridional jra55 wind velocity (m/s, hourly) - real :: uo_hr(:,:,:) !< Zonal OSCAR ocean surface velocity (m/s, hourly) - real :: vo_hr(:,:,:) !< Meridional OSCAR ocean surface velocity (m/s, hourly) - real :: ssh_hr(:,:,:) !< Effective sea surface height (m, hourly) - ! Local variables - character(len=1000) :: filename - - filename=trim(trim(data_dir)//'a68_experiment_wind_vel_ncep_10m_dec2020_HOURLY_ll_p125.nc') - call read_a68_data_3d(mpp_domain,filename,'ua',tauxa_hr) - call read_a68_data_3d(mpp_domain,filename,'va',tauya_hr) - - filename=trim(trim(data_dir)//'a68_experiment_ocean_surf_vel_oscar_dec2020_HOURLY_ll_p125.nc') - call read_a68_data_3d(mpp_domain,filename,'uo',uo_hr) - call read_a68_data_3d(mpp_domain,filename,'vo',vo_hr) - - filename=trim(trim(data_dir)//'a68_experiment_ssh_duacs_dec2020_HOURLY_ll_p125.nc') - call read_a68_data_3d(mpp_domain,filename,'SSH',ssh_hr) - end subroutine a68_prep_3d - - !> Read in time-varying A68 data - subroutine read_a68_data_3d(domain,filename,var,field) - ! Arguments - type(domain2D) :: domain !< Parallel decomposition - character(len=*), intent(in) :: filename, var - real :: field(:,:,:) - - call mpp_sync() - if (mpp_pe().eq.mpp_root_pe()) write(*,'(2a)') 'icebergs_driver: Reading ',var - if (file_exist(filename)) then - if (field_exist(filename, var)) then - call read_data(filename, var, field, domain) - else - if (mpp_pe().eq.mpp_root_pe()) then - call error_mesg('icebergs_driver:','Variable missing in file!', FATAL) - endif - endif - else - if (mpp_pe().eq.mpp_root_pe()) then - call error_mesg('icebergs_driver:','File for variable not found!', FATAL) - endif - endif - call mpp_sync() - end subroutine read_a68_data_3d +#ifndef USE_FMS2_IO + use driver_data_fms +#else + use driver_data_fms2 +#endif end module driver_data diff --git a/driver/driver_data_fms.F90 b/driver/driver_data_fms.F90 new file mode 100644 index 0000000..b0a09bc --- /dev/null +++ b/driver/driver_data_fms.F90 @@ -0,0 +1,185 @@ +module driver_data_fms +#ifndef USE_FMS2_IO + +use fms_mod, only : read_data +use fms_mod, only : file_exist +use fms_mod, only : field_exist +use fms_mod, only : error_mesg +use fms_mod, only : FATAL +use mpp_domains_mod, only : domain2D +use mpp_mod, only : mpp_pe +use mpp_mod, only : mpp_root_pe +use mpp_mod, only: mpp_sync +use constants_mod, only: pi +use netcdf +!include 'netcdf.inc' + +implicit none + +contains + + !> Determine dimensions of A68 background grid + subroutine a68_dims(data_dir,ni,nj) + ! Arguments + character(len=*) :: data_dir + integer, intent(out):: ni,nj + ! Local variables + character(len=1000):: infile + integer(kind=4) :: ncid + integer :: status + + infile=trim(trim(data_dir)//'a68_experiment_ll_p125_grid.nc') + + !Open netCDF file + call handle_err(nf90_open(infile, nf90_nowrite, ncid)) + !Inquire about the dimensions + call handle_err(nf90_inquire_dimension(ncid,1,len=nj)) + call handle_err(nf90_inquire_dimension(ncid,2,len=ni)) + !Close netCDF file + call handle_err(nf90_close(ncid)) + end subroutine a68_dims + + !> Handle for errors + subroutine handle_err(status) + integer, intent (in) :: status + + if(status /= nf90_noerr) then + print *, trim(nf90_strerror(status)) + stop "Stopped" + end if + end subroutine handle_err + + !>A68 grid setup + subroutine a68_prep(data_dir,mpp_domain,lon,lat,dx,dy,area,REarth) + ! Arguments + character(len=*) :: data_dir + type(domain2D) :: mpp_domain + real :: lon(:,:) !< Longitude or position in x + real :: lat(:,:) !< Latitude or position in y + real :: dx(:,:) !< Length of northern edge of cell (m) + real :: dy(:,:) !< Length of eatern edge of cell (m) + real :: area(:,:) !< Area of cell (m2) + real :: REarth + ! Local variables + character(len=1000) :: filename + real, parameter :: gres=0.125 + + !Grid specs: lon,lat,dx,dy,area + filename=trim(trim(data_dir)//'a68_experiment_ll_p125_grid.nc') + call read_a68_data(mpp_domain,filename,'longitude',lon) + lon=lon+360 + call read_a68_data(mpp_domain,filename,'latitude',lat) + + call haversine_dist_and_area(REarth,gres,lon,lat,dx,dy,area) + end subroutine a68_prep + + !> Haversine distance and area + subroutine haversine_dist_and_area(REarth,gres,lon1,lat1,dx,dy,area) + ! Arguments + real, intent(in) :: REarth !Earth radius (m) + real, intent(in) :: gres !grid resolution (degrees) + real, intent(in) :: lon1(:,:) !grid longitude + real, intent(in) :: lat1(:,:) !grid latitude + real, intent(out) :: dx(:,:) !grid dx (m) + real, intent(out) :: dy(:,:) !grid dy (m) + real, intent(out) :: area(:,:) !grid area (m^2) + ! Local variables + real, parameter :: pi_180=pi/180. + real, dimension(size(lon1,1),size(lon1,2)) :: lon2,lat2,p1,p2,dp,dm,a,c + + !calculate dx + lon2=lon1-gres + lat2=lat1 + p1 = lat1 * pi_180; p2 = lat2 * pi_180; dp = (lat2-lat1) * pi_180; dm = (lon2-lon1) * pi_180 + a = sin(0.5*dp)**2. + cos(p1) * cos(p2) * sin(0.5*dm)**2.; c = 2. * atan2(sqrt(a), sqrt(1-a)) + dx = REarth * c + + !calculate dy + lat2=lat1-gres + lon2=lon1 + p1 = lat1 * pi_180; p2 = lat2 * pi_180; dp = (lat2-lat1) * pi_180; dm = (lon2-lon1) * pi_180 + a = sin(0.5*dp)**2. + cos(p1) * cos(p2) * sin(0.5*dm)**2.; c = 2. * atan2(sqrt(a), sqrt(1-a)) + dy = REarth * c + + !calculate area + area=pi_180*(REarth**2.)*abs(sin(lat1*pi_180)-sin(lat2*pi_180))*abs(gres) + end subroutine haversine_dist_and_area + + + !> Read in static (not time-varying) data for A68 experiment from file + subroutine read_a68_data(domain,filename,var,field) + ! Arguments + type(domain2D) :: domain !< Parallel decomposition + character(len=*), intent(in) :: filename, var + real :: field(:,:) + + call mpp_sync() + if (mpp_pe().eq.mpp_root_pe()) write(*,'(2a)') 'icebergs_driver: Reading ',var + if (file_exist(filename)) then + if (field_exist(filename, var)) then + call read_data(filename, var, field, domain) + else + if (mpp_pe().eq.mpp_root_pe()) then + call error_mesg('icebergs_driver:','Variable missing in file!', FATAL) + endif + endif + else + if (mpp_pe().eq.mpp_root_pe()) then + call error_mesg('icebergs_driver:','File for variable not found!', FATAL) + endif + endif + call mpp_sync() + end subroutine read_a68_data + + !> Calls read_a68_data for wind velocity, ocean velocity, and ssh + subroutine a68_prep_3d(data_dir,mpp_domain,tauxa_hr,tauya_hr,uo_hr,vo_hr,ssh_hr) + ! Arguments + character(len=*) :: data_dir + + type(domain2D) :: mpp_domain + real :: tauxa_hr(:,:,:) !< Zonal jra55 wind velocity (m/s, hourly) + real :: tauya_hr(:,:,:) !< Meridional jra55 wind velocity (m/s, hourly) + real :: uo_hr(:,:,:) !< Zonal OSCAR ocean surface velocity (m/s, hourly) + real :: vo_hr(:,:,:) !< Meridional OSCAR ocean surface velocity (m/s, hourly) + real :: ssh_hr(:,:,:) !< Effective sea surface height (m, hourly) + ! Local variables + character(len=1000) :: filename + + filename=trim(trim(data_dir)//'a68_experiment_wind_vel_ncep_10m_dec2020_HOURLY_ll_p125.nc') + call read_a68_data_3d(mpp_domain,filename,'ua',tauxa_hr) + call read_a68_data_3d(mpp_domain,filename,'va',tauya_hr) + + filename=trim(trim(data_dir)//'a68_experiment_ocean_surf_vel_oscar_dec2020_HOURLY_ll_p125.nc') + call read_a68_data_3d(mpp_domain,filename,'uo',uo_hr) + call read_a68_data_3d(mpp_domain,filename,'vo',vo_hr) + + filename=trim(trim(data_dir)//'a68_experiment_ssh_duacs_dec2020_HOURLY_ll_p125.nc') + call read_a68_data_3d(mpp_domain,filename,'SSH',ssh_hr) + end subroutine a68_prep_3d + + !> Read in time-varying A68 data + subroutine read_a68_data_3d(domain,filename,var,field) + ! Arguments + type(domain2D) :: domain !< Parallel decomposition + character(len=*), intent(in) :: filename, var + real :: field(:,:,:) + + call mpp_sync() + if (mpp_pe().eq.mpp_root_pe()) write(*,'(2a)') 'icebergs_driver: Reading ',var + if (file_exist(filename)) then + if (field_exist(filename, var)) then + call read_data(filename, var, field, domain) + else + if (mpp_pe().eq.mpp_root_pe()) then + call error_mesg('icebergs_driver:','Variable missing in file!', FATAL) + endif + endif + else + if (mpp_pe().eq.mpp_root_pe()) then + call error_mesg('icebergs_driver:','File for variable not found!', FATAL) + endif + endif + call mpp_sync() + end subroutine read_a68_data_3d +#endif +end module driver_data_fms diff --git a/driver/driver_data_fms2.F90 b/driver/driver_data_fms2.F90 new file mode 100644 index 0000000..462cb8b --- /dev/null +++ b/driver/driver_data_fms2.F90 @@ -0,0 +1,229 @@ +module driver_data_fms2 +#ifdef USE_FMS2_IO + +use fms2_io_mod, only : read_data +use fms2_io_mod, only : variable_exists +use fms2_io_mod, only: get_instance_filename +use fms2_io_mod, only: open_file +use fms2_io_mod, only: get_num_dimensions, register_axis, get_dimension_names, get_dimension_size +use fms2_io_mod, only: FmsNetcdfDomainFile_t +use fms_mod, only : error_mesg, lowercase +use fms_mod, only : FATAL +use mpp_domains_mod, only : domain2D +use mpp_mod, only : mpp_pe +use mpp_mod, only : mpp_root_pe +use mpp_mod, only: mpp_sync +use constants_mod, only: pi +use netcdf +!include 'netcdf.inc' + +implicit none + +contains + + !> Determine dimensions of A68 background grid + subroutine a68_dims(data_dir,ni,nj) + ! Arguments + character(len=*) :: data_dir + integer, intent(out):: ni,nj + ! Local variables + character(len=1000):: infile + integer(kind=4) :: ncid + integer :: status + + infile=trim(trim(data_dir)//'a68_experiment_ll_p125_grid.nc') + + !Open netCDF file + call handle_err(nf90_open(infile, nf90_nowrite, ncid)) + !Inquire about the dimensions + call handle_err(nf90_inquire_dimension(ncid,1,len=nj)) + call handle_err(nf90_inquire_dimension(ncid,2,len=ni)) + !Close netCDF file + call handle_err(nf90_close(ncid)) + end subroutine a68_dims + + !> Handle for errors + subroutine handle_err(status) + integer, intent (in) :: status + + if(status /= nf90_noerr) then + print *, trim(nf90_strerror(status)) + stop "Stopped" + end if + end subroutine handle_err + + !>A68 grid setup + subroutine a68_prep(data_dir,mpp_domain,lon,lat,dx,dy,area,REarth) + ! Arguments + character(len=*) :: data_dir + type(domain2D) :: mpp_domain + real :: lon(:,:) !< Longitude or position in x + real :: lat(:,:) !< Latitude or position in y + real :: dx(:,:) !< Length of northern edge of cell (m) + real :: dy(:,:) !< Length of eatern edge of cell (m) + real :: area(:,:) !< Area of cell (m2) + real :: REarth + ! Local variables + character(len=1000) :: filename + real, parameter :: gres=0.125 + + !Grid specs: lon,lat,dx,dy,area + filename=trim(trim(data_dir)//'a68_experiment_ll_p125_grid.nc') + call read_a68_data(mpp_domain,filename,'longitude',lon) + lon=lon+360 + call read_a68_data(mpp_domain,filename,'latitude',lat) + + call haversine_dist_and_area(REarth,gres,lon,lat,dx,dy,area) + end subroutine a68_prep + + !> Haversine distance and area + subroutine haversine_dist_and_area(REarth,gres,lon1,lat1,dx,dy,area) + ! Arguments + real, intent(in) :: REarth !Earth radius (m) + real, intent(in) :: gres !grid resolution (degrees) + real, intent(in) :: lon1(:,:) !grid longitude + real, intent(in) :: lat1(:,:) !grid latitude + real, intent(out) :: dx(:,:) !grid dx (m) + real, intent(out) :: dy(:,:) !grid dy (m) + real, intent(out) :: area(:,:) !grid area (m^2) + ! Local variables + real, parameter :: pi_180=pi/180. + real, dimension(size(lon1,1),size(lon1,2)) :: lon2,lat2,p1,p2,dp,dm,a,c + + !calculate dx + lon2=lon1-gres + lat2=lat1 + p1 = lat1 * pi_180; p2 = lat2 * pi_180; dp = (lat2-lat1) * pi_180; dm = (lon2-lon1) * pi_180 + a = sin(0.5*dp)**2. + cos(p1) * cos(p2) * sin(0.5*dm)**2.; c = 2. * atan2(sqrt(a), sqrt(1-a)) + dx = REarth * c + + !calculate dy + lat2=lat1-gres + lon2=lon1 + p1 = lat1 * pi_180; p2 = lat2 * pi_180; dp = (lat2-lat1) * pi_180; dm = (lon2-lon1) * pi_180 + a = sin(0.5*dp)**2. + cos(p1) * cos(p2) * sin(0.5*dm)**2.; c = 2. * atan2(sqrt(a), sqrt(1-a)) + dy = REarth * c + + !calculate area + area=pi_180*(REarth**2.)*abs(sin(lat1*pi_180)-sin(lat2*pi_180))*abs(gres) + end subroutine haversine_dist_and_area + + + !> Read in static (not time-varying) data for A68 experiment from file + subroutine read_a68_data(domain,filename,var,field) + ! Arguments + type(domain2D) :: domain !< Parallel decomposition + character(len=*), intent(in) :: filename, var + character(len=120) :: actual_filename + real :: field(:,:) + type(FmsNetcdfDomainFile_t) :: fileobj !< Fms2_io fileobj + + call mpp_sync() + + if (mpp_pe().eq.mpp_root_pe()) write(*,'(2a)') 'icebergs_driver: Reading ',var + call get_instance_filename(filename, actual_filename) + if (open_file(fileobj, trim(actual_filename), "read", domain)) then + call register_axis_wrapper_a68(fileobj) + if (variable_exists(fileobj, var)) then + call read_data(fileobj, var, field) + !call read_data(actual_filename, var, field, domain) + else + if (mpp_pe().eq.mpp_root_pe()) then + call error_mesg('icebergs_driver:','Variable missing in file!', FATAL) + endif + endif + else + if (mpp_pe().eq.mpp_root_pe()) then + call error_mesg('icebergs_driver:','File for variable not found!', FATAL) + endif + endif + call mpp_sync() + end subroutine read_a68_data + + !> Calls read_a68_data for wind velocity, ocean velocity, and ssh + subroutine a68_prep_3d(data_dir,mpp_domain,tauxa_hr,tauya_hr,uo_hr,vo_hr,ssh_hr) + ! Arguments + character(len=*) :: data_dir + + type(domain2D) :: mpp_domain + real :: tauxa_hr(:,:,:) !< Zonal jra55 wind velocity (m/s, hourly) + real :: tauya_hr(:,:,:) !< Meridional jra55 wind velocity (m/s, hourly) + real :: uo_hr(:,:,:) !< Zonal OSCAR ocean surface velocity (m/s, hourly) + real :: vo_hr(:,:,:) !< Meridional OSCAR ocean surface velocity (m/s, hourly) + real :: ssh_hr(:,:,:) !< Effective sea surface height (m, hourly) + ! Local variables + character(len=1000) :: filename + + filename=trim(trim(data_dir)//'a68_experiment_wind_vel_ncep_10m_dec2020_HOURLY_ll_p125.nc') + call read_a68_data_3d(mpp_domain,filename,'ua',tauxa_hr) + call read_a68_data_3d(mpp_domain,filename,'va',tauya_hr) + + filename=trim(trim(data_dir)//'a68_experiment_ocean_surf_vel_oscar_dec2020_HOURLY_ll_p125.nc') + call read_a68_data_3d(mpp_domain,filename,'uo',uo_hr) + call read_a68_data_3d(mpp_domain,filename,'vo',vo_hr) + + filename=trim(trim(data_dir)//'a68_experiment_ssh_duacs_dec2020_HOURLY_ll_p125.nc') + call read_a68_data_3d(mpp_domain,filename,'SSH',ssh_hr) + end subroutine a68_prep_3d + + !> Read in time-varying A68 data + subroutine read_a68_data_3d(domain,filename,var,field) + ! Arguments + type(domain2D) :: domain !< Parallel decomposition + character(len=*), intent(in) :: filename, var + character(len=120) :: actual_filename + real :: field(:,:,:) + type(FmsNetcdfDomainFile_t) :: fileobj !< Fms2_io fileobj + + call mpp_sync() + if (mpp_pe().eq.mpp_root_pe()) write(*,'(2a)') 'icebergs_driver: Reading ',var + call get_instance_filename(filename, actual_filename) + if (open_file(fileobj, trim(actual_filename), "read", domain)) then + call register_axis_wrapper_a68(fileobj) + if (variable_exists(fileobj, var)) then + call read_data(fileobj, var, field) + !call read_data(actual_filename, var, field, domain) + else + if (mpp_pe().eq.mpp_root_pe()) then + call error_mesg('icebergs_driver:','Variable missing in file!', FATAL) + endif + endif + else + if (mpp_pe().eq.mpp_root_pe()) then + call error_mesg('icebergs_driver:','File for variable not found!', FATAL) + endif + endif + call mpp_sync() + end subroutine read_a68_data_3d + + !> Wrapper for fms2_io register_axis calls when reading domain decomposed files +subroutine register_axis_wrapper_a68(fileobj) + type(FmsNetcdfDomainFile_t), intent(inout) :: fileobj !< Fms2_io fileobj + + character(len=20), dimension(:), allocatable :: file_dim_names !< Array of dimension names + integer :: i !< For do loops + integer :: dim_size !< Size of the dimension + integer :: ndims !< Number of dimensions in the file + + ndims = get_num_dimensions(fileobj) + allocate(file_dim_names(ndims)) + + call get_dimension_names(fileobj, file_dim_names) + + !< When reading fms2_io requires that the domain decomposed dimensions ("x" and "y") are + !! registered so that fms2_io knows that they are domain decomposed + do i = 1, ndims + if (lowercase(file_dim_names(i)) .eq. 'lon') then + call register_axis(fileobj, file_dim_names(i), "x") + elseif (lowercase(file_dim_names(i)) .eq. 'lat') then + call register_axis(fileobj, file_dim_names(i), "y") + elseif (lowercase(file_dim_names(i)) .eq. 'time') then + call get_dimension_size(fileobj, file_dim_names(i), dim_size) + call register_axis(fileobj, file_dim_names(i), dim_size) + endif + enddo + + deallocate(file_dim_names) +end subroutine register_axis_wrapper_a68 +#endif +end module driver_data_fms2 diff --git a/driver/icebergs_driver.F90 b/driver/icebergs_driver.F90 index 125f3b8..2e17027 100644 --- a/driver/icebergs_driver.F90 +++ b/driver/icebergs_driver.F90 @@ -3,8 +3,7 @@ program icebergs_driver use constants_mod, only: pi use fms_mod, only : fms_init use fms_mod, only : fms_end -use fms_mod, only : open_namelist_file -use fms_mod, only : close_file +use fms_mod, only : check_nml_error use fms_mod, only : error_mesg use fms_mod, only : FATAL use mpp_domains_mod, only : mpp_define_layout @@ -19,6 +18,7 @@ program icebergs_driver use mpp_mod, only : mpp_root_pe use mpp_mod, only: mpp_sync use mpp_mod, only: mpp_exit +use mpp_mod, only: input_nml_file use time_manager_mod, only : time_type use time_manager_mod, only : set_date use time_manager_mod, only : set_calendar_type @@ -76,6 +76,7 @@ program icebergs_driver real :: bump_depth=0 real :: sst=-2 real :: REarth=6.378e6 +integer :: ierr integer :: ibhrs=2 integer :: nmax = 2000000000 !bond_trajectory @@ -5585,7 +5589,8 @@ subroutine append_bond_posn(bond_trajectory, bond_posn_vals) allocate(new_bond_posn) - if (dem) allocate(new_bond_posn%tangd1,new_bond_posn%tangd2,new_bond_posn%nstress,new_bond_posn%sstress,new_bond_posn%rel_rotation) + if (dem) allocate(new_bond_posn%tangd1,new_bond_posn%tangd2,new_bond_posn%nstress,& + new_bond_posn%sstress,new_bond_posn%rel_rotation,new_bond_posn%broken) new_bond_posn=bond_posn_vals new_bond_posn%next=>null() diff --git a/tests/a68_test/long_run.nml b/tests/a68_test/long_run.nml index 5566d06..2aa3ff7 100644 --- a/tests/a68_test/long_run.nml +++ b/tests/a68_test/long_run.nml @@ -4,7 +4,7 @@ &icebergs_driver_nml a68_test=.true. !specify your data dir (ocean currents, winds, etc) - data_dir='/lustre/f2/dev/gfdl/Alexander.Huth/MOM6-examples/src/icebergs/tests/a68_test/data/' + data_dir='./data/' Rearth= ibdt=1800.0 !seconds ibuo=0.0 !ocean x velocity @@ -148,8 +148,8 @@ verbose=.false. !(F) Be verbose to stderr verbose_hrs=1! (24) Period between verbose messages - traj_sample_hrs=0!1 !(24) sampling period for trajectory storage - traj_write_hrs=0!6 !(480) Period for writing sampled trajectories to disk + traj_sample_hrs=1 !(24) sampling period for trajectory storage + traj_write_hrs=6 !(480) Period for writing sampled trajectories to disk save_bond_traj=.false. !(F) If T, saves bond trajectories in a separate file budget=.false. !(T) Calculate budgets debug=.false. !(F) Turn on debugging