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 ! 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. diff --git a/src/icebergs_fms2io.F90 b/src/icebergs_fms2io.F90 index 844b660..ddc97f2 100644 --- a/src/icebergs_fms2io.F90 +++ b/src/icebergs_fms2io.F90 @@ -520,7 +520,7 @@ subroutine write_restart_bergs(bergs, time_stamp) enddo; enddo !End of loop over grid filename = "RESTART/"//trim("bonds_iceberg.res.nc") - if (open_file(fileobj, filename, "overwrite", bergs%grd%domain, is_restart=.true.)) & + if (.not. open_file(fileobj, filename, "overwrite", bergs%grd%domain, is_restart=.true.)) & call error_mesg('write_icebegrs_restart', "Error opening the bonds icebergs restart file to write", FATAL) call register_unlimited_compressed_axis(fileobj,'i',nbonds) diff --git a/src/icebergs_framework.F90 b/src/icebergs_framework.F90 index c9b5e0e..c327135 100644 --- a/src/icebergs_framework.F90 +++ b/src/icebergs_framework.F90 @@ -3974,7 +3974,8 @@ subroutine unpack_bond_traj_from_buffer2(first, buff, n) integer :: counter ! Position in stack integer :: cnt1, ij1, cnt2, ij2 - if (dem) allocate(bond_traj%tangd1,bond_traj%tangd2,bond_traj%nstress,bond_traj%sstress,bond_traj%rel_rotation) + if (dem) allocate(bond_traj%tangd1,bond_traj%tangd2,bond_traj%nstress,& + bond_traj%sstress,bond_traj%rel_rotation,bond_traj%broken) counter = 0 call pull_buffer_value(buff%data(:,n),counter,bond_traj%lon) @@ -5355,7 +5356,8 @@ subroutine record_posn(bergs) if (dem) allocate(posn%ang_vel,posn%ang_accel,posn%rot) - if (save_bond_traj .and. dem) allocate(bond_posn%tangd1,bond_posn%tangd2,bond_posn%nstress,bond_posn%sstress,bond_posn%rel_rotation) + if (save_bond_traj .and. dem) allocate(bond_posn%tangd1,bond_posn%tangd2,bond_posn%nstress,& + bond_posn%sstress,bond_posn%rel_rotation,bond_posn%broken) area_thres=bergs%traj_area_thres*1.e6 ! convert from km^2 to m^2 area_thres2=bergs%traj_area_thres_sntbc*1.e6 @@ -5479,6 +5481,7 @@ subroutine record_posn(bergs) bond_posn%nstress=current_bond%nstress bond_posn%sstress=current_bond%sstress bond_posn%rel_rotation=current_bond%rel_rotation + bond_posn%broken=current_bond%broken endif call push_bond_posn(current_bond%bond_trajectory, bond_posn) @@ -5530,7 +5533,8 @@ subroutine push_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=>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