Skip to content

Commit

Permalink
Bug fixes in ice_forcing.F90 to properly read and rotate ocean currents
Browse files Browse the repository at this point in the history 
from 3D netcdf data (gx1, oceanmixed_ice_depth.nc).  Thanks to Adrian 
Turner for spotting these.
  • Loading branch information
@eclare108213
eclare108213 committed Apr 1, 2015
1 parent 0b786dd commit 223e60f
Show file tree
Hide file tree
Showing 2 changed files with 150 additions and 7 deletions.
20 changes: 13 additions & 7 deletions cice/source/ice_forcing.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3202,13 +3202,15 @@ subroutine ocn_data_ncar_init_3D
field_loc_center, field_type_scalar
use ice_domain_size, only: max_blocks
use ice_grid, only: to_ugrid, ANGLET
use ice_read_write, only: ice_read_nc_uv
#ifdef ncdf
use netcdf
#endif

integer (kind=int_kind) :: &
n , & ! field index
m ! month index
m , & ! month index
nzlev ! z level of currents

character(char_len) :: &
vname(nfld) ! variable names to search for in file
Expand Down Expand Up @@ -3276,8 +3278,14 @@ subroutine ocn_data_ncar_init_3D
do m=1,12

! Note: netCDF does single to double conversion if necessary
call ice_read_nc(fid, m, vname(n), work1, dbug, &
field_loc_center, field_type_scalar)
if (n == 4 .or. n == 5) then ! 3D currents
nzlev = 1 ! surface currents
call ice_read_nc_uv(fid, m, nzlev, vname(n), work1, dbug, &
field_loc_center, field_type_scalar)
else
call ice_read_nc(fid, m, vname(n), work1, dbug, &
field_loc_center, field_type_scalar)
endif

! the land mask used in ocean_mixed_depth.nc does not
! match our gx1v3 mask (hm)
Expand All @@ -3301,10 +3309,8 @@ subroutine ocn_data_ncar_init_3D
ocn_frc_m(:,:,:,n+1,m) = work2(:,:,:)*cos(ANGLET(:,:,:)) &
- work1(:,:,:)*sin(ANGLET(:,:,:))

work1(:,:,:) = ocn_frc_m(:,:,:,n ,m)*cos(ANGLET(:,:,:)) &
+ ocn_frc_m(:,:,:,n+1,m)*sin(ANGLET(:,:,:))
work2(:,:,:) = ocn_frc_m(:,:,:,n+1,m)*cos(ANGLET(:,:,:)) &
- ocn_frc_m(:,:,:,n ,m)*sin(ANGLET(:,:,:))
work1(:,:,:) = ocn_frc_m(:,:,:,n ,m)
work2(:,:,:) = ocn_frc_m(:,:,:,n+1,m)
call to_ugrid(work1,ocn_frc_m(:,:,:,n ,m))
call to_ugrid(work2,ocn_frc_m(:,:,:,n+1,m))

Expand Down
137 changes: 137 additions & 0 deletions cice/source/ice_read_write.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ module ice_read_write
ice_read_nc, &
ice_read_global, &
ice_read_global_nc, &
ice_read_nc_uv, &
ice_write, &
ice_write_nc, &
ice_write_ext, &
Expand Down Expand Up @@ -1936,6 +1937,142 @@ subroutine ice_close_nc(fid)

end subroutine ice_close_nc

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

! Read a netCDF file and scatter to processors.
! If the optional variables field_loc and field_type are present,
! the ghost cells are filled using values from the global array.
! This prevents them from being filled with zeroes in land cells
! (subroutine ice_HaloUpdate need not be called).
!
! Adapted by Elizabeth Hunke for reading 3D ocean currents

subroutine ice_read_nc_uv(fid, nrec, nzlev, varname, work, diag, &
field_loc, field_type, restart_ext)

use ice_fileunits, only: nu_diag
use ice_gather_scatter, only: scatter_global, scatter_global_ext

integer (kind=int_kind), intent(in) :: &
fid , & ! file id
nrec , & ! record number
nzlev ! z level

logical (kind=log_kind), intent(in) :: &
diag ! if true, write diagnostic output

character (len=*), intent(in) :: &
varname ! field name in netcdf file

real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks), &
intent(out) :: &
work ! output array (real, 8-byte)

logical (kind=log_kind), optional, intent(in) :: &
restart_ext ! if true, read extended grid

integer (kind=int_kind), optional, intent(in) :: &
field_loc, & ! location of field on staggered grid
field_type ! type of field (scalar, vector, angle)

! local variables

#ifdef ncdf
! netCDF file diagnostics:
integer (kind=int_kind) :: &
varid , & ! variable id
status, & ! status output from netcdf routines
ndim, nvar, & ! sizes of netcdf file
id, & ! dimension index
dimlen ! size of dimension

real (kind=dbl_kind) :: &
amin, amax, asum ! min, max values and sum of input array

character (char_len) :: &
dimname ! dimension name

real (kind=dbl_kind), dimension(:,:), allocatable :: &
work_g1

integer (kind=int_kind) :: nx, ny

nx = nx_global
ny = ny_global

if (present(restart_ext)) then
if (restart_ext) then
nx = nx_global + 2*nghost
ny = ny_global + 2*nghost
endif
endif

if (my_task == master_task) then
allocate(work_g1(nx,ny))
else
allocate(work_g1(1,1)) ! to save memory
endif

if (my_task == master_task) then

!-------------------------------------------------------------
! Find out ID of required variable
!-------------------------------------------------------------

status = nf90_inq_varid(fid, trim(varname), varid)

if (status /= nf90_noerr) then
call abort_ice ( &
'ice_read_nc_xy: Cannot find variable '//trim(varname) )
endif

!--------------------------------------------------------------
! Read global array
!--------------------------------------------------------------

status = nf90_get_var( fid, varid, work_g1, &
start=(/1,1,nzlev,nrec/), &
count=(/nx,ny,1,1/) )

endif ! my_task = master_task

!-------------------------------------------------------------------
! optional diagnostics
!-------------------------------------------------------------------

if (my_task==master_task .and. diag) then
amin = minval(work_g1)
amax = maxval(work_g1, mask = work_g1 /= spval_dbl)
asum = sum (work_g1, mask = work_g1 /= spval_dbl)
write(nu_diag,*) ' min, max, sum =', amin, amax, asum
endif

!-------------------------------------------------------------------
! Scatter data to individual processors.
! NOTE: Ghost cells are not updated unless field_loc is present.
!-------------------------------------------------------------------

if (present(restart_ext)) then
if (restart_ext) then
call scatter_global_ext(work, work_g1, master_task, distrb_info)
endif
else
if (present(field_loc)) then
call scatter_global(work, work_g1, master_task, distrb_info, &
field_loc, field_type)
else
call scatter_global(work, work_g1, master_task, distrb_info, &
field_loc_noupdate, field_type_noupdate)
endif
endif

deallocate(work_g1)

#else
work = c0 ! to satisfy intent(out) attribute
#endif
end subroutine ice_read_nc_uv

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

end module ice_read_write
Expand Down

0 comments on commit 223e60f

Please sign in to comment.