From 3151d66d77603befd7b994473eca7cbb593ac2f0 Mon Sep 17 00:00:00 2001 From: "lin.gan" Date: Mon, 1 Mar 2021 19:31:34 +0000 Subject: [PATCH 1/5] Part of Doxygen updates for sorc/grid_tools.fd #366 doxygen update for pesg.f90 --- .../regional_esg_grid.fd/pesg.f90 | 1008 +++++++++-------- 1 file changed, 547 insertions(+), 461 deletions(-) diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 index 0e9b54bef..c92a5bd4a 100644 --- a/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 +++ b/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 @@ -62,59 +62,57 @@ module pesg contains -!============================================================================= +!> Inverse of xstoxc. I.e., cartesians to stereographic +!! @author R. J. Purser +!! @param xc xstoxc +!! @param xs Inverse of xstoxc subroutine xctoxs(xc,xs)! [xctoxs] -!============================================================================= -! Inverse of xstoxc. I.e., cartesians to stereographic -!============================================================================= implicit none real(dp),dimension(3),intent(in ):: xc real(dp),dimension(2),intent(out):: xs -!----------------------------------------------------------------------------- real(dp):: zp -!============================================================================= zp=u1+xc(3); xs=xc(1:2)/zp end subroutine xctoxs -!============================================================================= +!> Standard transformation from polar stereographic map coordinates, xs, to +!! cartesian, xc, assuming the projection axis is polar. +!! xcd=d(xc)/d(xs) is the jacobian matrix, encoding distortion and metric. +!! @author R. J. Purser +!! @param xs polar stereographic map coordinates +!! @param xc cartesian +!! @param xcd value of jacobian matrix, encoding distortion and metric subroutine xstoxc(xs,xc,xcd)! [xstoxc] -!============================================================================= -! Standard transformation from polar stereographic map coordinates, xs, to -! cartesian, xc, assuming the projection axis is polar. -! xcd=d(xc)/d(xs) is the jacobian matrix, encoding distortion and metric. -!============================================================================= use pmat4, only: outer_product implicit none real(dp),dimension(2), intent(in ):: xs real(dp),dimension(3), intent(out):: xc real(dp),dimension(3,2),intent(out):: xcd -!----------------------------------------------------------------------------- real(dp):: zp -!============================================================================= zp=u2/(u1+dot_product(xs,xs)); xc(1:2)=xs*zp; xc(3)=zp xcd=-outer_product(xc,xs)*zp; xcd(1,1)=xcd(1,1)+zp; xcd(2,2)=xcd(2,2)+zp xc(3)=xc(3)-u1 end subroutine xstoxc -!============================================================================= + +!> Standard transformation from polar stereographic map coordinates, xs, to +!! cartesian, xc, assuming the projection axis is polar. +!! xcd=d(xc)/d(xs) is the jacobian matrix, encoding distortion and metric. +!! xcdd is the further derivative, wrt xs, of xcd. +!! @author R. J. Purser +!! @param xs polar stereographic map coordinates +!! @param xc cartesian +!! @param xcd jacobian matrix, encoding distortion and metric +!! @param xcdd further derivative, wrt xs, of xcd subroutine xstoxc1(xs,xc,xcd,xcdd)! [xstoxc] -!============================================================================= -! Standard transformation from polar stereographic map coordinates, xs, to -! cartesian, xc, assuming the projection axis is polar. -! xcd=d(xc)/d(xs) is the jacobian matrix, encoding distortion and metric. -! xcdd is the further derivative, wrt xs, of xcd. -!============================================================================= use pmat4, only: outer_product implicit none real(dp),dimension(2), intent(in ):: xs real(dp),dimension(3), intent(out):: xc real(dp),dimension(3,2), intent(out):: xcd real(dp),dimension(3,2,2),intent(out):: xcdd -!----------------------------------------------------------------------------- real(dp),dimension(3,2):: zpxcdxs real(dp),dimension(3) :: zpxc real(dp) :: zp integer(spi) :: i -!============================================================================= zp=u2/(u1+dot_product(xs,xs)); xc(1:2)=xs*zp; xc(3)=zp xcd=-outer_product(xc,xs)*zp zpxc=zp*xc; xc(3)=xc(3)-u1; xcdd=u0 @@ -129,29 +127,25 @@ subroutine xstoxc1(xs,xc,xcd,xcdd)! [xstoxc] do i=1,2; xcd(i,i)=xcd(i,i)+zp; enddo end subroutine xstoxc1 -!============================================================================= +!> Inverse of xttoxs +!! @author R. J. Purser subroutine xstoxt(k,xs,xt,ff)! [xstoxt] -!============================================================================= -! Inverse of xttoxs. -!============================================================================= implicit none real(dp), intent(in ):: k real(dp),dimension(2),intent(in ):: xs real(dp),dimension(2),intent(out):: xt logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp):: s,sc -!============================================================================= s=k*(xs(1)*xs(1)+xs(2)*xs(2)); sc=u1-s ff=abs(s)>=u1; if(ff)return xt=u2*xs/sc end subroutine xstoxt -!============================================================================= +!> Scaled gnomonic plane xt to standard stereographic plane xs +!! @author R. J. Purser +!! @param xt Scaled gnomonic plane +!! @param xs standard stereographic plane subroutine xttoxs(k,xt,xs,xsd,ff)! [xttoxs -!============================================================================= -! Scaled gnomonic plane xt to standard stereographic plane xs -!============================================================================= use pmat4, only: outer_product implicit none real(dp), intent(in ):: k @@ -159,11 +153,9 @@ subroutine xttoxs(k,xt,xs,xsd,ff)! [xttoxs real(dp),dimension(2), intent(out):: xs real(dp),dimension(2,2),intent(out):: xsd logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp),dimension(2):: rspd real(dp) :: s,sp,rsp,rsppi,rsppis integer(spi) :: i -!============================================================================= s=k*dot_product(xt,xt); sp=u1+s ff=(sp<=u0); if(ff)return rsp=sqrt(sp) @@ -174,11 +166,13 @@ subroutine xttoxs(k,xt,xs,xsd,ff)! [xttoxs xsd=-outer_product(xt,rspd)*rsppis do i=1,2; xsd(i,i)=xsd(i,i)+rsppi; enddo end subroutine xttoxs -!============================================================================= + +!> Like xttoxs, but also, return the derivatives, wrt K, of xs and xsd +!! @author R. J. Purser +!! @param xt Scaled gnomonic plane +!! @param xs standard stereographic plane +!! @param K derivatives, wrt K, of xs and xsd subroutine xttoxs1(k,xt,xs,xsd,xsdd,xs1,xsd1,ff)! [xttoxs] -!============================================================================= -! Like xttoxs, but also, return the derivatives, wrt K, of xs and xsd -!============================================================================= use pmat4, only: outer_product implicit none real(dp), intent(in ):: k @@ -187,12 +181,10 @@ subroutine xttoxs1(k,xt,xs,xsd,xsdd,xs1,xsd1,ff)! [xttoxs] real(dp),dimension(2,2), intent(out):: xsd,xsd1 real(dp),dimension(2,2,2),intent(out):: xsdd logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp),dimension(2,2):: rspdd real(dp),dimension(2) :: rspd,rspd1,rsppid real(dp) :: s,sp,rsp,rsppi,rsppis,s1,rsp1 integer(spi) :: i -!============================================================================= s1=dot_product(xt,xt); s=k*s1; sp=u1+s ff=(sp<=u0); if(ff)return rsp=sqrt(sp); rsp1=o2*s1/rsp @@ -215,53 +207,49 @@ subroutine xttoxs1(k,xt,xs,xsd,xsdd,xs1,xsd1,ff)! [xttoxs] do i=1,2; xsdd(i,:,:)=xsdd(i,:,:)-xt(i)*rspdd*rsppi*k/sp; enddo end subroutine xttoxs1 -!============================================================================= +!> Inverse of xmtoxt +!! @author R. J. Purser subroutine xttoxm(a,xt,xm,ff)! [xttoxm] -!============================================================================= -! Inverse of xmtoxt -!============================================================================= implicit none real(dp), intent(in ):: a real(dp),dimension(2),intent(in ):: xt real(dp),dimension(2),intent(out):: xm logical ,intent(out):: ff -!----------------------------------------------------------------------------- integer(spi):: i -!============================================================================= do i=1,2; call zttozm(a,xt(i),xm(i),ff); if(ff)return; enddo end subroutine xttoxm -!============================================================================= +!> Like zmtozt, but for 2-vector xm and xt, and 2*2 diagonal Jacobian xtd +!! @author R. J. Purser +!! @param xm vector value +!! @param xt vector value +!! @param xtd 2*2 diagonal Jacobian subroutine xmtoxt(a,xm,xt,xtd,ff)! [xmtoxt] -!============================================================================= -! Like zmtozt, but for 2-vector xm and xt, and 2*2 diagonal Jacobian xtd -!============================================================================= implicit none real(dp), intent(in ):: a real(dp),dimension(2), intent(in ):: xm real(dp),dimension(2), intent(out):: xt real(dp),dimension(2,2),intent(out):: xtd logical, intent(out):: ff -!----------------------------------------------------------------------------- integer(spi):: i -!============================================================================= xtd=u0; do i=1,2; call zmtozt(a,xm(i),xt(i),xtd(i,i),ff); if(ff)return; enddo end subroutine xmtoxt -!============================================================================= + +!> Like zmtozt1, but for 2-vector xm and xt, and 2*2 diagonal Jacobian xtd +!! Also, the derivatives, wrt a, of these quantities. +!! @author R. J. Purser +!! @param xm vector value +!! @param xt vector value +!! @param xtd 2*2 diagonal Jacobian +!! @param a the derivatives of these quantities subroutine xmtoxt1(a,xm,xt,xtd,xt1,xtd1,ff)! [xmtoxt] -!============================================================================= -! Like zmtozt1, but for 2-vector xm and xt, and 2*2 diagonal Jacobian xtd -! Also, the derivatives, wrt a, of these quantities. -!============================================================================= implicit none real(dp), intent(in ):: a real(dp),dimension(2), intent(in ):: xm real(dp),dimension(2), intent(out):: xt,xt1 real(dp),dimension(2,2), intent(out):: xtd,xtd1 logical, intent(out):: ff -!----------------------------------------------------------------------------- integer(spi):: i -!============================================================================= xtd=u0 xtd1=u0 do i=1,2 @@ -270,18 +258,14 @@ subroutine xmtoxt1(a,xm,xt,xtd,xt1,xtd1,ff)! [xmtoxt] enddo end subroutine xmtoxt1 -!============================================================================= +!> Inverse of zmtozt +!! @author R. J. Purser subroutine zttozm(a,zt,zm,ff)! [zttozm] -!============================================================================= -! Inverse of zmtozt -!============================================================================= implicit none real(dp),intent(in ):: a,zt real(dp),intent(out):: zm logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp):: ra,razt -!============================================================================= ff=F if (a>u0)then; ra=sqrt( a); razt=ra*zt; zm=atan (razt)/ra elseif(a=u1; if(ff)return @@ -290,19 +274,17 @@ subroutine zttozm(a,zt,zm,ff)! [zttozm] endif end subroutine zttozm -!============================================================================= +!> Evaluate the function, zt = tan(sqrt(A)*z)/sqrt(A), and its derivative, ztd, +!! for positive and negative A and for the limiting case, A --> 0 +!! @author R. J. Purser +!! @param zt function to be evaluated +!! @param ztd derivative of the source function subroutine zmtozt(a,zm,zt,ztd,ff)! [zmtozt] -!============================================================================= -! Evaluate the function, zt = tan(sqrt(A)*z)/sqrt(A), and its derivative, ztd, -! for positive and negative A and for the limiting case, A --> 0 -!============================================================================= implicit none real(dp),intent(in ):: a,zm real(dp),intent(out):: zt,ztd logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp):: ra -!============================================================================= ff=f if (a>u0)then; ra=sqrt( a); zt=tan (ra*zm)/ra; ff=abs(ra*zm)>=pih elseif(a Like zmtozt, but +!! also, get the derivative with respect to a, zt1 of zt, and ztd1 of ztd. +!! @author R. J. Purser subroutine zmtozt1(a,zm,zt,ztd,zt1,ztd1,ff)! [zmtozt] -!============================================================================= -! Like zmtozt, but -! also, get the derivative with respect to a, zt1 of zt, and ztd1 of ztd. -!============================================================================= use pietc, only: o3 use pfun, only: sinoxm,sinox,sinhoxm,sinhox implicit none real(dp),intent(in ):: a,zm real(dp),intent(out):: zt,ztd,zt1,ztd1 logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp):: ra,rad,razm -!============================================================================= ff=f if (a>u0)then;ra=sqrt( a);razm=ra*zm; zt=tan(razm)/ra; ff=abs(razm)>=pih rad=o2/ra @@ -338,29 +317,29 @@ subroutine zmtozt1(a,zm,zt,ztd,zt1,ztd1,ff)! [zmtozt] ztd1=zt*zt +u2*a*zt*zt1 end subroutine zmtozt1 -!============================================================================= +!> Assuming the vector AK parameterization of the Extended Schmidt-transformed +!! Gnomonic (ESG) mapping with parameter vector, and given a map-space +!! 2-vector, xm, find the corresponding cartesian unit 3-vector and its +!! derivative wrt xm, the Jacobian matrix, xcd. +!! If for any reason the mapping cannot be done, return a raised failure flag,z +!! FF. +!! @author R. J. Purser +!! @param ak parameterization of the Extended ESG +!! @param xm vector +!! @param ff raised failure error code subroutine xmtoxc_vak(ak,xm,xc,xcd,ff)! [xmtoxc_ak] -!============================================================================= -! Assuming the vector AK parameterization of the Extended Schmidt-transformed -! Gnomonic (ESG) mapping with parameter vector, and given a map-space -! 2-vector, xm, find the corresponding cartesian unit 3-vector and its -! derivative wrt xm, the Jacobian matrix, xcd. -! If for any reason the mapping cannot be done, return a raised failure flag,z -! FF. -!============================================================================= implicit none real(dp),dimension(2), intent(in ):: ak,xm real(dp),dimension(3), intent(out):: xc real(dp),dimension(3,2),intent(out):: xcd logical, intent(out):: ff -!============================================================================= call xmtoxc_ak(ak(1),ak(2),xm,xc,xcd,ff) end subroutine xmtoxc_vak -!============================================================================= + +!> Like xmtoxc_vak, _ak, but also return derivatives wrt ak. +!! @author R. J. Purser +!! @param ak derivatives wrt ak subroutine xmtoxc_vak1(ak,xm,xc,xcd,xc1,xcd1,ff)! [xmtoxc_ak] -!============================================================================= -! Like xmtoxc_vak, _ak, but also return derivatives wrt ak. -!============================================================================= implicit none real(dp),dimension(2), intent(in ):: ak,xm real(dp),dimension(3), intent(out):: xc @@ -368,13 +347,11 @@ subroutine xmtoxc_vak1(ak,xm,xc,xcd,xc1,xcd1,ff)! [xmtoxc_ak] real(dp),dimension(3,2), intent(out):: xc1 real(dp),dimension(3,2,2),intent(out):: xcd1 logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp),dimension(3,2,2):: xcdd real(dp),dimension(2,2,2):: xsd1,xsdd real(dp),dimension(2,2) :: xtd,xsd,xs1,xtd1,xsdk real(dp),dimension(2) :: xt,xt1,xs,xsk integer(spi) :: i -!============================================================================= call xmtoxt1(ak(1),xm,xt,xtd,xt1,xtd1,ff); if(ff)return call xttoxs1(ak(2),xt,xs,xsd,xsdd,xsk,xsdk,ff); if(ff)return xs1(:,2)=xsk; xs1(:,1)=matmul(xsd,xt1) @@ -387,25 +364,28 @@ subroutine xmtoxc_vak1(ak,xm,xc,xcd,xc1,xcd1,ff)! [xmtoxc_ak] do i=1,2; xcd1(:,:,i)=xcd1(:,:,i)+matmul(xcd,xsd1(:,:,i)); enddo xcd=matmul(xcd,xsd) end subroutine xmtoxc_vak1 -!============================================================================= + +!> Assuming the A-K parameterization of the Extended Schmidt-transformed +!! Gnomonic (ESG) mapping, and given a map-space 2-vector, xm, find the +!! corresponding cartesian unit 3-vector and its derivative wrt xm, jacobian +!! matrix, xcd. If for any reason the mapping cannot be done, return a +!! raised failure flag, FF. +!! @author R. J. Purser +!! @param a ESG mapping parameterization +!! @param k ESG mapping parameterization +!! @param xm map-space vector +!! @param xm derivative +!! @param xcd jacobian matrix +!! @param ff error flag subroutine xmtoxc_ak(a,k,xm,xc,xcd,ff)! [xmtoxc_ak] -!============================================================================= -! Assuming the A-K parameterization of the Extended Schmidt-transformed -! Gnomonic (ESG) mapping, and given a map-space 2-vector, xm, find the -! corresponding cartesian unit 3-vector and its derivative wrt xm, jacobian -! matrix, xcd. If for any reason the mapping cannot be done, return a -! raised failure flag, FF. -!============================================================================= implicit none real(dp), intent(in ):: a,k real(dp),dimension(2), intent(in ):: xm real(dp),dimension(3), intent(out):: xc real(dp),dimension(3,2),intent(out):: xcd logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp),dimension(2,2):: xtd,xsd real(dp),dimension(2) :: xt,xs -!============================================================================= call xmtoxt(a,xm,xt,xtd,ff); if(ff)return call xttoxs(k,xt,xs,xsd,ff); if(ff)return xsd=matmul(xsd,xtd) @@ -413,63 +393,63 @@ subroutine xmtoxc_ak(a,k,xm,xc,xcd,ff)! [xmtoxc_ak] xcd=matmul(xcd,xsd) end subroutine xmtoxc_ak -!============================================================================= +!> Inverse mapping of xmtoxc_ak. That is, go from given cartesian unit +!! 3-vector, xc, to map coordinate 2-vector xm (or return a raised failure +!! flag, FF, if the attempt fails). +!! @author R. J. Purser +!! @param a Inverse mapping of xmtoxc +!! @param k Inverse mapping of xmtoxc +!! @param xc cartesian unit 3-vector +!! @param xm 2-vector map coordinate +!! @param ff error flag subroutine xctoxm_ak(a,k,xc,xm,ff)! [xctoxm_ak] -!============================================================================= -! Inverse mapping of xmtoxc_ak. That is, go from given cartesian unit -! 3-vector, xc, to map coordinate 2-vector xm (or return a raised failure -! flag, FF, if the attempt fails). -!============================================================================= implicit none real(dp), intent(in ):: a,k real(dp),dimension(3),intent(in ):: xc real(dp),dimension(2),intent(out):: xm logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp),dimension(2):: xs,xt -!============================================================================= ff=F call xctoxs(xc,xs) call xstoxt(k,xs,xt,ff); if(ff)return call xttoxm(a,xt,xm,ff) end subroutine xctoxm_ak -!============================================================================= +!> For angles (degrees) of the arcs spanning the halfwidths between the +!! region's center and its x and y edges, get the two cartesian vectors +!! that represent the locations of these edge midpoints in the positive x and y +!! directions. +!! @author R. J. Purser +!! @param edgex region's x edges +!! @param edgey region's y edges subroutine get_edges(arcx,arcy,edgex,edgey)! [get_edges] -!============================================================================= -! For angles (degrees) of the arcs spanning the halfwidths between the -! region's center and its x and y edges, get the two cartesian vectors -! that represent the locations of these edge midpoints in the positive x and y -! directions. -!============================================================================= implicit none real(dp), intent(in ):: arcx,arcy real(dp),dimension(3),intent(out):: edgex,edgey -!----------------------------------------------------------------------------- real(dp):: cx,sx,cy,sy,rarcx,rarcy -!============================================================================= rarcx=arcx*dtor; rarcy=arcy*dtor cx=cos(rarcx); sx=sin(rarcx) cy=cos(rarcy); sy=sin(rarcy) edgex=(/sx,u0,cx/); edgey=(/u0,sy,cy/) end subroutine get_edges -!============================================================================= +!> From a jacobian matrix, j0, get a sufficient set of v.. diagnostics such +!! that, from averages of these v, we can later compute the collective variance +!! of Q(lam) that they imply for any choice of the "lambda" parameter, lam. +!! Note that v1 and v4 are quadratic diagnostics of EL, while v2 and v3 are +!! linear. +!! @author R. J. Purser +!! @param j0 jacobian matrix +!! @param v1 quadratic diagnostics of EL +!! @param v4 quadratic diagnostics of EL +!! @param v2 linear diagnostics of EL +!! @param v3 linear diagnostics of EL subroutine get_qx(j0, v1,v2,v3,v4)! [get_qx] -!============================================================================= -! From a jacobian matrix, j0, get a sufficient set of v.. diagnostics such -! that, from averages of these v, we can later compute the collective variance -! of Q(lam) that they imply for any choice of the "lambda" parameter, lam. -! Note that v1 and v4 are quadratic diagnostics of EL, while v2 and v3 are -! linear. -!============================================================================= use psym2, only: logsym2 implicit none real(dp),dimension(3,2),intent(in ):: j0 real(dp), intent(out):: v1,v2,v3,v4 -!----------------------------------------------------------------------------- real(dp),dimension(2,2):: el,g -!============================================================================= g=matmul(transpose(j0),j0) call logsym2(g,el); el=el*o2 v1=el(1,1)**2+u2*el(1,2)**2+el(2,2)**2 @@ -477,25 +457,24 @@ subroutine get_qx(j0, v1,v2,v3,v4)! [get_qx] v3=el(2,2) v4=(el(1,1)+el(2,2))**2 end subroutine get_qx -!============================================================================= + +!> From a jacobian matrix, j0, and its derivative, j0d, get a sufficient set +!! of v.. diagnostics such that, from average of these diagnostics, we can +!! later compute the collective variance of Q and its derivative. +!! @author R. J. Purser +!! @param j0 jacobian matrix +!! @param j0d derivative of j0 subroutine get_qxd(j0,j0d, v1,v2,v3,v4,v1d,v2d,v3d,v4d)! [get_qx] -!============================================================================= -! From a jacobian matrix, j0, and its derivative, j0d, get a sufficient set -! of v.. diagnostics such that, from average of these diagnostics, we can -! later compute the collective variance of Q and its derivative. -!============================================================================= use psym2, only: logsym2 implicit none real(dp),dimension(3,2), intent(in ):: j0 real(dp),dimension(3,2,2),intent(in ):: j0d real(dp), intent(out):: v1,v2,v3,v4 real(dp),dimension(2), intent(out):: v1d,v2d,v3d,v4d -!----------------------------------------------------------------------------- real(dp),dimension(2,2,2,2):: deldg real(dp),dimension(2,2,2) :: eld,gd real(dp),dimension(2,2) :: el,g integer(spi) :: i,j,k -!============================================================================= g=matmul(transpose(j0),j0) do i=1,2 gd(:,:,i)=matmul(transpose(j0d(:,:,i)),j0)+matmul(transpose(j0),j0d(:,:,i)) @@ -515,25 +494,31 @@ subroutine get_qxd(j0,j0d, v1,v2,v3,v4,v1d,v2d,v3d,v4d)! [get_qx] v4d=u2*(el(1,1)+el(2,2))*(eld(1,1,:)+eld(2,2,:)) end subroutine get_qxd -!============================================================================= +!> For a parameter vector, ak and a map-space domain-parameter vector, ma, +!! return the lambda-parameterized quality diagnostic, Q, and the geographic +!! domain-parameter vector ga. Lambda is given by lam <1. Also, return +!! the derivatives, qdak and qdma, of Q wrt ak and ma, and the derivatives +!! gadak and gadma, of ga wrt ak and ma. +!! The domain averages of Q are accurately computed by bi-Gauss-Legendre +!! quadrature over the positive quadrant of the domain (exploiting the +!! symmetry) of the four constituent terms, v1, v2, v3, v4, from which +!! the mean Q is computed using a quadratic formula of these constituents. +!! The number of Gauss points in eaxh half-interval is ngh, and the +!! nodes themselves are, in proportion to the half-interval, at xg. +!! the normalized gauss weights are wg. +!! If a failure occurs, colmputations cease immediately and a failure +!! flag, FF, is raised on return. +!! @author R. J. Purser +!! @param ak parameter vector +!! @param ma map-space domain-parameter vector +!! @param q lambda-parameterized quality diagnostic +!! @param ga geographic domain-parameter vector +!! @param lam Lambda +!! @param qdak derivatives value +!! @param qdma derivatives value +!! @param ff error flag subroutine get_meanqd(ngh,lam,xg,wg,ak,ma, q,qdak,qdma, & ! [get_meanq] ga,gadak,gadma, ff) -!============================================================================= -! For a parameter vector, ak and a map-space domain-parameter vector, ma, -! return the lambda-parameterized quality diagnostic, Q, and the geographic -! domain-parameter vector ga. Lambda is given by lam <1. Also, return -! the derivatives, qdak and qdma, of Q wrt ak and ma, and the derivatives -! gadak and gadma, of ga wrt ak and ma. -! The domain averages of Q are accurately computed by bi-Gauss-Legendre -! quadrature over the positive quadrant of the domain (exploiting the -! symmetry) of the four constituent terms, v1, v2, v3, v4, from which -! the mean Q is computed using a quadratic formula of these constituents. -! The number of Gauss points in eaxh half-interval is ngh, and the -! nodes themselves are, in proportion to the half-interval, at xg. -! the normalized gauss weights are wg. -! If a failure occurs, colmputations cease immediately and a failure -! flag, FF, is raised on return. -!============================================================================= implicit none integer(spi), intent(in ):: ngh real(dp), intent(in ):: lam @@ -544,7 +529,6 @@ subroutine get_meanqd(ngh,lam,xg,wg,ak,ma, q,qdak,qdma, & ! [get_meanq] real(dp),dimension(2), intent(out):: ga real(dp),dimension(2,2),intent(out):: gadak,gadma logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp),dimension(3,2,2):: xcd1 real(dp),dimension(3,2) :: xcd,xc1 real(dp),dimension(3) :: xc @@ -553,7 +537,6 @@ subroutine get_meanqd(ngh,lam,xg,wg,ak,ma, q,qdak,qdma, & ! [get_meanq] real(dp) :: wx,wy, & v1xy,v2xy,v3xy,v4xy, v1L,v2L,v3L,v4L, v1,v2,v3,v4 integer(spi) :: i,ic,ix,iy -!============================================================================= v1 =u0; v2 =u0; v3 =u0; v4 =u0 v1d=u0; v2d=u0; v3d=u0; v4d=u0 do iy=1,ngh @@ -576,7 +559,6 @@ subroutine get_meanqd(ngh,lam,xg,wg,ak,ma, q,qdak,qdma, & ! [get_meanq] enddo call get_qofv(lam,v1,v2,v3,v4, q)! <- Q(lam) based on the v1,v2,v3,v4 call get_qofv(lam,v2,v3, v1d,v2d,v3d,v4d, qdak)! <- Derivative of Q wrt ak -!------ ! Derivatives of ga wrt ak, and of q and ga wrt ma: gadma=u0! <- needed because only diagonal elements are filled do i=1,2 @@ -600,12 +582,11 @@ subroutine get_meanqd(ngh,lam,xg,wg,ak,ma, q,qdak,qdma, & ! [get_meanq] enddo call get_qofv(lam,v2,v3, v1d,v2d,v3d,v4d, qdma) end subroutine get_meanqd -!============================================================================= + +!> Like getmeanqd, except for n different values, aks, of ak and n different +!! values, mas of ma, and without any of the derivatives. +!! @author R. J. Purser subroutine get_meanqs(n,ngh,lam,xg,wg,aks,mas, qs,ff)! [get_meanq] -!============================================================================= -! Like getmeanqd, except for n different values, aks, of ak and n different -! values, mas of ma, and without any of the derivatives. -!============================================================================= implicit none integer(spi), intent(in ):: n,ngh real(dp),dimension(ngh),intent(in ):: xg,wg @@ -613,7 +594,6 @@ subroutine get_meanqs(n,ngh,lam,xg,wg,aks,mas, qs,ff)! [get_meanq] real(dp),dimension(2,n),intent(in ):: aks,mas real(dp),dimension(n), intent(out):: qs logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp),dimension(n) :: v1s,v2s,v3s,v4s real(dp),dimension(n) :: v1sL,v2sL,v3sL,v4sL real(dp),dimension(3,2):: xcd @@ -621,7 +601,6 @@ subroutine get_meanqs(n,ngh,lam,xg,wg,aks,mas, qs,ff)! [get_meanq] real(dp),dimension(2) :: xm,xgs real(dp) :: wx,wy, v1xy,v2xy,v3xy,v4xy integer(spi) :: i,ix,iy -!============================================================================= v1s=u0; v2s=u0; v3s=u0; v4s=u0 do iy=1,ngh wy=wg(iy) @@ -642,82 +621,76 @@ subroutine get_meanqs(n,ngh,lam,xg,wg,aks,mas, qs,ff)! [get_meanq] call get_qofv(n,lam,v1s,v2s,v3s,v4s, qs) end subroutine get_meanqs -!============================================================================= +!> The quadratic quantity Q depends linearly on v1 and v4 (which are already +!! quadratic diagnostics of EL) and quadratically on v2 and v3 (which are +!! linear diagnostics of EL). EL = (1/2)log(G), where G=J^T.J, J the jacobian. +!! @author R. J. Purser +!! @param q quadratic quantity +!! @param v1 quadratic diagnostics of EL +!! @param v4 quadratic diagnostics of EL +!! @param v2 linear diagnostics of EL +!! @param v3 linear diagnostics of EL subroutine get_qofv(lam,v1,v2,v3,v4, q)! [get_qofv] -!============================================================================= -! The quadratic quantity Q depends linearly on v1 and v4 (which are already -! quadratic diagnostics of EL) and quadratically on v2 and v3 (which are -! linear diagnostics of EL). EL = (1/2)log(G), where G=J^T.J, J the jacobian. -!============================================================================= implicit none real(dp),intent(in ):: lam,v1,v2,v3,v4 real(dp),intent(out):: q -!----------------------------------------------------------------------------- real(dp):: lamc -!============================================================================= lamc=u1-lam q=lamc*(v1-(v2**2+v3**2)) +lam*(v4 -(v2+v3)**2) end subroutine get_qofv -!============================================================================= + +!> Like get_qofv, but for (only) the 2-vector derivatives of Q. Note that +!! the quadratic diagnostics v1 and v4 do not participate in this formula. +!! @author R. J. Purser subroutine get_qofvd(lam, v2,v3, v1d,v2d,v3d,v4d, qd)! [get_qofv] -!============================================================================= -! Like get_qofv, but for (only) the 2-vector derivatives of Q. Note that -! the quadratic diagnostics v1 and v4 do not participate in this formula. -!============================================================================= implicit none real(dp), intent(in ):: lam,v2,v3 real(dp),dimension(2),intent(in ):: v1d,v2d,v3d,v4d real(dp),dimension(2),intent(out):: qd -!----------------------------------------------------------------------------- real(dp):: lamc -!============================================================================= lamc=u1-lam qd=lamc*(v1d-u2*(v2d*v2+v3d*v3))+lam*(v4d-u2*(v2d+v3d)*(v2+v3)) end subroutine get_qofvd -!============================================================================= + +!> General util to convert value +!! @author R. J. Purser subroutine get_qsofvs(n,lam,v1s,v2s,v3s,v4s, qs)! [get_qofv] -!============================================================================= implicit none integer(spi), intent(in ):: n real(dp), intent(in ):: lam real(dp),dimension(n),intent(in ):: v1s,v2s,v3s,v4s real(dp),dimension(n),intent(out):: qs -!----------------------------------------------------------------------------- real(dp):: lamc -!============================================================================= lamc=u1-lam qs=lamc*(v1s-(v2s**2+v3s**2)) +lam*(v4s -(v2s+v3s)**2) end subroutine get_qsofvs -!============================================================================= +!> Given an aspect ratio, asp<=1, and major semi-axis, arc, in map-space +!! nondimensional units, return a first guess for the parameter vector, ak, +!! approximately optimal for the domain of the given dimensions. +!! @author R. J. Purser +!! @param asp aspect ratio +!! @param ak first guess for the parameter vector subroutine guessak_map(asp,tmarcx,ak)! [guessak_map] -!============================================================================= -! Given an aspect ratio, asp<=1, and major semi-axis, arc, in map-space -! nondimensional units, return a first guess for the parameter vector, ak, -! approximately optimal for the domain of the given dimensions. -!============================================================================= implicit none real(dp), intent(in ):: asp,tmarcx real(dp),dimension(2),intent(out):: ak -!----------------------------------------------------------------------------- real(dp):: gmarcx -!============================================================================= gmarcx=tmarcx*rtod call guessak_geo(asp,gmarcx,ak) end subroutine guessak_map -!============================================================================= +!> Given an aspect ratio, asp<=1, and major semi-axis, arc, in geographical +!! (degree) units measured along the rectangle's median, return a first guess +!! for the parameter vector, ak, approximately optimal for the domain of the +!! given dimensions. +!! @author R. J. Purser +!! @param asp aspect ratio +!! @param ak first guess or the parameter vector subroutine guessak_geo(asp,arc,ak)! [guessak_geo] -!============================================================================= -! Given an aspect ratio, asp<=1, and major semi-axis, arc, in geographical -! (degree) units measured along the rectangle's median, return a first guess -! for the parameter vector, ak, approximately optimal for the domain of the -! given dimensions. -!============================================================================= implicit none real(dp), intent(in ):: asp,arc real(dp),dimension(2),intent(out):: ak -!----------------------------------------------------------------------------- integer(spi),parameter :: narc=11,nasp=10! <- Table index bounds real(dp), parameter :: eps=1.e-7_dp,darc=10._dp+eps,dasp=.1_dp+eps real(dp),dimension(nasp,0:narc):: adarc,kdarc @@ -775,32 +748,35 @@ subroutine guessak_geo(asp,arc,ak)! [guessak_geo] wx1*(wa0*kdarc(iasp0,iarc1)+wa1*kdarc(iasp1,iarc1))/) end subroutine guessak_geo -!============================================================================= +!> Get the best Extended Schmidt Gnomonic parameter, (a,k), for the given +!! geographical half-spans, garcx and garcy, as well as the corresponding +!! map-space half-spans, garcx and garcy (in degrees) and the quality +!! diagnostic, Q(lam) for this optimal parameter choice. If this process +!! fails for any reason, the failure is alerted by a raised flag, FF, and +!! the other output arguments must then be taken to be meaningless. +!! +!! The diagnostic Q measures the variance over the domain of a local measure +!! of grid distortion. A logarithmic measure of local grid deformation is +!! give by L=log(J^t.J)/2, where J is the mapping Jacobian matrix, dX/dx, +!! where X is the cartesian unit 3-vector representation of the image of the +!! mapping of the map-coordinate 2-vector, x. +!! The Frobenius squared-norm, Trace(L*L), of L is the basis for the simplest +!! (lam=0) definition of the variance of L, but (Trace(L))**2 is another. +!! Here, we weight both contributions, by lam and (1-lam) respectively, with +!! 0 <= lam <1, to compute the variance Q(lam,a,k), and search for the (a,k) +!! that minimizes this Q. +!! +!! The domain averages are computed by double Gauss-Legendre quadrature (i.e., +!! in both the x and y directions), but restricted to a mere quadrant of the +!! domain (since bilateral symmetry pertains across both domain medians, +!! yielding a domain mean L that is strictly diagonal. +!! @author R. J. Purser +!! @param a Extended Schmidt Gnomonic parameter +!! @param k Extended Schmidt Gnomonic parameter +!! @param ff failure flag +!! @param garcx map-space half-spans +!! @param garcy map-space half-spans subroutine bestesg_geo(lam,garcx,garcy, a,k,marcx,marcy,q,ff)! [bestesg_geo] -!============================================================================= -! Get the best Extended Schmidt Gnomonic parameter, (a,k), for the given -! geographical half-spans, garcx and garcy, as well as the corresponding -! map-space half-spans, garcx and garcy (in degrees) and the quality -! diagnostic, Q(lam) for this optimal parameter choice. If this process -! fails for any reason, the failure is alerted by a raised flag, FF, and -! the other output arguments must then be taken to be meaningless. -! -! The diagnostic Q measures the variance over the domain of a local measure -! of grid distortion. A logarithmic measure of local grid deformation is -! give by L=log(J^t.J)/2, where J is the mapping Jacobian matrix, dX/dx, -! where X is the cartesian unit 3-vector representation of the image of the -! mapping of the map-coordinate 2-vector, x. -! The Frobenius squared-norm, Trace(L*L), of L is the basis for the simplest -! (lam=0) definition of the variance of L, but (Trace(L))**2 is another. -! Here, we weight both contributions, by lam and (1-lam) respectively, with -! 0 <= lam <1, to compute the variance Q(lam,a,k), and search for the (a,k) -! that minimizes this Q. -! -! The domain averages are computed by double Gauss-Legendre quadrature (i.e., -! in both the x and y directions), but restricted to a mere quadrant of the -! domain (since bilateral symmetry pertains across both domain medians, -! yielding a domain mean L that is strictly diagonal. -!============================================================================= use pietc, only: u5,o5,s18,s36,s54,s72,ms18,ms36,ms54,ms72 use pmat, only: inv use psym2, only: chol2 @@ -808,7 +784,6 @@ subroutine bestesg_geo(lam,garcx,garcy, a,k,marcx,marcy,q,ff)! [bestesg_geo] real(dp),intent(in ):: lam,garcx,garcy real(dp),intent(out):: a,k,marcx,marcy,q logical ,intent(out):: FF -!----------------------------------------------------------------------------- integer(spi),parameter :: nit=200,mit=20,ngh=25 real(dp) ,parameter :: u2o5=u2*o5,& f18=u2o5*s18,f36=u2o5*s36,f54=u2o5*s54,& @@ -832,7 +807,6 @@ subroutine bestesg_geo(lam,garcx,garcy, a,k,marcx,marcy,q,ff)! [bestesg_geo] ! First guess upper-triangular basis regularizing the samples used to ! estimate the Hessian of q by finite differencing: data basis0/0.1_dp,u0, 0.3_dp,0.3_dp/ -!============================================================================= ff=lam=u1 if(ff)then; print'("In bestesg_geo; lam out of range")';return; endif ff= garcx<=u0 .or. garcy<=u0 @@ -936,39 +910,42 @@ subroutine bestesg_geo(lam,garcx,garcy, a,k,marcx,marcy,q,ff)! [bestesg_geo] endif end subroutine bestesg_geo -!============================================================================= +!> Get the best Extended Schmidt Gnomonic parameter, (a,k), for the given +!! map-coordinate half-spans, marcx and marcy, as well as the corresponding +!! geographical half-spans, garcx and garcy (in degrees) and the quality +!! diagnostic, Q(lam) for this optimal parameter choice. If this process +!! fails for any reason, the failure is alerted by a raised flag, FF, and +!! the other output arguments must then be taken to be meaningless. +!! +!! The diagnostic Q measures the variance over the domain of a local measure +!! of grid distortion. A logarithmic measure of local grid deformation is +!! give by L=log(J^t.J)/2, where J is the mapping Jacobian matrix, dX/dx, +!! where X is the cartesian unit 3-vector representation of the image of the +!! mapping of the map-coordinate 2-vector, x. +!! The Frobenius squared-norm, Trace(L*L), of L is the basis for the simplest +!! (lam=0) definition of the variance of L, but (Trace(L))**2 is another. +!! Here, we weight both contributions, by lam and (1-lam) respectively, with +!! 0 <= lam <1, to compute the variance Q(lam,a,k), and search for the (a,k) +!! that minimizes this Q. +!! +!! The domain averages are computed by double Gauss-Legendre quadrature (i.e., +!! in both the x and y directions), but restricted to a mere quadrant of the +!! domain (since bilateral symmetry pertains across both domain medians, +!! yielding a domain mean L that is strictly diagonal. +!! @author R. J. Purser +!! @param a Extended Schmidt Gnomonic parameter +!! @param k Extended Schmidt Gnomonic parameter +!! @param marcx map-coordinate half-spans +!! @param marcy map-coordinate half-spans +!! @param garcx geographical half-spans +!! @param garcy geographical half-spans subroutine bestesg_map(lam,marcx,marcy, a,k,garcx,garcy,q,ff) ![bestesg_map] -!============================================================================= -! Get the best Extended Schmidt Gnomonic parameter, (a,k), for the given -! map-coordinate half-spans, marcx and marcy, as well as the corresponding -! geographical half-spans, garcx and garcy (in degrees) and the quality -! diagnostic, Q(lam) for this optimal parameter choice. If this process -! fails for any reason, the failure is alerted by a raised flag, FF, and -! the other output arguments must then be taken to be meaningless. -! -! The diagnostic Q measures the variance over the domain of a local measure -! of grid distortion. A logarithmic measure of local grid deformation is -! give by L=log(J^t.J)/2, where J is the mapping Jacobian matrix, dX/dx, -! where X is the cartesian unit 3-vector representation of the image of the -! mapping of the map-coordinate 2-vector, x. -! The Frobenius squared-norm, Trace(L*L), of L is the basis for the simplest -! (lam=0) definition of the variance of L, but (Trace(L))**2 is another. -! Here, we weight both contributions, by lam and (1-lam) respectively, with -! 0 <= lam <1, to compute the variance Q(lam,a,k), and search for the (a,k) -! that minimizes this Q. -! -! The domain averages are computed by double Gauss-Legendre quadrature (i.e., -! in both the x and y directions), but restricted to a mere quadrant of the -! domain (since bilateral symmetry pertains across both domain medians, -! yielding a domain mean L that is strictly diagonal. -!============================================================================= use pietc, only: u5,o5,s18,s36,s54,s72,ms18,ms36,ms54,ms72 use psym2, only: chol2 implicit none real(dp),intent(in ):: lam,marcx,marcy real(dp),intent(out):: a,k,garcx,garcy,q logical ,intent(out):: FF -!----------------------------------------------------------------------------- integer(spi),parameter :: nit=25,mit=7,ngh=25 real(dp),parameter :: u2o5=u2*o5, & f18=u2o5*s18,f36=u2o5*s36,f54=u2o5*s54, & @@ -992,7 +969,6 @@ subroutine bestesg_map(lam,marcx,marcy, a,k,garcx,garcy,q,ff) ![bestesg_map] ! First guess upper-triangular basis regularizing the samples used to ! estimate the Hessian of q by finite differencing: data basis0/0.1_dp,u0, 0.3_dp,0.3_dp/ -!============================================================================= ff=lam=u1 if(ff)then; print'("In bestesg_map; lam out of range")';return; endif ff= marcx<=u0 .or. marcy<=u0 @@ -1073,30 +1049,40 @@ subroutine bestesg_map(lam,marcx,marcy, a,k,garcx,garcy,q,ff) ![bestesg_map] endif end subroutine bestesg_map -!============================================================================= +!> Use a and k as the parameters of an Extended Schmidt-transformed +!! Gnomonic (ESG) mapping centered at (plat,plon) and twisted about this center +!! by an azimuth angle of pazi counterclockwise (these angles in radians). +!! +!! Assume the radius of the earth is unity, and using the central mapping +!! point as the coordinate origin, set up the grid with central x-spacing delx +!! and y-spacing dely. The grid index location of the left-lower +!! corner of the domain is (lx,ly) (typically both NEGATIVE). +!! The numbers of the grid spaces in x and y directions are nx and ny. +!! (Note that, for a centered rectangular grid lx and ly are negative and, in +!! magnitude, half the values of nx and ny respectively.) +!! Return the latitude and longitude, in radians again, of the grid points thus +!! 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 +!! +!! If all goes well, return a lowered failure flag, ff=.false. . +!! But if, for some reason, it is not possible to complete this task, +!! return the raised failure flag, ff=.TRUE. . +!! @author R. J. Purser +!! @param a parameters of an ESG mapping centered at (plat,plon) +!! @param k parameters of an ESG mapping centered at (plat,plon) +!! @param plat latitude center point +!! @param plon longitude center point +!! @param lx grid index location +!! @param ly grid index location +!! @param nx grid spaces +!! @param ny grid spaces +!! @param glat grid points for latitude +!! @param glon grid points for longitude +!! @param garea rectangular array +!! @param ff failure flag subroutine hgrid_ak_rr(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rr] delx,dely, glat,glon,garea, ff) -!============================================================================= -! Use a and k as the parameters of an Extended Schmidt-transformed -! Gnomonic (ESG) mapping centered at (plat,plon) and twisted about this center -! by an azimuth angle of pazi counterclockwise (these angles in radians). -! -! Assume the radius of the earth is unity, and using the central mapping -! point as the coordinate origin, set up the grid with central x-spacing delx -! and y-spacing dely. The grid index location of the left-lower -! corner of the domain is (lx,ly) (typically both NEGATIVE). -! The numbers of the grid spaces in x and y directions are nx and ny. -! (Note that, for a centered rectangular grid lx and ly are negative and, in -! magnitude, half the values of nx and ny respectively.) -! Return the latitude and longitude, in radians again, of the grid points thus -! 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 -! -! If all goes well, return a lowered failure flag, ff=.false. . -! But if, for some reason, it is not possible to complete this task, -! return the raised failure flag, ff=.TRUE. . -!============================================================================= use pmat4, only: sarea use pmat5, only: ctogr implicit none @@ -1106,7 +1092,6 @@ subroutine hgrid_ak_rr(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rr] 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 logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp),dimension(3,3):: prot,azirot real(dp),dimension(3,2):: xcd real(dp),dimension(3) :: xc @@ -1115,7 +1100,6 @@ subroutine hgrid_ak_rr(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rr] rlat,drlata,drlatb,drlatc, & rlon,drlona,drlonb,drlonc integer(spi) :: ix,iy,mx,my -!============================================================================= clat=cos(plat); slat=sin(plat) clon=cos(plon); slon=sin(plon) cazi=cos(pazi); sazi=sin(pazi) @@ -1162,37 +1146,51 @@ subroutine hgrid_ak_rr(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rr] enddo enddo end subroutine hgrid_ak_rr -!============================================================================= + +!> Use a and k as the parameters of an extended Schmidt-transformed +!! gnomonic (ESG) mapping centered at (plat,plon) and twisted about this center +!! by an azimuth angle of pazi counterclockwise (these angles in radians). +!! +!! Using the central mapping point as the coordinate origin, set up the grid +!! with central x-spacing delx and y-spacing dely in nondimensional units, +!! (i.e., as if the earth had unit radius) and with the location of the left- +!! lower corner of the grid at center-relative grid index pair, (lx,ly) and +!! with the number of the grid spaces in x and y directions given by nx and ny. +!! (Note that, for a centered rectangular grid lx and ly are negative and, in +!! magnitude, half the values of nx and ny respectively.) +!! Return the latitude and longitude, again, in radians, of the grid pts thus +!! defined in the arrays, glat and glon; return a rectangular array, garea, +!! of dimensions nx-1 by ny-1, that contains the areas of each of the grid +!! cells in nondimensional "steradian" units. +!! +!! 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 and 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 radians 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. +!! @author R. J. Purser +!! @param a Extended Schmidt Gnomonic parameter +!! @param k Extended Schmidt Gnomonic parameter +!! @param plat latitude define centered mapping +!! @param plon longitude define centered mapping +!! @param delx central x-spacing grid point +!! @param dely central y-spacing grid point +!! @param lx x grid index for left-lower corner of the grid at center +!! @param ly y grid index for left-lower corner of the grid at center +!! @param nx numbers of the grid spaces in x +!! @param ny numbers of the grid spaces in y +!! @param dx estimated grid steps x grid cell edges +!! @param dy estimated grid steps y grid cell edges +!! @param angle_dx x angles relative to local east and north +!! @param angle_dy y angles relative to local east and north +!! @param ff failure flag subroutine hgrid_ak_rr_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak_rr] delx,dely, glat,glon,garea,dx,dy,angle_dx,angle_dy, ff) -!============================================================================= -! Use a and k as the parameters of an extended Schmidt-transformed -! gnomonic (ESG) mapping centered at (plat,plon) and twisted about this center -! by an azimuth angle of pazi counterclockwise (these angles in radians). -! -! Using the central mapping point as the coordinate origin, set up the grid -! with central x-spacing delx and y-spacing dely in nondimensional units, -! (i.e., as if the earth had unit radius) and with the location of the left- -! lower corner of the grid at center-relative grid index pair, (lx,ly) and -! with the number of the grid spaces in x and y directions given by nx and ny. -! (Note that, for a centered rectangular grid lx and ly are negative and, in -! magnitude, half the values of nx and ny respectively.) -! Return the latitude and longitude, again, in radians, of the grid pts thus -! defined in the arrays, glat and glon; return a rectangular array, garea, -! of dimensions nx-1 by ny-1, that contains the areas of each of the grid -! cells in nondimensional "steradian" units. -! -! 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 and 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 radians 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 pmat4, only: cross_product,triple_product use pmat5, only: ctogr implicit none @@ -1205,7 +1203,6 @@ subroutine hgrid_ak_rr_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak_rr] 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 @@ -1216,7 +1213,6 @@ subroutine hgrid_ak_rr_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak_rr] real(dp),dimension(2) :: xm real(dp) :: clat,slat,clon,slon,cazi,sazi,delxy integer(spi) :: ix,iy,mx,my,lxm,lym,mxp,myp -!============================================================================= delxy=delx*dely clat=cos(plat); slat=sin(plat) clon=cos(plon); slon=sin(plon) @@ -1359,30 +1355,43 @@ subroutine hgrid_ak_rr_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak_rr] -(gat(lx:mx-1,lym:my-2)+gat(lx:mx-1,ly+2:myp)))/24_dp end subroutine hgrid_ak_rr_c -!============================================================================= +!> Use a and k as the parameters of an Extended Schmidt-transformed +!! Gnomonic (ESG) mapping centered at (plat,plon) and twisted about this center +!! by an azimuth angle of pazi counterclockwise (these angles in radians). +!! +!! Assume the radius of the earth is unity, and using the central mapping +!! point as the coordinate origin, set up the grid with central x-spacing delx +!! and y-spacing dely. The grid index location of the left-lower +!! corner of the domain is (lx,ly) (typically both NEGATIVE). +!! The numbers of the grid spaces in x and y directions are nx and ny. +!! (Note that, for a centered rectangular grid lx and ly are negative and, in +!! magnitude, half the values of nx and ny respectively.) +!! Return the unit cartesian vectors xc of the grid points and their jacobian +!! matrices xcd wrt the map coordinates, and return a rectangular array, garea, +!! of dimensions nx-1 by ny-1, that contains the areas of each of the grid +!! cells +!! +!! If all goes well, return a lowered failure flag, ff=.false. . +!! But if, for some reason, it is not possible to complete this task, +!! return the raised failure flag, ff=.TRUE. . +!! @author R. J. Purser +!! @param a parameters of an ESG mapping centered at (plat,plon) +!! @param k parameters of an ESG mapping centered at (plat,plon) +!! @param plat latitude define centered mapping +!! @param plon longitude define centered mapping +!! @param delx central x-spacing grid point +!! @param dely central y-spacing grid point +!! @param lx x grid index for left-lower corner of the grid at center +!! @param ly y grid index for left-lower corner of the grid at center +!! @param nx numbers of the grid spaces in x +!! @param ny numbers of the grid spaces in y +!! @param xc unit cartesian vectors +!! @param garea rectangular array +!! @param nx-1 x dimensions for garea +!! @param ny-1 y dimensions for garea +!! @param ff failure flag subroutine hgrid_ak_rc(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rc] delx,dely, xc,xcd,garea, ff) -!============================================================================= -! Use a and k as the parameters of an Extended Schmidt-transformed -! Gnomonic (ESG) mapping centered at (plat,plon) and twisted about this center -! by an azimuth angle of pazi counterclockwise (these angles in radians). -! -! Assume the radius of the earth is unity, and using the central mapping -! point as the coordinate origin, set up the grid with central x-spacing delx -! and y-spacing dely. The grid index location of the left-lower -! corner of the domain is (lx,ly) (typically both NEGATIVE). -! The numbers of the grid spaces in x and y directions are nx and ny. -! (Note that, for a centered rectangular grid lx and ly are negative and, in -! magnitude, half the values of nx and ny respectively.) -! Return the unit cartesian vectors xc of the grid points and their jacobian -! matrices xcd wrt the map coordinates, and return a rectangular array, garea, -! of dimensions nx-1 by ny-1, that contains the areas of each of the grid -! cells -! -! If all goes well, return a lowered failure flag, ff=.false. . -! But if, for some reason, it is not possible to complete this task, -! return the raised failure flag, ff=.TRUE. . -!============================================================================= use pmat4, only: sarea use pmat5, only: ctogr implicit none @@ -1392,14 +1401,12 @@ subroutine hgrid_ak_rc(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rc] real(dp),dimension(3,2,lx:lx+nx ,ly:ly+ny ),intent(out):: xcd real(dp),dimension( lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp),dimension(3,3):: prot,azirot real(dp),dimension(2) :: xm real(dp) :: clat,slat,clon,slon,cazi,sazi, & rlat,rlata,rlatb,rlatc,drlata,drlatb,drlatc, & rlon,rlona,rlonb,rlonc,drlona,drlonb,drlonc integer(spi) :: ix,iy,mx,my -!============================================================================= clat=cos(plat); slat=sin(plat) clon=cos(plon); slon=sin(plon) cazi=cos(pazi); sazi=sin(pazi) @@ -1446,19 +1453,28 @@ subroutine hgrid_ak_rc(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rc] enddo end subroutine hgrid_ak_rc -!============================================================================= +!> Use a and k as the parameters of an Extended Schmidt-transformed +!! Gnomonic (ESG) mapping centered at (pdlat,pdlon) and twisted about this +!! center. +!! by an azimuth angle of pdazi counterclockwise (these angles in degrees). +!! Like hgrid_ak_rr, return the grid points' lats and lons, except that here +!! the angles are returned in degrees. Garea, the area of each grid cell, is +!! returned as in hgrid_ak_rr, and a failure flag, ff, raised when a failure +!! occurs anywhere in these calculations. +!! @author R. J. Purser +!! @param a parameters of an ESG mapping +!! @param k parameters of an ESG mapping +!! @param pdlat latitude define centered mapping +!! @param pdlon longitude define centered mapping +!! @param lx x grid index for left-lower corner of the grid at center +!! @param ly y grid index for left-lower corner of the grid at center +!! @param nx numbers of the grid spaces in x +!! @param ny numbers of the grid spaces in y +!! @param delx central x-spacing grid point +!! @param dely central y-spacing grid point +!! @param ff failure flag subroutine hgrid_ak_dd(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dd] delx,dely, gdlat,gdlon,garea, ff) -!============================================================================= -! Use a and k as the parameters of an Extended Schmidt-transformed -! Gnomonic (ESG) mapping centered at (pdlat,pdlon) and twisted about this -! center. -! by an azimuth angle of pdazi counterclockwise (these angles in degrees). -! Like hgrid_ak_rr, return the grid points' lats and lons, except that here -! the angles are returned in degrees. Garea, the area of each grid cell, is -! returned as in hgrid_ak_rr, and a failure flag, ff, raised when a failure -! occurs anywhere in these calculations. -!============================================================================ implicit none integer(spi), intent(in ):: lx,ly,nx,ny real(dp), intent(in ):: A,K,pdlat,pdlon,& @@ -1466,9 +1482,7 @@ subroutine hgrid_ak_dd(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dd] real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: gdlat,gdlon real(dp),dimension(lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp):: plat,plon,pazi -!============================================================================= plat=pdlat*dtor ! Convert these angles from degrees to radians plon=pdlon*dtor ! pazi=pdazi*dtor ! @@ -1478,13 +1492,25 @@ subroutine hgrid_ak_dd(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dd] gdlat=gdlat*rtod ! Convert these angles from radians to degrees gdlon=gdlon*rtod ! end subroutine hgrid_ak_dd -!============================================================================= + +!> Like hgrid_ak_rr_c, except all the angle arguments (but not delx,dely) +!! are in degrees instead of radians. +!! @author R. J. Purser +!! @param a parameters of an ESG mapping +!! @param k parameters of an ESG mapping +!! @param pdlat latitude define centered mapping +!! @param pdlon longitude define centered mapping +!! @param lx x grid index for left-lower corner of the grid at center +!! @param ly y grid index for left-lower corner of the grid at center +!! @param nx numbers of the grid spaces in x +!! @param ny numbers of the grid spaces in y +!! @param delx central x-spacing grid point +!! @param dely central y-spacing grid point +!! @param dangle_dx x rotations of the steps +!! @param dangle_dy y rotations of the steps +!! @param ff failure flag subroutine hgrid_ak_dd_c(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, &! [hgrid_ak_dd] delx,dely, gdlat,gdlon,garea,dx,dy,dangle_dx,dangle_dy, ff) -!============================================================================= -! Like hgrid_ak_rr_c, except all the angle arguments (but not delx,dely) -! are in degrees instead of radians. -!============================================================================= implicit none integer(spi), intent(in ):: lx,ly,nx,ny real(dp), intent(in ):: a,k,pdlat,pdlon,& @@ -1495,9 +1521,7 @@ subroutine hgrid_ak_dd_c(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, &! [hgrid_ak_dd] real(dp),dimension(lx:lx+nx ,ly:ly+ny-1),intent(out):: dy real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: dangle_dx,dangle_dy logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp):: plat,plon,pazi -!============================================================================= plat=pdlat*dtor ! Convert these angles from degrees to radians plon=pdlon*dtor ! pazi=pdazi*dtor ! @@ -1510,19 +1534,31 @@ subroutine hgrid_ak_dd_c(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, &! [hgrid_ak_dd] dangle_dy=dangle_dy*rtod ! end subroutine hgrid_ak_dd_c -!============================================================================= +!> Use a and k as the parameters of an Extended Schmidt-transformed +!! Gnomonic (ESG) mapping centered at (pdlat,pdlon) and twisted about this +!! center by an azimuth angle of pdazi counterclockwise (these angles in +!! degrees). +!! Like hgrid_ak_rx, return the grid points' cartesians xc and Jacobian +!! matrices, xcd. Garea, the area of each grid cell, is also +!! returned as in hgrid_ak_rx, and a failure flag, ff, raised when a failure +!! occurs anywhere in these calculations. +!! @author R. J. Purser +!! @param a parameters of an ESG mapping +!! @param k parameters of an ESG mapping +!! @param pdlat latitude define centered mapping +!! @param pdlon longitude define centered mapping +!! @param lx x grid index for left-lower corner of the grid at center +!! @param ly y grid index for left-lower corner of the grid at center +!! @param nx numbers of the grid spaces in x +!! @param ny numbers of the grid spaces in y +!! @param delx central x-spacing grid point +!! @param dely central y-spacing grid point +!! @param xc grid points' cartesians +!! @param xcd Jacobian matrices +!! @param garea Jacobian matrices +!! @param ff failure flag subroutine hgrid_ak_dc(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dc] delx,dely, xc,xcd,garea, ff) -!============================================================================= -! Use a and k as the parameters of an Extended Schmidt-transformed -! Gnomonic (ESG) mapping centered at (pdlat,pdlon) and twisted about this -! center by an azimuth angle of pdazi counterclockwise (these angles in -! degrees). -! Like hgrid_ak_rx, return the grid points' cartesians xc and Jacobian -! matrices, xcd. Garea, the area of each grid cell, is also -! returned as in hgrid_ak_rx, and a failure flag, ff, raised when a failure -! occurs anywhere in these calculations. -!============================================================================ implicit none integer(spi),intent(in ):: lx,ly,nx,ny real(dp), intent(in ):: A,K,pdlat,pdlon,pdazi,delx,dely @@ -1530,9 +1566,7 @@ subroutine hgrid_ak_dc(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dc] real(dp),dimension(3,2,lx:lx+nx ,ly:ly+ny ),intent(out):: xcd real(dp),dimension( lx:lx+nx-1,ly:ly+ny-1),intent(out):: garea logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp):: plat,plon,pazi -!============================================================================= plat=pdlat*dtor plon=pdlon*dtor pazi=pdazi*dtor @@ -1540,16 +1574,26 @@ subroutine hgrid_ak_dc(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dc] delx,dely, xc,xcd,garea, ff) end subroutine hgrid_ak_dc -!============================================================================= +!> Like hgrid_ak_rr_c except the argument list includes the earth radius, re, +!! and this is used to express the map-space grid increments in the dimensional +!! units, delxre, delyre on entry, and to express the grid cell areas, garea, +!! in dimensional units upon return. +!! The gridded lats and lons, glat and glon, remain in radians. +!! @author R. J. Purser +!! @param a parameters of an ESG mapping +!! @param k parameters of an ESG mapping +!! @param pdlat latitude define centered mapping +!! @param pdlon longitude define centered mapping +!! @param lx x grid index for left-lower corner of the grid at center +!! @param ly y grid index for left-lower corner of the grid at center +!! @param nx numbers of the grid spaces in x +!! @param ny numbers of the grid spaces in y +!! @param re earth radius +!! @param delxre map-space grid increments in the dimensional units +!! @param delyre map-space grid increments in the dimensional units +!! @param ff failure flag subroutine hgrid_ak(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak] re,delxre,delyre, glat,glon,garea, ff) -!============================================================================= -! Like hgrid_ak_rr_c except the argument list includes the earth radius, re, -! and this is used to express the map-space grid increments in the dimensional -! units, delxre, delyre on entry, and to express the grid cell areas, garea, -! in dimensional units upon return. -! The gridded lats and lons, glat and glon, remain in radians. -!============================================================================ implicit none integer(spi), intent(in ):: lx,ly,nx,ny real(dp), intent(in ):: a,k,plat,plon,pazi, & @@ -1557,9 +1601,7 @@ subroutine hgrid_ak(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak] 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 logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp):: delx,dely,rere -!============================================================================= delx=delxre/re ! <- set nondimensional grid step delx dely=delyre/re ! <- set nondimensional grid step dely call hgrid_ak_rr(lx,ly,nx,ny,a,k,plat,plon,pazi, & @@ -1569,20 +1611,32 @@ subroutine hgrid_ak(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak] garea=garea*rere ! <- Convert from steradians to physical area units. end subroutine hgrid_ak -!============================================================================= +!> Like hgrid_ak_rr_c except the argument list includes the earth radius, re, +!! and this is used to express the map-space grid increments in the dimensional +!! units, delxre, delyre on entry, and to express the grid cell areas, garea, +!! and the x- and y- grid steps, dx and dy, in dimensional units upon return. +!! The gridded lats and lons, glat and glon, remain in radians. +!! Also, in order for the argument list to remain compatible with an earlier +!! version of this routine, the relative rotations of the steps, dangle_dx +!! and dangle_dy, are returned as degrees instead of radians (all other angles +!! in the argument list, i.e., plat,plon,pazi,glat,glon, remain radians). +!! @author R. J. Purser +!! @param a Extended Schmidt Gnomonic parameter +!! @param k Extended Schmidt Gnomonic parameter +!! @param plat latitude define centered mapping +!! @param plon longitude define centered mapping +!! @param re earth radius +!! @param delxre map-space grid increments in the dimensional units +!! @param delyre map-space grid increments in the dimensional units +!! @param garea grid cell areas +!! @param dx x- grid steps +!! @param dy y- grid steps +!! @param glat gridded lats +!! @param glon gridded lons +!! @param dangle_dx x rotations of the steps +!! @param dangle_dy y rotations of the steps subroutine hgrid_ak_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak] re,delxre,delyre, glat,glon,garea,dx,dy,dangle_dx,dangle_dy, ff) -!============================================================================= -! Like hgrid_ak_rr_c except the argument list includes the earth radius, re, -! and this is used to express the map-space grid increments in the dimensional -! units, delxre, delyre on entry, and to express the grid cell areas, garea, -! and the x- and y- grid steps, dx and dy, in dimensional units upon return. -! The gridded lats and lons, glat and glon, remain in radians. -! Also, in order for the argument list to remain compatible with an earlier -! version of this routine, the relative rotations of the steps, dangle_dx -! and dangle_dy, are returned as degrees instead of radians (all other angles -! in the argument list, i.e., plat,plon,pazi,glat,glon, remain radians). -!============================================================================= implicit none integer(spi), intent(in ):: lx,ly,nx,ny real(dp), intent(in ):: a,k,plat,plon,pazi, & @@ -1593,9 +1647,7 @@ subroutine hgrid_ak_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak] real(dp),dimension(lx:lx+nx ,ly:ly+ny-1),intent(out):: dy real(dp),dimension(lx:lx+nx ,ly:ly+ny ),intent(out):: dangle_dx,dangle_dy logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp):: delx,dely,rere -!============================================================================= delx=delxre/re ! <- set nondimensional grid step delx dely=delyre/re ! <- set nondimensional grid step dely call hgrid_ak_rr_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & @@ -1609,22 +1661,20 @@ subroutine hgrid_ak_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak] dangle_dy=dangle_dy*rtod ! <-Convert from radians to degrees end subroutine hgrid_ak_c -!============================================================================= +!> This Gauss-Legendre quadrature integrates exactly any even polynomial +!! up to degree m*4-2 in the half-interval [0,1]. This code is liberally +!! adapted from the algorithm given in Press et al., Numerical Recipes. +!! @param m number of nodes in half-interval +!! @param x nodes and weights +!! @param w nodes and weights subroutine gaulegh(m,x,w)! [gaulegh] -!============================================================================= -! This Gauss-Legendre quadrature integrates exactly any even polynomial -! up to degree m*4-2 in the half-interval [0,1]. This code is liberally -! adapted from the algorithm given in Press et al., Numerical Recipes. -!============================================================================= implicit none integer(spi), intent(IN ):: m ! <- number of nodes in half-interval real(dp),dimension(m),intent(OUT):: x,w ! <- nodes and weights -!----------------------------------------------------------------------------- integer(spi),parameter:: nit=8 real(dp), parameter:: eps=3.e-14_dp integer(spi) :: i,ic,j,jm,it,m2,m4p,m4p3 real(dp) :: z,zzm,p1,p2,p3,pp,z1 -!============================================================================= m2=m*2; m4p=m*4+1; m4p3=m4p+2 do i=1,m; ic=m4p3-4*i z=cos(pih*ic/m4p); zzm=z*z-u1 @@ -1638,14 +1688,19 @@ subroutine gaulegh(m,x,w)! [gaulegh] enddo end subroutine gaulegh -!============================================================================= +!> Given the map specification (angles in radians), the grid spacing (in +!! map-space units) and the sample lat-lon (in radian), return the the +!! image in map space in a 2-vector in grid units. If the transformation +!! is invalid, return a .true. failure flag. +!! @author R. J. Purser +!! @param a parameters of an ESG mapping +!! @param k parameters of an ESG mapping +!! @param plat latitude define centered mapping +!! @param plon longitude define centered mapping +!! @param lat radian latitude +!! @param lon radian longitude +!! @param ff failure flag subroutine gtoxm_ak_rr_m(A,K,plat,plon,pazi,lat,lon,xm,ff)! [gtoxm_ak_rr] -!============================================================================= -! Given the map specification (angles in radians), the grid spacing (in -! map-space units) and the sample lat-lon (in radian), return the the -! image in map space in a 2-vector in grid units. If the transformation -! is invalid, return a .true. failure flag. -!============================================================================= use pmat5, only: grtoc implicit none real(dp), intent(in ):: a,k,plat,plon,pazi,lat,lon @@ -1654,7 +1709,6 @@ subroutine gtoxm_ak_rr_m(A,K,plat,plon,pazi,lat,lon,xm,ff)! [gtoxm_ak_rr] real(dp),dimension(3,3):: prot,azirot real(dp) :: clat,slat,clon,slon,cazi,sazi real(dp),dimension(3) :: xc -!============================================================================= clat=cos(plat); slat=sin(plat) clon=cos(plon); slon=sin(plon) cazi=cos(pazi); sazi=sin(pazi) @@ -1672,37 +1726,47 @@ subroutine gtoxm_ak_rr_m(A,K,plat,plon,pazi,lat,lon,xm,ff)! [gtoxm_ak_rr] xc=matmul(transpose(prot),xc) call xctoxm_ak(a,k,xc,xm,ff) end subroutine gtoxm_ak_rr_m -!============================================================================= + +!> Given the map specification (angles in radians), the grid spacing (in +!! map-space units) and the sample lat-lon (in radian), return the the +!! image in map space in a 2-vector in grid units. If the transformation +!! is invalid, return a .true. failure flag. +!! @author R. J. Purser +!! @param a parameters of an ESG mapping +!! @param k parameters of an ESG mapping +!! @param plat latitude define centered mapping +!! @param plon longitude define centered mapping +!! @param delx central x-spacing grid point +!! @param dely central y-spacing grid point +!! @param lat radian latitude +!! @param lon radian longitude +!! @param ff failure flag subroutine gtoxm_ak_rr_g(A,K,plat,plon,pazi,delx,dely,lat,lon,&! [gtoxm_ak_rr] xm,ff) -!============================================================================= -! Given the map specification (angles in radians), the grid spacing (in -! map-space units) and the sample lat-lon (in radian), return the the -! image in map space in a 2-vector in grid units. If the transformation -! is invalid, return a .true. failure flag. -!============================================================================= implicit none real(dp), intent(in ):: a,k,plat,plon,pazi,delx,dely,lat,lon real(dp),dimension(2),intent(out):: xm logical, intent(out):: ff -!============================================================================= call gtoxm_ak_rr_m(A,K,plat,plon,pazi,lat,lon,xm,ff); if(ff)return xm(1)=xm(1)/delx; xm(2)=xm(2)/dely end subroutine gtoxm_ak_rr_g -!============================================================================= +!> Like gtoxm_ak_rr_m, except input angles are expressed in degrees +!! @author R. J. Purser +!! @param a parameters of an ESG mapping +!! @param k parameters of an ESG mapping +!! @param pdlat latitude define centered mapping +!! @param pdlon longitude define centered mapping +!! @param dlat radian latitude +!! @param dlon radian longitude +!! @param ff failure flag subroutine gtoxm_ak_dd_m(A,K,pdlat,pdlon,pdazi,dlat,dlon,&! [gtoxm_ak_dd] xm,ff) -!============================================================================= -! Like gtoxm_ak_rr_m, except input angles are expressed in degrees -!============================================================================= implicit none real(dp), intent(in ):: a,k,pdlat,pdlon,pdazi,dlat,dlon real(dp),dimension(2),intent(out):: xm logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp):: plat,plon,pazi,lat,lon -!============================================================================= plat=pdlat*dtor ! Convert these angles from degrees to radians plon=pdlon*dtor ! pazi=pdazi*dtor ! @@ -1710,19 +1774,23 @@ subroutine gtoxm_ak_dd_m(A,K,pdlat,pdlon,pdazi,dlat,dlon,&! [gtoxm_ak_dd] lon=dlon*dtor call gtoxm_ak_rr_m(A,K,plat,plon,pazi,lat,lon,xm,ff) end subroutine gtoxm_ak_dd_m -!============================================================================= + +!> Like gtoxm_ak_rr_g, except input angles are expressed in degrees +!! @author R. J. Purser +!! @param a parameters of an ESG mapping +!! @param k parameters of an ESG mapping +!! @param pdlat latitude define centered mapping +!! @param pdlon longitude define centered mapping +!! @param delx central x-spacing grid point +!! @param dely central y-spacing grid point +!! @param ff failure flag subroutine gtoxm_ak_dd_g(A,K,pdlat,pdlon,pdazi,delx,dely,&! [gtoxm_ak_dd] dlat,dlon, xm,ff) -!============================================================================= -! Like gtoxm_ak_rr_g, except input angles are expressed in degrees -!============================================================================= implicit none real(dp), intent(in ):: a,k,pdlat,pdlon,pdazi,delx,dely,dlat,dlon real(dp),dimension(2),intent(out):: xm logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp):: plat,plon,pazi,lat,lon -!============================================================================= plat=pdlat*dtor ! Convert these angles from degrees to radians plon=pdlon*dtor ! pazi=pdazi*dtor ! @@ -1731,27 +1799,30 @@ subroutine gtoxm_ak_dd_g(A,K,pdlat,pdlon,pdazi,delx,dely,&! [gtoxm_ak_dd] call gtoxm_ak_rr_g(A,K,plat,plon,pazi,delx,dely,lat,lon,xm,ff) end subroutine gtoxm_ak_dd_g -!============================================================================= +!> Given the ESG map specified by parameters (A,K) and geographical +!! orientation, plat,plon,pazi (radians), and a position, in map-space +!! coordinates given by the 2-vector, xm, return the geographical +!! coordinates, lat and lon (radians). If the transformation is invalid for +!! any reason, return instead with a raised failure flag, FF= .true. +!! @author R. J. Purser +!! @param a parameters of an ESG mapping +!! @param k parameters of an ESG mapping +!! @param plat latitude of geographical orientation +!! @param plon longitude of geographical orientation +!! @param lat latitude of geographical coordinates +!! @param lon longitude of geographical coordinates +!! @param ff failure flag subroutine xmtog_ak_rr_m(A,K,plat,plon,pazi,xm,lat,lon,ff)! [xmtog_ak_rr] -!============================================================================= -! Given the ESG map specified by parameters (A,K) and geographical -! orientation, plat,plon,pazi (radians), and a position, in map-space -! coordinates given by the 2-vector, xm, return the geographical -! coordinates, lat and lon (radians). If the transformation is invalid for -! any reason, return instead with a raised failure flag, FF= .true. -!============================================================================= use pmat5, only: ctogr implicit none real(dp), intent(in ):: a,k,plat,plon,pazi real(dp),dimension(2),intent(in ):: xm real(dp), intent(out):: lat,lon logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp),dimension(3,2):: xcd real(dp),dimension(3,3):: prot,azirot real(dp) :: clat,slat,clon,slon,cazi,sazi real(dp),dimension(3) :: xc -!============================================================================= clat=cos(plat); slat=sin(plat) clon=cos(plon); slon=sin(plon) cazi=cos(pazi); sazi=sin(pazi) @@ -1768,43 +1839,51 @@ subroutine xmtog_ak_rr_m(A,K,plat,plon,pazi,xm,lat,lon,ff)! [xmtog_ak_rr] xc=matmul(prot,xc) call ctogr(xc,lat,lon) end subroutine xmtog_ak_rr_m -!============================================================================= + +!> For an ESG map with parameters, (A,K), and geographical orientation, +!! given by plon,plat,pazi (radians), and given a point in grid-space units +!! as the 2-vector, xm, return the geographical coordinates, lat, lon, (radians) +!! of this point. If instead the transformation is invalid for any reason, then +!! return the raised failure flag, FF=.true. +!! @author R. J. Purser +!! @param a parameters of an ESG mapping +!! @param k parameters of an ESG mapping +!! @param plon longitude geographical orientation +!! @param plat latitude geographical orientation +!! @param xm grid-space units as the 2-vector +!! @param lat latitude geographical coordinates +!! @param lon longitude geographical coordinates +!! @param delx central x-spacing grid point +!! @param dely central y-spacing grid point +!! @param ff failure flag subroutine xmtog_ak_rr_g(A,K,plat,plon,pazi,delx,dely,xm,&! [xmtog_ak_rr] lat,lon,ff) -!============================================================================= -! For an ESG map with parameters, (A,K), and geographical orientation, -! given by plon,plat,pazi (radians), and given a point in grid-space units -! as the 2-vector, xm, return the geographical coordinates, lat, lon, (radians) -! of this point. If instead the transformation is invalid for any reason, then -! return the raised failure flag, FF=.true. -!============================================================================= implicit none real(dp), intent(in ):: a,k,plat,plon,pazi,delx,dely real(dp),dimension(2),intent(in ):: xm real(dp), intent(out):: lat,lon logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp),dimension(2):: xmt -!============================================================================= xmt(1)=xm(1)*delx ! Convert from grid units to intrinsic map-space units xmt(2)=xm(2)*dely ! call xmtog_ak_rr_m(A,K,plat,plon,pazi,xmt,lat,lon,ff) end subroutine xmtog_ak_rr_g -!============================================================================= +!> Like xmtog_ak_rr_m, except angles are expressed in degrees +!! @author R. J. Purser +!! @param a parameters of an ESG mapping +!! @param k parameters of an ESG mapping +!! @param plat latitude of geographical orientation +!! @param plon longitude of geographical orientation +!! @param ff failure flag subroutine xmtog_ak_dd_m(A,K,pdlat,pdlon,pdazi,xm,dlat,dlon,ff)! [xmtog_ak_dd] -!============================================================================= -! Like xmtog_ak_rr_m, except angles are expressed in degrees -!============================================================================= use pmat5, only: ctogr implicit none real(dp), intent(in ):: a,k,pdlat,pdlon,pdazi real(dp),dimension(2),intent(in ):: xm real(dp), intent(out):: dlat,dlon logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp):: plat,plon,pazi,lat,lon -!============================================================================= plat=pdlat*dtor ! Convert these angles from degrees to radians plon=pdlon*dtor ! pazi=pdazi*dtor ! @@ -1812,21 +1891,28 @@ subroutine xmtog_ak_dd_m(A,K,pdlat,pdlon,pdazi,xm,dlat,dlon,ff)! [xmtog_ak_dd] dlat=lat*rtod dlon=lon*rtod end subroutine xmtog_ak_dd_m -!============================================================================= + +!> Like xmtog_ak_rr_g, except angles are expressed in degrees +!! @author R. J. Purser +!! @param a parameters of an ESG mapping +!! @param k parameters of an ESG mapping +!! @param pdlat latitude define centered mapping +!! @param pdlon longitude define centered mapping +!! @param dlat latitude radians angle +!! @param dlon longitude radians angle +!! @param xm grid-space +!! @param delx central x-spacing grid point +!! @param dely central y-spacing grid point +!! @param ff failure flag subroutine xmtog_ak_dd_g(A,K,pdlat,pdlon,pdazi,delx,dely,xm,&! [xmtog_ak_dd] dlat,dlon,ff) -!============================================================================= -! Like xmtog_ak_rr_g, except angles are expressed in degrees -!============================================================================= implicit none real(dp), intent(in ):: a,k,pdlat,pdlon,pdazi,delx,dely real(dp),dimension(2),intent(in ):: xm real(dp), intent(out):: dlat,dlon logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp),dimension(2):: xmt real(dp) :: plat,plon,pazi,lat,lon -!============================================================================= xmt(1)=xm(1)*delx ! Convert from grid units to intrinsic map-space units xmt(2)=xm(2)*dely ! plat=pdlat*dtor ! Convert these angles from degrees to radians From c67d3ca480a40ea76001b7cf68023c356701978b Mon Sep 17 00:00:00 2001 From: "Lin.Gan" Date: Tue, 2 Mar 2021 15:16:18 +0000 Subject: [PATCH 2/5] Part of Doxygen updates for sorc/grid_tools.fd #366 doxygen update for sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 --- .../regional_esg_grid.fd/pesg.f90 | 84 +++++++++---------- 1 file changed, 42 insertions(+), 42 deletions(-) diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 index c92a5bd4a..f4be4550f 100644 --- a/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 +++ b/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 @@ -1,7 +1,4 @@ !> @file -!! @author R. J. Purser -!! @date May 2020 -!! !! Suite of routines to perform the 2-parameter family of Extended !! Schmidt Gnomonic (ESG) regional grid mappings, and to optimize the !! the two parameters, A and K, of those mappings for a given rectangular @@ -10,11 +7,13 @@ !! lam (for "lambda" in [0,1) ) which gives weight to additional weight !! areal inhomogeneities instead of treating all distortion components !! equally. -!! !! DEPENDENCIES !! Libraries: pmat, psym2, pfun !! Modules: pkind, pietc, pietc_s -!! +!! @author R. J. Purser @date May 2020 + +!> Module contain subroutines +!! @author R. J. Purser module pesg !============================================================================= use pkind, only: spi,dp @@ -63,9 +62,9 @@ module pesg contains !> Inverse of xstoxc. I.e., cartesians to stereographic -!! @author R. J. Purser !! @param xc xstoxc !! @param xs Inverse of xstoxc +!! @author R. J. Purser subroutine xctoxs(xc,xs)! [xctoxs] implicit none real(dp),dimension(3),intent(in ):: xc @@ -77,10 +76,10 @@ end subroutine xctoxs !> Standard transformation from polar stereographic map coordinates, xs, to !! cartesian, xc, assuming the projection axis is polar. !! xcd=d(xc)/d(xs) is the jacobian matrix, encoding distortion and metric. -!! @author R. J. Purser !! @param xs polar stereographic map coordinates !! @param xc cartesian !! @param xcd value of jacobian matrix, encoding distortion and metric +!! @author R. J. Purser subroutine xstoxc(xs,xc,xcd)! [xstoxc] use pmat4, only: outer_product implicit none @@ -97,11 +96,11 @@ end subroutine xstoxc !! cartesian, xc, assuming the projection axis is polar. !! xcd=d(xc)/d(xs) is the jacobian matrix, encoding distortion and metric. !! xcdd is the further derivative, wrt xs, of xcd. -!! @author R. J. Purser !! @param xs polar stereographic map coordinates !! @param xc cartesian !! @param xcd jacobian matrix, encoding distortion and metric !! @param xcdd further derivative, wrt xs, of xcd +!! @author R. J. Purser subroutine xstoxc1(xs,xc,xcd,xcdd)! [xstoxc] use pmat4, only: outer_product implicit none @@ -142,9 +141,9 @@ subroutine xstoxt(k,xs,xt,ff)! [xstoxt] end subroutine xstoxt !> Scaled gnomonic plane xt to standard stereographic plane xs -!! @author R. J. Purser !! @param xt Scaled gnomonic plane !! @param xs standard stereographic plane +!! @author R. J. Purser subroutine xttoxs(k,xt,xs,xsd,ff)! [xttoxs use pmat4, only: outer_product implicit none @@ -168,10 +167,10 @@ subroutine xttoxs(k,xt,xs,xsd,ff)! [xttoxs end subroutine xttoxs !> Like xttoxs, but also, return the derivatives, wrt K, of xs and xsd -!! @author R. J. Purser !! @param xt Scaled gnomonic plane !! @param xs standard stereographic plane !! @param K derivatives, wrt K, of xs and xsd +!! @author R. J. Purser subroutine xttoxs1(k,xt,xs,xsd,xsdd,xs1,xsd1,ff)! [xttoxs] use pmat4, only: outer_product implicit none @@ -220,10 +219,10 @@ subroutine xttoxm(a,xt,xm,ff)! [xttoxm] end subroutine xttoxm !> Like zmtozt, but for 2-vector xm and xt, and 2*2 diagonal Jacobian xtd -!! @author R. J. Purser !! @param xm vector value !! @param xt vector value !! @param xtd 2*2 diagonal Jacobian +!! @author R. J. Purser subroutine xmtoxt(a,xm,xt,xtd,ff)! [xmtoxt] implicit none real(dp), intent(in ):: a @@ -237,11 +236,11 @@ end subroutine xmtoxt !> Like zmtozt1, but for 2-vector xm and xt, and 2*2 diagonal Jacobian xtd !! Also, the derivatives, wrt a, of these quantities. -!! @author R. J. Purser !! @param xm vector value !! @param xt vector value !! @param xtd 2*2 diagonal Jacobian !! @param a the derivatives of these quantities +!! @author R. J. Purser subroutine xmtoxt1(a,xm,xt,xtd,xt1,xtd1,ff)! [xmtoxt] implicit none real(dp), intent(in ):: a @@ -276,9 +275,9 @@ end subroutine zttozm !> Evaluate the function, zt = tan(sqrt(A)*z)/sqrt(A), and its derivative, ztd, !! for positive and negative A and for the limiting case, A --> 0 -!! @author R. J. Purser !! @param zt function to be evaluated !! @param ztd derivative of the source function +!! @author R. J. Purser subroutine zmtozt(a,zm,zt,ztd,ff)! [zmtozt] implicit none real(dp),intent(in ):: a,zm @@ -323,10 +322,10 @@ end subroutine zmtozt1 !! derivative wrt xm, the Jacobian matrix, xcd. !! If for any reason the mapping cannot be done, return a raised failure flag,z !! FF. -!! @author R. J. Purser !! @param ak parameterization of the Extended ESG !! @param xm vector !! @param ff raised failure error code +!! @author R. J. Purser subroutine xmtoxc_vak(ak,xm,xc,xcd,ff)! [xmtoxc_ak] implicit none real(dp),dimension(2), intent(in ):: ak,xm @@ -337,8 +336,8 @@ subroutine xmtoxc_vak(ak,xm,xc,xcd,ff)! [xmtoxc_ak] end subroutine xmtoxc_vak !> Like xmtoxc_vak, _ak, but also return derivatives wrt ak. -!! @author R. J. Purser !! @param ak derivatives wrt ak +!! @author R. J. Purser subroutine xmtoxc_vak1(ak,xm,xc,xcd,xc1,xcd1,ff)! [xmtoxc_ak] implicit none real(dp),dimension(2), intent(in ):: ak,xm @@ -370,13 +369,13 @@ end subroutine xmtoxc_vak1 !! corresponding cartesian unit 3-vector and its derivative wrt xm, jacobian !! matrix, xcd. If for any reason the mapping cannot be done, return a !! raised failure flag, FF. -!! @author R. J. Purser !! @param a ESG mapping parameterization !! @param k ESG mapping parameterization !! @param xm map-space vector !! @param xm derivative !! @param xcd jacobian matrix !! @param ff error flag +!! @author R. J. Purser subroutine xmtoxc_ak(a,k,xm,xc,xcd,ff)! [xmtoxc_ak] implicit none real(dp), intent(in ):: a,k @@ -396,12 +395,12 @@ end subroutine xmtoxc_ak !> Inverse mapping of xmtoxc_ak. That is, go from given cartesian unit !! 3-vector, xc, to map coordinate 2-vector xm (or return a raised failure !! flag, FF, if the attempt fails). -!! @author R. J. Purser !! @param a Inverse mapping of xmtoxc !! @param k Inverse mapping of xmtoxc !! @param xc cartesian unit 3-vector !! @param xm 2-vector map coordinate !! @param ff error flag +!! @author R. J. Purser subroutine xctoxm_ak(a,k,xc,xm,ff)! [xctoxm_ak] implicit none real(dp), intent(in ):: a,k @@ -419,9 +418,9 @@ end subroutine xctoxm_ak !! region's center and its x and y edges, get the two cartesian vectors !! that represent the locations of these edge midpoints in the positive x and y !! directions. -!! @author R. J. Purser !! @param edgex region's x edges !! @param edgey region's y edges +!! @author R. J. Purser subroutine get_edges(arcx,arcy,edgex,edgey)! [get_edges] implicit none real(dp), intent(in ):: arcx,arcy @@ -438,12 +437,12 @@ end subroutine get_edges !! of Q(lam) that they imply for any choice of the "lambda" parameter, lam. !! Note that v1 and v4 are quadratic diagnostics of EL, while v2 and v3 are !! linear. -!! @author R. J. Purser !! @param j0 jacobian matrix !! @param v1 quadratic diagnostics of EL !! @param v4 quadratic diagnostics of EL !! @param v2 linear diagnostics of EL !! @param v3 linear diagnostics of EL +!! @author R. J. Purser subroutine get_qx(j0, v1,v2,v3,v4)! [get_qx] use psym2, only: logsym2 implicit none @@ -461,9 +460,9 @@ end subroutine get_qx !> From a jacobian matrix, j0, and its derivative, j0d, get a sufficient set !! of v.. diagnostics such that, from average of these diagnostics, we can !! later compute the collective variance of Q and its derivative. -!! @author R. J. Purser !! @param j0 jacobian matrix !! @param j0d derivative of j0 +!! @author R. J. Purser subroutine get_qxd(j0,j0d, v1,v2,v3,v4,v1d,v2d,v3d,v4d)! [get_qx] use psym2, only: logsym2 implicit none @@ -508,7 +507,6 @@ end subroutine get_qxd !! the normalized gauss weights are wg. !! If a failure occurs, colmputations cease immediately and a failure !! flag, FF, is raised on return. -!! @author R. J. Purser !! @param ak parameter vector !! @param ma map-space domain-parameter vector !! @param q lambda-parameterized quality diagnostic @@ -517,6 +515,7 @@ end subroutine get_qxd !! @param qdak derivatives value !! @param qdma derivatives value !! @param ff error flag +!! @author R. J. Purser subroutine get_meanqd(ngh,lam,xg,wg,ak,ma, q,qdak,qdma, & ! [get_meanq] ga,gadak,gadma, ff) implicit none @@ -624,12 +623,12 @@ end subroutine get_meanqs !> The quadratic quantity Q depends linearly on v1 and v4 (which are already !! quadratic diagnostics of EL) and quadratically on v2 and v3 (which are !! linear diagnostics of EL). EL = (1/2)log(G), where G=J^T.J, J the jacobian. -!! @author R. J. Purser !! @param q quadratic quantity !! @param v1 quadratic diagnostics of EL !! @param v4 quadratic diagnostics of EL !! @param v2 linear diagnostics of EL !! @param v3 linear diagnostics of EL +!! @author R. J. Purser subroutine get_qofv(lam,v1,v2,v3,v4, q)! [get_qofv] implicit none real(dp),intent(in ):: lam,v1,v2,v3,v4 @@ -668,9 +667,9 @@ end subroutine get_qsofvs !> Given an aspect ratio, asp<=1, and major semi-axis, arc, in map-space !! nondimensional units, return a first guess for the parameter vector, ak, !! approximately optimal for the domain of the given dimensions. -!! @author R. J. Purser !! @param asp aspect ratio !! @param ak first guess for the parameter vector +!! @author R. J. Purser subroutine guessak_map(asp,tmarcx,ak)! [guessak_map] implicit none real(dp), intent(in ):: asp,tmarcx @@ -684,9 +683,9 @@ end subroutine guessak_map !! (degree) units measured along the rectangle's median, return a first guess !! for the parameter vector, ak, approximately optimal for the domain of the !! given dimensions. -!! @author R. J. Purser !! @param asp aspect ratio !! @param ak first guess or the parameter vector +!! @author R. J. Purser subroutine guessak_geo(asp,arc,ak)! [guessak_geo] implicit none real(dp), intent(in ):: asp,arc @@ -770,12 +769,12 @@ end subroutine guessak_geo !! in both the x and y directions), but restricted to a mere quadrant of the !! domain (since bilateral symmetry pertains across both domain medians, !! yielding a domain mean L that is strictly diagonal. -!! @author R. J. Purser !! @param a Extended Schmidt Gnomonic parameter !! @param k Extended Schmidt Gnomonic parameter !! @param ff failure flag !! @param garcx map-space half-spans !! @param garcy map-space half-spans +!! @author R. J. Purser subroutine bestesg_geo(lam,garcx,garcy, a,k,marcx,marcy,q,ff)! [bestesg_geo] use pietc, only: u5,o5,s18,s36,s54,s72,ms18,ms36,ms54,ms72 use pmat, only: inv @@ -932,13 +931,13 @@ end subroutine bestesg_geo !! in both the x and y directions), but restricted to a mere quadrant of the !! domain (since bilateral symmetry pertains across both domain medians, !! yielding a domain mean L that is strictly diagonal. -!! @author R. J. Purser !! @param a Extended Schmidt Gnomonic parameter !! @param k Extended Schmidt Gnomonic parameter !! @param marcx map-coordinate half-spans !! @param marcy map-coordinate half-spans !! @param garcx geographical half-spans !! @param garcy geographical half-spans +!! @author R. J. Purser subroutine bestesg_map(lam,marcx,marcy, a,k,garcx,garcy,q,ff) ![bestesg_map] use pietc, only: u5,o5,s18,s36,s54,s72,ms18,ms36,ms54,ms72 use psym2, only: chol2 @@ -1068,7 +1067,6 @@ end subroutine bestesg_map !! If all goes well, return a lowered failure flag, ff=.false. . !! But if, for some reason, it is not possible to complete this task, !! return the raised failure flag, ff=.TRUE. . -!! @author R. J. Purser !! @param a parameters of an ESG mapping centered at (plat,plon) !! @param k parameters of an ESG mapping centered at (plat,plon) !! @param plat latitude center point @@ -1081,6 +1079,7 @@ end subroutine bestesg_map !! @param glon grid points for longitude !! @param garea rectangular array !! @param ff failure flag +!! @author R. J. Purser subroutine hgrid_ak_rr(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rr] delx,dely, glat,glon,garea, ff) use pmat4, only: sarea @@ -1173,7 +1172,6 @@ end subroutine hgrid_ak_rr !! 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. -!! @author R. J. Purser !! @param a Extended Schmidt Gnomonic parameter !! @param k Extended Schmidt Gnomonic parameter !! @param plat latitude define centered mapping @@ -1189,6 +1187,7 @@ end subroutine hgrid_ak_rr !! @param angle_dx x angles relative to local east and north !! @param angle_dy y angles relative to local east and north !! @param ff failure flag +!! @author R. J. Purser subroutine hgrid_ak_rr_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak_rr] delx,dely, glat,glon,garea,dx,dy,angle_dx,angle_dy, ff) use pmat4, only: cross_product,triple_product @@ -1374,7 +1373,6 @@ end subroutine hgrid_ak_rr_c !! If all goes well, return a lowered failure flag, ff=.false. . !! But if, for some reason, it is not possible to complete this task, !! return the raised failure flag, ff=.TRUE. . -!! @author R. J. Purser !! @param a parameters of an ESG mapping centered at (plat,plon) !! @param k parameters of an ESG mapping centered at (plat,plon) !! @param plat latitude define centered mapping @@ -1390,6 +1388,7 @@ end subroutine hgrid_ak_rr_c !! @param nx-1 x dimensions for garea !! @param ny-1 y dimensions for garea !! @param ff failure flag +!! @author R. J. Purser subroutine hgrid_ak_rc(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rc] delx,dely, xc,xcd,garea, ff) use pmat4, only: sarea @@ -1461,7 +1460,6 @@ end subroutine hgrid_ak_rc !! the angles are returned in degrees. Garea, the area of each grid cell, is !! returned as in hgrid_ak_rr, and a failure flag, ff, raised when a failure !! occurs anywhere in these calculations. -!! @author R. J. Purser !! @param a parameters of an ESG mapping !! @param k parameters of an ESG mapping !! @param pdlat latitude define centered mapping @@ -1473,6 +1471,7 @@ end subroutine hgrid_ak_rc !! @param delx central x-spacing grid point !! @param dely central y-spacing grid point !! @param ff failure flag +!! @author R. J. Purser subroutine hgrid_ak_dd(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dd] delx,dely, gdlat,gdlon,garea, ff) implicit none @@ -1495,7 +1494,6 @@ end subroutine hgrid_ak_dd !> Like hgrid_ak_rr_c, except all the angle arguments (but not delx,dely) !! are in degrees instead of radians. -!! @author R. J. Purser !! @param a parameters of an ESG mapping !! @param k parameters of an ESG mapping !! @param pdlat latitude define centered mapping @@ -1509,6 +1507,7 @@ end subroutine hgrid_ak_dd !! @param dangle_dx x rotations of the steps !! @param dangle_dy y rotations of the steps !! @param ff failure flag +!! @author R. J. Purser subroutine hgrid_ak_dd_c(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, &! [hgrid_ak_dd] delx,dely, gdlat,gdlon,garea,dx,dy,dangle_dx,dangle_dy, ff) implicit none @@ -1542,7 +1541,6 @@ end subroutine hgrid_ak_dd_c !! matrices, xcd. Garea, the area of each grid cell, is also !! returned as in hgrid_ak_rx, and a failure flag, ff, raised when a failure !! occurs anywhere in these calculations. -!! @author R. J. Purser !! @param a parameters of an ESG mapping !! @param k parameters of an ESG mapping !! @param pdlat latitude define centered mapping @@ -1557,6 +1555,7 @@ end subroutine hgrid_ak_dd_c !! @param xcd Jacobian matrices !! @param garea Jacobian matrices !! @param ff failure flag +!! @author R. J. Purser subroutine hgrid_ak_dc(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dc] delx,dely, xc,xcd,garea, ff) implicit none @@ -1579,7 +1578,6 @@ end subroutine hgrid_ak_dc !! units, delxre, delyre on entry, and to express the grid cell areas, garea, !! in dimensional units upon return. !! The gridded lats and lons, glat and glon, remain in radians. -!! @author R. J. Purser !! @param a parameters of an ESG mapping !! @param k parameters of an ESG mapping !! @param pdlat latitude define centered mapping @@ -1592,6 +1590,7 @@ end subroutine hgrid_ak_dc !! @param delxre map-space grid increments in the dimensional units !! @param delyre map-space grid increments in the dimensional units !! @param ff failure flag +!! @author R. J. Purser subroutine hgrid_ak(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak] re,delxre,delyre, glat,glon,garea, ff) implicit none @@ -1620,7 +1619,6 @@ end subroutine hgrid_ak !! version of this routine, the relative rotations of the steps, dangle_dx !! and dangle_dy, are returned as degrees instead of radians (all other angles !! in the argument list, i.e., plat,plon,pazi,glat,glon, remain radians). -!! @author R. J. Purser !! @param a Extended Schmidt Gnomonic parameter !! @param k Extended Schmidt Gnomonic parameter !! @param plat latitude define centered mapping @@ -1635,6 +1633,7 @@ end subroutine hgrid_ak !! @param glon gridded lons !! @param dangle_dx x rotations of the steps !! @param dangle_dy y rotations of the steps +!! @author R. J. Purser subroutine hgrid_ak_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak] re,delxre,delyre, glat,glon,garea,dx,dy,dangle_dx,dangle_dy, ff) implicit none @@ -1667,6 +1666,7 @@ end subroutine hgrid_ak_c !! @param m number of nodes in half-interval !! @param x nodes and weights !! @param w nodes and weights +!! @author R. J. Purser subroutine gaulegh(m,x,w)! [gaulegh] implicit none integer(spi), intent(IN ):: m ! <- number of nodes in half-interval @@ -1692,7 +1692,6 @@ end subroutine gaulegh !! map-space units) and the sample lat-lon (in radian), return the the !! image in map space in a 2-vector in grid units. If the transformation !! is invalid, return a .true. failure flag. -!! @author R. J. Purser !! @param a parameters of an ESG mapping !! @param k parameters of an ESG mapping !! @param plat latitude define centered mapping @@ -1700,6 +1699,7 @@ end subroutine gaulegh !! @param lat radian latitude !! @param lon radian longitude !! @param ff failure flag +!! @author R. J. Purser subroutine gtoxm_ak_rr_m(A,K,plat,plon,pazi,lat,lon,xm,ff)! [gtoxm_ak_rr] use pmat5, only: grtoc implicit none @@ -1731,7 +1731,6 @@ end subroutine gtoxm_ak_rr_m !! map-space units) and the sample lat-lon (in radian), return the the !! image in map space in a 2-vector in grid units. If the transformation !! is invalid, return a .true. failure flag. -!! @author R. J. Purser !! @param a parameters of an ESG mapping !! @param k parameters of an ESG mapping !! @param plat latitude define centered mapping @@ -1741,6 +1740,7 @@ end subroutine gtoxm_ak_rr_m !! @param lat radian latitude !! @param lon radian longitude !! @param ff failure flag +!! @author R. J. Purser subroutine gtoxm_ak_rr_g(A,K,plat,plon,pazi,delx,dely,lat,lon,&! [gtoxm_ak_rr] xm,ff) implicit none @@ -1752,7 +1752,6 @@ subroutine gtoxm_ak_rr_g(A,K,plat,plon,pazi,delx,dely,lat,lon,&! [gtoxm_ak_rr] end subroutine gtoxm_ak_rr_g !> Like gtoxm_ak_rr_m, except input angles are expressed in degrees -!! @author R. J. Purser !! @param a parameters of an ESG mapping !! @param k parameters of an ESG mapping !! @param pdlat latitude define centered mapping @@ -1760,6 +1759,7 @@ end subroutine gtoxm_ak_rr_g !! @param dlat radian latitude !! @param dlon radian longitude !! @param ff failure flag +!! @author R. J. Purser subroutine gtoxm_ak_dd_m(A,K,pdlat,pdlon,pdazi,dlat,dlon,&! [gtoxm_ak_dd] xm,ff) implicit none @@ -1776,7 +1776,6 @@ subroutine gtoxm_ak_dd_m(A,K,pdlat,pdlon,pdazi,dlat,dlon,&! [gtoxm_ak_dd] end subroutine gtoxm_ak_dd_m !> Like gtoxm_ak_rr_g, except input angles are expressed in degrees -!! @author R. J. Purser !! @param a parameters of an ESG mapping !! @param k parameters of an ESG mapping !! @param pdlat latitude define centered mapping @@ -1784,6 +1783,7 @@ end subroutine gtoxm_ak_dd_m !! @param delx central x-spacing grid point !! @param dely central y-spacing grid point !! @param ff failure flag +!! @author R. J. Purser subroutine gtoxm_ak_dd_g(A,K,pdlat,pdlon,pdazi,delx,dely,&! [gtoxm_ak_dd] dlat,dlon, xm,ff) implicit none @@ -1804,7 +1804,6 @@ end subroutine gtoxm_ak_dd_g !! coordinates given by the 2-vector, xm, return the geographical !! coordinates, lat and lon (radians). If the transformation is invalid for !! any reason, return instead with a raised failure flag, FF= .true. -!! @author R. J. Purser !! @param a parameters of an ESG mapping !! @param k parameters of an ESG mapping !! @param plat latitude of geographical orientation @@ -1812,6 +1811,7 @@ end subroutine gtoxm_ak_dd_g !! @param lat latitude of geographical coordinates !! @param lon longitude of geographical coordinates !! @param ff failure flag +!! @author R. J. Purser subroutine xmtog_ak_rr_m(A,K,plat,plon,pazi,xm,lat,lon,ff)! [xmtog_ak_rr] use pmat5, only: ctogr implicit none @@ -1845,7 +1845,6 @@ end subroutine xmtog_ak_rr_m !! as the 2-vector, xm, return the geographical coordinates, lat, lon, (radians) !! of this point. If instead the transformation is invalid for any reason, then !! return the raised failure flag, FF=.true. -!! @author R. J. Purser !! @param a parameters of an ESG mapping !! @param k parameters of an ESG mapping !! @param plon longitude geographical orientation @@ -1856,6 +1855,7 @@ end subroutine xmtog_ak_rr_m !! @param delx central x-spacing grid point !! @param dely central y-spacing grid point !! @param ff failure flag +!! @author R. J. Purser subroutine xmtog_ak_rr_g(A,K,plat,plon,pazi,delx,dely,xm,&! [xmtog_ak_rr] lat,lon,ff) implicit none @@ -1870,12 +1870,12 @@ subroutine xmtog_ak_rr_g(A,K,plat,plon,pazi,delx,dely,xm,&! [xmtog_ak_rr] end subroutine xmtog_ak_rr_g !> Like xmtog_ak_rr_m, except angles are expressed in degrees -!! @author R. J. Purser !! @param a parameters of an ESG mapping !! @param k parameters of an ESG mapping !! @param plat latitude of geographical orientation !! @param plon longitude of geographical orientation !! @param ff failure flag +!! @author R. J. Purser subroutine xmtog_ak_dd_m(A,K,pdlat,pdlon,pdazi,xm,dlat,dlon,ff)! [xmtog_ak_dd] use pmat5, only: ctogr implicit none @@ -1893,7 +1893,6 @@ subroutine xmtog_ak_dd_m(A,K,pdlat,pdlon,pdazi,xm,dlat,dlon,ff)! [xmtog_ak_dd] end subroutine xmtog_ak_dd_m !> Like xmtog_ak_rr_g, except angles are expressed in degrees -!! @author R. J. Purser !! @param a parameters of an ESG mapping !! @param k parameters of an ESG mapping !! @param pdlat latitude define centered mapping @@ -1904,6 +1903,7 @@ end subroutine xmtog_ak_dd_m !! @param delx central x-spacing grid point !! @param dely central y-spacing grid point !! @param ff failure flag +!! @author R. J. Purser subroutine xmtog_ak_dd_g(A,K,pdlat,pdlon,pdazi,delx,dely,xm,&! [xmtog_ak_dd] dlat,dlon,ff) implicit none From 69870531756e3409dc1c203abf25fee43367688d Mon Sep 17 00:00:00 2001 From: Edward Hartnett Date: Tue, 2 Mar 2021 13:07:05 -0700 Subject: [PATCH 3/5] doxygen fixes --- .../regional_esg_grid.fd/pesg.f90 | 145 +++++++++++------- 1 file changed, 93 insertions(+), 52 deletions(-) diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 index f4be4550f..a779356fe 100644 --- a/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 +++ b/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 @@ -1,18 +1,16 @@ !> @file -!! Suite of routines to perform the 2-parameter family of Extended -!! Schmidt Gnomonic (ESG) regional grid mappings, and to optimize the -!! the two parameters, A and K, of those mappings for a given rectangular -!! domain's principal (median) semi-arcs with respect to a domain-averaged -!! measure of distortion. This criterion is itself endowed with a parameter, -!! lam (for "lambda" in [0,1) ) which gives weight to additional weight -!! areal inhomogeneities instead of treating all distortion components -!! equally. -!! DEPENDENCIES -!! Libraries: pmat, psym2, pfun -!! Modules: pkind, pietc, pietc_s +!! @brief Routines to perform ESG regional grid mappings. !! @author R. J. Purser @date May 2020 -!> Module contain subroutines +!> Suite of routines to perform the 2-parameter family of Extended +!! Schmidt Gnomonic (ESG) regional grid mappings, and to optimize the +!! the two parameters, A and K, of those mappings for a given +!! rectangular domain's principal (median) semi-arcs with respect to a +!! domain-averaged measure of distortion. This criterion is itself +!! endowed with a parameter, lam (for "lambda" in [0,1) ) which gives +!! weight to additional weight areal inhomogeneities instead of +!! treating all distortion components equally. +!! !! @author R. J. Purser module pesg !============================================================================= @@ -61,9 +59,10 @@ module pesg contains -!> Inverse of xstoxc. I.e., cartesians to stereographic -!! @param xc xstoxc -!! @param xs Inverse of xstoxc +!> Inverse of xstoxc. I.e., cartesians to stereographic. +!! +!! @param[in] xc xstoxc +!! @param[out] xs Inverse of xstoxc !! @author R. J. Purser subroutine xctoxs(xc,xs)! [xctoxs] implicit none @@ -73,12 +72,14 @@ subroutine xctoxs(xc,xs)! [xctoxs] zp=u1+xc(3); xs=xc(1:2)/zp end subroutine xctoxs -!> Standard transformation from polar stereographic map coordinates, xs, to -!! cartesian, xc, assuming the projection axis is polar. -!! xcd=d(xc)/d(xs) is the jacobian matrix, encoding distortion and metric. -!! @param xs polar stereographic map coordinates -!! @param xc cartesian -!! @param xcd value of jacobian matrix, encoding distortion and metric +!> Standard transformation from polar stereographic map coordinates, +!! xs, to cartesian, xc, assuming the projection axis is polar. +!! xcd=d(xc)/d(xs) is the jacobian matrix, encoding distortion and +!! metric. +!! +!! @param[in] xs polar stereographic map coordinates +!! @param[out] xc cartesian +!! @param[out] xcd value of jacobian matrix, encoding distortion and metric !! @author R. J. Purser subroutine xstoxc(xs,xc,xcd)! [xstoxc] use pmat4, only: outer_product @@ -92,14 +93,15 @@ subroutine xstoxc(xs,xc,xcd)! [xstoxc] xc(3)=xc(3)-u1 end subroutine xstoxc -!> Standard transformation from polar stereographic map coordinates, xs, to -!! cartesian, xc, assuming the projection axis is polar. -!! xcd=d(xc)/d(xs) is the jacobian matrix, encoding distortion and metric. -!! xcdd is the further derivative, wrt xs, of xcd. -!! @param xs polar stereographic map coordinates -!! @param xc cartesian -!! @param xcd jacobian matrix, encoding distortion and metric -!! @param xcdd further derivative, wrt xs, of xcd +!> Standard transformation from polar stereographic map coordinates, +!! xs, to cartesian, xc, assuming the projection axis is polar. +!! xcd=d(xc)/d(xs) is the jacobian matrix, encoding distortion and +!! metric. xcdd is the further derivative, wrt xs, of xcd. +!! +!! @param[in] xs polar stereographic map coordinates +!! @param[out] xc cartesian +!! @param[out] xcd jacobian matrix, encoding distortion and metric +!! @param[out] xcdd further derivative, wrt xs, of xcd !! @author R. J. Purser subroutine xstoxc1(xs,xc,xcd,xcdd)! [xstoxc] use pmat4, only: outer_product @@ -126,7 +128,12 @@ subroutine xstoxc1(xs,xc,xcd,xcdd)! [xstoxc] do i=1,2; xcd(i,i)=xcd(i,i)+zp; enddo end subroutine xstoxc1 -!> Inverse of xttoxs +!> Inverse of xttoxs. +!! +!! @param[in] k ??? +!! @param[in] xs ??? +!! @param[out] xt ??? +!! @param[out] ff ??? !! @author R. J. Purser subroutine xstoxt(k,xs,xt,ff)! [xstoxt] implicit none @@ -140,9 +147,13 @@ subroutine xstoxt(k,xs,xt,ff)! [xstoxt] xt=u2*xs/sc end subroutine xstoxt -!> Scaled gnomonic plane xt to standard stereographic plane xs -!! @param xt Scaled gnomonic plane -!! @param xs standard stereographic plane +!> Scaled gnomonic plane xt to standard stereographic plane xs. +!! +!! @param[in] k ??? +!! @param[in] xt Scaled gnomonic plane +!! @param[out] xs standard stereographic plane +!! @param[out] xsd ??? +!! @param[out] ff ??? !! @author R. J. Purser subroutine xttoxs(k,xt,xs,xsd,ff)! [xttoxs use pmat4, only: outer_product @@ -166,10 +177,17 @@ subroutine xttoxs(k,xt,xs,xsd,ff)! [xttoxs do i=1,2; xsd(i,i)=xsd(i,i)+rsppi; enddo end subroutine xttoxs -!> Like xttoxs, but also, return the derivatives, wrt K, of xs and xsd -!! @param xt Scaled gnomonic plane -!! @param xs standard stereographic plane -!! @param K derivatives, wrt K, of xs and xsd +!> Like xttoxs, but also, return the derivatives, wrt K, of xs and +!! xsd. +!! +!! @param[in] k derivatives, wrt K, of xs and xsd +!! @param[in] xt Scaled gnomonic plane +!! @param[out] xs standard stereographic plane +!! @param[out] xsd ??? +!! @param[out] xsdd ??? +!! @param[out] xs1 ??? +!! @param[out] xsd1 ??? +!! @param[out] ff ??? !! @author R. J. Purser subroutine xttoxs1(k,xt,xs,xsd,xsdd,xs1,xsd1,ff)! [xttoxs] use pmat4, only: outer_product @@ -206,7 +224,12 @@ subroutine xttoxs1(k,xt,xs,xsd,xsdd,xs1,xsd1,ff)! [xttoxs] do i=1,2; xsdd(i,:,:)=xsdd(i,:,:)-xt(i)*rspdd*rsppi*k/sp; enddo end subroutine xttoxs1 -!> Inverse of xmtoxt +!> Inverse of xmtoxt. +!! +!! @param[in] a ??? +!! @param[in] xt ??? +!! @param[out] xm ??? +!! @param[out] ff ??? !! @author R. J. Purser subroutine xttoxm(a,xt,xm,ff)! [xttoxm] implicit none @@ -218,10 +241,14 @@ subroutine xttoxm(a,xt,xm,ff)! [xttoxm] do i=1,2; call zttozm(a,xt(i),xm(i),ff); if(ff)return; enddo end subroutine xttoxm -!> Like zmtozt, but for 2-vector xm and xt, and 2*2 diagonal Jacobian xtd -!! @param xm vector value -!! @param xt vector value -!! @param xtd 2*2 diagonal Jacobian +!> Like zmtozt, but for 2-vector xm and xt, and 2*2 diagonal Jacobian +!! xtd. +!! +!! @param[in] a ??? +!! @param[in] xm vector value +!! @param[out] xt vector value +!! @param[out] xtd 2*2 diagonal Jacobian +!! @param[out] ff ??? !! @author R. J. Purser subroutine xmtoxt(a,xm,xt,xtd,ff)! [xmtoxt] implicit none @@ -234,12 +261,16 @@ subroutine xmtoxt(a,xm,xt,xtd,ff)! [xmtoxt] xtd=u0; do i=1,2; call zmtozt(a,xm(i),xt(i),xtd(i,i),ff); if(ff)return; enddo end subroutine xmtoxt -!> Like zmtozt1, but for 2-vector xm and xt, and 2*2 diagonal Jacobian xtd -!! Also, the derivatives, wrt a, of these quantities. -!! @param xm vector value -!! @param xt vector value -!! @param xtd 2*2 diagonal Jacobian -!! @param a the derivatives of these quantities +!> Like zmtozt1, but for 2-vector xm and xt, and 2*2 diagonal Jacobian +!! xtd Also, the derivatives, wrt a, of these quantities. +!! +!! @param[in] a the derivatives of these quantities +!! @param[in] xm vector value +!! @param[out] xt vector value +!! @param[out] xtd 2*2 diagonal Jacobian +!! @param[out] xt1 ??? +!! @param[out] xtd1 ??? +!! @param[out] ff ??? !! @author R. J. Purser subroutine xmtoxt1(a,xm,xt,xtd,xt1,xtd1,ff)! [xmtoxt] implicit none @@ -258,6 +289,11 @@ subroutine xmtoxt1(a,xm,xt,xtd,xt1,xtd1,ff)! [xmtoxt] end subroutine xmtoxt1 !> Inverse of zmtozt +!! +!! @param[in] a ??? +!! @param[in] zt ??? +!! @param[out] zm ??? +!! @param[out] ff ??? !! @author R. J. Purser subroutine zttozm(a,zt,zm,ff)! [zttozm] implicit none @@ -273,10 +309,15 @@ subroutine zttozm(a,zt,zm,ff)! [zttozm] endif end subroutine zttozm -!> Evaluate the function, zt = tan(sqrt(A)*z)/sqrt(A), and its derivative, ztd, -!! for positive and negative A and for the limiting case, A --> 0 -!! @param zt function to be evaluated -!! @param ztd derivative of the source function +!> Evaluate the function, zt = tan(sqrt(A)*z)/sqrt(A), and its +!! derivative, ztd, for positive and negative A and for the limiting +!! case, A --> 0. +!! +!! @param[in] a ??? +!! @param[in] zm ??? +!! @param[out] zt function to be evaluated +!! @param[out] ztd derivative of the source function +!! @param[out] ff ??? !! @author R. J. Purser subroutine zmtozt(a,zm,zt,ztd,ff)! [zmtozt] implicit none From 453a403ce2fbff79086411e7ab73e694815c1b24 Mon Sep 17 00:00:00 2001 From: Edward Hartnett Date: Wed, 3 Mar 2021 05:45:49 -0700 Subject: [PATCH 4/5] more doxygen work --- docs/Doxyfile.in | 4 ++-- .../regional_esg_grid.fd/regional_esg_grid.f90 | 9 +++++++-- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/docs/Doxyfile.in b/docs/Doxyfile.in index d596f9b16..bed0df074 100644 --- a/docs/Doxyfile.in +++ b/docs/Doxyfile.in @@ -913,13 +913,13 @@ STRIP_CODE_COMMENTS = YES # function all documented functions referencing it will be listed. # The default value is: NO. -REFERENCED_BY_RELATION = NO +REFERENCED_BY_RELATION = YES # If the REFERENCES_RELATION tag is set to YES then for each documented function # all documented entities called/used by that function will be listed. # The default value is: NO. -REFERENCES_RELATION = NO +REFERENCES_RELATION = YES # If the REFERENCES_LINK_SOURCE tag is set to YES and SOURCE_BROWSER tag is set # to YES, then the hyperlinks from functions in REFERENCES_RELATION and diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/regional_esg_grid.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/regional_esg_grid.f90 index 22a1f2e2b..b66fe5c1b 100644 --- a/sorc/grid_tools.fd/regional_esg_grid.fd/regional_esg_grid.f90 +++ b/sorc/grid_tools.fd/regional_esg_grid.fd/regional_esg_grid.f90 @@ -1,7 +1,12 @@ !> @file -!============================================================================= +!! @brief ??? +!! @author R. J. Purser + +!! Program which does ??? +!! +!! @author R. J. Purser +!! @return 0 for success, error code otherwise. program regional_grid -!============================================================================= use pkind, only: dp use pietc, only: dtor,rtod From 2622af8689295ca9225f1ab11952c4125f964328 Mon Sep 17 00:00:00 2001 From: Edward Hartnett Date: Wed, 3 Mar 2021 06:28:05 -0700 Subject: [PATCH 5/5] more doxygen --- sorc/grid_tools.fd/CMakeLists.txt | 3 + .../regional_esg_grid.fd/CMakeLists.txt | 4 + .../regional_esg_grid.fd/pesg.f90 | 907 +++++++++++------- 3 files changed, 543 insertions(+), 371 deletions(-) diff --git a/sorc/grid_tools.fd/CMakeLists.txt b/sorc/grid_tools.fd/CMakeLists.txt index a4e66910c..8befb8124 100644 --- a/sorc/grid_tools.fd/CMakeLists.txt +++ b/sorc/grid_tools.fd/CMakeLists.txt @@ -1,3 +1,6 @@ +# This is the CMake file for the grid_tools utilites of UFS_UTILS. +# +# George Gayno add_subdirectory(shave.fd) add_subdirectory(filter_topo.fd) add_subdirectory(regional_esg_grid.fd) diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/CMakeLists.txt b/sorc/grid_tools.fd/regional_esg_grid.fd/CMakeLists.txt index 2e3e93295..592290472 100644 --- a/sorc/grid_tools.fd/regional_esg_grid.fd/CMakeLists.txt +++ b/sorc/grid_tools.fd/regional_esg_grid.fd/CMakeLists.txt @@ -1,3 +1,7 @@ +# This is the CMake file for the regional_esf_grid utility in +# UFS_UTILS. +# +# George Gayno set(fortran_src pesg.f90 pfun.f90 diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 index a779356fe..3c5c1e4fb 100644 --- a/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 +++ b/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 @@ -333,8 +333,16 @@ subroutine zmtozt(a,zm,zt,ztd,ff)! [zmtozt] ztd=u1+a*zt*zt end subroutine zmtozt -!> Like zmtozt, but -!! also, get the derivative with respect to a, zt1 of zt, and ztd1 of ztd. +!> Like zmtozt, but also, get the derivative with respect to a, zt1 of +!! zt, and ztd1 of ztd. +!! +!! @param[in] a ??? +!! @param[in] zm ??? +!! @param[in] zt ??? +!! @param[in] ztd ??? +!! @param[in] zt1 ??? +!! @param[in] ztd1 ??? +!! @param[in] ff ??? !! @author R. J. Purser subroutine zmtozt1(a,zm,zt,ztd,zt1,ztd1,ff)! [zmtozt] use pietc, only: o3 @@ -365,6 +373,8 @@ end subroutine zmtozt1 !! FF. !! @param ak parameterization of the Extended ESG !! @param xm vector +!! @param xc +!! @param xcd !! @param ff raised failure error code !! @author R. J. Purser subroutine xmtoxc_vak(ak,xm,xc,xcd,ff)! [xmtoxc_ak] @@ -377,7 +387,14 @@ subroutine xmtoxc_vak(ak,xm,xc,xcd,ff)! [xmtoxc_ak] end subroutine xmtoxc_vak !> Like xmtoxc_vak, _ak, but also return derivatives wrt ak. +!! !! @param ak derivatives wrt ak +!! @param[in] xm +!! @param[out] xc +!! @param[out] xcd +!! @param[out] xc1 +!! @param[out] xcd1 +!! @param[out] ff !! @author R. J. Purser subroutine xmtoxc_vak1(ak,xm,xc,xcd,xc1,xcd1,ff)! [xmtoxc_ak] implicit none @@ -405,11 +422,12 @@ subroutine xmtoxc_vak1(ak,xm,xc,xcd,xc1,xcd1,ff)! [xmtoxc_ak] xcd=matmul(xcd,xsd) end subroutine xmtoxc_vak1 -!> Assuming the A-K parameterization of the Extended Schmidt-transformed -!! Gnomonic (ESG) mapping, and given a map-space 2-vector, xm, find the -!! corresponding cartesian unit 3-vector and its derivative wrt xm, jacobian -!! matrix, xcd. If for any reason the mapping cannot be done, return a -!! raised failure flag, FF. +!> Assuming the A-K parameterization of the Extended +!! Schmidt-transformed Gnomonic (ESG) mapping, and given a map-space +!! 2-vector, xm, find the corresponding cartesian unit 3-vector and +!! its derivative wrt xm, jacobian matrix, xcd. If for any reason the +!! mapping cannot be done, return a raised failure flag, FF. +!! !! @param a ESG mapping parameterization !! @param k ESG mapping parameterization !! @param xm map-space vector @@ -433,14 +451,15 @@ subroutine xmtoxc_ak(a,k,xm,xc,xcd,ff)! [xmtoxc_ak] xcd=matmul(xcd,xsd) end subroutine xmtoxc_ak -!> Inverse mapping of xmtoxc_ak. That is, go from given cartesian unit -!! 3-vector, xc, to map coordinate 2-vector xm (or return a raised failure -!! flag, FF, if the attempt fails). -!! @param a Inverse mapping of xmtoxc -!! @param k Inverse mapping of xmtoxc -!! @param xc cartesian unit 3-vector -!! @param xm 2-vector map coordinate -!! @param ff error flag +!! Inverse mapping of xmtoxc_ak. That is, go from given cartesian unit +!! 3-vector, xc, to map coordinate 2-vector xm (or return a raised +!! failure flag, FF, if the attempt fails). +!! +!! @param[in] a Inverse mapping of xmtoxc +!! @param[in] k Inverse mapping of xmtoxc +!! @param[in] xc cartesian unit 3-vector +!! @param[out] xm 2-vector map coordinate +!! @param[out] ff error flag !! @author R. J. Purser subroutine xctoxm_ak(a,k,xc,xm,ff)! [xctoxm_ak] implicit none @@ -455,12 +474,15 @@ subroutine xctoxm_ak(a,k,xc,xm,ff)! [xctoxm_ak] call xttoxm(a,xt,xm,ff) end subroutine xctoxm_ak -!> For angles (degrees) of the arcs spanning the halfwidths between the -!! region's center and its x and y edges, get the two cartesian vectors -!! that represent the locations of these edge midpoints in the positive x and y -!! directions. -!! @param edgex region's x edges -!! @param edgey region's y edges +!> For angles (degrees) of the arcs spanning the halfwidths between +!! the region's center and its x and y edges, get the two cartesian +!! vectors that represent the locations of these edge midpoints in the +!! positive x and y directions. +!! +!! @param[in] arcx +!! @param[in] arcy +!! @param[out] edgex region's x edges +!! @param[out] edgey region's y edges !! @author R. J. Purser subroutine get_edges(arcx,arcy,edgex,edgey)! [get_edges] implicit none @@ -473,16 +495,17 @@ subroutine get_edges(arcx,arcy,edgex,edgey)! [get_edges] edgex=(/sx,u0,cx/); edgey=(/u0,sy,cy/) end subroutine get_edges -!> From a jacobian matrix, j0, get a sufficient set of v.. diagnostics such -!! that, from averages of these v, we can later compute the collective variance -!! of Q(lam) that they imply for any choice of the "lambda" parameter, lam. -!! Note that v1 and v4 are quadratic diagnostics of EL, while v2 and v3 are -!! linear. -!! @param j0 jacobian matrix -!! @param v1 quadratic diagnostics of EL -!! @param v4 quadratic diagnostics of EL -!! @param v2 linear diagnostics of EL -!! @param v3 linear diagnostics of EL +!> From a jacobian matrix, j0, get a sufficient set of v.. diagnostics +!! such that, from averages of these v, we can later compute the +!! collective variance of Q(lam) that they imply for any choice of the +!! "lambda" parameter, lam. Note that v1 and v4 are quadratic +!! diagnostics of EL, while v2 and v3 are linear. +!! +!! @param[in] j0 jacobian matrix +!! @param[out] v1 quadratic diagnostics of EL +!! @param[out] v2 linear diagnostics of EL +!! @param[out] v3 linear diagnostics of EL +!! @param[out] v4 quadratic diagnostics of EL !! @author R. J. Purser subroutine get_qx(j0, v1,v2,v3,v4)! [get_qx] use psym2, only: logsym2 @@ -498,11 +521,21 @@ subroutine get_qx(j0, v1,v2,v3,v4)! [get_qx] v4=(el(1,1)+el(2,2))**2 end subroutine get_qx -!> From a jacobian matrix, j0, and its derivative, j0d, get a sufficient set -!! of v.. diagnostics such that, from average of these diagnostics, we can -!! later compute the collective variance of Q and its derivative. -!! @param j0 jacobian matrix -!! @param j0d derivative of j0 +!> From a jacobian matrix, j0, and its derivative, j0d, get a +!! sufficient set of v.. diagnostics such that, from average of these +!! diagnostics, we can later compute the collective variance of Q and +!! its derivative. +!! +!! @param[in] j0 jacobian matrix +!! @param[in] j0d derivative of j0 +!! @param[in] v1 +!! @param[in] v2 +!! @param[in] v3 +!! @param[in] v4 +!! @param[in] v1d +!! @param[in] v2d +!! @param[in] v3d +!! @param[in] v4d !! @author R. J. Purser subroutine get_qxd(j0,j0d, v1,v2,v3,v4,v1d,v2d,v3d,v4d)! [get_qx] use psym2, only: logsym2 @@ -534,28 +567,36 @@ subroutine get_qxd(j0,j0d, v1,v2,v3,v4,v1d,v2d,v3d,v4d)! [get_qx] v4d=u2*(el(1,1)+el(2,2))*(eld(1,1,:)+eld(2,2,:)) end subroutine get_qxd -!> For a parameter vector, ak and a map-space domain-parameter vector, ma, -!! return the lambda-parameterized quality diagnostic, Q, and the geographic -!! domain-parameter vector ga. Lambda is given by lam <1. Also, return -!! the derivatives, qdak and qdma, of Q wrt ak and ma, and the derivatives -!! gadak and gadma, of ga wrt ak and ma. -!! The domain averages of Q are accurately computed by bi-Gauss-Legendre -!! quadrature over the positive quadrant of the domain (exploiting the -!! symmetry) of the four constituent terms, v1, v2, v3, v4, from which -!! the mean Q is computed using a quadratic formula of these constituents. -!! The number of Gauss points in eaxh half-interval is ngh, and the -!! nodes themselves are, in proportion to the half-interval, at xg. -!! the normalized gauss weights are wg. +!> For a parameter vector, ak and a map-space domain-parameter vector, +!! ma, return the lambda-parameterized quality diagnostic, Q, and the +!! geographic domain-parameter vector ga. Lambda is given by lam +!! <1. Also, return the derivatives, qdak and qdma, of Q wrt ak and +!! ma, and the derivatives gadak and gadma, of ga wrt ak and ma. +!! +!! The domain averages of Q are accurately computed by +!! bi-Gauss-Legendre quadrature over the positive quadrant of the +!! domain (exploiting the symmetry) of the four constituent terms, v1, +!! v2, v3, v4, from which the mean Q is computed using a quadratic +!! formula of these constituents. The number of Gauss points in eaxh +!! half-interval is ngh, and the nodes themselves are, in proportion +!! to the half-interval, at xg. the normalized gauss weights are wg. +!! !! If a failure occurs, colmputations cease immediately and a failure !! flag, FF, is raised on return. -!! @param ak parameter vector -!! @param ma map-space domain-parameter vector -!! @param q lambda-parameterized quality diagnostic -!! @param ga geographic domain-parameter vector -!! @param lam Lambda -!! @param qdak derivatives value -!! @param qdma derivatives value -!! @param ff error flag +!! +!! @param[in] ngh +!! @param[in] lam Lambda +!! @param[in] xg +!! @param[in] wg +!! @param[in] ak parameter vector +!! @param[in] ma map-space domain-parameter vector +!! @param[out] q lambda-parameterized quality diagnostic +!! @param[out] qdak derivatives value +!! @param[out] qdma derivatives value +!! @param[out] ga geographic domain-parameter vector +!! @param[out] gadak +!! @param[out] gadma +!! @param[out] ff error flag !! @author R. J. Purser subroutine get_meanqd(ngh,lam,xg,wg,ak,ma, q,qdak,qdma, & ! [get_meanq] ga,gadak,gadma, ff) @@ -623,8 +664,18 @@ subroutine get_meanqd(ngh,lam,xg,wg,ak,ma, q,qdak,qdma, & ! [get_meanq] call get_qofv(lam,v2,v3, v1d,v2d,v3d,v4d, qdma) end subroutine get_meanqd -!> Like getmeanqd, except for n different values, aks, of ak and n different -!! values, mas of ma, and without any of the derivatives. +!> Like getmeanqd, except for n different values, aks, of ak and n +!! different values, mas of ma, and without any of the derivatives. +!! +!! @param[in] n +!! @param[in] ngh +!! @param[in] lam +!! @param[in] xg +!! @param[in] wg +!! @param[in] aks +!! @param[in] mas +!! @param[out] qs +!! @param[out] ff !! @author R. J. Purser subroutine get_meanqs(n,ngh,lam,xg,wg,aks,mas, qs,ff)! [get_meanq] implicit none @@ -661,14 +712,17 @@ subroutine get_meanqs(n,ngh,lam,xg,wg,aks,mas, qs,ff)! [get_meanq] call get_qofv(n,lam,v1s,v2s,v3s,v4s, qs) end subroutine get_meanqs -!> The quadratic quantity Q depends linearly on v1 and v4 (which are already -!! quadratic diagnostics of EL) and quadratically on v2 and v3 (which are -!! linear diagnostics of EL). EL = (1/2)log(G), where G=J^T.J, J the jacobian. -!! @param q quadratic quantity -!! @param v1 quadratic diagnostics of EL -!! @param v4 quadratic diagnostics of EL -!! @param v2 linear diagnostics of EL -!! @param v3 linear diagnostics of EL +!> The quadratic quantity Q depends linearly on v1 and v4 (which are +!! already quadratic diagnostics of EL) and quadratically on v2 and v3 +!! (which are linear diagnostics of EL). EL = (1/2)log(G), where +!! G=J^T.J, J the jacobian. +!! +!! @param[in] lam +!! @param[in] v1 quadratic diagnostics of EL +!! @param[in] v2 linear diagnostics of EL +!! @param[in] v3 linear diagnostics of EL +!! @param[in] v4 quadratic diagnostics of EL +!! @param[out] q quadratic quantity !! @author R. J. Purser subroutine get_qofv(lam,v1,v2,v3,v4, q)! [get_qofv] implicit none @@ -679,8 +733,18 @@ subroutine get_qofv(lam,v1,v2,v3,v4, q)! [get_qofv] q=lamc*(v1-(v2**2+v3**2)) +lam*(v4 -(v2+v3)**2) end subroutine get_qofv -!> Like get_qofv, but for (only) the 2-vector derivatives of Q. Note that -!! the quadratic diagnostics v1 and v4 do not participate in this formula. +!> Like get_qofv, but for (only) the 2-vector derivatives of Q. Note +!! that the quadratic diagnostics v1 and v4 do not participate in this +!! formula. +!! +!! @param[in] lam +!! @param[in] v2 +!! @param[in] v3 +!! @param[in] v1d +!! @param[in] v2d +!! @param[in] v3d +!! @param[in] v4d +!! @param[out] qd !! @author R. J. Purser subroutine get_qofvd(lam, v2,v3, v1d,v2d,v3d,v4d, qd)! [get_qofv] implicit none @@ -692,7 +756,15 @@ subroutine get_qofvd(lam, v2,v3, v1d,v2d,v3d,v4d, qd)! [get_qofv] qd=lamc*(v1d-u2*(v2d*v2+v3d*v3))+lam*(v4d-u2*(v2d+v3d)*(v2+v3)) end subroutine get_qofvd -!> General util to convert value +!> General util to convert value. +!! +!! @param[in] n +!! @param[in] lam +!! @param[in] v1s +!! @param[in] v2s +!! @param[in] v3s +!! @param[in] v4s +!! @param[out] qs !! @author R. J. Purser subroutine get_qsofvs(n,lam,v1s,v2s,v3s,v4s, qs)! [get_qofv] implicit none @@ -705,11 +777,14 @@ subroutine get_qsofvs(n,lam,v1s,v2s,v3s,v4s, qs)! [get_qofv] qs=lamc*(v1s-(v2s**2+v3s**2)) +lam*(v4s -(v2s+v3s)**2) end subroutine get_qsofvs -!> Given an aspect ratio, asp<=1, and major semi-axis, arc, in map-space -!! nondimensional units, return a first guess for the parameter vector, ak, -!! approximately optimal for the domain of the given dimensions. -!! @param asp aspect ratio -!! @param ak first guess for the parameter vector +!> Given an aspect ratio, asp<=1, and major semi-axis, arc, in +!! map-space nondimensional units, return a first guess for the +!! parameter vector, ak, approximately optimal for the domain of the +!! given dimensions. +!! +!! @param[in] asp aspect ratio +!! @param[in] tmarcx +!! @param[out] ak first guess for the parameter vector !! @author R. J. Purser subroutine guessak_map(asp,tmarcx,ak)! [guessak_map] implicit none @@ -720,11 +795,13 @@ subroutine guessak_map(asp,tmarcx,ak)! [guessak_map] call guessak_geo(asp,gmarcx,ak) end subroutine guessak_map -!> Given an aspect ratio, asp<=1, and major semi-axis, arc, in geographical -!! (degree) units measured along the rectangle's median, return a first guess -!! for the parameter vector, ak, approximately optimal for the domain of the -!! given dimensions. +!> Given an aspect ratio, asp<=1, and major semi-axis, arc, in +!! geographical (degree) units measured along the rectangle's median, +!! return a first guess for the parameter vector, ak, approximately +!! optimal for the domain of the given dimensions. +!! !! @param asp aspect ratio +!! @param arc ??? !! @param ak first guess or the parameter vector !! @author R. J. Purser subroutine guessak_geo(asp,arc,ak)! [guessak_geo] @@ -810,11 +887,16 @@ end subroutine guessak_geo !! in both the x and y directions), but restricted to a mere quadrant of the !! domain (since bilateral symmetry pertains across both domain medians, !! yielding a domain mean L that is strictly diagonal. -!! @param a Extended Schmidt Gnomonic parameter -!! @param k Extended Schmidt Gnomonic parameter -!! @param ff failure flag -!! @param garcx map-space half-spans -!! @param garcy map-space half-spans +!! +!! @param[in] lam +!! @param[in] garcx map-space half-spans +!! @param[in] garcy map-space half-spans +!! @param[out] a Extended Schmidt Gnomonic parameter +!! @param[out] k Extended Schmidt Gnomonic parameter +!! @param[out] marcx +!! @param[out] marcy +!! @param[out] q +!! @param[out] ff failure flag !! @author R. J. Purser subroutine bestesg_geo(lam,garcx,garcy, a,k,marcx,marcy,q,ff)! [bestesg_geo] use pietc, only: u5,o5,s18,s36,s54,s72,ms18,ms36,ms54,ms72 @@ -950,34 +1032,41 @@ subroutine bestesg_geo(lam,garcx,garcy, a,k,marcx,marcy,q,ff)! [bestesg_geo] endif end subroutine bestesg_geo -!> Get the best Extended Schmidt Gnomonic parameter, (a,k), for the given -!! map-coordinate half-spans, marcx and marcy, as well as the corresponding -!! geographical half-spans, garcx and garcy (in degrees) and the quality -!! diagnostic, Q(lam) for this optimal parameter choice. If this process -!! fails for any reason, the failure is alerted by a raised flag, FF, and -!! the other output arguments must then be taken to be meaningless. +!> Get the best Extended Schmidt Gnomonic parameter, (a,k), for the +!! given map-coordinate half-spans, marcx and marcy, as well as the +!! corresponding geographical half-spans, garcx and garcy (in degrees) +!! and the quality diagnostic, Q(lam) for this optimal parameter +!! choice. If this process fails for any reason, the failure is +!! alerted by a raised flag, FF, and the other output arguments must +!! then be taken to be meaningless. !! -!! The diagnostic Q measures the variance over the domain of a local measure -!! of grid distortion. A logarithmic measure of local grid deformation is -!! give by L=log(J^t.J)/2, where J is the mapping Jacobian matrix, dX/dx, -!! where X is the cartesian unit 3-vector representation of the image of the -!! mapping of the map-coordinate 2-vector, x. -!! The Frobenius squared-norm, Trace(L*L), of L is the basis for the simplest -!! (lam=0) definition of the variance of L, but (Trace(L))**2 is another. -!! Here, we weight both contributions, by lam and (1-lam) respectively, with -!! 0 <= lam <1, to compute the variance Q(lam,a,k), and search for the (a,k) -!! that minimizes this Q. +!! The diagnostic Q measures the variance over the domain of a local +!! measure of grid distortion. A logarithmic measure of local grid +!! deformation is give by L=log(J^t.J)/2, where J is the mapping +!! Jacobian matrix, dX/dx, where X is the cartesian unit 3-vector +!! representation of the image of the mapping of the map-coordinate +!! 2-vector, x. The Frobenius squared-norm, Trace(L*L), of L is the +!! basis for the simplest (lam=0) definition of the variance of L, but +!! (Trace(L))**2 is another. Here, we weight both contributions, by +!! lam and (1-lam) respectively, with 0 <= lam <1, to compute the +!! variance Q(lam,a,k), and search for the (a,k) that minimizes this +!! Q. !! -!! The domain averages are computed by double Gauss-Legendre quadrature (i.e., -!! in both the x and y directions), but restricted to a mere quadrant of the -!! domain (since bilateral symmetry pertains across both domain medians, -!! yielding a domain mean L that is strictly diagonal. -!! @param a Extended Schmidt Gnomonic parameter -!! @param k Extended Schmidt Gnomonic parameter -!! @param marcx map-coordinate half-spans -!! @param marcy map-coordinate half-spans -!! @param garcx geographical half-spans -!! @param garcy geographical half-spans +!! The domain averages are computed by double Gauss-Legendre +!! quadrature (i.e., in both the x and y directions), but restricted +!! to a mere quadrant of the domain (since bilateral symmetry pertains +!! across both domain medians, yielding a domain mean L that is +!! strictly diagonal. +!! +!! @param[in] lam +!! @param[in] marcx map-coordinate half-spans +!! @param[in] marcy map-coordinate half-spans +!! @param[out] a Extended Schmidt Gnomonic parameter +!! @param[out] k Extended Schmidt Gnomonic parameter +!! @param[out] garcx geographical half-spans +!! @param[out] garcy geographical half-spans +!! @param[out] q +!! @param[out] ff failure flag !! @author R. J. Purser subroutine bestesg_map(lam,marcx,marcy, a,k,garcx,garcy,q,ff) ![bestesg_map] use pietc, only: u5,o5,s18,s36,s54,s72,ms18,ms36,ms54,ms72 @@ -1090,36 +1179,41 @@ subroutine bestesg_map(lam,marcx,marcy, a,k,garcx,garcy,q,ff) ![bestesg_map] end subroutine bestesg_map !> Use a and k as the parameters of an Extended Schmidt-transformed -!! Gnomonic (ESG) mapping centered at (plat,plon) and twisted about this center -!! by an azimuth angle of pazi counterclockwise (these angles in radians). +!! Gnomonic (ESG) mapping centered at (plat,plon) and twisted about +!! this center by an azimuth angle of pazi counterclockwise (these +!! angles in radians). !! -!! Assume the radius of the earth is unity, and using the central mapping -!! point as the coordinate origin, set up the grid with central x-spacing delx -!! and y-spacing dely. The grid index location of the left-lower -!! corner of the domain is (lx,ly) (typically both NEGATIVE). -!! The numbers of the grid spaces in x and y directions are nx and ny. -!! (Note that, for a centered rectangular grid lx and ly are negative and, in -!! magnitude, half the values of nx and ny respectively.) -!! Return the latitude and longitude, in radians again, of the grid points thus -!! 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 +!! Assume the radius of the earth is unity, and using the central +!! mapping point as the coordinate origin, set up the grid with +!! central x-spacing delx and y-spacing dely. The grid index location +!! of the left-lower corner of the domain is (lx,ly) (typically both +!! NEGATIVE). The numbers of the grid spaces in x and y directions are +!! nx and ny. (Note that, for a centered rectangular grid lx and ly +!! are negative and, in magnitude, half the values of nx and ny +!! respectively.) Return the latitude and longitude, in radians +!! again, of the grid points thus 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 !! -!! If all goes well, return a lowered failure flag, ff=.false. . -!! But if, for some reason, it is not possible to complete this task, +!! If all goes well, return a lowered failure flag, ff=.false. . But +!! if, for some reason, it is not possible to complete this task, !! return the raised failure flag, ff=.TRUE. . -!! @param a parameters of an ESG mapping centered at (plat,plon) -!! @param k parameters of an ESG mapping centered at (plat,plon) -!! @param plat latitude center point -!! @param plon longitude center point -!! @param lx grid index location -!! @param ly grid index location -!! @param nx grid spaces -!! @param ny grid spaces -!! @param glat grid points for latitude -!! @param glon grid points for longitude -!! @param garea rectangular array -!! @param ff failure flag +!! +!! @param[in] lx grid index location +!! @param[in] ly grid index location +!! @param[in] nx grid spaces +!! @param[in] ny grid spaces +!! @param[in] A parameters of an ESG mapping centered at (plat,plon) +!! @param[in] K parameters of an ESG mapping centered at (plat,plon) +!! @param[in] plat latitude center point +!! @param[in] plon longitude center point +!! @param[in] pazi +!! @param[in] delx central x-spacing grid point +!! @param[in] dely central y-spacing grid point +!! @param[out] glat grid points for latitude +!! @param[out] glon grid points for longitude +!! @param[out] garea rectangular array +!! @param[out] ff failure flag !! @author R. J. Purser subroutine hgrid_ak_rr(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rr] delx,dely, glat,glon,garea, ff) @@ -1188,46 +1282,55 @@ subroutine hgrid_ak_rr(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rr] end subroutine hgrid_ak_rr !> Use a and k as the parameters of an extended Schmidt-transformed -!! gnomonic (ESG) mapping centered at (plat,plon) and twisted about this center -!! by an azimuth angle of pazi counterclockwise (these angles in radians). +!! gnomonic (ESG) mapping centered at (plat,plon) and twisted about +!! this center by an azimuth angle of pazi counterclockwise (these +!! angles in radians). !! -!! Using the central mapping point as the coordinate origin, set up the grid -!! with central x-spacing delx and y-spacing dely in nondimensional units, -!! (i.e., as if the earth had unit radius) and with the location of the left- -!! lower corner of the grid at center-relative grid index pair, (lx,ly) and -!! with the number of the grid spaces in x and y directions given by nx and ny. -!! (Note that, for a centered rectangular grid lx and ly are negative and, in -!! magnitude, half the values of nx and ny respectively.) -!! Return the latitude and longitude, again, in radians, of the grid pts thus -!! defined in the arrays, glat and glon; return a rectangular array, garea, -!! of dimensions nx-1 by ny-1, that contains the areas of each of the grid -!! cells in nondimensional "steradian" units. +!! Using the central mapping point as the coordinate origin, set up +!! the grid with central x-spacing delx and y-spacing dely in +!! nondimensional units, (i.e., as if the earth had unit radius) and +!! with the location of the left- lower corner of the grid at +!! center-relative grid index pair, (lx,ly) and with the number of the +!! grid spaces in x and y directions given by nx and ny. (Note that, +!! for a centered rectangular grid lx and ly are negative and, in +!! magnitude, half the values of nx and ny respectively.) Return the +!! latitude and longitude, again, in radians, of the grid pts thus +!! defined in the arrays, glat and glon; return a rectangular array, +!! garea, of dimensions nx-1 by ny-1, that contains the areas of each +!! of the grid cells in nondimensional "steradian" units. !! -!! 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 and 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 radians counterclockwise. +!! 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 and 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 radians +!! 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. -!! @param a Extended Schmidt Gnomonic parameter -!! @param k Extended Schmidt Gnomonic parameter -!! @param plat latitude define centered mapping -!! @param plon longitude define centered mapping -!! @param delx central x-spacing grid point -!! @param dely central y-spacing grid point -!! @param lx x grid index for left-lower corner of the grid at center -!! @param ly y grid index for left-lower corner of the grid at center -!! @param nx numbers of the grid spaces in x -!! @param ny numbers of the grid spaces in y -!! @param dx estimated grid steps x grid cell edges -!! @param dy estimated grid steps y grid cell edges -!! @param angle_dx x angles relative to local east and north -!! @param angle_dy y angles relative to local east and north -!! @param ff failure flag +!! reason, it is not possible to complete this task, return the +!! failure flag as .TRUE. +!! +!! @param[in] lx x grid index for left-lower corner of the grid at center +!! @param[in] ly y grid index for left-lower corner of the grid at center +!! @param[in] nx numbers of the grid spaces in x +!! @param[in] ny numbers of the grid spaces in y +!! @param[in] a Extended Schmidt Gnomonic parameter +!! @param[in] k Extended Schmidt Gnomonic parameter +!! @param[in] plat latitude define centered mapping +!! @param[in] plon longitude define centered mapping +!! @param[in] pazi +!! @param[in] delx central x-spacing grid point +!! @param[in] dely central y-spacing grid point +!! @param[out] glat grid points for latitude +!! @param[out] glon grid points for longitude +!! @param[out] garea rectangular array +!! @param[out] dx estimated grid steps x grid cell edges +!! @param[out] dy estimated grid steps y grid cell edges +!! @param[out] angle_dx x angles relative to local east and north +!! @param[out] angle_dy y angles relative to local east and north +!! @param[out] ff failure flag !! @author R. J. Purser subroutine hgrid_ak_rr_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak_rr] delx,dely, glat,glon,garea,dx,dy,angle_dx,angle_dy, ff) @@ -1396,38 +1499,40 @@ subroutine hgrid_ak_rr_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak_rr] end subroutine hgrid_ak_rr_c !> Use a and k as the parameters of an Extended Schmidt-transformed -!! Gnomonic (ESG) mapping centered at (plat,plon) and twisted about this center -!! by an azimuth angle of pazi counterclockwise (these angles in radians). +!! Gnomonic (ESG) mapping centered at (plat,plon) and twisted about +!! this center by an azimuth angle of pazi counterclockwise (these +!! angles in radians). !! -!! Assume the radius of the earth is unity, and using the central mapping -!! point as the coordinate origin, set up the grid with central x-spacing delx -!! and y-spacing dely. The grid index location of the left-lower -!! corner of the domain is (lx,ly) (typically both NEGATIVE). -!! The numbers of the grid spaces in x and y directions are nx and ny. -!! (Note that, for a centered rectangular grid lx and ly are negative and, in -!! magnitude, half the values of nx and ny respectively.) -!! Return the unit cartesian vectors xc of the grid points and their jacobian -!! matrices xcd wrt the map coordinates, and return a rectangular array, garea, -!! of dimensions nx-1 by ny-1, that contains the areas of each of the grid -!! cells +!! Assume the radius of the earth is unity, and using the central +!! mapping point as the coordinate origin, set up the grid with +!! central x-spacing delx and y-spacing dely. The grid index location +!! of the left-lower corner of the domain is (lx,ly) (typically both +!! NEGATIVE). The numbers of the grid spaces in x and y directions +!! are nx and ny. (Note that, for a centered rectangular grid lx and +!! ly are negative and, in magnitude, half the values of nx and ny +!! respectively.) Return the unit cartesian vectors xc of the grid +!! points and their jacobian matrices xcd wrt the map coordinates, and +!! return a rectangular array, garea, of dimensions nx-1 by ny-1, that +!! contains the areas of each of the grid cells !! -!! If all goes well, return a lowered failure flag, ff=.false. . -!! But if, for some reason, it is not possible to complete this task, +!! If all goes well, return a lowered failure flag, ff=.false. . But +!! if, for some reason, it is not possible to complete this task, !! return the raised failure flag, ff=.TRUE. . +!! +!! @param lx x grid index for left-lower corner of the grid at center +!! @param ly y grid index for left-lower corner of the grid at center +!! @param nx numbers of the grid spaces in x +!! @param ny numbers of the grid spaces in y !! @param a parameters of an ESG mapping centered at (plat,plon) !! @param k parameters of an ESG mapping centered at (plat,plon) !! @param plat latitude define centered mapping !! @param plon longitude define centered mapping +!! @param pazi ??? !! @param delx central x-spacing grid point !! @param dely central y-spacing grid point -!! @param lx x grid index for left-lower corner of the grid at center -!! @param ly y grid index for left-lower corner of the grid at center -!! @param nx numbers of the grid spaces in x -!! @param ny numbers of the grid spaces in y !! @param xc unit cartesian vectors +!! @param xcd ??? !! @param garea rectangular array -!! @param nx-1 x dimensions for garea -!! @param ny-1 y dimensions for garea !! @param ff failure flag !! @author R. J. Purser subroutine hgrid_ak_rc(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rc] @@ -1493,25 +1598,31 @@ subroutine hgrid_ak_rc(lx,ly,nx,ny,A,K,plat,plon,pazi, & ! [hgrid_ak_rc] enddo end subroutine hgrid_ak_rc -!> Use a and k as the parameters of an Extended Schmidt-transformed -!! Gnomonic (ESG) mapping centered at (pdlat,pdlon) and twisted about this -!! center. -!! by an azimuth angle of pdazi counterclockwise (these angles in degrees). -!! Like hgrid_ak_rr, return the grid points' lats and lons, except that here -!! the angles are returned in degrees. Garea, the area of each grid cell, is -!! returned as in hgrid_ak_rr, and a failure flag, ff, raised when a failure -!! occurs anywhere in these calculations. -!! @param a parameters of an ESG mapping -!! @param k parameters of an ESG mapping -!! @param pdlat latitude define centered mapping -!! @param pdlon longitude define centered mapping -!! @param lx x grid index for left-lower corner of the grid at center -!! @param ly y grid index for left-lower corner of the grid at center -!! @param nx numbers of the grid spaces in x -!! @param ny numbers of the grid spaces in y -!! @param delx central x-spacing grid point -!! @param dely central y-spacing grid point -!! @param ff failure flag +!! Use a and k as the parameters of an Extended Schmidt-transformed +!! Gnomonic (ESG) mapping centered at (pdlat,pdlon) and twisted about +!! this center by an azimuth angle of pdazi counterclockwise (these +!! angles in degrees). +!! +!! Like hgrid_ak_rr, return the grid points' lats and lons, except +!! that here the angles are returned in degrees. Garea, the area of +!! each grid cell, is returned as in hgrid_ak_rr, and a failure flag, +!! ff, raised when a failure occurs anywhere in these calculations. +!! +!! @param[in] lx x grid index for left-lower corner of the grid at center +!! @param[in] ly y grid index for left-lower corner of the grid at center +!! @param[in] nx numbers of the grid spaces in x +!! @param[in] ny numbers of the grid spaces in y +!! @param[in] a parameters of an ESG mapping +!! @param[in] k parameters of an ESG mapping +!! @param[in] pdlat latitude define centered mapping +!! @param[in] pdlon longitude define centered mapping +!! @param[in] pdazi ??? +!! @param[in] delx central x-spacing grid point +!! @param[in] dely central y-spacing grid point +!! @param[out] gdlat ??? +!! @param[out] gdlon ??? +!! @param[out] garea ??? +!! @param[out] ff failure flag !! @author R. J. Purser subroutine hgrid_ak_dd(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dd] delx,dely, gdlat,gdlon,garea, ff) @@ -1533,21 +1644,28 @@ subroutine hgrid_ak_dd(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dd] gdlon=gdlon*rtod ! end subroutine hgrid_ak_dd -!> Like hgrid_ak_rr_c, except all the angle arguments (but not delx,dely) -!! are in degrees instead of radians. -!! @param a parameters of an ESG mapping -!! @param k parameters of an ESG mapping -!! @param pdlat latitude define centered mapping -!! @param pdlon longitude define centered mapping -!! @param lx x grid index for left-lower corner of the grid at center -!! @param ly y grid index for left-lower corner of the grid at center -!! @param nx numbers of the grid spaces in x -!! @param ny numbers of the grid spaces in y -!! @param delx central x-spacing grid point -!! @param dely central y-spacing grid point -!! @param dangle_dx x rotations of the steps -!! @param dangle_dy y rotations of the steps -!! @param ff failure flag +!! Like hgrid_ak_rr_c, except all the angle arguments (but not +!! delx,dely) are in degrees instead of radians. +!! +!! @param[in] lx x grid index for left-lower corner of the grid at center +!! @param[in] ly y grid index for left-lower corner of the grid at center +!! @param[in] nx numbers of the grid spaces in x +!! @param[in] ny numbers of the grid spaces in y +!! @param[in] a parameters of an ESG mapping +!! @param[in] k parameters of an ESG mapping +!! @param[in] pdlat latitude define centered mapping +!! @param[in] pdlon longitude define centered mapping +!! @param[in] pdazi ??? +!! @param[in] delx central x-spacing grid point +!! @param[in] dely central y-spacing grid point +!! @param[out] gdlat ??? +!! @param[out] gdlon ??? +!! @param[out] garea ??? +!! @param[out] dx ??? +!! @param[out] dy ??? +!! @param[out] dangle_dx x rotations of the steps +!! @param[out] dangle_dy y rotations of the steps +!! @param[out] ff failure flag !! @author R. J. Purser subroutine hgrid_ak_dd_c(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, &! [hgrid_ak_dd] delx,dely, gdlat,gdlon,garea,dx,dy,dangle_dx,dangle_dy, ff) @@ -1575,27 +1693,30 @@ subroutine hgrid_ak_dd_c(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, &! [hgrid_ak_dd] end subroutine hgrid_ak_dd_c !> Use a and k as the parameters of an Extended Schmidt-transformed -!! Gnomonic (ESG) mapping centered at (pdlat,pdlon) and twisted about this -!! center by an azimuth angle of pdazi counterclockwise (these angles in -!! degrees). -!! Like hgrid_ak_rx, return the grid points' cartesians xc and Jacobian -!! matrices, xcd. Garea, the area of each grid cell, is also -!! returned as in hgrid_ak_rx, and a failure flag, ff, raised when a failure -!! occurs anywhere in these calculations. -!! @param a parameters of an ESG mapping -!! @param k parameters of an ESG mapping -!! @param pdlat latitude define centered mapping -!! @param pdlon longitude define centered mapping -!! @param lx x grid index for left-lower corner of the grid at center -!! @param ly y grid index for left-lower corner of the grid at center -!! @param nx numbers of the grid spaces in x -!! @param ny numbers of the grid spaces in y -!! @param delx central x-spacing grid point -!! @param dely central y-spacing grid point -!! @param xc grid points' cartesians -!! @param xcd Jacobian matrices -!! @param garea Jacobian matrices -!! @param ff failure flag +!! Gnomonic (ESG) mapping centered at (pdlat,pdlon) and twisted about +!! this center by an azimuth angle of pdazi counterclockwise (these +!! angles in degrees). +!! +!! Like hgrid_ak_rx, return the grid points' cartesians xc and +!! Jacobian matrices, xcd. Garea, the area of each grid cell, is also +!! returned as in hgrid_ak_rx, and a failure flag, ff, raised when a +!! failure occurs anywhere in these calculations. +!! +!! @param[in] lx x grid index for left-lower corner of the grid at center +!! @param[in] ly y grid index for left-lower corner of the grid at center +!! @param[in] nx numbers of the grid spaces in x +!! @param[in] ny numbers of the grid spaces in y +!! @param[in] a parameters of an ESG mapping +!! @param[in] k parameters of an ESG mapping +!! @param[in] pdlat latitude define centered mapping +!! @param[in] pdlon longitude define centered mapping +!! @param[in] pdazi ??? +!! @param[in] delx central x-spacing grid point +!! @param[in] dely central y-spacing grid point +!! @param[out] xc grid points' cartesians +!! @param[out] xcd Jacobian matrices +!! @param[out] garea Jacobian matrices +!! @param[out] ff failure flag !! @author R. J. Purser subroutine hgrid_ak_dc(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dc] delx,dely, xc,xcd,garea, ff) @@ -1614,23 +1735,30 @@ subroutine hgrid_ak_dc(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & ! [hgrid_ak_dc] delx,dely, xc,xcd,garea, ff) end subroutine hgrid_ak_dc -!> Like hgrid_ak_rr_c except the argument list includes the earth radius, re, -!! and this is used to express the map-space grid increments in the dimensional -!! units, delxre, delyre on entry, and to express the grid cell areas, garea, -!! in dimensional units upon return. +!> Like hgrid_ak_rr_c except the argument list includes the earth +!! radius, re, and this is used to express the map-space grid +!! increments in the dimensional units, delxre, delyre on entry, and +!! to express the grid cell areas, garea, in dimensional units upon +!! return. +!! !! The gridded lats and lons, glat and glon, remain in radians. -!! @param a parameters of an ESG mapping -!! @param k parameters of an ESG mapping -!! @param pdlat latitude define centered mapping -!! @param pdlon longitude define centered mapping -!! @param lx x grid index for left-lower corner of the grid at center -!! @param ly y grid index for left-lower corner of the grid at center -!! @param nx numbers of the grid spaces in x -!! @param ny numbers of the grid spaces in y -!! @param re earth radius -!! @param delxre map-space grid increments in the dimensional units -!! @param delyre map-space grid increments in the dimensional units -!! @param ff failure flag +!! +!! @param[in] lx x grid index for left-lower corner of the grid at center +!! @param[in] ly y grid index for left-lower corner of the grid at center +!! @param[in] nx numbers of the grid spaces in x +!! @param[in] ny numbers of the grid spaces in y +!! @param[in] a parameters of an ESG mapping +!! @param[in] k parameters of an ESG mapping +!! @param[in] plat latitude define centered mapping +!! @param[in] plon longitude define centered mapping +!! @param[in] pazi ??? +!! @param[in] re earth radius +!! @param[in] delxre map-space grid increments in the dimensional units +!! @param[in] delyre map-space grid increments in the dimensional units +!! @param[out] glat grid points for latitude +!! @param[out] glon grid points for longitude +!! @param[out] garea rectangular array +!! @param[out] ff failure flag !! @author R. J. Purser subroutine hgrid_ak(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak] re,delxre,delyre, glat,glon,garea, ff) @@ -1651,29 +1779,38 @@ subroutine hgrid_ak(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak] garea=garea*rere ! <- Convert from steradians to physical area units. end subroutine hgrid_ak -!> Like hgrid_ak_rr_c except the argument list includes the earth radius, re, -!! and this is used to express the map-space grid increments in the dimensional -!! units, delxre, delyre on entry, and to express the grid cell areas, garea, -!! and the x- and y- grid steps, dx and dy, in dimensional units upon return. -!! The gridded lats and lons, glat and glon, remain in radians. -!! Also, in order for the argument list to remain compatible with an earlier -!! version of this routine, the relative rotations of the steps, dangle_dx -!! and dangle_dy, are returned as degrees instead of radians (all other angles -!! in the argument list, i.e., plat,plon,pazi,glat,glon, remain radians). -!! @param a Extended Schmidt Gnomonic parameter -!! @param k Extended Schmidt Gnomonic parameter -!! @param plat latitude define centered mapping -!! @param plon longitude define centered mapping -!! @param re earth radius -!! @param delxre map-space grid increments in the dimensional units -!! @param delyre map-space grid increments in the dimensional units -!! @param garea grid cell areas -!! @param dx x- grid steps -!! @param dy y- grid steps -!! @param glat gridded lats -!! @param glon gridded lons -!! @param dangle_dx x rotations of the steps -!! @param dangle_dy y rotations of the steps +!> Like hgrid_ak_rr_c except the argument list includes the earth +!! radius, re, and this is used to express the map-space grid +!! increments in the dimensional units, delxre, delyre on entry, and +!! to express the grid cell areas, garea, and the x- and y- grid +!! steps, dx and dy, in dimensional units upon return. The gridded +!! lats and lons, glat and glon, remain in radians. Also, in order +!! for the argument list to remain compatible with an earlier version +!! of this routine, the relative rotations of the steps, dangle_dx and +!! dangle_dy, are returned as degrees instead of radians (all other +!! angles in the argument list, i.e., plat,plon,pazi,glat,glon, remain +!! radians). +!! +!! @param[in] lx x grid index for left-lower corner of the grid at center +!! @param[in] ly y grid index for left-lower corner of the grid at center +!! @param[in] nx ??? +!! @param[in] ny ??? +!! @param[in] a Extended Schmidt Gnomonic parameter +!! @param[in] k Extended Schmidt Gnomonic parameter +!! @param[in] plat latitude define centered mapping +!! @param[in] plon longitude define centered mapping +!! @param[in] pazi ??? +!! @param[in] re earth radius +!! @param[in] delxre map-space grid increments in the dimensional units +!! @param[in] delyre map-space grid increments in the dimensional units +!! @param[out] glat gridded lats +!! @param[out] glon gridded lons +!! @param[out] garea grid cell areas +!! @param[out] dx x- grid steps +!! @param[out] dy y- grid steps +!! @param[out] dangle_dx x rotations of the steps +!! @param[out] dangle_dy y rotations of the steps +!! @param[out] ff failure flag !! @author R. J. Purser subroutine hgrid_ak_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak] re,delxre,delyre, glat,glon,garea,dx,dy,dangle_dx,dangle_dy, ff) @@ -1701,9 +1838,11 @@ subroutine hgrid_ak_c(lx,ly,nx,ny,a,k,plat,plon,pazi, & ! [hgrid_ak] dangle_dy=dangle_dy*rtod ! <-Convert from radians to degrees end subroutine hgrid_ak_c -!> This Gauss-Legendre quadrature integrates exactly any even polynomial -!! up to degree m*4-2 in the half-interval [0,1]. This code is liberally -!! adapted from the algorithm given in Press et al., Numerical Recipes. +!> This Gauss-Legendre quadrature integrates exactly any even +!! polynomial up to degree m*4-2 in the half-interval [0,1]. This code +!! is liberally adapted from the algorithm given in Press et al., +!! Numerical Recipes. +!! !! @param m number of nodes in half-interval !! @param x nodes and weights !! @param w nodes and weights @@ -1729,17 +1868,20 @@ subroutine gaulegh(m,x,w)! [gaulegh] enddo end subroutine gaulegh -!> Given the map specification (angles in radians), the grid spacing (in -!! map-space units) and the sample lat-lon (in radian), return the the -!! image in map space in a 2-vector in grid units. If the transformation -!! is invalid, return a .true. failure flag. -!! @param a parameters of an ESG mapping -!! @param k parameters of an ESG mapping -!! @param plat latitude define centered mapping -!! @param plon longitude define centered mapping -!! @param lat radian latitude -!! @param lon radian longitude -!! @param ff failure flag +!> Given the map specification (angles in radians), the grid spacing +!! (in map-space units) and the sample lat-lon (in radian), return the +!! the image in map space in a 2-vector in grid units. If the +!! transformation is invalid, return a .true. failure flag. +!! +!! @param[in] a parameters of an ESG mapping +!! @param[in] k parameters of an ESG mapping +!! @param[in] plat latitude define centered mapping +!! @param[in] plon longitude define centered mapping +!! @param[in] pazi ??? +!! @param[in] lat radian latitude +!! @param[in] lon radian longitude +!! @param[out] xm ??? +!! @param[out] ff failure flag !! @author R. J. Purser subroutine gtoxm_ak_rr_m(A,K,plat,plon,pazi,lat,lon,xm,ff)! [gtoxm_ak_rr] use pmat5, only: grtoc @@ -1768,19 +1910,22 @@ subroutine gtoxm_ak_rr_m(A,K,plat,plon,pazi,lat,lon,xm,ff)! [gtoxm_ak_rr] call xctoxm_ak(a,k,xc,xm,ff) end subroutine gtoxm_ak_rr_m -!> Given the map specification (angles in radians), the grid spacing (in -!! map-space units) and the sample lat-lon (in radian), return the the -!! image in map space in a 2-vector in grid units. If the transformation -!! is invalid, return a .true. failure flag. -!! @param a parameters of an ESG mapping -!! @param k parameters of an ESG mapping -!! @param plat latitude define centered mapping -!! @param plon longitude define centered mapping -!! @param delx central x-spacing grid point -!! @param dely central y-spacing grid point -!! @param lat radian latitude -!! @param lon radian longitude -!! @param ff failure flag +!> Given the map specification (angles in radians), the grid spacing +!! (in map-space units) and the sample lat-lon (in radian), return the +!! the image in map space in a 2-vector in grid units. If the +!! transformation is invalid, return a .true. failure flag. +!! +!! @param[in] a parameters of an ESG mapping +!! @param[in] k parameters of an ESG mapping +!! @param[in] plat latitude define centered mapping +!! @param[in] plon longitude define centered mapping +!! @param[in] pazi ??? +!! @param[in] delx central x-spacing grid point +!! @param[in] dely central y-spacing grid point +!! @param[in] lat radian latitude +!! @param[in] lon radian longitude +!! @param[out] xm ??? +!! @param[out] ff failure flag !! @author R. J. Purser subroutine gtoxm_ak_rr_g(A,K,plat,plon,pazi,delx,dely,lat,lon,&! [gtoxm_ak_rr] xm,ff) @@ -1792,14 +1937,17 @@ subroutine gtoxm_ak_rr_g(A,K,plat,plon,pazi,delx,dely,lat,lon,&! [gtoxm_ak_rr] xm(1)=xm(1)/delx; xm(2)=xm(2)/dely end subroutine gtoxm_ak_rr_g -!> Like gtoxm_ak_rr_m, except input angles are expressed in degrees -!! @param a parameters of an ESG mapping -!! @param k parameters of an ESG mapping -!! @param pdlat latitude define centered mapping -!! @param pdlon longitude define centered mapping -!! @param dlat radian latitude -!! @param dlon radian longitude -!! @param ff failure flag +!> Like gtoxm_ak_rr_m, except input angles are expressed in degrees. +!! +!! @param[in] a parameters of an ESG mapping +!! @param[in] k parameters of an ESG mapping +!! @param[in] pdlat latitude define centered mapping +!! @param[in] pdlon longitude define centered mapping +!! @param[in] pdazi ??? +!! @param[in] dlat radian latitude +!! @param[in] dlon radian longitude +!! @param[out] xm ??? +!! @param[out] ff failure flag !! @author R. J. Purser subroutine gtoxm_ak_dd_m(A,K,pdlat,pdlon,pdazi,dlat,dlon,&! [gtoxm_ak_dd] xm,ff) @@ -1816,14 +1964,17 @@ subroutine gtoxm_ak_dd_m(A,K,pdlat,pdlon,pdazi,dlat,dlon,&! [gtoxm_ak_dd] call gtoxm_ak_rr_m(A,K,plat,plon,pazi,lat,lon,xm,ff) end subroutine gtoxm_ak_dd_m -!> Like gtoxm_ak_rr_g, except input angles are expressed in degrees -!! @param a parameters of an ESG mapping -!! @param k parameters of an ESG mapping -!! @param pdlat latitude define centered mapping -!! @param pdlon longitude define centered mapping -!! @param delx central x-spacing grid point -!! @param dely central y-spacing grid point -!! @param ff failure flag +!> Like gtoxm_ak_rr_g, except input angles are expressed in degrees. +!! +!! @param[in] a parameters of an ESG mapping +!! @param[in] k parameters of an ESG mapping +!! @param[in] pdlat latitude define centered mapping +!! @param[in] pdlon longitude define centered mapping +!! @param[in] pdazi ??? +!! @param[in] delx central x-spacing grid point +!! @param[in] dely central y-spacing grid point +!! @param[out] xm ??? +!! @param[out] ff failure flag !! @author R. J. Purser subroutine gtoxm_ak_dd_g(A,K,pdlat,pdlon,pdazi,delx,dely,&! [gtoxm_ak_dd] dlat,dlon, xm,ff) @@ -1843,15 +1994,19 @@ end subroutine gtoxm_ak_dd_g !> Given the ESG map specified by parameters (A,K) and geographical !! orientation, plat,plon,pazi (radians), and a position, in map-space !! coordinates given by the 2-vector, xm, return the geographical -!! coordinates, lat and lon (radians). If the transformation is invalid for -!! any reason, return instead with a raised failure flag, FF= .true. -!! @param a parameters of an ESG mapping -!! @param k parameters of an ESG mapping -!! @param plat latitude of geographical orientation -!! @param plon longitude of geographical orientation -!! @param lat latitude of geographical coordinates -!! @param lon longitude of geographical coordinates -!! @param ff failure flag +!! coordinates, lat and lon (radians). If the transformation is +!! invalid for any reason, return instead with a raised failure flag, +!! FF= .true. +!! +!! @param[in] a parameters of an ESG mapping +!! @param[in] k parameters of an ESG mapping +!! @param[in] plat latitude of geographical orientation +!! @param[in] plon longitude of geographical orientation +!! @param[in] pazi ??? +!! @param[in] xm ??? +!! @param[out] lat latitude of geographical coordinates +!! @param[out] lon longitude of geographical coordinates +!! @param[out] ff failure flag !! @author R. J. Purser subroutine xmtog_ak_rr_m(A,K,plat,plon,pazi,xm,lat,lon,ff)! [xmtog_ak_rr] use pmat5, only: ctogr @@ -1881,21 +2036,24 @@ subroutine xmtog_ak_rr_m(A,K,plat,plon,pazi,xm,lat,lon,ff)! [xmtog_ak_rr] call ctogr(xc,lat,lon) end subroutine xmtog_ak_rr_m -!> For an ESG map with parameters, (A,K), and geographical orientation, -!! given by plon,plat,pazi (radians), and given a point in grid-space units -!! as the 2-vector, xm, return the geographical coordinates, lat, lon, (radians) -!! of this point. If instead the transformation is invalid for any reason, then -!! return the raised failure flag, FF=.true. -!! @param a parameters of an ESG mapping -!! @param k parameters of an ESG mapping -!! @param plon longitude geographical orientation -!! @param plat latitude geographical orientation -!! @param xm grid-space units as the 2-vector -!! @param lat latitude geographical coordinates -!! @param lon longitude geographical coordinates -!! @param delx central x-spacing grid point -!! @param dely central y-spacing grid point -!! @param ff failure flag +!> For an ESG map with parameters, (A,K), and geographical +!! orientation, given by plon,plat,pazi (radians), and given a point +!! in grid-space units as the 2-vector, xm, return the geographical +!! coordinates, lat, lon, (radians) of this point. If instead the +!! transformation is invalid for any reason, then return the raised +!! failure flag, FF=.true. +!! +!! @param[in] a parameters of an ESG mapping +!! @param[in] k parameters of an ESG mapping +!! @param[in] plat latitude geographical orientation +!! @param[in] plon longitude geographical orientation +!! @param[in] pazi ??? +!! @param[in] delx central x-spacing grid point +!! @param[in] dely central y-spacing grid point +!! @param[in] xm grid-space units as the 2-vector +!! @param[out] lat latitude geographical coordinates +!! @param[out] lon longitude geographical coordinates +!! @param[out] ff failure flag !! @author R. J. Purser subroutine xmtog_ak_rr_g(A,K,plat,plon,pazi,delx,dely,xm,&! [xmtog_ak_rr] lat,lon,ff) @@ -1910,12 +2068,17 @@ subroutine xmtog_ak_rr_g(A,K,plat,plon,pazi,delx,dely,xm,&! [xmtog_ak_rr] call xmtog_ak_rr_m(A,K,plat,plon,pazi,xmt,lat,lon,ff) end subroutine xmtog_ak_rr_g -!> Like xmtog_ak_rr_m, except angles are expressed in degrees -!! @param a parameters of an ESG mapping -!! @param k parameters of an ESG mapping -!! @param plat latitude of geographical orientation -!! @param plon longitude of geographical orientation -!! @param ff failure flag +!> Like xmtog_ak_rr_m, except angles are expressed in degrees. +!! +!! @param[in] a parameters of an ESG mapping +!! @param[in] k parameters of an ESG mapping +!! @param[in] pdlat ??? +!! @param[in] pdlon ??? +!! @param[in] pdazi ??? +!! @param[in] xm ??? +!! @param[out] dlat latitude of geographical orientation +!! @param[out] dlon longitude of geographical orientation +!! @param[out] ff failure flag !! @author R. J. Purser subroutine xmtog_ak_dd_m(A,K,pdlat,pdlon,pdazi,xm,dlat,dlon,ff)! [xmtog_ak_dd] use pmat5, only: ctogr @@ -1933,17 +2096,19 @@ subroutine xmtog_ak_dd_m(A,K,pdlat,pdlon,pdazi,xm,dlat,dlon,ff)! [xmtog_ak_dd] dlon=lon*rtod end subroutine xmtog_ak_dd_m -!> Like xmtog_ak_rr_g, except angles are expressed in degrees -!! @param a parameters of an ESG mapping -!! @param k parameters of an ESG mapping -!! @param pdlat latitude define centered mapping -!! @param pdlon longitude define centered mapping -!! @param dlat latitude radians angle -!! @param dlon longitude radians angle -!! @param xm grid-space -!! @param delx central x-spacing grid point -!! @param dely central y-spacing grid point -!! @param ff failure flag +!> Like xmtog_ak_rr_g, except angles are expressed in degrees. +!! +!! @param[in] a parameters of an ESG mapping +!! @param[in] k parameters of an ESG mapping +!! @param[in] pdlat latitude define centered mapping +!! @param[in] pdlon longitude define centered mapping +!! @param[in] pdazi ??? +!! @param[in] delx central x-spacing grid point +!! @param[in] dely central y-spacing grid point +!! @param[in] xm grid-space +!! @param[out] dlat latitude radians angle +!! @param[out] dlon longitude radians angle +!! @param[out] ff failure flag !! @author R. J. Purser subroutine xmtog_ak_dd_g(A,K,pdlat,pdlon,pdazi,delx,dely,xm,&! [xmtog_ak_dd] dlat,dlon,ff)