| C384.mx025.tile[1-6].nc | the mapped ocean mask on the ATM tiles | used to create ATM ICs consistent with the
fractional grid
-
+
-
-* If the optional post-weights are generated, the following files will be produced in the output location:
-
+
+* If the optional post-weights are generated, the following files will be produced in the output location:
+
Optional post-weights files for 1/4deg
@@ -147,27 +150,29 @@ The exact list of files produced by the *cpld_gridgen.sh* script will vary depen
same tripole grid using bilinear mapping
| tripole.mx025.Ct.to.rect.[destination resolution].[bilinear][conserve].nc | the ESMF weights for mapping variables
on the center (Ct) stagger location on
the tripole grid to a rectilinear grid
with [destination resolution] using
either bilinear or conservative mapping
|
-
+
-* If a resolution other than 1/4 degree is used in *cpld_gridgen.sh*, the following file will be produced in the output location:
-
-
+* The following file will be produced in the output location to map 1/4 degree tripole values to a tripole grid of lower resolution.
+
-Output files for CICE6 IC creation at tripole destination resolution
+Output files for down-sampled IC creation at tripole destination resolution
| File name | Function
- |
|---|
| tripole.mx025.Ct.to.mx[destination resolution].Ct.neareststod.nc | the ESMF weights for mapping the 1/4 CICE ICs to
- a tripole [destination resolution] using nearest
source-to-destination mapping
+ | | tripole.mx025.Ct.to.mx[destination resolution].Ct.neareststod.nc | the ESMF weights for mapping 1/4 deg tripole ICs to
+ a lower tripole destination resolution using nearest
source-to-destination mapping
+ | | tripole.mx025.Ct.to.mx[destination resolution].Ct.bilinear.nc | the ESMF weights for mapping 1/4 deg tripole ICs to
+ a lower tripole destination resolution using
bilinear mapping
+ | | tripole.[destination resolution].Ct.to.[Cu][Cv][Bu].bilinear.nc | the ESMF weights for mapping downscaled IC values on a
+ tripole grid from Ct locations to the native stagger locations
|
-
+
-
+
* If run-time land mask changes for MOM6 are requested, the following file will be produced in the output location:
-
-
+
+
Output files for run-time modification of MOM6 land mask
| File name | Function
|
|---|
| ufs.[Default filename].nc | Topo-edits required for UFS application. These are appended to the existing default topo
edits file and implemented at run time with the parameter flag
``ALLOW_LANDMASK_CHANGES=true``. All files produced by the *cpld_gridgen.sh* will be
consistent with this run-time land mask.
|
-
diff --git a/sorc/cpld_gridgen.fd/gen_fixgrid.F90 b/sorc/cpld_gridgen.fd/gen_fixgrid.F90
index 36f118eed..f2fe0dc88 100644
--- a/sorc/cpld_gridgen.fd/gen_fixgrid.F90
+++ b/sorc/cpld_gridgen.fd/gen_fixgrid.F90
@@ -20,7 +20,7 @@ program gen_fixgrid
use grdvars
use inputnml
use gengrid_kinds, only: CL, CS, dbl_kind, real_kind, int_kind
- use angles, only: find_angq, find_ang
+ use angles, only: find_ang, find_angq, find_angchk
use vertices, only: fill_vertices, fill_bottom, fill_top
use mapped_mask, only: make_frac_land
use postwgts, only: make_postwgts
@@ -49,7 +49,7 @@ program gen_fixgrid
integer :: rc,ncid,id,xtype
integer :: i,j,k,i2,j2
- integer :: ii,jj
+ integer :: ii
integer :: localPet, nPet
logical :: fexist = .false.
@@ -147,9 +147,6 @@ program gen_fixgrid
if(xtype.eq. 6)wet4 = real(wet8,4)
- !print *,minval(wet8),maxval(wet8)
- !print *,minval(wet4),maxval(wet4)
-
!---------------------------------------------------------------------
! read the MOM6 depth file
!---------------------------------------------------------------------
@@ -201,17 +198,17 @@ program gen_fixgrid
! this modified topoedits file
!---------------------------------------------------------------------
- fsrc = trim(dirsrc)//'/'//trim(editsfile)
- if(editmask)fsrc = trim(dirout)//'/'//'ufs.'//trim(editsfile)
-
- if (trim(editsfile) /= 'none') then
- inquire(file=trim(fsrc),exist=fexist)
- if (.not. fexist) then
- print '(a)', 'Required topoedits file '//trim(fsrc)//' is missing '
- call abort()
- end if
+ fsrc = trim(dirsrc)//'/'//trim(editsfile)
+ if(editmask)fsrc = trim(dirout)//'/'//'ufs.'//trim(editsfile)
+
+ if (trim(editsfile) /= 'none') then
+ inquire(file=trim(fsrc),exist=fexist)
+ if (.not. fexist) then
+ print '(a)', 'Required topoedits file '//trim(fsrc)//' is missing '
+ call abort()
end if
- call apply_topoedits(fsrc)
+ end if
+ call apply_topoedits(fsrc)
!---------------------------------------------------------------------
! read MOM6 supergrid file
@@ -236,15 +233,9 @@ program gen_fixgrid
rc = nf90_get_var(ncid, id, dy)
rc = nf90_close(ncid)
- !print *,'super grid size ',size(y,1),size(y,2)
- !print *,'max lat in super grid ',maxval(y)
sg_maxlat = maxval(y)
-
- !---------------------------------------------------------------------
- ! find the angle on corners---this requires the supergrid
- !---------------------------------------------------------------------
-
- call find_angq
+ write(logmsg,'(a,f12.2)')'max lat in super grid ',maxval(y)
+ print '(a)',trim(logmsg)
!---------------------------------------------------------------------
! fill grid variables
@@ -256,8 +247,6 @@ program gen_fixgrid
!deg->rad
ulon(i,j) = x(i2,j2)*deg2rad
ulat(i,j) = y(i2,j2)*deg2rad
- !in rad already
- angle(i,j) = -angq(i2,j2)
!m->cm
htn(i,j) = (dx(i2-1,j2) + dx(i2,j2))*100._dbl_kind
hte(i,j) = (dy(i2,j2-1) + dy(i2,j2))*100._dbl_kind
@@ -280,12 +269,63 @@ program gen_fixgrid
enddo
!---------------------------------------------------------------------
- ! find the angle on centers---this does not requires the supergrid
+ ! locate the ith index of the two poles on j=nj
+ ! the corner points must lie on the pole
+ !---------------------------------------------------------------------
+
+ ipole = -1
+ j = nj
+ do i = 1,ni/2
+ if(latBu(i,j) .eq. sg_maxlat)ipole(1) = i
+ enddo
+ do i = ni/2+1,ni
+ if(latBu(i,j) .eq. sg_maxlat)ipole(2) = i
+ enddo
+ write(logmsg,'(a,2i6,2f12.2)')'poles found at i = ',ipole, latBu(ipole(1),nj), &
+ latBu(ipole(2),nj)
+ print '(a)',trim(logmsg)
+
+ !---------------------------------------------------------------------
+ ! find the angle on centers using the same procedure as MOM6
!---------------------------------------------------------------------
- call find_ang
- print *,'ANGLET ',minval(anglet),maxval(anglet)
- print *,'ANGLE ',minval(angle),maxval(angle)
+ call find_ang((/1,ni/),(/1,nj/),lonBu,latBu,lonCt,anglet)
+ write(logmsg,'(a,2f12.2)')'ANGLET min,max: ',minval(anglet),maxval(anglet)
+ print '(a)',trim(logmsg)
+ write(logmsg,'(a,2f12.2)')'ANGLET edges i=1,i=ni: ',anglet(1,nj),anglet(ni,nj)
+ print '(a)',trim(logmsg)
+
+ xangCt(:) = 0.0
+ do i = 1,ni
+ i2 = ipole(2)+(ipole(1)-i)+1
+ xangCt(i) = -anglet(i2,nj) ! angle changes sign across seam
+ end do
+
+ !---------------------------------------------------------------------
+ ! find the angle on corners using the same procedure as CICE6
+ !---------------------------------------------------------------------
+
+ call find_angq((/1,ni/),(/1,nj/),xangCt,anglet,angle)
+ angle(ni,:) = -angle(1,:)
+ ! reverse angle for CICE
+ angle = -angle
+ write(logmsg,'(a,2f12.2)')'ANGLE min,max: ',minval(angle),maxval(angle)
+ print '(a)',trim(logmsg)
+ write(logmsg,'(a,2f12.2)')'ANGLE edges i=1,i=ni: ',angle(1,nj),angle(ni,nj)
+ print '(a)',trim(logmsg)
+
+ !---------------------------------------------------------------------
+ ! check the Bu angle
+ !---------------------------------------------------------------------
+
+ call find_angchk((/1,ni/),(/1,nj/),angle,angchk)
+ angchk(1,:) = -angchk(ni,:)
+ ! reverse angle for MOM6
+ angchk = -angchk
+ write(logmsg,'(a,2f12.2)')'ANGCHK min,max: ',minval(angchk),maxval(angchk)
+ print '(a)',trim(logmsg)
+ write(logmsg,'(a,2f12.2)')'ANGCHK edges i=1,i=ni: ',angchk(1,nj),angchk(ni,nj)
+ print '(a)',trim(logmsg)
!---------------------------------------------------------------------
! For the 1/4deg grid, hte at j=720 and j = 1440 is identically=0.0 for
@@ -295,8 +335,7 @@ program gen_fixgrid
! hte < 1.0
!---------------------------------------------------------------------
- write(logmsg,'(a,2e12.5)')'min vals of hte at folds ', &
- minval(hte(ni/2,:)),minval(hte(ni,:))
+ write(logmsg,'(a,2e12.5)')'min vals of hte at folds ', minval(hte(ni/2,:)),minval(hte(ni,:))
print '(a)',trim(logmsg)
do j = 1,nj
ii = ni/2
@@ -304,37 +343,13 @@ program gen_fixgrid
ii = ni
if(hte(ii,j) .le. 1.0)hte(ii,j) = 0.5*(hte(ii-1,j) + hte( 1,j))
enddo
- write(logmsg,'(a,2e12.5)')'min vals of hte at folds ', &
- minval(hte(ni/2,:)),minval(hte(ni,:))
+ write(logmsg,'(a,2e12.5)')'min vals of hte at folds ', minval(hte(ni/2,:)),minval(hte(ni,:))
print '(a)',trim(logmsg)
!---------------------------------------------------------------------
- !
+ ! find required extended values for setting all vertices
!---------------------------------------------------------------------
- where(lonCt .lt. 0.0)lonCt = lonCt + 360._dbl_kind
- where(lonCu .lt. 0.0)lonCu = lonCu + 360._dbl_kind
- where(lonCv .lt. 0.0)lonCv = lonCv + 360._dbl_kind
- where(lonBu .lt. 0.0)lonBu = lonBu + 360._dbl_kind
-
- !---------------------------------------------------------------------
- ! some basic error checking
- ! find the i-th index of the poles at j= nj
- ! the corner points must lie on the pole
- !---------------------------------------------------------------------
-
- ipole = -1
- j = nj
- do i = 1,ni/2
- if(latBu(i,j) .eq. sg_maxlat)ipole(1) = i
- enddo
- do i = ni/2+1,ni
- if(latBu(i,j) .eq. sg_maxlat)ipole(2) = i
- enddo
- write(logmsg,'(a,2i6,2f12.2)')'poles found at i = ',ipole,latBu(ipole(1),nj), &
- latBu(ipole(2),nj)
- print '(a)',trim(logmsg)
-
if(debug)call checkseam
do i = 1,ni
@@ -435,7 +450,7 @@ program gen_fixgrid
write(form1,'(a)')'('//trim(cnx)//'f14.8)'
write(form2,'(a)')'('//trim(cnx)//'i2)'
- allocate(ww3mask(1:ni,1:nj)); ww3mask = wet4
+ allocate(ww3mask(1:ni,1:nj)); ww3mask = int(wet4)
allocate(ww3dpth(1:ni,1:nj)); ww3dpth = dp4
where(latCt .ge. maximum_lat)ww3mask = 3
@@ -504,7 +519,17 @@ program gen_fixgrid
fwgt = trim(dirout)//'/'//'tripole.mx025.Ct.to.mx'//trim(res)//'.Ct.neareststod.nc'
logmsg = 'creating weight file '//trim(fwgt)
print '(a)',trim(logmsg)
+ call ESMF_RegridWeightGen(srcFile=trim(fsrc),dstFile=trim(fdst), &
+ weightFile=trim(fwgt), regridmethod=method, &
+ ignoreDegenerate=.true., unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)
+ method=ESMF_REGRIDMETHOD_BILINEAR
+ fdst = trim(dirout)//'/'//'Ct.mx'//trim(res)//'_SCRIP.nc'
+ fwgt = trim(dirout)//'/'//'tripole.mx025.Ct.to.mx'//trim(res)//'.Ct.bilinear.nc'
+ logmsg = 'creating weight file '//trim(fwgt)
+ print '(a)',trim(logmsg)
call ESMF_RegridWeightGen(srcFile=trim(fsrc),dstFile=trim(fdst), &
weightFile=trim(fwgt), regridmethod=method, &
ignoreDegenerate=.true., unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, rc=rc)
@@ -517,6 +542,19 @@ program gen_fixgrid
end if
end if
+ ! tripole Ct->tripole Bu for CICE are only for CICE IC creation
+ fsrc = trim(dirout)//'/'//'Ct.mx'//trim(res)//'_SCRIP.nc'
+ fdst = trim(dirout)//'/'//'Bu.mx'//trim(res)//'_SCRIP.nc'
+ fwgt = trim(dirout)//'/'//'tripole.mx'//trim(res)//'.Ct.to.Bu.bilinear.nc'
+ logmsg = 'creating weight file '//trim(fwgt)
+ print '(a)',trim(logmsg)
+
+ call ESMF_RegridWeightGen(srcFile=trim(fsrc),dstFile=trim(fdst), &
+ weightFile=trim(fwgt), regridmethod=method, &
+ ignoreDegenerate=.true., unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)
+
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
@@ -527,8 +565,8 @@ program gen_fixgrid
! clean up
!---------------------------------------------------------------------
- deallocate(x,y, angq, dx, dy, xsgp1, ysgp1)
- deallocate(areaCt, anglet, angle)
+ deallocate(x, y, dx, dy)
+ deallocate(areaCt, anglet, angle, angchk)
deallocate(latCt, lonCt)
deallocate(latCv, lonCv)
deallocate(latCu, lonCu)
diff --git a/sorc/cpld_gridgen.fd/grdvars.F90 b/sorc/cpld_gridgen.fd/grdvars.F90
index 578ac269a..9c549a494 100644
--- a/sorc/cpld_gridgen.fd/grdvars.F90
+++ b/sorc/cpld_gridgen.fd/grdvars.F90
@@ -68,17 +68,10 @@ module grdvars
! Super-grid source grid variables
real(dbl_kind), allocatable, dimension(:,:) :: x !< The longitudes of the MOM6 supergrid
real(dbl_kind), allocatable, dimension(:,:) :: y !< The latitudes of the MOM6 supergrid
- real(dbl_kind), allocatable, dimension(:,:) :: angq !< The grid rotation angle at the Bu (or corner)
- !! grid points of the super grid
-
real(dbl_kind), allocatable, dimension(:,:) :: dx !< The grid cell width in meters of the supergrid
!! in the x-direction (i-dimension)
real(dbl_kind), allocatable, dimension(:,:) :: dy !< The grid cell width in meters of the supergrid
!! in the y-direction (j-dimension)
- real(dbl_kind), allocatable, dimension(:,:) :: xsgp1 !< The longitudes of the super-grid replicated
- !! across the tripole seam
- real(dbl_kind), allocatable, dimension(:,:) :: ysgp1 !< The latitudes of the super-grid replicated
- !! across the tripole seam
! Output grid variables
real(dbl_kind), allocatable, dimension(:,:) :: latCt !< The latitude of the center (tracer) grid points
@@ -103,6 +96,8 @@ module grdvars
real(dbl_kind), allocatable, dimension(:,:) :: anglet !< The rotation angle on Ct points (opposite sense
!! from angle)
real(dbl_kind), allocatable, dimension(:,:) :: angle !< The rotation angle on Bu points
+ real(dbl_kind), allocatable, dimension(:,:) :: angchk !< The rotation angle on Ct points, as calculated by
+ !! CICE internally using angle on Bu
real(dbl_kind), allocatable, dimension(:,:,:) :: latCt_vert !< The latitudes of the 4 vertices of each Ct grid
!! point
@@ -129,6 +124,9 @@ module grdvars
!! opposite side of the tripole seam
real(dbl_kind), allocatable, dimension(:) :: xlatCt !< The latitude of the Ct grid points on the
!! opposite side of the tripole seam
+ real(dbl_kind), allocatable, dimension(:) :: xangCt !< The rotation angle on the Ct grid points on the
+ !! opposite side of the tripole seam
+
real(dbl_kind), allocatable, dimension(:) :: xlonCu !< The longitude of the Cu grid points on the
!! opposite side of the tripole seam
real(dbl_kind), allocatable, dimension(:) :: xlatCu !< The latitude of the Cu grid points on the
@@ -172,24 +170,22 @@ module grdvars
subroutine allocate_all
- allocate( x(0:nx,0:ny), y(0:nx,0:ny), angq(0:nx,0:ny) )
+ allocate( x(0:nx,0:ny), y(0:nx,0:ny) )
allocate( dx(nx,0:ny), dy(0:nx,ny) )
- allocate( xsgp1(0:nx,0:ny+1), ysgp1(0:nx,0:ny+1) )
-
allocate( latCt(ni,nj), lonCt(ni,nj) )
allocate( latCv(ni,nj), lonCv(ni,nj) )
allocate( latCu(ni,nj), lonCu(ni,nj) )
allocate( latBu(ni,nj), lonBu(ni,nj) )
- allocate( areaCt(ni,nj), anglet(ni,nj), angle(ni,nj) )
+ allocate( areaCt(ni,nj), anglet(ni,nj), angle(ni,nj), angchk(ni,nj))
allocate( latCt_vert(ni,nj,nv), lonCt_vert(ni,nj,nv) )
allocate( latCv_vert(ni,nj,nv), lonCv_vert(ni,nj,nv) )
allocate( latCu_vert(ni,nj,nv), lonCu_vert(ni,nj,nv) )
allocate( latBu_vert(ni,nj,nv), lonBu_vert(ni,nj,nv) )
- allocate( xlonCt(ni), xlatCt(ni) )
+ allocate( xlonCt(ni), xlatCt(ni), xangCt(ni) )
allocate( xlonCu(ni), xlatCu(ni) )
allocate( dlatBu(ni), dlatCv(ni) )
diff --git a/sorc/cpld_gridgen.fd/inputnml.F90 b/sorc/cpld_gridgen.fd/inputnml.F90
index 6e363657b..96d86f4c3 100644
--- a/sorc/cpld_gridgen.fd/inputnml.F90
+++ b/sorc/cpld_gridgen.fd/inputnml.F90
@@ -26,7 +26,8 @@ subroutine read_inputnml(fname)
character(len=*), intent(in) :: fname
! local variables
- integer :: stderr, iounit, rc
+ integer :: iounit, rc
+ character(len=200) :: tmpstr
namelist /grid_nml/ ni, nj, dirsrc, dirout, fv3dir, topofile, editsfile, &
res, atmres, npx, editmask, debug, &
@@ -34,24 +35,25 @@ subroutine read_inputnml(fname)
! Check whether file exists.
inquire (file=trim(fname), iostat=rc)
-
if (rc /= 0) then
- write (stderr, '(3a)') 'Error: input file "', trim(fname), '" does not exist.'
- return
+ write (0, '(3a)') 'Error: input file "', trim(fname), '" does not exist.'
+ stop 1
end if
! Open and read Namelist file.
open (action='read', file=trim(fname), iostat=rc, newunit=iounit)
read (nml=grid_nml, iostat=rc, unit=iounit)
+ if (rc /= 0) then
+ backspace(iounit)
+ read(iounit,'(a)')tmpstr
+ write (6, '(a)') 'Error: invalid Namelist format '//trim(tmpstr)
+ stop 1
+ end if
+ close(iounit)
! set supergrid dimensions
nx = ni*2
ny = nj*2
- if (rc /= 0) then
- write (stderr, '(a)') 'Error: invalid Namelist format.'
- end if
-
- close (iounit)
end subroutine read_inputnml
end module inputnml
diff --git a/sorc/cpld_gridgen.fd/topoedits.F90 b/sorc/cpld_gridgen.fd/topoedits.F90
index 6e1889474..e0b617b01 100644
--- a/sorc/cpld_gridgen.fd/topoedits.F90
+++ b/sorc/cpld_gridgen.fd/topoedits.F90
@@ -52,8 +52,8 @@ subroutine add_topoedits(fsrc,fdst)
! return the existing values
allocate(ieds1(cnt1)); ieds1 = 0
- allocate(jeds1(cnt1)); jeds2 = 0
- allocate(zeds1(cnt1)); zeds2 = 0.0
+ allocate(jeds1(cnt1)); jeds1 = 0
+ allocate(zeds1(cnt1)); zeds1 = 0.0
rc = nf90_open(fsrc, nf90_nowrite, ncid)
rc = nf90_inq_varid(ncid, 'iEdit', id)
diff --git a/sorc/cpld_gridgen.fd/tripolegrid.F90 b/sorc/cpld_gridgen.fd/tripolegrid.F90
index feb287101..c707833ca 100644
--- a/sorc/cpld_gridgen.fd/tripolegrid.F90
+++ b/sorc/cpld_gridgen.fd/tripolegrid.F90
@@ -13,7 +13,7 @@ module tripolegrid
use grdvars, only: lonCu,latCu,lonCu_vert,latCu_vert
use grdvars, only: lonCv,latCv,lonCv_vert,latCv_vert
use grdvars, only: lonBu,latBu,lonBu_vert,latBu_vert
- use grdvars, only: wet4,areaCt,angleT,dp4
+ use grdvars, only: wet4,areaCt,angleT,dp4,angle,angchk
use charstrings, only: logmsg,history
use vartypedefs, only: maxvars, fixvars, fixvars_typedefine
use netcdf
@@ -67,6 +67,12 @@ subroutine write_tripolegrid(fname)
!angleT
rc = nf90_def_var(ncid, 'anglet', nf90_double, dim2, id)
rc = nf90_put_att(ncid, id, 'units', 'radians')
+ !angle (angBu)
+ rc = nf90_def_var(ncid, 'angle', nf90_double, dim2, id)
+ rc = nf90_put_att(ncid, id, 'units', 'radians')
+ !angchk
+ rc = nf90_def_var(ncid, 'angchk', nf90_double, dim2, id)
+ rc = nf90_put_att(ncid, id, 'units', 'radians')
!bathymetry
rc = nf90_def_var(ncid, 'depth', nf90_float, dim2, id)
rc = nf90_put_att(ncid, id, 'units', 'm')
@@ -102,6 +108,12 @@ subroutine write_tripolegrid(fname)
rc = nf90_inq_varid(ncid,'anglet', id)
rc = nf90_put_var(ncid, id, anglet)
+ rc = nf90_inq_varid(ncid, 'angle', id)
+ rc = nf90_put_var(ncid, id, angle)
+
+ rc = nf90_inq_varid(ncid,'angchk', id)
+ rc = nf90_put_var(ncid, id, angchk)
+
rc = nf90_inq_varid(ncid, 'depth', id)
rc = nf90_put_var(ncid, id, dp4)
diff --git a/tests/cpld_gridgen/ftst_find_angq.F90 b/tests/cpld_gridgen/ftst_find_angq.F90
index bcddfb8eb..b24c6965e 100644
--- a/tests/cpld_gridgen/ftst_find_angq.F90
+++ b/tests/cpld_gridgen/ftst_find_angq.F90
@@ -1,108 +1,335 @@
-! Unit test for cpld_gridgen routine "find_angq".
-!
-! Reads a sample MOM6 supergrid and calculates the
-! rotation angle on corner points
-!
-! Author Denise Worthen 2/08/2022
+!> @file
+!! @brief unit test for angles
+!! @author Denise.Worthen@noaa.gov
+!!
+!! Given a 5x5 block of coords, calculate the angle, anglet and angchk
+!! and compare against values from MOM6 and CICE6 history output. The
+!! 3 5x5 blocks are for the 1deg tripole grid. The first two blocks contain
+!! the two polar points. The final block is located in the lower right quad
+!! of the tripole region (i=300:304,j=306:310). Because the angle calculations
+!! involve points outside of the 5x5 blocks (eg, i+1), not all 5x5 points are
+!! checked.
+!!
+!! @author Denise.Worthen@noaa.gov
program ftst_find_angq
-
- use netcdf
- use grdvars, only : ni,nj,nx,ny
- use grdvars, only : x,y,xsgp1,ysgp1,sg_maxlat
- use grdvars, only : angq
- use angles, only : find_angq
-
- implicit none
-
- integer :: i,j,i1, i2
- integer :: rc, ncid, id
-
- logical :: mastertask = .false.
- logical :: debug = .false.
-
- ! pole locations on SG
- integer :: ipolesg(2)
- ! unit test values
- real(kind=8) :: puny = 1.0e-12
- real(kind=8) :: delta(15)
- real(kind=8) :: sumdelta
-
- print *,"Starting test of cpld_gridgen routine find_angq"
-
- ! 1deg MOM6 dimensions
- ni = 360
- nj = 320
- ! super grid dimensions
- nx = 2*ni
- ny = 2*nj
-
- ! supergrid x,y and angles on corners
- allocate (x(0:nx,0:ny), y(0:nx,0:ny), angq(0:nx,0:ny))
- ! supergrid "plus 1" arrays
- allocate (xsgp1(0:nx,0:ny+1), ysgp1(0:nx,0:ny+1))
-
- !open the supergrid file and read the x,y coords
- rc = nf90_open('./data/ocean_hgrid.nc', nf90_nowrite, ncid)
- rc = nf90_inq_varid(ncid, 'x', id) !lon
- rc = nf90_get_var(ncid, id, x)
-
- rc = nf90_inq_varid(ncid, 'y', id) !lat
- rc = nf90_get_var(ncid, id, y)
- rc = nf90_close(ncid)
-
- ! max lat on supergrid
- sg_maxlat = maxval(y)
-
- !pole index on supergrid
- ipolesg = -1
- j = ny
- do i = 1,nx/2
- if(y(i,j) .eq. sg_maxlat)ipolesg(1) = i
- enddo
- do i = nx/2+1,nx
- if(y(i,j) .eq. sg_maxlat)ipolesg(2) = i
- enddo
-
- ! test angleq calculation
- call find_angq
-
- ! required for checking longitudes across seam
- where(xsgp1 .lt. 0.0)xsgp1 = xsgp1 + 360.0
-
- j = ny+1
- i1 = ipolesg(1); i2 = ipolesg(2)-(ipolesg(1)-i1)
- delta = 0.0
- ! check lons match across seam
- delta( 1) = xsgp1(i1-2,j)-xsgp1(i2+2,j)
- delta( 2) = xsgp1(i1-1,j)-xsgp1(i2+1,j)
- delta( 3) = xsgp1(i1, j)-xsgp1(i2, j)
- delta( 4) = xsgp1(i1+1,j)-xsgp1(i2-1,j)
- delta( 5) = xsgp1(i1+2,j)-xsgp1(i2-2,j)
- ! check lats match across seam
- delta( 6) = ysgp1(i1-2,j)-ysgp1(i2+2,j)
- delta( 7) = ysgp1(i1-1,j)-ysgp1(i2+1,j)
- delta( 8) = ysgp1(i1, j)-ysgp1(i2, j)
- delta( 9) = ysgp1(i1+1,j)-ysgp1(i2-1,j)
- delta(10) = ysgp1(i1+2,j)-ysgp1(i2-2,j)
- ! check angq match across seam
- j = ny
- delta(11)=angq(i1-2,j)-angq(i2-2,j)
- delta(12)=angq(i1-1,j)-angq(i2-1,j)
- delta(13)=angq(i1, j)-angq(i2, j)
- delta(14)=angq(i1+1,j)-angq(i2+1,j)
- delta(15)=angq(i1+2,j)-angq(i2+2,j)
-
- sumdelta = 0.0
- sumdelta = sum(delta)
-
- if (sumdelta >= puny) then
- print *,'OK'
- print *,'SUCCESS!'
- deallocate(x,y,xsgp1,ysgp1,angq)
- else
- print *,'ftst_find_angq failed'
- stop 1
- endif
-
- end program ftst_find_angq
+
+ use gengrid_kinds, only : dbl_kind, real_kind
+ use angles , only : find_angq, find_ang, find_angchk
+
+ implicit none
+
+ integer , parameter :: ni = 5, nj = 5, nblocks=3
+
+ ! verification values from model history files, mx100
+ ! blocks 1 and 2 center on the left and right tripole
+ ! block 3 is in lower right quad of the arctic at (300:304,306:310)
+ real(dbl_kind), dimension(ni,nj,nblocks) :: cice6_angle
+ real(dbl_kind), dimension(ni,nj,nblocks) :: cice6_anglet
+ real(dbl_kind), dimension(ni,nj,nblocks) :: mom6_sinrot
+
+ ! block grid point values
+ real(dbl_kind), dimension(ni,nj,nblocks) :: lonct, lonbu, latbu
+
+ ! test values
+ real(dbl_kind), dimension(ni,nj,nblocks) :: anglet, angle, angchk
+ real(dbl_kind), dimension(ni,nblocks) :: xangCt
+
+ real(dbl_kind) :: dmax, diff
+ character(len=1) :: ck
+
+ integer :: i,i2,j,k
+ integer :: ipole(nblocks)
+ logical :: debug = .false.
+
+ !CICE6
+ data ((cice6_angle(i,j,1), i=1,ni), j=1,nj) / &
+ 0.55923, 0.29966, 0.00000, -0.29966, -0.55923, &
+ 0.70064, 0.38859, 0.00000, -0.38859, -0.70064, &
+ 0.90943, 0.53371, 0.00000, -0.53371, -0.90943, &
+ 1.20822, 0.75580, 0.00000, -0.75580, -1.20822, &
+ 1.38117, 0.89328, 0.00000, -0.89328, -1.38117/
+ data ((cice6_angle(i,j,2), i=1,ni), j=1,nj) / &
+ 0.55923, 0.29966, 0.00000, -0.29966, -0.55923, &
+ 0.70064, 0.38859, 0.00000, -0.38859, -0.70064, &
+ 0.90943, 0.53371, 0.00000, -0.53371, -0.90943, &
+ 1.20822, 0.75580, 0.00000, -0.75580, -1.20822, &
+ 1.38117, 0.89328, 0.00000, -0.89328, -1.38117/
+ data ((cice6_angle(i,j,3), i=1,ni), j=1,nj) / &
+ -1.14949, -1.16044, -1.17076, -1.18049, -1.18968, &
+ -1.18188, -1.19217, -1.20185, -1.21097, -1.21957, &
+ -1.21425, -1.22383, -1.23284, -1.24131, -1.24930, &
+ -1.24649, -1.25534, -1.26364, -1.27145, -1.27880, &
+ -1.27854, -1.28661, -1.29419, -1.30130, -1.30799 /
+
+ data ((cice6_anglet(i,j,1), i=1,ni), j=1,nj) / &
+ 0.60483, 0.38996, 0.13521, -0.13521, -0.38996, &
+ 0.73204, 0.48692, 0.17198, -0.17198, -0.48692, &
+ 0.90492, 0.63285, 0.23026, -0.23026, -0.63285, &
+ 1.13597, 0.85126, 0.32134, -0.32134, -0.85126, &
+ 1.34026, 1.05945, 0.41175, -0.41175, -1.05945/
+ data ((cice6_anglet(i,j,2), i=1,ni), j=1,nj) / &
+ 0.60483, 0.38996, 0.13521, -0.13521, -0.38996, &
+ 0.73204, 0.48692, 0.17198, -0.17198, -0.48692, &
+ 0.90492, 0.63285, 0.23026, -0.23026, -0.63285, &
+ 1.13597, 0.85126, 0.32134, -0.32134, -0.85126, &
+ 1.34026, 1.05945, 0.41175, -0.41175, -1.05945/
+ data ((cice6_anglet(i,j,3), i=1,ni), j=1,nj) / &
+ -1.12735, -1.13895, -1.14989, -1.16020, -1.16995, &
+ -1.16004, -1.17100, -1.18130, -1.19102, -1.20018, &
+ -1.19278, -1.20303, -1.21267, -1.22174, -1.23029, &
+ -1.22546, -1.23498, -1.24391, -1.25231, -1.26021, &
+ -1.25800, -1.26675, -1.27495, -1.28264, -1.28988 /
+
+ ! MOM6 sin_rot
+ data ((mom6_sinrot(i,j,1), i=1,ni), j=1,nj) / &
+ -0.57231, -0.38355, -0.13608, 0.13608, 0.38355, &
+ -0.67404, -0.47460, -0.17342, 0.17342, 0.47460, &
+ -0.79354, -0.60678, -0.23239, 0.23239, 0.60678, &
+ -0.91250, -0.79512, -0.32692, 0.32692, 0.79512, &
+ -0.98949, -0.97272, -0.43486, 0.43486, 0.97272/
+ data ((mom6_sinrot(i,j,2), i=1,ni), j=1,nj) / &
+ -0.57231, -0.38355, -0.13608, 0.13608, 0.38355, &
+ -0.67404, -0.47460, -0.17342, 0.17342, 0.47460, &
+ -0.79354, -0.60678, -0.23239, 0.23239, 0.60678, &
+ -0.91250, -0.79512, -0.32692, 0.32692, 0.79512, &
+ -0.98949, -0.97272, -0.43486, 0.43486, 0.97272/
+ data ((mom6_sinrot(i,j,3), i=1,ni), j=1,nj) / &
+ 0.90334, 0.90826, 0.91278, 0.91694, 0.92079, &
+ 0.91689, 0.92120, 0.92516, 0.92881, 0.93216, &
+ 0.92946, 0.93320, 0.93662, 0.93976, 0.94264, &
+ 0.94103, 0.94420, 0.94711, 0.94977, 0.95221, &
+ 0.95154, 0.95419, 0.95661, 0.95882, 0.96085 /
+
+ ! lon Ct,Bu and latBu
+ data ((lonct(i,j,1), i=1,ni), j=1,nj) / &
+ -245.26937, -232.98519, -218.04609, -201.95391, -187.01481, &
+ -252.74561, -239.00064, -220.46665, -199.53335, -180.99936, &
+ -262.76769, -248.28139, -224.73626, -195.26374, -171.71861, &
+ -275.86250, -263.23478, -234.03925, -185.96075, -156.76522, &
+ -291.65329, -286.25383, -263.72015, -156.27985, -133.74617 /
+ data ((lonbu(i,j,1), i=1,ni), j=1,nj) / &
+ -242.68153, -227.77943, -210.00000, -192.22057, -177.31847, &
+ -251.01234, -233.49458, -210.00000, -186.50542, -168.98766, &
+ -262.99894, -243.55611, -210.00000, -176.44389, -157.00106, &
+ -279.68531, -263.47519, -210.00000, -156.52481, -140.31469, &
+ -300.00000, -300.00000, -210.00000, -120.00000, -120.00000 /
+ data ((latbu(i,j,1), i=1,ni), j=1,nj) / &
+ 88.99001, 89.10706, 89.14964, 89.10706, 88.99001, &
+ 89.16921, 89.31628, 89.37292, 89.31628, 89.16921, &
+ 89.31750, 89.50699, 89.58912, 89.50699, 89.31750, &
+ 89.41886, 89.66093, 89.79818, 89.66093, 89.41886, &
+ 89.45503, 89.72754, 90.00000, 89.72754, 89.45503 /
+
+ data ((lonct(i,j,2), i=1,ni), j=1,nj) / &
+ -65.26937, -52.98519, -38.04609, -21.95391, -7.01481, &
+ -72.74561, -59.00064, -40.46665, -19.53335, -0.99936, &
+ -82.76769, -68.28139, -44.73626, -15.26374, 8.28139, &
+ -95.86250, -83.23478, -54.03925, -5.96075, 23.23478, &
+ -111.65329, -106.25383, -83.72015, 23.72015, 46.25383 /
+ data ((lonbu(i,j,2), i=1,ni), j=1,nj) / &
+ -62.68153, -47.77943, -30.00000, -12.22057, 2.68153, &
+ -71.01234, -53.49458, -30.00000, -6.50542, 11.01234, &
+ -82.99894, -63.55611, -30.00000, 3.55611, 22.99894, &
+ -99.68531, -83.47519, -30.00000, 23.47519, 39.68531, &
+ -120.00000, -120.00000, -30.00000, 60.00000, 60.00000 /
+ data ((latbu(i,j,2), i=1,ni), j=1,nj) / &
+ 88.99001, 89.10706, 89.14964, 89.10706, 88.99001, &
+ 89.16921, 89.31628, 89.37292, 89.31628, 89.16921, &
+ 89.31750, 89.50699, 89.58912, 89.50699, 89.31750, &
+ 89.41886, 89.66093, 89.79818, 89.66093, 89.41886, &
+ 89.45503, 89.72754, 90.00000, 89.72754, 89.45503 /
+
+ data ((lonct(i,j,3), i=1,ni), j=1,nj) / &
+ 38.08144, 38.86965, 39.62025, 40.33605, 41.01961, &
+ 39.67340, 40.41426, 41.11877, 41.78975, 42.42974, &
+ 41.27339, 41.96446, 42.62075, 43.24505, 43.83985, &
+ 42.87655, 43.51566, 44.12188, 44.69789, 45.24611, &
+ 44.47789, 45.06321, 45.61778, 46.14417, 46.64470 /
+ data ((lonbu(i,j,3), i=1,ni), j=1,nj) / &
+ 39.26336, 40.00933, 40.71974, 41.39725, 42.04431, &
+ 40.83545, 41.53354, 42.19746, 42.82988, 43.43321, &
+ 42.41197, 43.06002, 43.67560, 44.26131, 44.81949, &
+ 43.98813, 44.58428, 45.14993, 45.68757, 46.19946, &
+ 45.55910, 46.10182, 46.61624, 47.10473, 47.56942 /
+ data ((latbu(i,j,3), i=1,ni), j=1,nj) / &
+ 80.97195, 80.70286, 80.43134, 80.15742, 79.88114, &
+ 81.07710, 80.80479, 80.53019, 80.25334, 79.97425, &
+ 81.17218, 80.89690, 80.61947, 80.33992, 80.05826, &
+ 81.25741, 80.97942, 80.69942, 80.41741, 80.13341, &
+ 81.33309, 81.05265, 80.77032, 80.48610, 80.19999 /
+
+ print *,"Starting test of cpld_gridgen routine find_angq"
+
+ angle = 0.0
+ anglet = 0.0
+ angchk = 0.0
+ ipole = 0
+
+ j = nj
+ do i = 1,ni
+ if(latBu(i,j,1) .eq. 90.00)ipole(1) = i
+ if(latBu(i,j,2) .eq. 90.00)ipole(2) = i
+ enddo
+
+ !------------------------------------------
+ ! find anglet and test against mom6 values
+ !------------------------------------------
+
+ call find_ang((/2,5/),(/2,5/),lonBu(:,:,1),latBu(:,:,1),lonCt(:,:,1),anglet(:,:,1))
+ call find_ang((/2,5/),(/2,5/),lonBu(:,:,2),latBu(:,:,2),lonCt(:,:,2),anglet(:,:,2))
+ call find_ang((/2,5/),(/2,5/),lonBu(:,:,3),latBu(:,:,3),lonCt(:,:,3),anglet(:,:,3))
+
+ do k = 1,nblocks
+ write(ck,'(i1.1)')k
+ dmax = 1.0e-30
+ do j = 1,nj
+ do i = 1,ni
+ if (abs(anglet(i,j,k)) .gt. 0.0)then
+ diff = abs(anglet(i,j,k) - asin(mom6_sinrot(i,j,k)))
+ dmax = max(diff,dmax)
+ end if
+ end do
+ end do
+ call passfail(dmax, 1.0e-4, 'MOM6 anglet, block '//trim(ck))
+ end do
+
+ if (debug) then
+ print *,'anglet '
+ print *,'left'
+ do j = 1,nj
+ print '(5f15.6)',(anglet(i,j,1), i = 1,ni)
+ end do
+ print *,'right'
+ do j = 1,nj
+ print '(5f15.6)',(anglet(i,j,2), i = 1,ni)
+ end do
+ print *,'quad'
+ do j = 1,nj
+ print '(5f15.6)',(anglet(i,j,3), i = 1,ni)
+ end do
+ end if
+
+ !--------------------------------------------------------
+ !
+ !---------------------------------------------------------
+
+ xangCt = 0.0
+ do i = 2,ni
+ i2 = ipole(2)+(ipole(1)-i)+1
+ xangCt(i,1) = -anglet(i2,nj,1) ! angle changes sign across seam
+ xangCt(i,2) = -anglet(i2,nj,2)
+ xangCt(i,3) = anglet(i,nj,3)
+ end do
+
+ !-------------------------------------------
+ ! find angle and test against cice6 values
+ !-------------------------------------------
+
+ call find_angq((/2,5/),(/2,5/),xangCt(:,1),anglet(:,:,1),angle(:,:,1))
+ call find_angq((/2,5/),(/2,5/),xangCt(:,2),anglet(:,:,2),angle(:,:,2))
+ call find_angq((/2,5/),(/2,4/),xangCt(:,3),anglet(:,:,3),angle(:,:,3))
+ ! reverse angle for CICE
+ angle = -angle
+
+ do k = 1,nblocks
+ write(ck,'(i1.1)')k
+ dmax = 1.0e-30
+ do j = 1,nj
+ do i = 1,ni
+ if (abs(angle(i,j,k)) .gt. 0.0) then
+ diff = abs(angle(i,j,k) - cice6_angle(i,j,k))
+ dmax = max(diff,dmax)
+ end if
+ end do
+ end do
+ call passfail(dmax, 1.0e-4, 'CICE6 angle, block '//trim(ck))
+ end do
+
+ if (debug) then
+ print *,'angle'
+ print *,'left'
+ do j = 1,nj
+ print '(5f15.6)',(angle(i,j,1), i = 1,ni)
+ end do
+ print *,'right'
+ do j = 1,nj
+ print '(5f15.6)',(angle(i,j,2), i = 1,ni)
+ end do
+ print *,'quad'
+ do j = 1,nj
+ print '(5f15.6)',(angle(i,j,3), i = 1,ni)
+ end do
+ end if
+
+ !-----------------------------------------------------------------------
+ ! find anglet calculated by CICE and test against cice6 and mom6 values
+ !-----------------------------------------------------------------------
+ call find_angchk((/2,4/),(/3,4/),angle(:,:,1),angchk(:,:,1))
+ call find_angchk((/2,4/),(/3,4/),angle(:,:,2),angchk(:,:,2))
+ call find_angchk((/2,4/),(/3,4/),angle(:,:,3),angchk(:,:,3))
+ ! reverse angle for MOM6
+ angchk = -angchk
+
+ do k = 1,nblocks
+ write(ck,'(i1.1)')k
+ dmax = 1.0e-30
+ do j = 1,nj
+ do i = 1,ni
+ if (abs(angchk(i,j,k)) .gt. 0.0) then
+ diff = abs(-angchk(i,j,k) - cice6_anglet(i,j,k))
+ dmax = max(diff,dmax)
+ end if
+ end do
+ end do
+ call passfail(dmax, 1.0e-4, 'angchk vs CICE6 anglet, block '//trim(ck))
+ end do
+
+ do k = 1,nblocks
+ write(ck,'(i1.1)')k
+ dmax = 1.0e-30
+ do j = 1,nj
+ do i = 1,ni
+ if (abs(angchk(i,j,k)) .gt. 0.0) then
+ diff = abs(angchk(i,j,k) - asin(mom6_sinrot(i,j,k)))
+ dmax = max(diff,dmax)
+ end if
+ end do
+ end do
+ call passfail(dmax, 1.2e-2, 'angchk vs MOM6 anglet, block '//trim(ck))
+ end do
+
+ if (debug) then
+ print *,'angchk'
+ print *,'left'
+ do j = 1,nj
+ print '(5f15.6)',(angchk(i,j,1), i = 1,ni)
+ end do
+ print *,'right'
+ do j = 1,nj
+ print '(5f15.6)',(angchk(i,j,2), i = 1,ni)
+ end do
+ print *,'quad'
+ do j = 1,nj
+ print '(5f15.6)',(angchk(i,j,3), i = 1,ni)
+ end do
+ end if
+ end program
+
+ subroutine passfail(dmax,tolerance,msg)
+
+ use gengrid_kinds, only : dbl_kind
+
+ implicit none
+
+ real(dbl_kind), intent(in) :: dmax
+ real(dbl_kind), intent(in) :: tolerance
+ character(len=*), intent(in) :: msg
+
+ if (dmax .le. tolerance) then
+ print '(a,2f12.8)','SUCCESS! '//trim(msg)//' ',dmax,tolerance
+ else
+ print '(a,2f12.8)','FAIL! '//trim(msg)//' ',dmax,tolerance
+ stop 1
+ endif
+ end subroutine passfail
|