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

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

Merged
merged 3 commits into from
Oct 27, 2023
Merged

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

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
ee03042
Split N/E grid computation out of Tlonlat, create NElonlat subroutine.
apcraig Oct 26, 2023
80e88ae
Add gx3 grid/kmt files with TLON, TLAT, ANGLET netcdf grid test.
apcraig Oct 27, 2023
a219ed9
Merge branch 'main' of https://github.com/cice-consortium/cice into n…
apcraig Oct 27, 2023
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
191 changes: 124 additions & 67 deletions cicecore/cicedyn/infrastructure/ice_grid.F90
View file
Open in desktop
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

Comment on lines +2369 to +2372
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

did we want to keep this? it looks like a debugging print...

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
Comment on lines -2584 to +2648
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't these prints logically be in Tlatlon instead of NElatlon ?

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
Open in desktop
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we already had NetCDF versions of the gx3 grid, so this file could have been named set_nml.gx3t to be in line with the name of the new files.

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
Open in desktop
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