Skip to content

Commit

Permalink
Split N/E grid computation out of Tlonlat, create NElonlat subroutine. (
Browse files Browse the repository at this point in the history 
CICE-Consortium#899)

* Split N/E grid computation out of Tlonlat, create NElonlat subroutine.

See CICE-Consortium#897

When TLON, TLAT, ANGLET are on the CICE grid, Tlonlat is NOT called.  This
meant N and E grid info was never computed.  This would fail during history
writing with invalid values in N and E grid arrays.  And it would also
cause problem if the C-grid were run with this type of CICE grid.

There are no test grids that have TLON, TLAT, ANGLET on them, so this
error was not found in standard test suites.  This was detected by
users.

* Add gx3 grid/kmt files with TLON, TLAT, ANGLET netcdf grid test.

The grid and kmt files were produced from a gx3 history file.  Results
are not bit-for-bit with the standard gx3 runs, but seem to be roundoff
different initially (as expected).
  • Loading branch information
@apcraig
apcraig authored Oct 27, 2023
1 parent 0b5ca09 commit 0484dcd
Show file tree
Hide file tree
Showing 3 changed files with 129 additions and 67 deletions.
191 changes: 124 additions & 67 deletions cicecore/cicedyn/infrastructure/ice_grid.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -675,36 +675,37 @@ subroutine init_grid2
if (trim(grid_type) == 'cpom_grid') then
ANGLET(:,:,:) = ANGLE(:,:,:)
else if (.not. (l_readCenter)) then
ANGLET = c0
ANGLET = c0

!$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, &
!$OMP angle_0,angle_w,angle_s,angle_sw)
do iblk = 1, nblocks
this_block = get_block(blocks_ice(iblk),iblk)
ilo = this_block%ilo
ihi = this_block%ihi
jlo = this_block%jlo
jhi = this_block%jhi
!$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, &
!$OMP angle_0,angle_w,angle_s,angle_sw)
do iblk = 1, nblocks
this_block = get_block(blocks_ice(iblk),iblk)
ilo = this_block%ilo
ihi = this_block%ihi
jlo = this_block%jlo
jhi = this_block%jhi

do j = jlo, jhi
do i = ilo, ihi
angle_0 = ANGLE(i ,j ,iblk) ! w----0
angle_w = ANGLE(i-1,j ,iblk) ! | |
angle_s = ANGLE(i, j-1,iblk) ! | |
angle_sw = ANGLE(i-1,j-1,iblk) ! sw---s
ANGLET(i,j,iblk) = atan2(p25*(sin(angle_0)+ &
sin(angle_w)+ &
sin(angle_s)+ &
sin(angle_sw)),&
p25*(cos(angle_0)+ &
cos(angle_w)+ &
cos(angle_s)+ &
cos(angle_sw)))
enddo
do j = jlo, jhi
do i = ilo, ihi
angle_0 = ANGLE(i ,j ,iblk) ! w----0
angle_w = ANGLE(i-1,j ,iblk) ! | |
angle_s = ANGLE(i, j-1,iblk) ! | |
angle_sw = ANGLE(i-1,j-1,iblk) ! sw---s
ANGLET(i,j,iblk) = atan2(p25*(sin(angle_0)+ &
sin(angle_w)+ &
sin(angle_s)+ &
sin(angle_sw)),&
p25*(cos(angle_0)+ &
cos(angle_w)+ &
cos(angle_s)+ &
cos(angle_sw)))
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
!$OMP END PARALLEL DO
endif ! cpom_grid

if (trim(grid_type) == 'regional' .and. &
(.not. (l_readCenter))) then
! for W boundary extrapolate from interior
Expand Down Expand Up @@ -734,8 +735,10 @@ subroutine init_grid2

call makemask ! velocity mask, hemisphere masks
if (.not. (l_readCenter)) then
call Tlatlon ! get lat, lon on the T grid
call Tlatlon ! get lat, lon on the T grid
endif
call NElatlon ! get lat, lon on the N, E grid

!-----------------------------------------------------------------
! bathymetry
!-----------------------------------------------------------------
Expand Down Expand Up @@ -1961,8 +1964,8 @@ subroutine cpomgrid
close (nu_kmt)
endif

write(nu_diag,*) "min/max HTN: ", minval(HTN), maxval(HTN)
write(nu_diag,*) "min/max HTE: ", minval(HTE), maxval(HTE)
write(nu_diag,*) subname," min/max HTN: ", minval(HTN), maxval(HTN)
write(nu_diag,*) subname," min/max HTE: ", minval(HTE), maxval(HTE)

end subroutine cpomgrid

Expand Down Expand Up @@ -2363,17 +2366,17 @@ subroutine Tlatlon

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

if (my_task==master_task) then
write(nu_diag,*) subname,' called'
endif

call icepack_query_parameters(rad_to_deg_out=rad_to_deg)
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
file=__FILE__, line=__LINE__)

TLAT(:,:,:) = c0
TLON(:,:,:) = c0
NLAT(:,:,:) = c0
NLON(:,:,:) = c0
ELAT(:,:,:) = c0
ELON(:,:,:) = c0

!$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, &
!$OMP x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4, &
Expand Down Expand Up @@ -2426,15 +2429,87 @@ subroutine Tlatlon
! TLAT in radians North
TLAT(i,j,iblk) = asin(tz)

! these two loops should be merged to save cos/sin calculations,
! but atan2 is not bit-for-bit. This suggests the result for atan2 depends on
! the prior atan2 call ??? not sure what's going on.
#if (1 == 1)
enddo ! i
enddo ! j
enddo ! iblk
!$OMP END PARALLEL DO

if (trim(grid_type) == 'regional') then
! for W boundary extrapolate from interior
!$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
do iblk = 1, nblocks
this_block = get_block(blocks_ice(iblk),iblk)
ilo = this_block%ilo
ihi = this_block%ihi
jlo = this_block%jlo
jhi = this_block%jhi

i = ilo
if (this_block%i_glob(i) == 1) then
do j = jlo, jhi
TLON(i,j,iblk) = c2*TLON(i+1,j,iblk) - &
TLON(i+2,j,iblk)
TLAT(i,j,iblk) = c2*TLAT(i+1,j,iblk) - &
TLAT(i+2,j,iblk)
enddo
endif
enddo
!$OMP END PARALLEL DO
endif ! regional

call ice_timer_start(timer_bound)
call ice_HaloUpdate (TLON, halo_info, &
field_loc_center, field_type_scalar, &
fillValue=c1)
call ice_HaloUpdate (TLAT, halo_info, &
field_loc_center, field_type_scalar, &
fillValue=c1)
call ice_HaloExtrapolate(TLON, distrb_info, &
ew_boundary_type, ns_boundary_type)
call ice_HaloExtrapolate(TLAT, distrb_info, &
ew_boundary_type, ns_boundary_type)

end subroutine Tlatlon

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

! Initializes latitude and longitude on N, E grid
!
! author: T. Craig from Tlatlon

subroutine NElatlon

use ice_constants, only: c0, c1, c1p5, c2, c4, p5, &
field_loc_center, field_loc_Nface, field_loc_Eface, &
field_type_scalar

integer (kind=int_kind) :: &
i, j, iblk , & ! horizontal indices
ilo,ihi,jlo,jhi ! beginning and end of physical domain

real (kind=dbl_kind) :: &
z1,x1,y1,z2,x2,y2,z3,x3,y3,z4,x4,y4,tx,ty,tz,da, &
rad_to_deg

type (block) :: &
this_block ! block information for current block

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

if (my_task==master_task) then
write(nu_diag,*) subname,' called'
endif

call icepack_query_parameters(rad_to_deg_out=rad_to_deg)
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
file=__FILE__, line=__LINE__)

NLAT(:,:,:) = c0
NLON(:,:,:) = c0
ELAT(:,:,:) = c0
ELON(:,:,:) = c0

!$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, &
!$OMP x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4, &
!$OMP tx,ty,tz,da)
Expand Down Expand Up @@ -2467,7 +2542,7 @@ subroutine Tlatlon
x4 = cos(ULON(i,j,iblk))*z4
y4 = sin(ULON(i,j,iblk))*z4
z4 = sin(ULAT(i,j,iblk))
#endif

! ---------
! NLON/NLAT 2 pt computation (pts 3, 4)
! ---------
Expand Down Expand Up @@ -2522,10 +2597,6 @@ subroutine Tlatlon
i = ilo
if (this_block%i_glob(i) == 1) then
do j = jlo, jhi
TLON(i,j,iblk) = c2*TLON(i+1,j,iblk) - &
TLON(i+2,j,iblk)
TLAT(i,j,iblk) = c2*TLAT(i+1,j,iblk) - &
TLAT(i+2,j,iblk)
NLON(i,j,iblk) = c1p5*TLON(i+1,j,iblk) - &
p5*TLON(i+2,j,iblk)
NLAT(i,j,iblk) = c1p5*TLAT(i+1,j,iblk) - &
Expand All @@ -2537,12 +2608,6 @@ subroutine Tlatlon
endif ! regional

call ice_timer_start(timer_bound)
call ice_HaloUpdate (TLON, halo_info, &
field_loc_center, field_type_scalar, &
fillValue=c1)
call ice_HaloUpdate (TLAT, halo_info, &
field_loc_center, field_type_scalar, &
fillValue=c1)
call ice_HaloUpdate (NLON, halo_info, &
field_loc_Nface, field_type_scalar, &
fillValue=c1)
Expand All @@ -2555,10 +2620,6 @@ subroutine Tlatlon
call ice_HaloUpdate (ELAT, halo_info, &
field_loc_Eface, field_type_scalar, &
fillValue=c1)
call ice_HaloExtrapolate(TLON, distrb_info, &
ew_boundary_type, ns_boundary_type)
call ice_HaloExtrapolate(TLAT, distrb_info, &
ew_boundary_type, ns_boundary_type)
call ice_HaloExtrapolate(NLON, distrb_info, &
ew_boundary_type, ns_boundary_type)
call ice_HaloExtrapolate(NLAT, distrb_info, &
Expand All @@ -2581,12 +2642,10 @@ subroutine Tlatlon

if (my_task==master_task) then
write(nu_diag,*) ' '
! if (nx_block > 5+2*nghost .and. ny_block > 5+2*nghost) then
write(nu_diag,*) 'min/max ULON:', y1*rad_to_deg, y2*rad_to_deg
write(nu_diag,*) 'min/max ULAT:', y3*rad_to_deg, y4*rad_to_deg
! endif
write(nu_diag,*) 'min/max TLON:', x1*rad_to_deg, x2*rad_to_deg
write(nu_diag,*) 'min/max TLAT:', x3*rad_to_deg, x4*rad_to_deg
write(nu_diag,*) subname,' min/max ULON:', y1*rad_to_deg, y2*rad_to_deg
write(nu_diag,*) subname,' min/max ULAT:', y3*rad_to_deg, y4*rad_to_deg
write(nu_diag,*) subname,' min/max TLON:', x1*rad_to_deg, x2*rad_to_deg
write(nu_diag,*) subname,' min/max TLAT:', x3*rad_to_deg, x4*rad_to_deg
endif ! my_task

x1 = global_minval(NLON, distrb_info, nmask)
Expand All @@ -2601,15 +2660,13 @@ subroutine Tlatlon

if (my_task==master_task) then
write(nu_diag,*) ' '
! if (nx_block > 5+2*nghost .and. ny_block > 5+2*nghost) then
write(nu_diag,*) 'min/max NLON:', x1*rad_to_deg, x2*rad_to_deg
write(nu_diag,*) 'min/max NLAT:', x3*rad_to_deg, x4*rad_to_deg
write(nu_diag,*) 'min/max ELON:', y1*rad_to_deg, y2*rad_to_deg
write(nu_diag,*) 'min/max ELAT:', y3*rad_to_deg, y4*rad_to_deg
! endif
write(nu_diag,*) subname,' min/max NLON:', x1*rad_to_deg, x2*rad_to_deg
write(nu_diag,*) subname,' min/max NLAT:', x3*rad_to_deg, x4*rad_to_deg
write(nu_diag,*) subname,' min/max ELON:', y1*rad_to_deg, y2*rad_to_deg
write(nu_diag,*) subname,' min/max ELAT:', y3*rad_to_deg, y4*rad_to_deg
endif ! my_task

end subroutine Tlatlon
end subroutine NElatlon

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

Expand Down Expand Up @@ -4677,7 +4734,7 @@ subroutine read_seabedstress_bathy
fieldname='Bathymetry'

if (my_task == master_task) then
write(nu_diag,*) 'reading ',TRIM(fieldname)
write(nu_diag,*) subname,' reading ',TRIM(fieldname)
call icepack_warnings_flush(nu_diag)
endif
call ice_read_nc(fid_init,1,fieldname,bathymetry,diag, &
Expand All @@ -4687,7 +4744,7 @@ subroutine read_seabedstress_bathy
call ice_close_nc(fid_init)

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

Expand Down
3 changes: 3 additions & 0 deletions configuration/scripts/options/set_nml.gx3nc
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
grid_format = 'nc'
grid_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/grid_gx3t.nc'
kmt_file = 'ICE_MACHINE_INPUTDATA/CICE_data/grid/gx3/kmt_gx3t.nc'
2 changes: 2 additions & 0 deletions configuration/scripts/tests/base_suite.ts
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ smoke gx3 1x1 debug,diag1,run2day
smoke gx3 1x4 debug,diag1,run2day
smoke gx3 4x1 debug,diag1,run5day
restart gx3 8x2 debug
restart gx3 8x2 debug,gx3nc
smoke gx3 8x2 diag24,run1year,medium
smoke gx3 7x2 diag1,bigdiag,run1day,diagpt1
decomp gx3 4x2x25x29x5 none
Expand All @@ -14,6 +15,7 @@ restart gx1 40x4 droundrobin,medium
restart tx1 40x4 dsectrobin,medium
restart tx1 40x4 dsectrobin,medium,jra55do
restart gx3 4x4 none
restart gx3 4x4 gx3nc
restart gx3 10x4 maskhalo
restart gx3 6x2 alt01
restart gx3 8x2 alt02
Expand Down

0 comments on commit 0484dcd

Please sign in to comment.