Skip to content

Commit

Permalink
feature/regional_grid This commit references NOAA-EMC#4.
Browse files Browse the repository at this point in the history 
Update Jim's regional code to the latest version from
Dusan's regional workflow branch.  This update adds
new records 'dx', 'dy', 'angle_dx', and 'angle_dy' to
the grid files.
  • Loading branch information
@GeorgeGayno-NOAA
GeorgeGayno-NOAA committed Jun 17, 2020
1 parent d9ba1d5 commit c9ef8bd
Show file tree
Hide file tree
Showing 2 changed files with 176 additions and 41 deletions.
180 changes: 141 additions & 39 deletions sorc/fre-nctools.fd/tools/regional_grid.fd/hgrid_ak.f90
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
!=============================================================================
subroutine hgrid_ak(lx,ly,nx,ny,a,k,plat,plon,pazi, &
re,delx,dely, glat,glon,garea, ff)
re,delx,dely, glat,glon,garea,dx,dy,angle_dx,angle_dy, ff)
!=============================================================================
! Use a and k as the parameters of a generalized Schmidt-transformed
! gnomonic mapping centered at (plat,plon) and twisted about this center
Expand All @@ -17,33 +17,49 @@ subroutine hgrid_ak(lx,ly,nx,ny,a,k,plat,plon,pazi, &
! defined in the arrays, glat and glon, and return a rectangular array, garea,
! of dimensions nx-1 by ny-1, that contains the areas of each of the grid cells
! in the SQUARE of the same physical length unit that was employed to define
! the radius of the earth, re (and the central grid steps, delx and dely).
! the radius of the earth, re, and the central grid steps, delx and dely.
! In this version, these grid cell areas are computed by 2D integrating the
! scalar jacobian of the transformation, using a 4th-order centered scheme.
! The estimated grid steps, dx snd dy, are returned at the grid cell edges,
! using the same 4th-order scheme to integrate the 1D projected jacobian.
! The angles, relative to local east and north, are returned respectively
! as angle_dx and angle_dy at grid cell corners, in degrees counterclockwise.
!
! if all goes well, return a .FALSE. failure flag, ff. If, for some
! reason, it is not possible to complete this task, return the failure flag
! as .TRUE.
!=============================================================================
use pkind, only: dp
use pietc, only: u0,u1,dtor
use pmat4, only: sarea
use pietc, only: u0,u1,dtor,rtod
use pmat4, only: cross_product,triple_product
use pmat5, only: ctog
implicit none
integer, intent(in ):: lx,ly,nx,ny
real(dp), intent(in ):: a,k,plat,plon,pazi, &
re,delx,dely
real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: glat,glon
real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea
real(dp),dimension(lx:lx+nx-1,ly:ly+ny ),intent(out):: dx
real(dp),dimension(lx:lx+nx ,ly:ly+ny-1),intent(out):: dy
real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: angle_dx,angle_dy
logical, intent(out):: ff
!-----------------------------------------------------------------------------
real(dp),dimension(lx-1:lx+nx+1,ly-1:ly+ny+1):: gat ! Temporary area array
real(dp),dimension(lx-1:lx+nx+1,ly :ly+ny ):: dxt ! Temporary dx array
real(dp),dimension(lx :lx+nx ,ly-1:ly+ny+1):: dyt ! Temporary dy array
real(dp),dimension(3,3):: prot,azirot
real(dp),dimension(3,2):: xcd
real(dp),dimension(3) :: xc
real(dp),dimension(3,2):: xcd,eano
real(dp),dimension(2,2):: xcd2
real(dp),dimension(3) :: xc,east,north
real(dp),dimension(2) :: xm
real(dp) :: clat,slat,clon,slon,cazi,sazi,&
rlat,drlata,drlatb,drlatc, &
rlon,drlona,drlonb,drlonc, rre
integer :: ix,iy,mx,my
rlon,drlona,drlonb,drlonc, delxy,delxore,delyore
integer :: ix,iy,mx,my,lxm,lym,mxp,myp
!=============================================================================
delxore=delx/re
delyore=dely/re
delxy=delx*dely
clat=cos(plat); slat=sin(plat)
clon=cos(plon); slon=sin(plon)
cazi=cos(pazi); sazi=sin(pazi)
Expand All @@ -56,49 +72,135 @@ subroutine hgrid_ak(lx,ly,nx,ny,a,k,plat,plon,pazi, &
prot(:,2)=(/-slat*clon, -slat*slon, clat/)
prot(:,3)=(/ clat*clon, clat*slon, slat/)
prot=matmul(prot,azirot)

mx=lx+nx ! Index of the 'right' edge of the rectangular grid
my=ly+ny ! Index of the 'top' edge of the rectangular grid
!This code assumes symmetry about the grid center
lxm=lx-1; mxp=mx+1 ! Indices of extra left and right edges
lym=ly-1; myp=my+1 ! Indices of extra bottom and top edges

!-- main body of horizontal grid:
do iy=ly,my
!xm(2)=iy*dely/re
xm(2)=-iy*dely/re
xm(2)=iy*delyore
do ix=lx,mx
!xm(1)=ix*delx/re
xm(1)=-ix*delx/re
call xmtoxc_ak(a,k,xm,xc,xcd,ff)
if(ff)return
xm(1)=ix*delxore
call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
xcd=matmul(prot,xcd)
xc =matmul(prot,xc )
call ctog(xc,glat(ix,iy),glon(ix,iy))
east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
north=cross_product(xc,east)
eano(:,1)=east; eano(:,2)=north
xcd2=matmul(transpose(eano),xcd)
angle_dx(ix,iy)=atan2( xcd2(2,1),xcd2(1,1))*rtod
angle_dy(ix,iy)=atan2(-xcd2(1,2),xcd2(2,2))*rtod
dxt(ix,iy)=sqrt(dot_product(xcd2(:,1),xcd2(:,1)))*delx
dyt(ix,iy)=sqrt(dot_product(xcd2(:,2),xcd2(:,2)))*dely
gat(ix,iy)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
enddo
enddo

!-- extra left edge, gat, dxt only:
xm(1)=lxm*delxore
do iy=ly,my
xm(2)=iy*delyore
call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
xcd=matmul(prot,xcd)
xc =matmul(prot,xc )
east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
north=cross_product(xc,east)
eano(:,1)=east; eano(:,2)=north
xcd2=matmul(transpose(eano),xcd)
dxt(lxm,iy)=sqrt(dot_product(xcd2(:,1),xcd2(:,1)))*delx
gat(lxm,iy)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
enddo

!-- extra right edge, gat, dxt only:
xm(1)=mxp*delxore
do iy=ly,my
xm(2)=iy*delyore
call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
xcd=matmul(prot,xcd)
xc =matmul(prot,xc )
east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
north=cross_product(xc,east)
eano(:,1)=east; eano(:,2)=north
xcd2=matmul(transpose(eano),xcd)
dxt(mxp,iy)=sqrt(dot_product(xcd2(:,1),xcd2(:,1)))*delx
gat(mxp,iy)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
enddo

!-- extra bottom edge, gat, dyt only:
xm(2)=lym*delyore
do ix=lx,mx
xm(1)=ix*delxore
call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
xcd=matmul(prot,xcd)
xc =matmul(prot,xc )
east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
north=cross_product(xc,east)
eano(:,1)=east; eano(:,2)=north
xcd2=matmul(transpose(eano),xcd)
dyt(ix,lym)=sqrt(dot_product(xcd2(:,2),xcd2(:,2)))*dely
gat(ix,lym)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
enddo

!-- extra top edge, gat, dyt only:
xm(2)=myp*delyore
do ix=lx,mx
xm(1)=ix*delxore
call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
xcd=matmul(prot,xcd)
xc =matmul(prot,xc )
east=(/-xc(2),xc(1),u0/); east=east/sqrt(dot_product(east,east))
north=cross_product(xc,east)
eano(:,1)=east; eano(:,2)=north
xcd2=matmul(transpose(eano),xcd)
dyt(ix,myp)=sqrt(dot_product(xcd2(:,2),xcd2(:,2)))*dely
gat(ix,myp)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy
enddo

! Extra four corners, gat only:
xm(2)=lym*delyore
!-- extra bottom left corner:
xm(1)=lxm*delxore
call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
xcd=matmul(prot,xcd)
xc =matmul(prot,xc )
gat(lxm,lym)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy

!-- extra bottom right corner:
xm(1)=mxp*delxore
call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
xcd=matmul(prot,xcd)
xc =matmul(prot,xc )
gat(mxp,lym)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy

xm(2)=myp*delyore
!-- extra top left corner:
xm(1)=lxm*delxore
call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
xcd=matmul(prot,xcd)
xc =matmul(prot,xc )
gat(lxm,myp)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy

!-- extra top right corner:
xm(1)=mxp*delxore
call xmtoxc_ak(a,k,xm,xc,xcd,ff); if(ff)return
xcd=matmul(prot,xcd)
xc =matmul(prot,xc )
gat(mxp,myp)=triple_product(xc,xcd(:,1),xcd(:,2))*delxy

!-- 4th-order averaging over each central interval using 4-pt. stencils:
dx =(13*(dxt(lx :mx-1,:)+dxt(lx+1:mx ,:)) &
-(dxt(lxm:mx-2,:)+dxt(lx+2:mxp,:)))/24
dy =(13*(dyt(:,ly :my-1)+dyt(:,ly+1:my )) &
-(dyt(:,lym:my-2)+dyt(:,ly+2:myp)))/24
gat(lx:mx-1,:)=(13*(gat(lx :mx-1,:)+gat(lx+1:mx ,:)) &
-(gat(lxm:mx-2,:)+gat(lx+2:mxp,:)))/24
garea =(13*(gat(lx:mx-1,ly :my-1)+gat(lx:mx-1,ly+1:my )) &
-(gat(lx:mx-1,lym:my-2)+gat(lx:mx-1,ly+2:myp)))/24
! Convert degrees to radians in the glat and glon arrays:
glat=glat*dtor
glon=glon*dtor

! Compute the areas of the quadrilateral grid cells:
do iy=ly,my-1
do ix=lx,mx-1
rlat =glat(ix ,iy )
drlata=glat(ix+1,iy )-rlat
drlatb=glat(ix+1,iy+1)-rlat
drlatc=glat(ix ,iy+1)-rlat
rlon =glon(ix ,iy )
drlona=glon(ix+1,iy )-rlon
drlonb=glon(ix+1,iy+1)-rlon
drlonc=glon(ix ,iy+1)-rlon
! If 'I' is the grid point (ix,iy), 'A' is (ix+1,iy); 'B' is (ix+1,iy+1)
! and 'C' is (ix,iy+1), then the area of the grid cell IABC is the sum of
! the areas of the traingles, IAB and IBC (the latter being the negative
! of the signed area of triangle, ICB):
garea(ix,iy)=sarea(rlat, drlata,drlona, drlatb,drlonb) &
-sarea(rlat, drlatc,drlonc, drlatb,drlonb)
enddo
enddo
! Convert the areas to area units consistent with the length units used for
! the radius, re, of the sphere:
rre=re*re
garea=garea*rre

end subroutine hgrid_ak
37 changes: 35 additions & 2 deletions sorc/fre-nctools.fd/tools/regional_grid.fd/regional_grid.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,13 +23,15 @@ program regional_grid

real(dp),dimension(:,:),allocatable:: glat,glon
real(dp),dimension(:,:),allocatable:: garea
real(dp),dimension(:,:),allocatable:: dx,dy,angle_dx,angle_dy

character(len=256) :: nml_fname

! netcdf
integer :: ncid
integer :: string_dimid, nxp_dimid, nyp_dimid, nx_dimid, ny_dimid
integer :: tile_varid, x_varid, y_varid, area_varid
integer :: dx_varid, dy_varid, angle_dx_varid, angle_dy_varid
integer, dimension(2) :: dimids

!=============================================================================
Expand All @@ -56,9 +58,14 @@ program regional_grid
allocate(glon(0:nx,0:ny))
allocate(garea(0:nxm,0:nym))

allocate(dx(0:nxm,0:ny))
allocate(dy(0:nx,0:nym))
allocate(angle_dx(0:nx,0:ny))
allocate(angle_dy(0:nx,0:ny))

call hgrid_ak(lx,ly,nx,ny,a,k,plat*dtor,plon*dtor,pazi*dtor, &
re,redelx,redely, glat,glon,garea, ff)
if(ff)stop 'Failure flag raised in hgrid routine'
re,redelx,redely, glat,glon,garea,dx,dy,angle_dx,angle_dy, ff)
if(ff)stop 'Failure flag raised in hgrid_ak routine'

glon = glon*rtod
glat = glat*rtod
Expand Down Expand Up @@ -89,6 +96,28 @@ program regional_grid
call check( nf90_put_att(ncid, area_varid, "units", "m2") )
call check( nf90_put_att(ncid, area_varid, "hstagger", "H") )

dimids = (/ nx_dimid, nyp_dimid /)
call check( nf90_def_var(ncid, "dx", NF90_DOUBLE, dimids, dx_varid) )
call check( nf90_put_att(ncid, dx_varid, "standard_name", "dx") )
call check( nf90_put_att(ncid, dx_varid, "units", "m") )
call check( nf90_put_att(ncid, dx_varid, "hstagger", "H") )

dimids = (/ nxp_dimid, ny_dimid /)
call check( nf90_def_var(ncid, "dy", NF90_DOUBLE, dimids, dy_varid) )
call check( nf90_put_att(ncid, dy_varid, "standard_name", "dy") )
call check( nf90_put_att(ncid, dy_varid, "units", "m") )
call check( nf90_put_att(ncid, dy_varid, "hstagger", "H") )

dimids = (/ nxp_dimid, nyp_dimid /)
call check( nf90_def_var(ncid, "angle_dx", NF90_DOUBLE, dimids, angle_dx_varid) )
call check( nf90_put_att(ncid, angle_dx_varid, "standard_name", "angle_dx") )
call check( nf90_put_att(ncid, angle_dx_varid, "units", "deg") )
call check( nf90_put_att(ncid, angle_dx_varid, "hstagger", "C") )
call check( nf90_def_var(ncid, "angle_dy", NF90_DOUBLE, dimids, angle_dy_varid) )
call check( nf90_put_att(ncid, angle_dy_varid, "standard_name", "angle_dy") )
call check( nf90_put_att(ncid, angle_dy_varid, "units", "deg") )
call check( nf90_put_att(ncid, angle_dy_varid, "hstagger", "C") )

call check( nf90_put_att(ncid, NF90_GLOBAL, "history", "gnomonic_ed") )
call check( nf90_put_att(ncid, NF90_GLOBAL, "source", "FV3GFS") )
call check( nf90_put_att(ncid, NF90_GLOBAL, "grid", "akappa") )
Expand All @@ -108,6 +137,10 @@ program regional_grid
call check( nf90_put_var(ncid, x_varid, glon) )
call check( nf90_put_var(ncid, y_varid, glat) )
call check( nf90_put_var(ncid, area_varid, garea) )
call check( nf90_put_var(ncid, dx_varid, dx) )
call check( nf90_put_var(ncid, dy_varid, dy) )
call check( nf90_put_var(ncid, angle_dx_varid, angle_dx) )
call check( nf90_put_var(ncid, angle_dy_varid, angle_dy) )

call check( nf90_close(ncid) )

Expand Down

0 comments on commit c9ef8bd

Please sign in to comment.