Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Split bathymetry subroutines into seperate file #978

Draft
wants to merge 3 commits into
base: main
Choose a base branch
from
Draft

Split bathymetry subroutines into seperate file #978

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e6a5c51
Split read_seabedstress_bathy into seperate ice_grid_bathy.f90
anton-seaice Oct 2, 2024
13381f9
move get_bathymetry to ice_grid_bathy
anton-seaice Oct 2, 2024
bad36f7
move get_bathymetry_popfile to ice_grid_bathy.f90
anton-seaice Oct 2, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
243 changes: 3 additions & 240 deletions cicecore/cicedyn/infrastructure/ice_grid.F90 → .../cicedyn/infrastructure/grid/ice_grid.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@ module ice_grid
use ice_global_reductions, only: global_minval, global_maxval
use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
use icepack_intfc, only: icepack_query_parameters, icepack_init_parameters
use ice_grid_bathy

implicit none
private
Expand Down Expand Up @@ -682,9 +683,9 @@ subroutine init_grid2
!-----------------------------------------------------------------

if (trim(bathymetry_format) == 'default') then
call get_bathymetry
call get_bathymetry(use_bathymetry, bathymetry_file, bathymetry, kmt)
elseif (trim(bathymetry_format) == 'pop') then
call get_bathymetry_popfile
call get_bathymetry_popfile(use_bathymetry, bathymetry_file, bathymetry, kmt)
else
call abort_ice(subname//' ERROR: bathymetry_format value must be default or pop', &
file=__FILE__, line=__LINE__)
Expand Down Expand Up @@ -4432,244 +4433,6 @@ subroutine gridbox_verts(work_g,vbounds)

end subroutine gridbox_verts

!=======================================================================
! ocean bathymetry for grounded sea ice (seabed stress) or icebergs
! currently hardwired for 40 levels (gx3, gx1 grids)
! should be read from a file instead (see subroutine read_seabedstress_bathy)

subroutine get_bathymetry

integer (kind=int_kind) :: &
i, j, k, iblk ! loop indices

integer (kind=int_kind), parameter :: &
nlevel = 40 ! number of layers (gx3 grid)

real (kind=dbl_kind), dimension(nlevel) :: &
depth ! total depth, m

real (kind=dbl_kind) :: &
puny

logical (kind=log_kind) :: &
calc_dragio

real (kind=dbl_kind), dimension(nlevel), parameter :: &
thick = (/ & ! ocean layer thickness, m
10.01244_dbl_kind, 10.11258_dbl_kind, 10.31682_dbl_kind, &
10.63330_dbl_kind, 11.07512_dbl_kind, 11.66145_dbl_kind, &
12.41928_dbl_kind, 13.38612_dbl_kind, 14.61401_dbl_kind, &
16.17561_dbl_kind, 18.17368_dbl_kind, 20.75558_dbl_kind, &
24.13680_dbl_kind, 28.63821_dbl_kind, 34.74644_dbl_kind, &
43.20857_dbl_kind, 55.16812_dbl_kind, 72.30458_dbl_kind, &
96.74901_dbl_kind, 130.0392_dbl_kind, 170.0489_dbl_kind, &
207.9933_dbl_kind, 233.5694_dbl_kind, 245.2719_dbl_kind, &
248.9804_dbl_kind, 249.8322_dbl_kind, 249.9787_dbl_kind, &
249.9979_dbl_kind, 249.9998_dbl_kind, 250.0000_dbl_kind, &
250.0000_dbl_kind, 250.0000_dbl_kind, 250.0000_dbl_kind, &
250.0000_dbl_kind, 250.0000_dbl_kind, 250.0000_dbl_kind, &
250.0000_dbl_kind, 250.0000_dbl_kind, 250.0000_dbl_kind, &
250.0000_dbl_kind /)

character(len=*), parameter :: subname = '(get_bathymetry)'

call icepack_query_parameters(puny_out=puny, calc_dragio_out=calc_dragio)
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
file=__FILE__, line=__LINE__)

if (use_bathymetry) then

call read_seabedstress_bathy

else

! convert to total depth
depth(1) = thick(1)
do k = 2, nlevel
depth(k) = depth(k-1) + thick(k)
enddo

bathymetry = c0
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
k = min(nint(kmt(i,j,iblk)),nlevel)
if (k > nlevel) call abort_ice(subname//' kmt gt nlevel error', &
file=__FILE__, line=__LINE__)
if (k > 0) bathymetry(i,j,iblk) = depth(k)
enddo
enddo
enddo

! For consistency, set thickness_ocn_layer1 in Icepack if 'calc_dragio' is active
if (calc_dragio) then
call icepack_init_parameters(thickness_ocn_layer1_in=thick(1))
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
file=__FILE__, line=__LINE__)
endif

endif ! bathymetry_file

end subroutine get_bathymetry

!=======================================================================
! with use_bathymetry = false, vertical depth profile generated for max KMT
! with use_bathymetry = true, expects to read in pop vert_grid file

subroutine get_bathymetry_popfile

integer (kind=int_kind) :: &
i, j, k, iblk ! loop indices

integer (kind=int_kind) :: &
ntmp, nlevel , & ! number of levels (max KMT)
k1 , & ! levels
ierr , & ! error tag
fid ! fid unit number

real (kind=dbl_kind), dimension(:),allocatable :: &
depth , & ! total depth, m
thick ! layer thickness, cm -> m

logical (kind=log_kind) :: &
calc_dragio

character(len=*), parameter :: subname = '(get_bathymetry_popfile)'

ntmp = maxval(nint(KMT))
nlevel = global_maxval(ntmp,distrb_info)

if (my_task==master_task) then
write(nu_diag,*) subname,' KMT max = ',nlevel
endif

allocate(depth(nlevel),thick(nlevel))
thick = -999999.
depth = -999999.

if (use_bathymetry) then

write (nu_diag,*) subname,' Bathymetry file = ', trim(bathymetry_file)
if (my_task == master_task) then
call get_fileunit(fid)
open(fid,file=bathymetry_file,form='formatted',iostat=ierr)
if (ierr/=0) call abort_ice(subname//' open error', file=__FILE__, line=__LINE__)
do k = 1,nlevel
read(fid,*,iostat=ierr) thick(k)
if (ierr/=0) call abort_ice(subname//' read error', file=__FILE__, line=__LINE__)
enddo
call release_fileunit(fid)
endif

call broadcast_array(thick,master_task)

else

! create thickness profile
k1 = min(5,nlevel)
do k = 1,k1
thick(k) = max(10000._dbl_kind/float(nlevel),500._dbl_kind)
enddo
do k = k1+1,nlevel
thick(k) = min(thick(k-1)*1.2_dbl_kind,20000._dbl_kind)
enddo

endif

! convert thick from cm to m
thick = thick / 100._dbl_kind

! convert to total depth
depth(1) = thick(1)
do k = 2, nlevel
depth(k) = depth(k-1) + thick(k)
if (depth(k) < 0.) call abort_ice(subname//' negative depth error', file=__FILE__, line=__LINE__)
enddo

if (my_task==master_task) then
do k = 1,nlevel
write(nu_diag,'(2a,i6,2f13.7)') subname,' k, thick(m), depth(m) = ',k,thick(k),depth(k)
enddo
endif

bathymetry = 0._dbl_kind
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
k = nint(kmt(i,j,iblk))
if (k > nlevel) call abort_ice(subname//' kmt gt nlevel error', file=__FILE__, line=__LINE__)
if (k > 0) bathymetry(i,j,iblk) = depth(k)
enddo
enddo
enddo

! For consistency, set thickness_ocn_layer1 in Icepack if 'calc_dragio' is active
call icepack_query_parameters(calc_dragio_out=calc_dragio)
if (calc_dragio) then
call icepack_init_parameters(thickness_ocn_layer1_in=thick(1))
endif
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
file=__FILE__, line=__LINE__)

deallocate(depth,thick)

end subroutine get_bathymetry_popfile

!=======================================================================

! Read bathymetry data for seabed stress calculation (grounding scheme for
! landfast ice) in CICE stand-alone mode. When CICE is in coupled mode
! (e.g. CICE-NEMO), hwater should be uptated at each time level so that
! it varies with ocean dynamics.
!
! author: Fred Dupont, CMC

subroutine read_seabedstress_bathy

! use module
use ice_read_write

! local variables
integer (kind=int_kind) :: &
fid_init ! file id for netCDF init file

character (char_len_long) :: & ! input data file names
fieldname

logical (kind=log_kind) :: diag=.true.

character(len=*), parameter :: subname = '(read_seabedstress_bathy)'

if (my_task == master_task) then
write (nu_diag,*) ' '
write (nu_diag,*) 'Bathymetry file: ', trim(bathymetry_file)
call icepack_warnings_flush(nu_diag)
endif

call ice_open_nc(bathymetry_file,fid_init)

fieldname='Bathymetry'

if (my_task == master_task) then
write(nu_diag,*) subname,' reading ',TRIM(fieldname)
call icepack_warnings_flush(nu_diag)
endif
call ice_read_nc(fid_init,1,fieldname,bathymetry,diag, &
field_loc=field_loc_center, &
field_type=field_type_scalar)

call ice_close_nc(fid_init)

if (my_task == master_task) then
write(nu_diag,*) subname,' closing file ',TRIM(bathymetry_file)
call icepack_warnings_flush(nu_diag)
endif

end subroutine read_seabedstress_bathy

!=======================================================================

end module ice_grid
Expand Down
Loading
Loading