Code to read new bathymetry files for basal stress parameterization. (C…

…ICE-Consortium#222)

* Read bathymetry file

* Add use_bathymetry flag

* change default ice_in

* Update basalstress documentation

* Remove use_bathymetry from gx1 and gx3 test

* Remove bathymetry_file from alt01 and alt03

cicecore/cicedynB/general/ice_init.F90

@@ -86,7 +86,9 @@ subroutine input_data
           oceanmixed_file, restore_ocn,   trestore
       use ice_arrays_column, only: bgc_data_dir, &
           sil_data_type, nit_data_type, fe_data_type
-      use ice_grid, only: grid_file, gridcpl_file, kmt_file, grid_type, grid_format, &
+      use ice_grid, only: grid_file, gridcpl_file, kmt_file, &
+                          bathymetry_file, use_bathymetry, &
+                          grid_type, grid_format, &
                           dxrect, dyrect
       use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, &
                                 basalstress, k1, Ktens, e_ratio, coriolis, &
@@ -153,6 +155,7 @@ subroutine input_data
 
       namelist /grid_nml/ &
         grid_format,    grid_type,       grid_file,     kmt_file,       &
+        bathymetry_file, use_bathymetry,                                &
         ncat,           nilyr,           nslyr,         nblyr,          &
         kcatbound,      gridcpl_file,    dxrect,        dyrect,         &
         close_boundaries
@@ -264,6 +267,8 @@ subroutine input_data
       grid_type    = 'rectangular'  ! define rectangular grid internally
       grid_file    = 'unknown_grid_file'
       gridcpl_file = 'unknown_gridcpl_file'
+      bathymetry_file    = 'unknown_bathymetry_file'
+      use_bathymetry = .false.
       kmt_file     = 'unknown_kmt_file'
       version_name = 'unknown_version_name'
       ncat  = 0
@@ -538,6 +543,8 @@ subroutine input_data
       call broadcast_scalar(grid_type,          master_task)
       call broadcast_scalar(grid_file,          master_task)
       call broadcast_scalar(gridcpl_file,       master_task)
+      call broadcast_scalar(bathymetry_file,    master_task)
+      call broadcast_scalar(use_bathymetry,     master_task)
       call broadcast_scalar(kmt_file,           master_task)
       call broadcast_scalar(kitd,               master_task)
       call broadcast_scalar(kcatbound,          master_task)
@@ -966,6 +973,10 @@ subroutine input_data
                                trim(grid_file)
             write(nu_diag,*) ' gridcpl_file              = ', &
                                trim(gridcpl_file)
+            write(nu_diag,*) ' bathymetry_file           = ', &
+                               trim(bathymetry_file)
+            write(nu_diag,*) ' use_bathymetry            = ', &
+                               use_bathymetry
             write(nu_diag,*) ' kmt_file                  = ', &
                                trim(kmt_file)
          endif

cicecore/cicedynB/infrastructure/ice_grid.F90

@@ -41,6 +41,7 @@ module ice_grid
          gridcpl_file , & !  input file for POP coupling grid info
          grid_file    , & !  input file for POP grid info
          kmt_file     , & !  input file for POP grid info
+         bathymetry_file, & !  input bathymetry for basalstress
          grid_spacing , & !  default of 30.e3m or set by user in namelist 
          grid_type        !  current options are rectangular (default),
                           !  displaced_pole, tripole, regional
@@ -114,6 +115,9 @@ module ice_grid
          uvm    , & ! land/boundary mask, velocity (U-cell)
          kmt        ! ocean topography mask for bathymetry (T-cell)
 
+      logical (kind=log_kind), public :: &
+         use_bathymetry     ! flag for reading in bathymetry_file
+
       logical (kind=log_kind), &
          dimension (:,:,:), allocatable, public :: &
          tmask  , & ! land/boundary mask, thickness (T-cell)
@@ -2335,20 +2339,28 @@ subroutine get_bathymetry
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
          file=__FILE__, line=__LINE__)
 
-      ! convert to total depth
-      depth(1) = thick(1)
-      do k = 2, nlevel
-         depth(k) = depth(k-1) + thick(k)
-      enddo
+      if (use_bathymetry) then
 
-      do iblk = 1, nblocks
-         do j = 1, ny_block
-         do i = 1, nx_block
-            k = kmt(i,j,iblk)
-            if (k > puny) bathymetry(i,j,iblk) = depth(k)
+         call read_basalstress_bathy
+
+      else
+
+         ! convert to total depth
+         depth(1) = thick(1)
+         do k = 2, nlevel
+            depth(k) = depth(k-1) + thick(k)
          enddo
+
+         do iblk = 1, nblocks
+            do j = 1, ny_block
+            do i = 1, nx_block
+               k = kmt(i,j,iblk)
+               if (k > puny) bathymetry(i,j,iblk) = depth(k)
+            enddo
+            enddo
          enddo
-      enddo
+
+      endif ! bathymetry_file
 
       end subroutine get_bathymetry
 
@@ -2373,25 +2385,22 @@ subroutine read_basalstress_bathy
          fid_init        ! file id for netCDF init file
 
       character (char_len_long) :: &        ! input data file names
-         init_file, &
          fieldname
 
       logical (kind=log_kind) :: diag=.true.
 
       character(len=*), parameter :: subname = '(read_basalstress_bathy)'
 
-      init_file='bathymetry.nc'
-
       if (my_task == master_task) then
 
           write (nu_diag,*) ' '
-          write (nu_diag,*) 'Initial ice file: ', trim(init_file)
-          write (*,*) 'Initial ice file: ', trim(init_file)
+          write (nu_diag,*) 'Bathymetry file: ', trim(bathymetry_file)
+          write (*,*) 'Bathymetry file: ', trim(bathymetry_file)
           call icepack_warnings_flush(nu_diag)
 
       endif
 
-      call ice_open_nc(init_file,fid_init)
+      call ice_open_nc(bathymetry_file,fid_init)
 
       fieldname='Bathymetry'
 
@@ -2407,7 +2416,7 @@ subroutine read_basalstress_bathy
       call ice_close_nc(fid_init)
 
       if (my_task == master_task) then
-         write(nu_diag,*) 'closing file ',TRIM(init_file)
+         write(nu_diag,*) 'closing file ',TRIM(bathymetry_file)
          call icepack_warnings_flush(nu_diag)
       endif
 

configuration/scripts/ice_in

@@ -48,6 +48,8 @@
     grid_type    = 'displaced_pole'
     grid_file    = 'grid'
     kmt_file     = 'kmt'
+    bathymetry_file = 'unknown_bathymetry_file'
+    use_bathymetry = .false.
     gridcpl_file = 'unknown_gridcpl_file'
     kcatbound    = 0
     dxrect         = 30.e5

configuration/scripts/options/set_nml.alt01

@@ -16,6 +16,7 @@ ktherm = 0
 conduct = 'bubbly'
 kdyn  =  0
 basalstress = .true.
+use_bathymetry = .true.
 shortwave = 'ccsm3'
 albedo_type = 'constant'
 formdrag = .true.

configuration/scripts/options/set_nml.alt03

@@ -14,3 +14,4 @@ revised_evp = .false.
 Ktens          = 0.
 e_ratio        = 2.
 basalstress    = .true.
+use_bathymetry = .true.

configuration/scripts/options/set_nml.gx1

@@ -5,6 +5,7 @@ grid_format = 'bin'
 grid_type   = 'displaced_pole'
 grid_file   = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/grid_gx1.bin'
 kmt_file    = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/kmt_gx1.bin'
+bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx1/global_gx1.bathy.nc'
 maskhalo_dyn   = .true.
 maskhalo_remap = .true.
 maskhalo_bound = .true.

configuration/scripts/options/set_nml.gx3

@@ -5,6 +5,7 @@ grid_format = 'bin'
 grid_type   = 'displaced_pole'
 grid_file   = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/grid_gx3.bin'
 kmt_file    = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/kmt_gx3.bin'
+bathymetry_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/global_gx3.bathy.nc'
 fyear_init      = 1997
 atm_data_format = 'bin'
 atm_data_type   = 'ncar'

doc/source/user_guide/ug_case_settings.rst

@@ -151,6 +151,8 @@ Table of namelist options
    "", "``tripole``", "read from file in *popgrid*", ""
    "", "``regional``", "read from file in *popgrid*", ""
    "``grid_file``", "filename", "name of grid file to be read", "‘grid’"
+   "``bathymetry_file``", "filename", "name of bathymetry file to be read", "‘grid’"
+   "``use_bathymetry``", "true/false", "use read in bathymetry file for basalstress option", ""
    "``kmt_file``", "filename", "name of land mask file to be read", "‘kmt’"
    "``gridcpl_file``", "filename", "input file for coupling grid info", ""
    "``kcatbound``", "``0``", "original category boundary formula", "0"
@@ -255,6 +257,7 @@ Table of namelist options
    "", "``2``", "EAP dynamics", ""
    "", "``1``", "dynamics ON", ""
    "``revised_evp``", "true/false", "use revised EVP formulation", ""
+   "``basalstress``", "true/false", "use bottom basal stress", ""
    "``ndte``", "integer", "number of EVP subcycles", "120"
    "``advection``", "``remap``", "linear remapping advection", "‘remap’"
    "", "``upwind``", "donor cell advection", ""