From 2e58280e97304a3ffb147619673583c4c52b3a10 Mon Sep 17 00:00:00 2001 From: "Lin.Gan" Date: Tue, 2 Mar 2021 18:54:28 +0000 Subject: [PATCH 1/3] Part of #366 Doxygen improvements for sorc/grid_tools.fd/regional_esg_grid.fd/pmat4.f90 --- .../regional_esg_grid.fd/pmat4.f90 | 824 +++++++++--------- 1 file changed, 411 insertions(+), 413 deletions(-) diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/pmat4.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/pmat4.f90 index 454d6b677..8eb24628a 100644 --- a/sorc/grid_tools.fd/regional_esg_grid.fd/pmat4.f90 +++ b/sorc/grid_tools.fd/regional_esg_grid.fd/pmat4.f90 @@ -1,8 +1,7 @@ !> @file -!! @author R. J. Purser @date Oct 2005 -!! -!! Euclidean geometry, geometric (stereographic) projections, +!! @brief Euclidean geometry, geometric (stereographic) projections, !! related transformations (Mobius). +!! !! Package for handy vector and matrix operations in Euclidean geometry. !! This package is primarily intended for 3D operations and three of the !! functions (Cross_product, Triple_product and Axial) do not possess simple @@ -16,7 +15,6 @@ !! exponentials of matrices (without resort to eigen methods). Also added !! Quaternion and spinor representations of 3D rotations, and their !! conversion routines. -!! !! FUNCTION: !!- absv: Absolute magnitude of vector as its euclidean length !!- Normalized: Normalized version of given real vector @@ -45,7 +43,10 @@ !! DIRECT DEPENDENCIES !! Libraries[their Modules]: pmat[pmat] !! Additional Modules : pkind, pietc -!! +!! @author R. J. Purser @date Oct 2005 + +!> Module for handy vector and matrix operations in Euclidean geometry +!! @author R. J. Purser module pmat4 !============================================================================ use pkind, only: spi,sp,dp,dpc @@ -114,26 +115,27 @@ module pmat4 contains -!============================================================================= +!> sqrt calculation for absv_s function +!! @author R. J. Purser function absv_s(a)result(s)! [absv] -!============================================================================= implicit none real(sp),dimension(:),intent(in):: a real(sp) :: s s=sqrt(dot_product(a,a)) end function absv_s -!============================================================================= + +!> sqrt calculation for absv_d function +!! @author R. J. Purser function absv_d(a)result(s)! [absv] -!============================================================================= implicit none real(dp),dimension(:),intent(in):: a real(dp) :: s s=sqrt(dot_product(a,a)) end function absv_d -!============================================================================= +!> calculation for normalized_s function +!! @author R. J. Purser function normalized_s(a)result(b)! [normalized] -!============================================================================= use pietc_s, only: u0 implicit none real(sp),dimension(:),intent(IN):: a @@ -141,9 +143,10 @@ function normalized_s(a)result(b)! [normalized] real(sp) :: s s=absv_s(a); if(s==u0)then; b=u0;else;b=a/s;endif end function normalized_s -!============================================================================= + +!> calculation for normalized_d function +!! @author R. J. Purser function normalized_d(a)result(b)! [normalized] -!============================================================================= use pietc, only: u0 implicit none real(dp),dimension(:),intent(IN):: a @@ -152,9 +155,9 @@ function normalized_d(a)result(b)! [normalized] s=absv_d(a); if(s==u0)then; b=u0;else;b=a/s;endif end function normalized_d -!============================================================================= +!> calculation for orthogonalized_s function +!! @author R. J. Purser function orthogonalized_s(u,a)result(b)! [orthogonalized] -!============================================================================= implicit none real(sp),dimension(:),intent(in):: u,a real(sp),dimension(size(u)) :: b @@ -162,9 +165,10 @@ function orthogonalized_s(u,a)result(b)! [orthogonalized] ! Note: this routine assumes u is already normalized s=dot_product(u,a); b=a-u*s end function orthogonalized_s -!============================================================================= + +!> calculation for orthogonalized_d function +!! @author R. J. Purser function orthogonalized_d(u,a)result(b)! [orthogonalized] -!============================================================================= implicit none real(dp),dimension(:),intent(in):: u,a real(dp),dimension(size(u)) :: b @@ -173,36 +177,38 @@ function orthogonalized_d(u,a)result(b)! [orthogonalized] s=dot_product(u,a); b=a-u*s end function orthogonalized_d -!============================================================================= +!> calculation for cross_product_s function +!! @author R. J. Purser function cross_product_s(a,b)result(c)! [cross_product] -!============================================================================= implicit none real(sp),dimension(3),intent(in):: a,b real(sp),dimension(3) :: c c(1)=a(2)*b(3)-a(3)*b(2); c(2)=a(3)*b(1)-a(1)*b(3); c(3)=a(1)*b(2)-a(2)*b(1) end function cross_product_s -!============================================================================= + +!> calculation for cross_product_d function +!! @author R. J. Purser function cross_product_d(a,b)result(c)! [cross_product] -!============================================================================= implicit none real(dp),dimension(3),intent(in):: a,b real(dp),dimension(3) :: c c(1)=a(2)*b(3)-a(3)*b(2); c(2)=a(3)*b(1)-a(1)*b(3); c(3)=a(1)*b(2)-a(2)*b(1) end function cross_product_d -!============================================================================= + +!> Deliver the triple-cross-product, x, of the +!! three 4-vectors, u, v, w, with the sign convention +!! that ordered, {u,v,w,x} form a right-handed quartet +!! in the generic case (determinant >= 0). +!! @param x triple-cross-product vector +!! @param u vector +!! @param v vector +!! @param w vector +!! @author R. J. Purser function triple_cross_product_s(u,v,w)result(x)! [cross_product] -!============================================================================= -! Deliver the triple-cross-product, x, of the -! three 4-vectors, u, v, w, with the sign convention -! that ordered, {u,v,w,x} form a right-handed quartet -! in the generic case (determinant >= 0). -!============================================================================= implicit none real(sp),dimension(4),intent(in ):: u,v,w real(sp),dimension(4) :: x -!----------------------------------------------------------------------------- real(sp):: uv12,uv13,uv14,uv23,uv24,uv34 -!============================================================================= uv12=u(1)*v(2)-u(2)*v(1); uv13=u(1)*v(3)-u(3)*v(1); uv14=u(1)*v(4)-u(4)*v(1) uv23=u(2)*v(3)-u(3)*v(2); uv24=u(2)*v(4)-u(4)*v(2) uv34=u(3)*v(4)-u(4)*v(3) @@ -211,15 +217,14 @@ function triple_cross_product_s(u,v,w)result(x)! [cross_product] x(3)=-uv12*w(4) +uv14*w(2)-uv24*w(1) x(4)= uv12*w(3)-uv13*w(2)+uv23*w(1) end function triple_cross_product_s -!============================================================================= + +!> calculation for triple_cross_product_d function +!! @author R. J. Purser function triple_cross_product_d(u,v,w)result(x)! [cross_product] -!============================================================================= implicit none real(dp),dimension(4),intent(in ):: u,v,w real(dp),dimension(4) :: x -!----------------------------------------------------------------------------- real(dp):: uv12,uv13,uv14,uv23,uv24,uv34 -!============================================================================= uv12=u(1)*v(2)-u(2)*v(1); uv13=u(1)*v(3)-u(3)*v(1); uv14=u(1)*v(4)-u(4)*v(1) uv23=u(2)*v(3)-u(3)*v(2); uv24=u(2)*v(4)-u(4)*v(2) uv34=u(3)*v(4)-u(4)*v(3) @@ -229,7 +234,8 @@ function triple_cross_product_d(u,v,w)result(x)! [cross_product] x(4)= uv12*w(3)-uv13*w(2)+uv23*w(1) end function triple_cross_product_d -!============================================================================= +!> calculation for outer_product_s function +!! @author R. J. Purser function outer_product_s(a,b)result(c)! [outer_product] !============================================================================= implicit none @@ -240,9 +246,10 @@ function outer_product_s(a,b)result(c)! [outer_product] nb=size(b) do i=1,nb; c(:,i)=a*b(i); enddo end function outer_product_s -!============================================================================= + +!> Calculation for outer_product_d function +!! @author R. J. Purser function outer_product_d(a,b)result(c)! [outer_product] -!============================================================================= implicit none real(dp),dimension(:), intent(in ):: a real(dp),dimension(:), intent(in ):: b @@ -251,9 +258,10 @@ function outer_product_d(a,b)result(c)! [outer_product] nb=size(b) do i=1,nb; c(:,i)=a*b(i); enddo end function outer_product_d -!============================================================================= + +!> Calculation for outer_product_i function +!! @author R. J. Purser function outer_product_i(a,b)result(c)! [outer_product] -!============================================================================= implicit none integer(spi),dimension(:), intent(in ):: a integer(spi),dimension(:), intent(in ):: b @@ -263,26 +271,27 @@ function outer_product_i(a,b)result(c)! [outer_product] do i=1,nb; c(:,i)=a*b(i); enddo end function outer_product_i -!============================================================================= +!> Calculation for triple_product_s function +!! @author R. J. Purser function triple_product_s(a,b,c)result(tripleproduct)! [triple_product] -!============================================================================= implicit none real(sp),dimension(3),intent(IN ):: a,b,c real(sp) :: tripleproduct tripleproduct=dot_product( cross_product(a,b),c ) end function triple_product_s -!============================================================================= + +!> Calculation for triple_product_d function +!! @author R. J. Purser function triple_product_d(a,b,c)result(tripleproduct)! [triple_product] -!============================================================================= implicit none real(dp),dimension(3),intent(IN ):: a,b,c real(dp) :: tripleproduct tripleproduct=dot_product( cross_product(a,b),c ) end function triple_product_d -!============================================================================= +!> Calculation for det_s function +!! @author R. J. Purser function det_s(a)result(det)! [det] -!============================================================================= use pietc_s, only: u0 implicit none real(sp),dimension(:,:),intent(IN ) :: a @@ -297,9 +306,10 @@ function det_s(a)result(det)! [det] if(nrank Calculation for det_d function +!! @author R. J. Purser function det_d(a)result(det)! [det] -!============================================================================= use pietc, only: u0 implicit none real(dp),dimension(:,:),intent(IN ) ::a @@ -314,9 +324,10 @@ function det_d(a)result(det)! [det] if(nrank Calculation for det_i function +!! @author R. J. Purser function det_i(a)result(idet)! [det] -!============================================================================= implicit none integer(spi), dimension(:,:),intent(IN ) :: a integer(spi) :: idet @@ -325,9 +336,9 @@ function det_i(a)result(idet)! [det] b=a; bdet=det(b); idet=nint(bdet) end function det_i -!============================================================================= +!> Calculation for det_id function +!! @author R. J. Purser function det_id(a)result(idet)! [det] -!============================================================================= use pkind, only: dp,dpi implicit none integer(dpi), dimension(:,:),intent(IN ):: a @@ -337,36 +348,39 @@ function det_id(a)result(idet)! [det] b=a; bdet=det(b); idet=nint(bdet) end function det_id -!============================================================================= +!> Calculation for axial3_s function +!! @author R. J. Purser function axial3_s(a)result(b)! [axial] -!============================================================================= use pietc_s, only: u0 implicit none real(sp),dimension(3),intent(IN ):: a real(sp),dimension(3,3) :: b b=u0;b(3,2)=a(1);b(1,3)=a(2);b(2,1)=a(3);b(2,3)=-a(1);b(3,1)=-a(2);b(1,2)=-a(3) end function axial3_s -!============================================================================= + +!> Calculation for axial3_d function +!! @author R. J. Purser function axial3_d(a)result(b)! [axial] -!============================================================================= use pietc, only: u0 implicit none real(dp),dimension(3),intent(IN ):: a real(dp),dimension(3,3) :: b b=u0;b(3,2)=a(1);b(1,3)=a(2);b(2,1)=a(3);b(2,3)=-a(1);b(3,1)=-a(2);b(1,2)=-a(3) end function axial3_d -!============================================================================= + +!> Calculation for axial33_s function +!! @author R. J. Purser function axial33_s(b)result(a)! [axial] -!============================================================================= use pietc_s, only: o2 implicit none real(sp),dimension(3,3),intent(IN ):: b real(sp),dimension(3) :: a a(1)=(b(3,2)-b(2,3))*o2; a(2)=(b(1,3)-b(3,1))*o2; a(3)=(b(2,1)-b(1,2))*o2 end function axial33_s -!============================================================================= + +!> Calculation for axial33_d function +!! @author R. J. Purser function axial33_d(b)result(a)! [axial] -!============================================================================= use pietc, only: o2 implicit none real(dp),dimension(3,3),intent(IN ):: b @@ -374,9 +388,9 @@ function axial33_d(b)result(a)! [axial] a(1)=(b(3,2)-b(2,3))*o2; a(2)=(b(1,3)-b(3,1))*o2; a(3)=(b(2,1)-b(1,2))*o2 end function axial33_d -!============================================================================= +!> Calculation for diagn_s function +!! @author R. J. Purser function diagn_s(a)result(b)! [diag] -!============================================================================= use pietc, only: u0 implicit none real(sp),dimension(:),intent(IN ) :: a @@ -385,9 +399,10 @@ function diagn_s(a)result(b)! [diag] n=size(a) b=u0; do i=1,n; b(i,i)=a(i); enddo end function diagn_s -!============================================================================= + +!> Calculation for diagn_d function +!! @author R. J. Purser function diagn_d(a)result(b)! [diag] -!============================================================================= use pietc, only: u0 implicit none real(dp),dimension(:),intent(IN ) :: a @@ -396,9 +411,10 @@ function diagn_d(a)result(b)! [diag] n=size(a) b=u0; do i=1,n; b(i,i)=a(i); enddo end function diagn_d -!============================================================================= + +!> Calculation for diagn_i function +!! @author R. J. Purser function diagn_i(a)result(b)! [diag] -!============================================================================= implicit none integer(spi),dimension(:),intent(IN ) :: a integer(spi),dimension(size(a),size(a)):: b @@ -406,9 +422,10 @@ function diagn_i(a)result(b)! [diag] n=size(a) b=0; do i=1,n; b(i,i)=a(i); enddo end function diagn_i -!============================================================================= + +!> Calculation for diagnn_s function +!! @author R. J. Purser function diagnn_s(b)result(a)! [diag] -!============================================================================= implicit none real(sp),dimension(:,:),intent(IN ):: b real(sp),dimension(size(b,1)) :: a @@ -416,9 +433,10 @@ function diagnn_s(b)result(a)! [diag] n=size(b,1) do i=1,n; a(i)=b(i,i); enddo end function diagnn_s -!============================================================================= + +!> Calculation for diagnn_d function +!! @author R. J. Purser function diagnn_d(b)result(a)! [diag] -!============================================================================= implicit none real(dp),dimension(:,:),intent(IN ):: b real(dp),dimension(size(b,1)) :: a @@ -426,9 +444,10 @@ function diagnn_d(b)result(a)! [diag] n=size(b,1) do i=1,n; a(i)=b(i,i); enddo end function diagnn_d -!============================================================================= + +!> Calculation for diagnn_i function +!! @author R. J. Purser function diagnn_i(b)result(a)! [diag] -!============================================================================= implicit none integer(spi),dimension(:,:),intent(IN ):: b integer(spi),dimension(size(b,1)) :: a @@ -437,93 +456,92 @@ function diagnn_i(b)result(a)! [diag] do i=1,n; a(i)=b(i,i); enddo end function diagnn_i -!============================================================================= +!> Calculation for trace_s function +!! @author R. J. Purser function trace_s(b)result(s)! [trace] -!============================================================================= implicit none real(sp),dimension(:,:),intent(IN ):: b real(sp) :: s s=sum(diag(b)) end function trace_s -!============================================================================= + +!> Calculation for trace_d function +!! @author R. J. Purser function trace_d(b)result(s)! [trace] -!============================================================================= implicit none real(dp),dimension(:,:),intent(IN ):: b real(dp) :: s s=sum(diag(b)) end function trace_d -!============================================================================= + +!> Calculation for trace_i function +!! @author R. J. Purser function trace_i(b)result(s)! [trace] -!============================================================================= implicit none integer(spi),dimension(:,:),intent(IN ):: b integer(spi) :: s s=sum(diag(b)) end function trace_i -!============================================================================= +!> Calculation for identity_i function +!! @author R. J. Purser function identity_i(n)result(a)! [identity] -!============================================================================= implicit none integer(spi),intent(IN ) :: n integer(spi),dimension(n,n):: a integer(spi) :: i a=0; do i=1,n; a(i,i)=1; enddo end function identity_i -!============================================================================= + +!> Calculation for identity3_i function +!! @author R. J. Purser function identity3_i()result(a)! [identity] -!============================================================================= implicit none integer(spi),dimension(3,3):: a integer(spi) :: i a=0; do i=1,3; a(i,i)=1; enddo end function identity3_i -!============================================================================= +!> Spherical area of right-angle triangle whose orthogonal sides have +!! orthographic projection dimensions, sa and sb. +!! @author R. J. Purser function huarea_s(sa,sb)result(area)! [huarea] -!============================================================================= -! Spherical area of right-angle triangle whose orthogonal sides have -! orthographic projection dimensions, sa and sb. -!============================================================================= implicit none real(sp),intent(IN ):: sa,sb real(sp) :: area real(sp) :: ca,cb -!----------------------------------------------------------------------------- ca=sqrt(1-sa**2) cb=sqrt(1-sb**2) area=asin(sa*sb/(1+ca*cb)) end function huarea_s -!============================================================================= + +!> Calculation for huarea_d function +!! @author R. J. Purser function huarea_d(sa,sb)result(area)! [huarea] -!============================================================================= implicit none real(dp),intent(IN ):: sa,sb real(dp) :: area real(dp) :: ca,cb -!----------------------------------------------------------------------------- ca=sqrt(1-sa**2) cb=sqrt(1-sb**2) area=asin(sa*sb/(1+ca*cb)) end function huarea_d -!============================================================================= +!> Compute the area of the spherical triangle, {v1,v2,v3}, measured in the +!! right-handed sense, by dropping a perpendicular to u0 on the longest side so +!! that the area becomes the sum of areas of the two simpler right-angled +!! triangles. +!! @param v1 area of the spherical triangle +!! @param v2 area of the spherical triangle +!! @param v3 area of the spherical triangle +!! @author R. J. Purser function sarea_s(v1,v2,v3)result(area)! [sarea] -!============================================================================= -! Compute the area of the spherical triangle, {v1,v2,v3}, measured in the -! right-handed sense, by dropping a perpendicular to u0 on the longest side so -! that the area becomes the sum of areas of the two simpler right-angled -! triangles. -!============================================================================= use pietc_s, only: zero=>u0 implicit none real(sp),dimension(3),intent(IN ):: v1,v2,v3 real(sp) :: area -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - real(sp) :: s123,a1,a2,b,d1,d2,d3 real(sp),dimension(3) :: u0,u1,u2,u3,x,y -!============================================================================= area=zero u1=normalized(v1); u2=normalized(v2); u3=normalized(v3) s123=triple_product(u1,u2,u3) @@ -544,9 +562,16 @@ function sarea_s(v1,v2,v3)result(area)! [sarea] area=huarea(a1,b)+huarea(a2,b) contains -!----------------------------------------------------------------------------- + +!> Shift switch variable +!! @param u1 real variable to be shifted +!! @param u2 real variable to be shifted +!! @param u3 real variable to be shifted +!! @param d1 real variable to be shifted +!! @param d2 real variable to be shifted +!! @param d3 real variable to be shifted +!! @author R. J. Purser subroutine cyclic(u1,u2,u3,d1,d2,d3) -!----------------------------------------------------------------------------- implicit none real(sp),dimension(3),intent(INOUT):: u1,u2,u3 real(sp), intent(INOUT):: d1,d2,d3 @@ -556,17 +581,19 @@ subroutine cyclic(u1,u2,u3,d1,d2,d3) ut=u1; u1=u2; u2=u3; u3=ut end subroutine cyclic end function sarea_s -!============================================================================= + +!> Compute the area of sarea_d, {v1,v2,v3} +!! @param v1 area of the spherical triangle +!! @param v2 area of the spherical triangle +!! @param v3 area of the spherical triangle +!! @author R. J. Purser function sarea_d(v1,v2,v3)result(area)! [sarea] -!============================================================================= use pietc, only: zero=>u0 implicit none real(dp),dimension(3),intent(IN ):: v1,v2,v3 real(dp) :: area -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - real(dp) :: s123,a1,a2,b,d1,d2,d3 real(dp),dimension(3) :: u0,u1,u2,u3,x,y -!============================================================================= area=zero u1=normalized(v1); u2=normalized(v2); u3=normalized(v3) s123=triple_product(u1,u2,u3) @@ -587,9 +614,16 @@ function sarea_d(v1,v2,v3)result(area)! [sarea] area=huarea(a1,b)+huarea(a2,b) contains -!----------------------------------------------------------------------------- + +!> Shift switch variable +!! @param u1 real variable to be shifted +!! @param u2 real variable to be shifted +!! @param u3 real variable to be shifted +!! @param d1 real variable to be shifted +!! @param d2 real variable to be shifted +!! @param d3 real variable to be shifted +!! @author R. J. Purser subroutine cyclic(u1,u2,u3,d1,d2,d3) -!----------------------------------------------------------------------------- implicit none real(dp),dimension(3),intent(INOUT):: u1,u2,u3 real(dp), intent(INOUT):: d1,d2,d3 @@ -599,24 +633,26 @@ subroutine cyclic(u1,u2,u3,d1,d2,d3) ut=u1; u1=u2; u2=u3; u3=ut end subroutine cyclic end function sarea_d -!============================================================================= + +!> Compute the area of the spherical triangle with a vertex at latitude +!! rlat, and two other vertices, A and B, whose incremented latitudes +!! and longitudes are drlata,drlona (for A) and drlatb,drlonb (for B). +!! The computations are designed to give a proportionately accurate area +!! estimate even when the triangle is very small, provided the B-increment +!! is not disproportionately small compared to the other two sides. +!! @param rlat latitude +!! @param drlata latitudes +!! @param drlona longitudes +!! @param drlatb latitudes +!! @param drlonb longitudes +!! @author R. J. Purser function dtarea_s(rlat,drlata,drlona,drlatb,drlonb) result(area)! [sarea] -!============================================================================= -! Compute the area of the spherical triangle with a vertex at latitude -! rlat, and two other vertices, A and B, whose incremented latitudes -! and longitudes are drlata,drlona (for A) and drlatb,drlonb (for B). -! The computations are designed to give a proportionately accurate area -! estimate even when the triangle is very small, provided the B-increment -! is not disproportionately small compared to the other two sides. -!============================================================================= use pietc_s, only: u0,u1 implicit none real(sp),intent(in ):: rlat,drlata,drlona,drlatb,drlonb real(sp) :: area -!----------------------------------------------------------------------------- real(sp),dimension(2):: x2a,x2b,xb,yb real(sp) :: sb,ssb,cb,xa,sa,ca,sc,cc -!============================================================================= call dlltoxy(rlat,drlata,drlona,x2a) call dlltoxy(rlat,drlatb,drlonb,x2b) ssb=dot_product(x2b,x2b); sb=sqrt(ssb) @@ -633,17 +669,21 @@ function dtarea_s(rlat,drlata,drlona,drlatb,drlonb) result(area)! [sarea] sb=sb*cc-cb*sc area=huarea(-sa,sb)+huarea(sc,-sa) end function dtarea_s -!============================================================================= + +!> Compute the area with dtarea_d +!! @param rlat latitude +!! @param drlata latitudes +!! @param drlona longitudes +!! @param drlatb latitudes +!! @param drlonb longitudes +!! @author R. J. Purser function dtarea_d(rlat,drlata,drlona,drlatb,drlonb) result(area)! [sarea] -!============================================================================= use pietc, only: u0,u1 implicit none real(dp),intent(in ):: rlat,drlata,drlona,drlatb,drlonb real(dp) :: area -!----------------------------------------------------------------------------- real(dp),dimension(2):: x2a,x2b,xb,yb real(dp) :: sb,ssb,cb,xa,sa,ca,sc,cc -!============================================================================= call dlltoxy(rlat,drlata,drlona,x2a) call dlltoxy(rlat,drlatb,drlonb,x2b) ssb=dot_product(x2b,x2b); sb=sqrt(ssb) @@ -660,78 +700,93 @@ function dtarea_d(rlat,drlata,drlona,drlatb,drlonb) result(area)! [sarea] sb=sb*cc-cb*sc area=huarea(-sa,sb)+huarea(sc,-sa) end function dtarea_d -!============================================================================= + +!> Compute the area of the spherical quadrilateral with a vertex at latitude +!! rlat, and three other vertices at A, B, and C inturn, +!! whose incremented latitudes and longitudes are drlata,drlona (for A), +!! drlatb,drlonb (for B), and drlatc,drlonc (for C). +!! The computations are designed to give a proportionately accurate area +!! estimate even when the quadrilateral is very small, provided the +!! diagonal making the B-increment is not disproportionately small compared to +!! the characteristic size of the quadrilateral. +!! @param rlat latitude +!! @param drlata latitudes +!! @param drlona longitudes +!! @param drlatb latitudes +!! @param drlonb longitudes +!! @param drlatc latitudes +!! @param drlonc longitudes +!! @author R. J. Purser function dqarea_s &! [sarea] (rlat,drlata,drlona,drlatb,drlonb,drlatc,drlonc) result(area) -!============================================================================= -! Compute the area of the spherical quadrilateral with a vertex at latitude -! rlat, and three other vertices at A, B, and C inturn, -! whose incremented latitudes and longitudes are drlata,drlona (for A), -! drlatb,drlonb (for B), and drlatc,drlonc (for C). -! The computations are designed to give a proportionately accurate area -! estimate even when the quadrilateral is very small, provided the -! diagonal making the B-increment is not disproportionately small compared to -! the characteristic size of the quadrilateral. -!============================================================================= implicit none real(sp),intent(in ):: rlat,drlata,drlona,drlatb,drlonb,drlatc,drlonc real(sp) :: area -!============================================================================= area=sarea(rlat,drlata,drlona,drlatb,drlonb)& -sarea(rlat,drlatc,drlonc,drlatb,drlonb) end function dqarea_s -!============================================================================= + +!> Compute the area using dqarea_d +!! @param rlat latitude +!! @param drlata latitudes +!! @param drlona longitudes +!! @param drlatb latitudes +!! @param drlonb longitudes +!! @param drlatc latitudes +!! @param drlonc longitudes +!! @author R. J. Purser function dqarea_d &! [sarea] (rlat,drlata,drlona,drlatb,drlonb,drlatc,drlonc) result(area) -!============================================================================= implicit none real(dp),intent(in ):: rlat,drlata,drlona,drlatb,drlonb,drlatc,drlonc real(dp) :: area -!============================================================================= area=sarea(rlat,drlata,drlona,drlatb,drlonb)& -sarea(rlat,drlatc,drlonc,drlatb,drlonb) end function dqarea_d -!============================================================================= +!> Calculate dlltoxy_s +!! @param rlat latitude +!! @param drlat latitude +!! @param drlon longitudes +!! @author R. J. Purser subroutine dlltoxy_s(rlat,drlat,drlon,x2)! [dlltoxy] -!============================================================================= use pietc_s, only: u2 implicit none real(sp), intent(in ):: rlat,drlat,drlon real(sp),dimension(2),intent(out):: x2 -!----------------------------------------------------------------------------- real(sp):: clata -!============================================================================= clata=cos(rlat+drlat) x2=(/clata*sin(drlon),sin(drlat)+u2*sin(rlat)*clata*hav(drlon)/) end subroutine dlltoxy_s -!============================================================================= + +!> Calculate dlltoxy_d +!! @param rlat latitude +!! @param drlat latitude +!! @param drlon longitudes +!! @author R. J. Purser subroutine dlltoxy_d(rlat,drlat,drlon,x2)! [dlltoxy] -!============================================================================= use pietc, only: u2 implicit none real(dp), intent(in ):: rlat,drlat,drlon real(dp),dimension(2),intent(out):: x2 -!----------------------------------------------------------------------------- real(dp):: clata -!============================================================================= clata=cos(rlat+drlat) x2=(/clata*sin(drlon),sin(drlat)+u2*sin(rlat)*clata*hav(drlon)/) end subroutine dlltoxy_d -!============================================================================= +!> Haversine function +!! @author R. J. Purser function hav_s(t) result(a)! [hav] -!============================================================================= -! Haversine function use pietc_s, only: o2 implicit none real(sp),intent(in ):: t real(sp) :: a a=(sin(t*o2))**2 end function hav_s -!============================================================================= + +!> hav_d function +!! @author R. J. Purser function hav_d(t) result(a)! [hav] -!============================================================================= use pietc, only: o2 implicit none real(dp),intent(in ):: t @@ -739,19 +794,21 @@ function hav_d(t) result(a)! [hav] a=(sin(t*o2))**2 end function hav_d -!============================================================================= +!> Normalize the given vector +!! @param v vector +!! @author R. J. Purser subroutine normalize_s(v)! [normalize] -!============================================================================= -! Normalize the given vector. use pietc_s, only: u0,u1 implicit none real(sp),dimension(:),intent(inout):: v real(sp) :: s s=absv(v); if(s==0)then; v=u0; v(1)=u1; else; v=v/s; endif end subroutine normalize_s -!============================================================================= + +!> normalize_d calculation for given vector +!! @param v vector +!! @author R. J. Purser subroutine normalize_d(v)! [normalize] -!============================================================================= use pietc, only: u0,u1 implicit none real(dp),dimension(:),intent(inout):: v @@ -759,16 +816,15 @@ subroutine normalize_d(v)! [normalize] s=absv(v); if(s==u0)then; v=0; v(1)=u1; else; v=v/s; endif end subroutine normalize_d -!============================================================================= +!> gram_s routine +!! @author R. J. Purser subroutine gram_s(as,b,nrank,det)! [gram] -!============================================================================= use pietc_s, only: u0,u1 implicit none real(sp),dimension(:,:),intent(IN ) :: as real(sp),dimension(:,:),intent(OUT) :: b integer(spi), intent(OUT) :: nrank real(sp), intent(OUT) :: det -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - real(sp),parameter :: crit=1.e-5_sp real(sp),dimension(size(as,1),size(as,2)):: a real(sp),dimension(size(as,2),size(as,1)):: ab @@ -776,7 +832,6 @@ subroutine gram_s(as,b,nrank,det)! [gram] real(sp) :: val,s,vcrit integer(spi) :: i,j,k,l,m,n integer(spi),dimension(2) :: ii -!============================================================================= n=size(as,1) m=size(as,2) if(n/=size(b,1) .or. n/=size(b,2))stop 'In gram; incompatible dimensions' @@ -822,17 +877,16 @@ subroutine gram_s(as,b,nrank,det)! [gram] enddo enddo end subroutine gram_s - -!============================================================================= + +!> gram_d routine +!! @author R. J. Purser subroutine gram_d(as,b,nrank,det)! [gram] -!============================================================================= use pietc, only: u0,u1 implicit none real(dp),dimension(:,:),intent(IN ) :: as real(dp),dimension(:,:),intent(OUT) :: b integer(spi), intent(OUT) :: nrank real(dp), intent(OUT) :: det -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - real(dp),parameter :: crit=1.e-9_dp real(dp),dimension(size(as,1),size(as,2)):: a real(dp),dimension(size(as,2),size(as,1)):: ab @@ -840,7 +894,6 @@ subroutine gram_d(as,b,nrank,det)! [gram] real(dp) :: val,s,vcrit integer(spi) :: i,j,k,l,m,n integer(spi),dimension(2) :: ii -!============================================================================= n=size(as,1) m=size(as,2) if(n/=size(b,1) .or. n/=size(b,2))stop 'In gram; incompatible dimensions' @@ -887,15 +940,14 @@ subroutine gram_d(as,b,nrank,det)! [gram] enddo end subroutine gram_d -!============================================================================= +!> A version of gram_d where the determinant information is returned in +!! logarithmic form (to avoid overflows for large matrices). When the +!! matrix is singular, the "sign" of the determinant, detsign, is returned +!! as zero (instead of either +1 or -1) and ldet is then just the log of +!! the nonzero factors found by the process. +!! @param detsign singular determinant +!! @author R. J. Purser subroutine graml_d(as,b,nrank,detsign,ldet)! [gram] -!============================================================================= -! A version of gram_d where the determinant information is returned in -! logarithmic form (to avoid overflows for large matrices). When the -! matrix is singular, the "sign" of the determinant, detsign, is returned -! as zero (instead of either +1 or -1) and ldet is then just the log of -! the nonzero factors found by the process. -!============================================================================= use pietc, only: u0 implicit none real(dp),dimension(:,:),intent(IN ) :: as @@ -903,7 +955,6 @@ subroutine graml_d(as,b,nrank,detsign,ldet)! [gram] integer(spi), intent(OUT) :: nrank integer(spi), intent(out) :: detsign real(dp), intent(OUT) :: ldet -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - real(dp),parameter :: crit=1.e-9_dp real(dp),dimension(size(as,1),size(as,2)):: a real(dp),dimension(size(as,2),size(as,1)):: ab @@ -911,14 +962,12 @@ subroutine graml_d(as,b,nrank,detsign,ldet)! [gram] real(dp) :: val,s,vcrit integer(spi) :: i,j,k,l,m,n integer(spi),dimension(2) :: ii -!============================================================================= detsign=1 n=size(as,1) m=size(as,2) if(n/=size(b,1) .or. n/=size(b,2))stop 'In gram; incompatible dimensions' a=as b=identity(n) - ldet=u0 val=maxval(abs(a)) if(val==u0)then @@ -967,20 +1016,18 @@ subroutine graml_d(as,b,nrank,detsign,ldet)! [gram] enddo enddo end subroutine graml_d - -!============================================================================= + + +!> A "plain" (unpivoted) version of Gram-Schmidt, for square matrices only. +!! @author R. J. Purser subroutine plaingram_s(b,nrank)! [gram] -!============================================================================= -! A "plain" (unpivoted) version of Gram-Schmidt, for square matrices only. use pietc_s, only: u0 implicit none real(sp),dimension(:,:),intent(INOUT) :: b integer(spi), intent( OUT) :: nrank -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - real(sp),parameter :: crit=1.e-5_sp real(sp) :: val,vcrit integer(spi) :: j,k,n -!============================================================================= n=size(b,1); if(n/=size(b,2))stop 'In gram; matrix needs to be square' val=maxval(abs(b)) nrank=0 @@ -1003,19 +1050,16 @@ subroutine plaingram_s(b,nrank)! [gram] enddo end subroutine plaingram_s -!============================================================================= +!> A "plain" (unpivoted) version of Gram-Schmidt, for square matrices only. +!! @author R. J. Purser subroutine plaingram_d(b,nrank)! [gram] -!============================================================================= -! A "plain" (unpivoted) version of Gram-Schmidt, for square matrices only. use pietc, only: u0 implicit none real(dp),dimension(:,:),intent(INOUT):: b integer(spi), intent( OUT):: nrank -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - real(dp),parameter:: crit=1.e-9_dp real(dp) :: val,vcrit integer(spi) :: j,k,n -!============================================================================= n=size(b,1); if(n/=size(b,2))stop 'In gram; matrix needs to be square' val=maxval(abs(b)) nrank=0 @@ -1038,23 +1082,25 @@ subroutine plaingram_d(b,nrank)! [gram] enddo end subroutine plaingram_d -!============================================================================= +!> Without changing (tall) rectangular input matrix a, perform pivoted gram- +!! Schmidt operations to orthogonalize the rows, until rows that remain become +!! negligible. Record the pivoting sequence in ipiv, and the row-normalization +!! in tt(j,j) and the row-orthogonalization in tt(i,j), for i>j. Note that +!! tt(i,j)=0 for ij. Note that -! tt(i,j)=0 for i=n please' nepss=n*epss rank=n aa=a tt=u0 do ii=1,n - ! At this stage, all rows less than ii are already orthonormalized and are ! orthogonal to all rows at and beyond ii. Find the norms of these lower ! rows and pivot the largest of them into position ii: @@ -1113,7 +1156,6 @@ subroutine rowgram(m,n,a,ipiv,tt,b,rank)! [gram] rank=ii-1 return endif - ipiv(ii)=maxi if(maxi/=ii)then rowv =aa(ii, :) @@ -1133,22 +1175,20 @@ subroutine rowgram(m,n,a,ipiv,tt,b,rank)! [gram] b=aa(1:n,:) end subroutine rowgram -!============================================================================= +!> Apply the row-operations, implied by ipiv and tt returned by rowgram, to +!! the single column vector, v, to produce the transformed vector vv. +!! @param v vector +!! @param vv vector +!! @author R. J. Purser subroutine rowops(m,n,ipiv,tt,v,vv)! [rowops] -!============================================================================= -! Apply the row-operations, implied by ipiv and tt returned by rowgram, to -! the single column vector, v, to produce the transformed vector vv. -!============================================================================= implicit none integer(spi), intent(in ):: m,n integer(spi),dimension(n),intent(in ):: ipiv real(dp),dimension(m,n), intent(in ):: tt real(dp),dimension(m), intent(in ):: v real(dp),dimension(m), intent(out):: vv -!----------------------------------------------------------------------------- integer(spi):: i,j,k real(dp) :: p -!============================================================================= vv=v do j=1,n k=ipiv(j) @@ -1163,19 +1203,21 @@ subroutine rowops(m,n,ipiv,tt,v,vv)! [rowops] enddo enddo end subroutine rowops - -!============================================================================= + +!> Find positive diagonals D and E and a Lagrange multiplier F that minimize +!! the row-sum +column-sum of masked terms, +!! (D_i +log(|A_ij|) +E_j)^2 +!! subject to the single constraint, sum_j E_j =0, where the mask permits +!! only nonnegligible A_ij to participate in the quadratic quantities. +!! Once a solution for D and E is found, return their exponentials, d and e, +!! together with the rescaled matrix aa such that a = d.aa.e when d and e are +!! interpreted as diagonal matrices. +!! @param d positive diagonals +!! @param e positive diagonals +!! @param a positive diagonals +!! @param f Lagrange multiplier +!! @author R. J. Purser subroutine corral(m,n,mask,a,d,aa,e)! [corral] -!============================================================================= -! Find positive diagonals D and E and a Lagrange multiplier F that minimize -! the row-sum +column-sum of masked terms, -! (D_i +log(|A_ij|) +E_j)^2 -! subject to the single constraint, sum_j E_j =0, where the mask permits -! only nonnegligible A_ij to participate in the quadratic quantities. -! Once a solution for D and E is found, return their exponentials, d and e, -! together with the rescaled matrix aa such that a = d.aa.e when d and e are -! interpreted as diagonal matrices. -!============================================================================= use pietc, only: u0,u1 use pmat, only: inv implicit none @@ -1185,11 +1227,9 @@ subroutine corral(m,n,mask,a,d,aa,e)! [corral] real(dp),dimension(m ),intent(out):: d real(dp),dimension(m,n),intent(out):: aa real(dp),dimension( n),intent(out):: e -!----------------------------------------------------------------------------- real(dp),dimension(0:m+n,0:m+n):: g real(dp),dimension(0:m+n) :: h integer(spi) :: i,j,k,nh -!============================================================================= nh=1+m+n aa=u0 do j=1,n @@ -1197,16 +1237,13 @@ subroutine corral(m,n,mask,a,d,aa,e)! [corral] if(mask(i,j))aa(i,j)=log(abs(a(i,j))) enddo enddo - h=u0 g=u0 - ! Equations on row 0 enforcing Lagrange multiplier F-constraint: do j=1,n k=m+j g(0,k)=u1 enddo - ! Equations on rows 1:m minimizing row sums of quadratic terms: do i=1,m do j=1,n @@ -1218,7 +1255,6 @@ subroutine corral(m,n,mask,a,d,aa,e)! [corral] endif enddo enddo - ! Equations on rows m+1:m+n minimizing col sums subject to constraint do j=1,n k=m+j @@ -1231,10 +1267,8 @@ subroutine corral(m,n,mask,a,d,aa,e)! [corral] endif enddo enddo - ! Invert the normal equations: call inv(g,h) - ! Exponentiate the parts that become final scaling diagnonal matrices d and e: do i=1,m d(i)=exp(h(i)) @@ -1243,7 +1277,6 @@ subroutine corral(m,n,mask,a,d,aa,e)! [corral] k=m+j e(j)=exp(h(k)) enddo - ! Compute the rescaled matrix directly: do j=1,n do i=1,m @@ -1252,25 +1285,23 @@ subroutine corral(m,n,mask,a,d,aa,e)! [corral] enddo end subroutine corral -!============================================================================= +!> Assuming that given orth33 is a 3*3 proper rotation matrix, derive an axial +!! 3-vector, ax3, such that orth33 is implied by ax3 when the latter is +!! interpreted as encoding a rotation (as in subroutine axtorot). Note that +!! such ax3 are not unique -- adding any multiple of 2*pi times the parallel +!! unit vector leads to the same orth33. +!! @param orth33 3*3 proper rotation matrix +!! @param ax3 axial 3-vector +!! @author R. J. Purser subroutine rottoax(orth33,ax3)! [rottoax] -!============================================================================= -! Assuming that given orth33 is a 3*3 proper rotation matrix, derive an axial -! 3-vector, ax3, such that orth33 is implied by ax3 when the latter is -! interpreted as encoding a rotation (as in subroutine axtorot). Note that -! such ax3 are not unique -- adding any multiple of 2*pi times the parallel -! unit vector leads to the same orth33. -!============================================================================= implicit none real(dp),dimension(3,3),intent(in ):: orth33 real(dp),dimension(3), intent(out):: ax3 -!----------------------------------------------------------------------------- real(dp),dimension(3,3) :: plane real(dp),dimension(3) :: x,y,z real(dp) :: s integer(spi),dimension(1):: ii integer(spi) :: i,j,k -!============================================================================= plane=orth33-identity()! Columns must be coplanar vectors do i=1,3; z(i)=dot_product(plane(:,i),plane(:,i)); enddo ii=minloc(z) @@ -1285,67 +1316,61 @@ subroutine rottoax(orth33,ax3)! [rottoax] ax3=ax3*atan2(dot_product(y,z),dot_product(x,z))! multiply ax3 by the angle end subroutine rottoax -!============================================================================= +!> Construct the 3*3 orthogonal matrix, orth33, that corresponds to the +!! proper rotation encoded by the 3-vector, ax3. The antisymmetric matrix +!! ax33 equivalent to the axial vector ax3 is exponentiated to obtain orth33. +!! @param orth33 3*3 orthogonal matrix +!! @param ax3 axial 3-vector +!! @author R. J. Purser subroutine axtorot(ax3,orth33)! [axtorot] -!============================================================================= -! Construct the 3*3 orthogonal matrix, orth33, that corresponds to the -! proper rotation encoded by the 3-vector, ax3. The antisymmetric matrix -! ax33 equivalent to the axial vector ax3 is exponentiated to obtain orth33. -!============================================================================= implicit none real(dp),dimension(3), intent(in ):: ax3 real(dp),dimension(3,3),intent(out):: orth33 -!----------------------------------------------------------------------------- real(dp),dimension(3,3):: ax33 real(dp) :: d -!============================================================================= ax33=axial(ax3); call expmat(3,ax33,orth33,d) end subroutine axtorot -!============================================================================= +!> Go from the spinor to the quaternion representation. +!! @param cspin spinor representation +!! @param q quaternion representation +!! @author R. J. Purser subroutine spintoq(cspin,q)! [spintoq] -!============================================================================= -! Go from the spinor to the quaternion representation -!============================================================================= implicit none complex(dpc),dimension(2,2),intent(IN ):: cspin real(dp), dimension(0:3),intent(OUT):: q -!------------------------------------------- q(0)=real(cspin(1,1)); q(3)=aimag(cspin(1,1)) q(2)=real(cspin(2,1)); q(1)=aimag(cspin(2,1)) end subroutine spintoq -!============================================================================== +!> Go from the quaternion to the spinor representation +!! @param cspin spinor representation +!! @param q quaternion representation +!! @author R. J. Purser subroutine qtospin(q,cspin)! [qtospin] -!============================================================================== -! Go from the quaternion to the spinor representation -!============================================================================== implicit none real(dp), dimension(0:3),intent(IN ):: q complex(dpc),dimension(2,2),intent(OUT):: cspin -!------------------------------------------- cspin(1,1)=cmplx( q(0), q(3)) cspin(2,1)=cmplx( q(2), q(1)) cspin(1,2)=cmplx(-q(2), q(1)) cspin(2,2)=cmplx( q(0),-q(3)) end subroutine qtospin -!============================================================================= +!> Go from rotation matrix to quaternion representation +!! @param rot rotation matrix +!! @param q quaternion representation +!! @author R. J. Purser subroutine rottoq(rot,q)! [rottoq] -!============================================================================= -! Go from rotation matrix to quaternion representation -!============================================================================= use pietc, only: zero=>u0,one=>u1,two=>u2 implicit none real(dp),dimension(3,3),intent(IN ):: rot real(dp),dimension(0:3),intent(OUT):: q -!------------------------------------------------------------------------------ real(dp),dimension(3,3) :: t1,t2 real(dp),dimension(3) :: u1,u2 real(dp) :: gamma,gammah,s,ss integer(spi) :: i,j integer(spi),dimension(1):: ii -!============================================================================== ! construct the orthogonal matrix, t1, whose third row is the rotation axis ! of rot: t1=rot; do i=1,3; t1(i,i)=t1(i,i)-1; u1(i)=dot_product(t1(i,:),t1(i,:)); enddo @@ -1372,18 +1397,14 @@ subroutine rottoq(rot,q)! [rottoq] if(ss==zero)stop 'In rotov; invalid rot' if(j/=2)t1(2,:)=t1(3,:) t1(2,:)=t1(2,:)/sqrt(ss) - ! Form t1(3,:) as the cross product of t1(1,:) and t1(2,:) t1(3,1)=t1(1,2)*t1(2,3)-t1(1,3)*t1(2,2) t1(3,2)=t1(1,3)*t1(2,1)-t1(1,1)*t1(2,3) t1(3,3)=t1(1,1)*t1(2,2)-t1(1,2)*t1(2,1) - ! Project rot into the frame whose axes are the rows of t1: t2=matmul(t1,matmul(rot,transpose(t1))) - ! Obtain the rotation angle, gamma, implied by rot, and gammah=gamma/2: gamma=atan2(t2(2,1),t2(1,1)); gammah=gamma/two - ! Hence deduce coefficients (in the form of a real 4-vector) of one of the two ! possible equivalent spinors: s=sin(gammah) @@ -1391,57 +1412,50 @@ subroutine rottoq(rot,q)! [rottoq] q(1:3)=t1(3,:)*s end subroutine rottoq -!============================================================================== +!> Go from quaternion to rotation matrix representations +!! @param q quaternion representation +!! @param rot rotation matrix representations +!! @author R. J. Purser subroutine qtorot(q,rot)! [qtorot] -!============================================================================== -! Go from quaternion to rotation matrix representations -!============================================================================== implicit none real(dp),dimension(0:3),intent(IN ):: q real(dp),dimension(3,3),intent(OUT):: rot -!============================================================================= call setem(q(0),q(1),q(2),q(3),rot) end subroutine qtorot -!============================================================================= +!> Go from an axial 3-vector to its equivalent quaternion +!! @param q quaternion +!! @param v axial 3-vector +!! @author R. J. Purser subroutine axtoq(v,q)! [axtoq] -!============================================================================= -! Go from an axial 3-vector to its equivalent quaternion -!============================================================================= implicit none real(dp),dimension(3), intent(in ):: v real(dp),dimension(0:3),intent(out):: q -!----------------------------------------------------------------------------- real(dp),dimension(3,3):: rot -!============================================================================= call axtorot(v,rot) call rottoq(rot,q) end subroutine axtoq -!============================================================================= +!> Go from quaternion to axial 3-vector +!! @param q quaternion +!! @param v axial 3-vector +!! @author R. J. Purser subroutine qtoax(q,v)! [qtoax] -!============================================================================= -! Go from quaternion to axial 3-vector -!============================================================================= implicit none real(dp),dimension(0:3),intent(in ):: q real(dp),dimension(3), intent(out):: v -!----------------------------------------------------------------------------- real(dp),dimension(3,3):: rot -!============================================================================= call qtorot(q,rot) call rottoax(rot,v) end subroutine qtoax -!============================================================================= +!> setem routine +!! @author R. J. Purser subroutine setem(c,d,e,g,r)! [setem] -!============================================================================= implicit none real(dp), intent(IN ):: c,d,e,g real(dp),dimension(3,3),intent(OUT):: r -!----------------------------------------------------------------------------- real(dp):: cc,dd,ee,gg,de,dg,eg,dc,ec,gc -!============================================================================= cc=c*c; dd=d*d; ee=e*e; gg=g*g de=d*e; dg=d*g; eg=e*g dc=d*c; ec=e*c; gc=g*c @@ -1450,40 +1464,39 @@ subroutine setem(c,d,e,g,r)! [setem] r(3,2)=2*(eg+dc); r(1,3)=2*(dg+ec); r(2,1)=2*(de+gc) end subroutine setem -!============================================================================= +!> Multiply quaternions, a*b, assuming operation performed from right to left +!! @param a real quaternion +!! @param b real quaternion +!! @author R. J. Purser function mulqq(a,b)result(c)! [mulqq] -!============================================================================= -! Multiply quaternions, a*b, assuming operation performed from right to left -!============================================================================= implicit none real(dp),dimension(0:3),intent(IN ):: a,b real(dp),dimension(0:3) :: c -!------------------------------------------- c(0)=a(0)*b(0) -a(1)*b(1) -a(2)*b(2) -a(3)*b(3) c(1)=a(0)*b(1) +a(1)*b(0) +a(2)*b(3) -a(3)*b(2) c(2)=a(0)*b(2) +a(2)*b(0) +a(3)*b(1) -a(1)*b(3) c(3)=a(0)*b(3) +a(3)*b(0) +a(1)*b(2) -a(2)*b(1) end function mulqq -!============================================================================= + +!> Evaluate the exponential, b, of a matrix, a, of degree n. +!! Apply the iterated squaring method, m times, to the approximation to +!! exp(a/(2**m)) obtained as a Taylor expansion of degree L +!! See Fung, T. C., 2004, Int. J. Numer. Meth. Engng, 59, 1273--1286. +!! @param b exponential +!! @param a matrix +!! @param n degree +!! @author R. J. Purser subroutine expmat(n,a,b,detb)! [expmat] -!============================================================================= -! Evaluate the exponential, b, of a matrix, a, of degree n. -! Apply the iterated squaring method, m times, to the approximation to -! exp(a/(2**m)) obtained as a Taylor expansion of degree L -! See Fung, T. C., 2004, Int. J. Numer. Meth. Engng, 59, 1273--1286. -!============================================================================= use pietc, only: u0,u1,u2,o2 implicit none integer(spi), intent(IN ):: n real(dp),dimension(n,n),intent(IN ):: a real(dp),dimension(n,n),intent(OUT):: b real(dp), intent(OUT):: detb -!----------------------------------------------------------------------------- integer(spi),parameter :: L=5 real(dp),dimension(n,n):: c,p real(dp) :: t integer(spi) :: i,m -!============================================================================= m=10+floor(log(u1+maxval(abs(a)))/log(u2)) t=o2**m c=a*t @@ -1502,11 +1515,12 @@ subroutine expmat(n,a,b,detb)! [expmat] detb=u0; do i=1,n; detb=detb+a(i,i); enddo; detb=exp(detb) end subroutine expmat -!============================================================================= +!> Like expmat, but for the 1st derivatives also. +!! @param b exponential +!! @param a matrix +!! @param n degree +!! @author R. J. Purser subroutine expmatd(n,a,b,bd,detb,detbd)! [expmat] -!============================================================================= -! Like expmat, but for the 1st derivatives also. -!============================================================================= use pietc, only: u0,u1,u2,o2 implicit none integer(spi), intent(IN ):: n @@ -1515,13 +1529,11 @@ subroutine expmatd(n,a,b,bd,detb,detbd)! [expmat] real(dp),dimension(n,n,(n*(n+1))/2),intent(OUT):: bd real(dp), intent(OUT):: detb real(dp),dimension((n*(n+1))/2), intent(OUT):: detbd -!----------------------------------------------------------------------------- integer(spi),parameter :: L=5 real(dp),dimension(n,n) :: c,p real(dp),dimension(n,n,(n*(n+1))/2):: pd,cd real(dp) :: t integer(spi) :: i,j,k,m,n1 -!============================================================================= n1=(n*(n+1))*o2 m=10+floor(log(u1+maxval(abs(a)))/log(u2)) t=o2**m @@ -1543,7 +1555,6 @@ subroutine expmatd(n,a,b,bd,detb,detbd)! [expmat] cd=pd b=p bd=pd - do i=2,L do k=1,n1 pd(:,:,k)=(matmul(cd(:,:,k),p)+matmul(c,pd(:,:,k)))/i @@ -1565,11 +1576,12 @@ subroutine expmatd(n,a,b,bd,detb,detbd)! [expmat] detbd=u0; do k=1,n; detbd(k)=detb; enddo end subroutine expmatd -!============================================================================= +!> Like expmat, but for the 1st and 2nd derivatives also. +!! @param b exponential +!! @param a matrix +!! @param n degree +!! @author R. J. Purser subroutine expmatdd(n,a,b,bd,bdd,detb,detbd,detbdd)! [expmat] -!============================================================================= -! Like expmat, but for the 1st and 2nd derivatives also. -!============================================================================= use pietc, only: u0,u1,u2,o2 implicit none integer(spi), intent(IN ):: n @@ -1580,14 +1592,12 @@ subroutine expmatdd(n,a,b,bd,bdd,detb,detbd,detbdd)! [expmat] real(dp), intent(OUT):: detb real(dp),dimension((n*(n+1))/2), intent(OUT):: detbd real(dp),dimension((n*(n+1))/2,(n*(n+1))/2), intent(OUT):: detbdd -!----------------------------------------------------------------------------- integer(spi),parameter :: L=5 real(dp),dimension(n,n) :: c,p real(dp),dimension(n,n,(n*(n+1))/2) :: pd,cd real(dp),dimension(n,n,(n*(n+1))/2,(n*(n+1))/2):: pdd,cdd real(dp) :: t integer(spi) :: i,j,k,ki,kj,m,n1 -!============================================================================= n1=(n*(n+1))/2 m=10+floor(log(u1+maxval(abs(a)))/log(u2)) t=o2**m @@ -1612,7 +1622,6 @@ subroutine expmatdd(n,a,b,bd,bdd,detb,detbd,detbdd)! [expmat] b=p bd=pd bdd=u0 - do i=2,L do ki=1,n1 do kj=1,n1 @@ -1652,20 +1661,18 @@ subroutine expmatdd(n,a,b,bd,bdd,detb,detbd,detbdd)! [expmat] detbdd=u0; do ki=1,n; do kj=1,n; detbdd(ki,kj)=detb; enddo; enddo end subroutine expmatdd -!============================================================================= +!> routine zntay +!! @author R. J. Purser subroutine zntay(n,z,zn)! [zntay] -!============================================================================= use pietc, only: u2 implicit none integer(spi), intent(IN ):: n real(dp), intent(IN ):: z real(dp), intent(OUT):: zn -!----------------------------------------------------------------------------- integer(spi),parameter:: ni=100 real(dp),parameter :: eps0=1.e-16_dp integer(spi) :: i,i2,n2 real(dp) :: t,eps,z2 -!============================================================================= z2=z*u2 n2=n*2 t=1 @@ -1684,18 +1691,16 @@ subroutine zntay(n,z,zn)! [zntay] print'("In zntay; full complement of iterations used")' end subroutine zntay -!============================================================================= +!> routine znfun +!! @author R. J. Purser subroutine znfun(n,z,zn,znd,zndd,znddd)! [znfun] -!============================================================================= use pietc, only: u0,u2,u3 implicit none integer(spi),intent(IN ):: n real(dp), intent(IN ):: z real(dp), intent(OUT):: zn,znd,zndd,znddd -!----------------------------------------------------------------------------- real(dp) :: z2,rz2 integer(spi):: i,i2p3 -!============================================================================= z2=abs(z*u2) rz2=sqrt(z2) if(z2zv by the transformatn -! with coefficients aa3,bb3,cc3,dd3, such that, as 2*2 complex matrices: -! -! [ aa3, bb3 ] [ aa2, bb2 ] [ aa1, bb1 ] -! [ ] = [ ] * [ ] -! [ cc3, dd3 ] [ cc2, dd2 ] [ cc1, dd1 ] . -! -! Note that the determinant of these matrices is always +1 -! -!============================================================================= +!> Utility code for various Mobius transformations. If aa1,bb1,cc1,dd1 are +!! the coefficients for one transformation, and aa2,bb2,cc2,dd2 are the +!! coefficients for a second one, then the coefficients for the mapping +!! of a test point, zz, by aa1 etc to zw, followed by a mapping of zw, by +!! aa2 etc to zv, is equivalent to a single mapping zz-->zv by the transformatn +!! with coefficients aa3,bb3,cc3,dd3, such that, as 2*2 complex matrices: +!! +!! [ aa3, bb3 ] [ aa2, bb2 ] [ aa1, bb1 ] +!! [ ] = [ ] * [ ] +!! [ cc3, dd3 ] [ cc2, dd2 ] [ cc1, dd1 ] . +!! +!! Note that the determinant of these matrices is always +1 +!! @author R. J. Purser subroutine ctoz(v, z,infz)! [ctoz] -!============================================================================= use pietc, only: u0,u1 implicit none real(dp),dimension(3),intent(IN ):: v complex(dpc), intent(OUT):: z logical, intent(OUT):: infz -!----------------------------------------------------------------------------- real(dp) :: rr,zzpi -!============================================================================= infz=.false. z=cmplx(v(1),v(2),dpc) if(v(3)>u0)then @@ -1769,17 +1769,15 @@ subroutine ctoz(v, z,infz)! [ctoz] z=z*zzpi end subroutine ctoz -!============================================================================= +!> routine ztoc +!! @author R. J. Purser subroutine ztoc(z,infz, v)! [ztoc] -!============================================================================= implicit none complex(dpc), intent(IN ):: z logical, intent(IN ):: infz real(dp),dimension(3),intent(OUT):: v -!----------------------------------------------------------------------------- real(dp),parameter:: zero=0_dp,one=1_dp,two=2_dp real(dp) :: r,q,rs,rsc,rsbi -!============================================================================= if(infz)then; v=(/zero,zero,-one/); return; endif r=real(z); q=aimag(z); rs=r*r+q*q rsc=one-rs @@ -1789,27 +1787,25 @@ subroutine ztoc(z,infz, v)! [ztoc] v(3)=rsc*rsbi end subroutine ztoc -!============================================================================= +!> The convention adopted for the complex derivative is that, for a complex +!! infinitesimal map displacement, delta_z, the corresponding infinitesimal +!! change of cartesian vector position is delta_v given by: +!! delta_v = Real(vd*delta_z). +!! Thus, by a kind of Cauchy-Riemann relation, Imag(vd)=v CROSS Real(vd). +!! THE DERIVATIVE FOR THE IDEAL POINT AT INFINITY HAS NOT BEEN CODED YET!!! +!! @param z delta_z +!! @param v delta_v +!! @author R. J. Purser subroutine ztocd(z,infz, v,vd)! [ztoc] -!============================================================================= -! The convention adopted for the complex derivative is that, for a complex -! infinitesimal map displacement, delta_z, the corresponding infinitesimal -! change of cartesian vector position is delta_v given by: -! delta_v = Real(vd*delta_z). -! Thus, by a kind of Cauchy-Riemann relation, Imag(vd)=v CROSS Real(vd). -! THE DERIVATIVE FOR THE IDEAL POINT AT INFINITY HAS NOT BEEN CODED YET!!! -!============================================================================= implicit none complex(dpc), intent(IN ):: z logical, intent(IN ):: infz real(dp),dimension(3), intent(OUT):: v complex(dpc),dimension(3),intent(OUT):: vd -!----------------------------------------------------------------------------- real(dp),parameter :: zero=0_dp,one=1_dp,two=2_dp,four=4_dp real(dp) :: r,q,rs,rsc,rsbi,rsbis real(dp),dimension(3):: u1,u2 integer(spi) :: i -!============================================================================= if(infz)then; v=(/zero,zero,-one/); return; endif r=real(z); q=aimag(z); rs=r*r+q*q rsc=one-rs @@ -1827,22 +1823,22 @@ subroutine ztocd(z,infz, v,vd)! [ztoc] enddo end subroutine ztocd -!============================================================================ +!> Find the Mobius transformation complex coefficients, aa,bb,cc,dd, +!! with aa*dd-bb*cc=1, for a standard (north-)polar stereographic transformation +!! that takes cartesian point, xc0 to the north pole, xc1 to (lat=0,lon=0), +!! xc2 to the south pole (=complex infinity). +!! @param aa Mobius transformation complex coefficients +!! @param bb Mobius transformation complex coefficients +!! @param cc Mobius transformation complex coefficients +!! @param dd Mobius transformation complex coefficients +!! @author R. J. Purser subroutine setmobius(xc0,xc1,xc2, aa,bb,cc,dd)! [setmobius] -!============================================================================ -! Find the Mobius transformation complex coefficients, aa,bb,cc,dd, -! with aa*dd-bb*cc=1, for a standard (north-)polar stereographic transformation -! that takes cartesian point, xc0 to the north pole, xc1 to (lat=0,lon=0), -! xc2 to the south pole (=complex infinity). -!============================================================================ implicit none real(dp),dimension(3),intent(IN ):: xc0,xc1,xc2 complex(dpc), intent(OUT):: aa,bb,cc,dd -!---------------------------------------------------------------------------- real(dp),parameter:: zero=0_dp,one=1_dp logical :: infz0,infz1,infz2 complex(dpc) :: z0,z1,z2,z02,z10,z21 -!============================================================================ call ctoz(xc0,z0,infz0) call ctoz(xc1,z1,infz1) call ctoz(xc2,z2,infz2) @@ -1854,7 +1850,6 @@ subroutine setmobius(xc0,xc1,xc2, aa,bb,cc,dd)! [setmobius] (z1==z2.and.infz1.eqv.infz2).or.& (z2==z0.and.infz2.eqv.infz0)) & stop 'In setmobius; anchor points must be distinct' - if(infz2 .or. (.not.infz0 .and. abs(z0) Find the Mobius transformation complex coefficients, aa,bb,cc,dd, +!! with aa*dd-bb*cc=1, +!! that takes polar stereographic point, z0 to the north pole, +!! z1 to (lat=0,lon=0), z2 to the south pole (=complex infinity). +!! Should any one of z0,z1,z2 be itself the "point at infinity" its +!! corresponding infz will be set "true" (and the z value itself not used). +!! This routine is like setmobius, except the three fixed points defining +!! the mapping are given in standard complex stereographic form, together +!! with the logical codes "infzn" that are TRUE if that point is itself +!! the projection pole (i.e., the South Pole for a north polar stereographic). +!! @param aa Mobius transformation complex coefficients +!! @param bb Mobius transformation complex coefficients +!! @param cc Mobius transformation complex coefficients +!! @param dd Mobius transformation complex coefficients +!! @author R. J. Purser subroutine zsetmobius(z0,infz0, z1,infz1, z2,infz2, aa,bb,cc,dd) -! [setmobius] -!============================================================================ -! Find the Mobius transformation complex coefficients, aa,bb,cc,dd, -! with aa*dd-bb*cc=1, -! that takes polar stereographic point, z0 to the north pole, -! z1 to (lat=0,lon=0), z2 to the south pole (=complex infinity). -! Should any one of z0,z1,z2 be itself the "point at infinity" its -! corresponding infz will be set "true" (and the z value itself not used). -! This routine is like setmobius, except the three fixed points defining -! the mapping are given in standard complex stereographic form, together -! with the logical codes "infzn" that are TRUE if that point is itself -! the projection pole (i.e., the South Pole for a north polar stereographic). -!============================================================================ implicit none complex(dp), intent(IN ):: z0,z1,z2 logical, intent(IN ):: infz0,infz1,infz2 complex(dpc), intent(OUT):: aa,bb,cc,dd -!---------------------------------------------------------------------------- real(dp),parameter:: zero=0_dp,one=1_dp complex(dpc) :: z02,z10,z21 -!============================================================================ z21=z2-z1 z02=z0-z2 z10=z1-z0 - if( (z0==z1.and.infz0.eqv.infz1).or.& (z1==z2.and.infz1.eqv.infz2).or.& (z2==z0.and.infz2.eqv.infz0)) & stop 'In setmobius; anchor points must be distinct' - if(infz2 .or. (.not.infz0 .and. abs(z0) Perform a complex Mobius transformation from (z,infz) to (w,infw) +!! where the transformation coefficients are the standard aa,bb,cc,dd. +!! Infz is .TRUE. only when z is at complex infinity; likewise infw and w. +!! For these infinite cases, it is important that numerical z==(0,0). +!! @param z Mobius transformation +!! @param infz Mobius transformation +!! @param w Mobius transformation +!! @param infw Mobius transformation +!! @param aa Mobius transformation complex coefficients +!! @param bb Mobius transformation complex coefficients +!! @param cc Mobius transformation complex coefficients +!! @param dd Mobius transformation complex coefficients +!! @author R. J. Purser subroutine zmobius(aa,bb,cc,dd, z,infz, w,infw)! [mobius] -!============================================================================= -! Perform a complex Mobius transformation from (z,infz) to (w,infw) -! where the transformation coefficients are the standard aa,bb,cc,dd. -! Infz is .TRUE. only when z is at complex infinity; likewise infw and w. -! For these infinite cases, it is important that numerical z==(0,0). -!============================================================================= implicit none complex(dpc),intent(IN ):: aa,bb,cc,dd,z logical, intent(IN ):: infz complex(dpc),intent(OUT):: w logical, intent(OUT):: infw -!----------------------------------------------------------------------------- real(dp),parameter:: zero=0_dp complex(dpc) :: top,bot -!============================================================================= w=0 infw=.false. if(infz)then @@ -2007,40 +2003,42 @@ subroutine zmobius(aa,bb,cc,dd, z,infz, w,infw)! [mobius] endif end subroutine zmobius -!============================================================================= +!> Perform a complex Mobius transformation from cartesian vz to cartesian vw +!! where the transformation coefficients are the standard aa,bb,cc,dd. +!! @param aa Mobius transformation coefficients +!! @param bb Mobius transformation coefficients +!! @param cc Mobius transformation coefficients +!! @param dd Mobius transformation coefficients +!! @param vz cartesian +!! @param vw cartesian +!! @author R. J. Purser subroutine cmobius(aa,bb,cc,dd, vz,vw)! [mobius] -!============================================================================= -! Perform a complex Mobius transformation from cartesian vz to cartesian vw -! where the transformation coefficients are the standard aa,bb,cc,dd. -!============================================================================= implicit none complex(dpc), intent(IN ):: aa,bb,cc,dd real(dp),dimension(3),intent(IN ):: vz real(dp),dimension(3),intent(OUT):: vw -!----------------------------------------------------------------------------- complex(dpc):: z,w logical :: infz,infw -!============================================================================= call ctoz(vz, z,infz) call zmobius(aa,bb,cc,dd, z,infz, w,infw) call ztoc(w,infw, vw) end subroutine cmobius -!============================================================================ +!> Perform the inverse of the mobius transformation with coefficients, +!! {aa,bb,cc,dd}. +!! @param aa Mobius transformation coefficients +!! @param bb Mobius transformation coefficients +!! @param cc Mobius transformation coefficients +!! @param dd Mobius transformation coefficients +!! @author R. J. Purser subroutine zmobiusi(aa,bb,cc,dd, zz,infz, zw,infw) ! [mobiusi] -!============================================================================ -! Perform the inverse of the mobius transformation with coefficients, -! {aa,bb,cc,dd}. -!============================================================================ implicit none complex(dpc),intent(IN ):: aa,bb,cc,dd,zz logical, intent(IN ):: infz complex(dpc),intent(OUT):: zw logical, intent(OUT):: infw -!---------------------------------------------------------------------------- real(dp),parameter:: one=1_dp complex(dpc) :: aai,bbi,cci,ddi,d -!============================================================================ d=one/(aa*dd-bb*cc) aai=dd*d ddi=aa*d From 0c52d57fbb6322957094e270d88bfe0e3593f784 Mon Sep 17 00:00:00 2001 From: "Lin.Gan" Date: Wed, 3 Mar 2021 19:52:28 +0000 Subject: [PATCH 2/3] Doxygen improvements for sorc/grid_tools.fd/regional_esg_grid.fd/pmat4.f90 Part of #366 --- .../regional_esg_grid.fd/pmat4.f90 | 547 ++++++++++++------ 1 file changed, 366 insertions(+), 181 deletions(-) diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/pmat4.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/pmat4.f90 index 8eb24628a..506d7d99a 100644 --- a/sorc/grid_tools.fd/regional_esg_grid.fd/pmat4.f90 +++ b/sorc/grid_tools.fd/regional_esg_grid.fd/pmat4.f90 @@ -45,7 +45,7 @@ !! Additional Modules : pkind, pietc !! @author R. J. Purser @date Oct 2005 -!> Module for handy vector and matrix operations in Euclidean geometry +!> Module for handy vector and matrix operations in Euclidean geometry. !! @author R. J. Purser module pmat4 !============================================================================ @@ -115,7 +115,10 @@ module pmat4 contains -!> sqrt calculation for absv_s function +!> Doing sqrt calculation for absv_s function. +!! +!! @param[in] a real type input value +!! @param[out] s result !! @author R. J. Purser function absv_s(a)result(s)! [absv] implicit none @@ -124,7 +127,10 @@ function absv_s(a)result(s)! [absv] s=sqrt(dot_product(a,a)) end function absv_s -!> sqrt calculation for absv_d function +!> Doing sqrt calculation for absv_d function. +!! +!! @param[in] a real type input value +!! @param[out] s result !! @author R. J. Purser function absv_d(a)result(s)! [absv] implicit none @@ -133,7 +139,10 @@ function absv_d(a)result(s)! [absv] s=sqrt(dot_product(a,a)) end function absv_d -!> calculation for normalized_s function +!> Doing calculation for normalized_s function. +!! +!! @param[in] a real type input value +!! @param[out] b result !! @author R. J. Purser function normalized_s(a)result(b)! [normalized] use pietc_s, only: u0 @@ -144,7 +153,10 @@ function normalized_s(a)result(b)! [normalized] s=absv_s(a); if(s==u0)then; b=u0;else;b=a/s;endif end function normalized_s -!> calculation for normalized_d function +!> Doing calculation for normalized_d function. +!! +!! @param[in] a real type input value +!! @param[out] b result !! @author R. J. Purser function normalized_d(a)result(b)! [normalized] use pietc, only: u0 @@ -155,7 +167,11 @@ function normalized_d(a)result(b)! [normalized] s=absv_d(a); if(s==u0)then; b=u0;else;b=a/s;endif end function normalized_d -!> calculation for orthogonalized_s function +!> Doing calculation for orthogonalized_s function. +!! +!! @param[in] u real type input value +!! @param[in] a real type input value +!! @param[out] b result !! @author R. J. Purser function orthogonalized_s(u,a)result(b)! [orthogonalized] implicit none @@ -166,7 +182,11 @@ function orthogonalized_s(u,a)result(b)! [orthogonalized] s=dot_product(u,a); b=a-u*s end function orthogonalized_s -!> calculation for orthogonalized_d function +!> Doing calculation for orthogonalized_d function. +!! +!! @param[in] u real type input value +!! @param[in] a real type input value +!! @param[out] b result !! @author R. J. Purser function orthogonalized_d(u,a)result(b)! [orthogonalized] implicit none @@ -177,7 +197,11 @@ function orthogonalized_d(u,a)result(b)! [orthogonalized] s=dot_product(u,a); b=a-u*s end function orthogonalized_d -!> calculation for cross_product_s function +!> Doing calculation for cross_product_s function. +!! +!! @param[in] a real type input value +!! @param[in] b real type input value +!! @param[out] c result !! @author R. J. Purser function cross_product_s(a,b)result(c)! [cross_product] implicit none @@ -186,7 +210,11 @@ function cross_product_s(a,b)result(c)! [cross_product] c(1)=a(2)*b(3)-a(3)*b(2); c(2)=a(3)*b(1)-a(1)*b(3); c(3)=a(1)*b(2)-a(2)*b(1) end function cross_product_s -!> calculation for cross_product_d function +!> Doing calculation for cross_product_d function. +!! +!! @param[in] a real type input value +!! @param[in] b real type input value +!! @param[out] c result !! @author R. J. Purser function cross_product_d(a,b)result(c)! [cross_product] implicit none @@ -199,10 +227,11 @@ end function cross_product_d !! three 4-vectors, u, v, w, with the sign convention !! that ordered, {u,v,w,x} form a right-handed quartet !! in the generic case (determinant >= 0). -!! @param x triple-cross-product vector -!! @param u vector -!! @param v vector -!! @param w vector +!! +!! @param[in] u vector +!! @param[in] v vector +!! @param[in] w vector +!! @param[out] x triple-cross-product vector !! @author R. J. Purser function triple_cross_product_s(u,v,w)result(x)! [cross_product] implicit none @@ -218,7 +247,12 @@ function triple_cross_product_s(u,v,w)result(x)! [cross_product] x(4)= uv12*w(3)-uv13*w(2)+uv23*w(1) end function triple_cross_product_s -!> calculation for triple_cross_product_d function +!> Doing calculation for triple_cross_product_d function. +!! +!! @param[in] u vector +!! @param[in] v vector +!! @param[in] w vector +!! @param[out] x result !! @author R. J. Purser function triple_cross_product_d(u,v,w)result(x)! [cross_product] implicit none @@ -234,10 +268,13 @@ function triple_cross_product_d(u,v,w)result(x)! [cross_product] x(4)= uv12*w(3)-uv13*w(2)+uv23*w(1) end function triple_cross_product_d -!> calculation for outer_product_s function +!> Doing calculation for outer_product_s function. +!! +!! @param[in] a real type input value +!! @param[in] b real type input value +!! @param[out] c result !! @author R. J. Purser function outer_product_s(a,b)result(c)! [outer_product] -!============================================================================= implicit none real(sp),dimension(:), intent(in ):: a real(sp),dimension(:), intent(in ):: b @@ -247,7 +284,11 @@ function outer_product_s(a,b)result(c)! [outer_product] do i=1,nb; c(:,i)=a*b(i); enddo end function outer_product_s -!> Calculation for outer_product_d function +!> Calculation for outer_product_d function. +!! +!! @param[in] a real type input value +!! @param[in] b real type input value +!! @param[out] c result !! @author R. J. Purser function outer_product_d(a,b)result(c)! [outer_product] implicit none @@ -259,7 +300,11 @@ function outer_product_d(a,b)result(c)! [outer_product] do i=1,nb; c(:,i)=a*b(i); enddo end function outer_product_d -!> Calculation for outer_product_i function +!> Calculation for outer_product_i function. +!! +!! @param[in] a input value +!! @param[in] b input value +!! @param[out] c result !! @author R. J. Purser function outer_product_i(a,b)result(c)! [outer_product] implicit none @@ -271,7 +316,12 @@ function outer_product_i(a,b)result(c)! [outer_product] do i=1,nb; c(:,i)=a*b(i); enddo end function outer_product_i -!> Calculation for triple_product_s function +!> Calculation for triple_product_s function. +!! +!! @param[in] a real type input value +!! @param[in] b real type input value +!! @param[in] c real type input value +!! @param[out] tripleproduct result !! @author R. J. Purser function triple_product_s(a,b,c)result(tripleproduct)! [triple_product] implicit none @@ -280,7 +330,12 @@ function triple_product_s(a,b,c)result(tripleproduct)! [triple_product] tripleproduct=dot_product( cross_product(a,b),c ) end function triple_product_s -!> Calculation for triple_product_d function +!> Calculation for triple_product_d function. +!! +!! @param[in] a real type input value +!! @param[in] b real type input value +!! @param[in] c real type input value +!! @param[out] tripleproduct result !! @author R. J. Purser function triple_product_d(a,b,c)result(tripleproduct)! [triple_product] implicit none @@ -289,7 +344,10 @@ function triple_product_d(a,b,c)result(tripleproduct)! [triple_product] tripleproduct=dot_product( cross_product(a,b),c ) end function triple_product_d -!> Calculation for det_s function +!> Calculation for det_s function. +!! +!! @param[in] a real type input value +!! @param[out] det result !! @author R. J. Purser function det_s(a)result(det)! [det] use pietc_s, only: u0 @@ -307,7 +365,10 @@ function det_s(a)result(det)! [det] endif end function det_s -!> Calculation for det_d function +!> Calculation for det_d function. +!! +!! @param[in] a real type input value +!! @param[out] det result !! @author R. J. Purser function det_d(a)result(det)! [det] use pietc, only: u0 @@ -325,7 +386,10 @@ function det_d(a)result(det)! [det] endif end function det_d -!> Calculation for det_i function +!> Calculation for det_i function. +!! +!! @param[in] a real type input value +!! @param[out] idet result !! @author R. J. Purser function det_i(a)result(idet)! [det] implicit none @@ -336,7 +400,10 @@ function det_i(a)result(idet)! [det] b=a; bdet=det(b); idet=nint(bdet) end function det_i -!> Calculation for det_id function +!> Calculation for det_id function. +!! +!! @param[in] a real type input value +!! @param[out] idet result !! @author R. J. Purser function det_id(a)result(idet)! [det] use pkind, only: dp,dpi @@ -348,7 +415,10 @@ function det_id(a)result(idet)! [det] b=a; bdet=det(b); idet=nint(bdet) end function det_id -!> Calculation for axial3_s function +!> Calculation for axial3_s function. +!! +!! @param[in] a real type input value +!! @param[out] b result !! @author R. J. Purser function axial3_s(a)result(b)! [axial] use pietc_s, only: u0 @@ -358,7 +428,10 @@ function axial3_s(a)result(b)! [axial] b=u0;b(3,2)=a(1);b(1,3)=a(2);b(2,1)=a(3);b(2,3)=-a(1);b(3,1)=-a(2);b(1,2)=-a(3) end function axial3_s -!> Calculation for axial3_d function +!> Calculation for axial3_d function. +!! +!! @param[in] a real type input value +!! @param[out] b result !! @author R. J. Purser function axial3_d(a)result(b)! [axial] use pietc, only: u0 @@ -368,7 +441,10 @@ function axial3_d(a)result(b)! [axial] b=u0;b(3,2)=a(1);b(1,3)=a(2);b(2,1)=a(3);b(2,3)=-a(1);b(3,1)=-a(2);b(1,2)=-a(3) end function axial3_d -!> Calculation for axial33_s function +!> Calculation for axial33_s function. +!! +!! @param[in] b real type input value +!! @param[out] a result !! @author R. J. Purser function axial33_s(b)result(a)! [axial] use pietc_s, only: o2 @@ -378,7 +454,10 @@ function axial33_s(b)result(a)! [axial] a(1)=(b(3,2)-b(2,3))*o2; a(2)=(b(1,3)-b(3,1))*o2; a(3)=(b(2,1)-b(1,2))*o2 end function axial33_s -!> Calculation for axial33_d function +!> Calculation for axial33_d function. +!! +!! @param[in] b real type input value +!! @param[out] a result !! @author R. J. Purser function axial33_d(b)result(a)! [axial] use pietc, only: o2 @@ -388,7 +467,10 @@ function axial33_d(b)result(a)! [axial] a(1)=(b(3,2)-b(2,3))*o2; a(2)=(b(1,3)-b(3,1))*o2; a(3)=(b(2,1)-b(1,2))*o2 end function axial33_d -!> Calculation for diagn_s function +!> Calculation for diagn_s function. +!! +!! @param[in] a real type input value +!! @param[out] b result !! @author R. J. Purser function diagn_s(a)result(b)! [diag] use pietc, only: u0 @@ -400,7 +482,10 @@ function diagn_s(a)result(b)! [diag] b=u0; do i=1,n; b(i,i)=a(i); enddo end function diagn_s -!> Calculation for diagn_d function +!> Calculation for diagn_d function. +!! +!! @param[in] a real type input value +!! @param[out] b result !! @author R. J. Purser function diagn_d(a)result(b)! [diag] use pietc, only: u0 @@ -412,7 +497,10 @@ function diagn_d(a)result(b)! [diag] b=u0; do i=1,n; b(i,i)=a(i); enddo end function diagn_d -!> Calculation for diagn_i function +!> Calculation for diagn_i function. +!! +!! @param[in] a input value +!! @param[out] b result !! @author R. J. Purser function diagn_i(a)result(b)! [diag] implicit none @@ -423,7 +511,10 @@ function diagn_i(a)result(b)! [diag] b=0; do i=1,n; b(i,i)=a(i); enddo end function diagn_i -!> Calculation for diagnn_s function +!> Calculation for diagnn_s function. +!! +!! @param[in] b real type input value +!! @param[out] a result !! @author R. J. Purser function diagnn_s(b)result(a)! [diag] implicit none @@ -434,7 +525,10 @@ function diagnn_s(b)result(a)! [diag] do i=1,n; a(i)=b(i,i); enddo end function diagnn_s -!> Calculation for diagnn_d function +!> Calculation for diagnn_d function. +!! +!! @param[in] b real type input value +!! @param[out] a result !! @author R. J. Purser function diagnn_d(b)result(a)! [diag] implicit none @@ -445,7 +539,10 @@ function diagnn_d(b)result(a)! [diag] do i=1,n; a(i)=b(i,i); enddo end function diagnn_d -!> Calculation for diagnn_i function +!> Calculation for diagnn_i function. +!! +!! @param[in] b integer type input value +!! @param[out] a result !! @author R. J. Purser function diagnn_i(b)result(a)! [diag] implicit none @@ -456,7 +553,10 @@ function diagnn_i(b)result(a)! [diag] do i=1,n; a(i)=b(i,i); enddo end function diagnn_i -!> Calculation for trace_s function +!> Calculation for trace_s function. +!! +!! @param[in] b real type input value +!! @param[out] s result !! @author R. J. Purser function trace_s(b)result(s)! [trace] implicit none @@ -465,7 +565,10 @@ function trace_s(b)result(s)! [trace] s=sum(diag(b)) end function trace_s -!> Calculation for trace_d function +!> Calculation for trace_d function. +!! +!! @param[in] b real type input value +!! @param[out] s result !! @author R. J. Purser function trace_d(b)result(s)! [trace] implicit none @@ -474,7 +577,10 @@ function trace_d(b)result(s)! [trace] s=sum(diag(b)) end function trace_d -!> Calculation for trace_i function +!> Calculation for trace_i function. +!! +!! @param[in] b input value +!! @param[out] s result !! @author R. J. Purser function trace_i(b)result(s)! [trace] implicit none @@ -483,7 +589,10 @@ function trace_i(b)result(s)! [trace] s=sum(diag(b)) end function trace_i -!> Calculation for identity_i function +!> Calculation for identity_i function. +!! +!! @param[in] n input value +!! @param[out] a result !! @author R. J. Purser function identity_i(n)result(a)! [identity] implicit none @@ -493,7 +602,9 @@ function identity_i(n)result(a)! [identity] a=0; do i=1,n; a(i,i)=1; enddo end function identity_i -!> Calculation for identity3_i function +!> Calculation for identity3_i function. +!! +!! @param[out] a result !! @author R. J. Purser function identity3_i()result(a)! [identity] implicit none @@ -504,6 +615,10 @@ end function identity3_i !> Spherical area of right-angle triangle whose orthogonal sides have !! orthographic projection dimensions, sa and sb. +!! +!! @param[in] sa result +!! @param[in] sb result +!! @param[out] area result !! @author R. J. Purser function huarea_s(sa,sb)result(area)! [huarea] implicit none @@ -515,7 +630,11 @@ function huarea_s(sa,sb)result(area)! [huarea] area=asin(sa*sb/(1+ca*cb)) end function huarea_s -!> Calculation for huarea_d function +!> Calculation for huarea_d function. +!! +!! @param[in] sa result +!! @param[in] sb result +!! @param[out] area result !! @author R. J. Purser function huarea_d(sa,sb)result(area)! [huarea] implicit none @@ -531,9 +650,11 @@ end function huarea_d !! right-handed sense, by dropping a perpendicular to u0 on the longest side so !! that the area becomes the sum of areas of the two simpler right-angled !! triangles. -!! @param v1 area of the spherical triangle -!! @param v2 area of the spherical triangle -!! @param v3 area of the spherical triangle +!! +!! @param[in] v1 area of the spherical triangle +!! @param[in] v2 area of the spherical triangle +!! @param[in] v3 area of the spherical triangle +!! @param[out] area result !! @author R. J. Purser function sarea_s(v1,v2,v3)result(area)! [sarea] use pietc_s, only: zero=>u0 @@ -563,15 +684,16 @@ function sarea_s(v1,v2,v3)result(area)! [sarea] contains -!> Shift switch variable -!! @param u1 real variable to be shifted -!! @param u2 real variable to be shifted -!! @param u3 real variable to be shifted -!! @param d1 real variable to be shifted -!! @param d2 real variable to be shifted -!! @param d3 real variable to be shifted +!> Shift switch variable. +!! +!! @param[inout] u1 real variable to be shifted +!! @param[inout] u2 real variable to be shifted +!! @param[inout] u3 real variable to be shifted +!! @param[inout] d1 real variable to be shifted +!! @param[inout] d2 real variable to be shifted +!! @param[inout] d3 real variable to be shifted !! @author R. J. Purser - subroutine cyclic(u1,u2,u3,d1,d2,d3) +subroutine cyclic(u1,u2,u3,d1,d2,d3) implicit none real(sp),dimension(3),intent(INOUT):: u1,u2,u3 real(sp), intent(INOUT):: d1,d2,d3 @@ -582,10 +704,12 @@ subroutine cyclic(u1,u2,u3,d1,d2,d3) end subroutine cyclic end function sarea_s -!> Compute the area of sarea_d, {v1,v2,v3} -!! @param v1 area of the spherical triangle -!! @param v2 area of the spherical triangle -!! @param v3 area of the spherical triangle +!> Compute the area of sarea_d, {v1,v2,v3}. +!! +!! @param[in] v1 area of the spherical triangle +!! @param[in] v2 area of the spherical triangle +!! @param[in] v3 area of the spherical triangle +!! @param[out] area result !! @author R. J. Purser function sarea_d(v1,v2,v3)result(area)! [sarea] use pietc, only: zero=>u0 @@ -615,15 +739,16 @@ function sarea_d(v1,v2,v3)result(area)! [sarea] contains -!> Shift switch variable -!! @param u1 real variable to be shifted -!! @param u2 real variable to be shifted -!! @param u3 real variable to be shifted -!! @param d1 real variable to be shifted -!! @param d2 real variable to be shifted -!! @param d3 real variable to be shifted +!> Shift switch variable. +!! +!! @param[inout] u1 real variable to be shifted +!! @param[inout] u2 real variable to be shifted +!! @param[inout] u3 real variable to be shifted +!! @param[inout] d1 real variable to be shifted +!! @param[inout] d2 real variable to be shifted +!! @param[inout] d3 real variable to be shifted !! @author R. J. Purser - subroutine cyclic(u1,u2,u3,d1,d2,d3) +subroutine cyclic(u1,u2,u3,d1,d2,d3) implicit none real(dp),dimension(3),intent(INOUT):: u1,u2,u3 real(dp), intent(INOUT):: d1,d2,d3 @@ -640,11 +765,13 @@ end function sarea_d !! The computations are designed to give a proportionately accurate area !! estimate even when the triangle is very small, provided the B-increment !! is not disproportionately small compared to the other two sides. -!! @param rlat latitude -!! @param drlata latitudes -!! @param drlona longitudes -!! @param drlatb latitudes -!! @param drlonb longitudes +!! +!! @param[in] rlat latitude +!! @param[in] drlata latitudes +!! @param[in] drlona longitudes +!! @param[in] drlatb latitudes +!! @param[in] drlonb longitudes +!! @param[out] area result !! @author R. J. Purser function dtarea_s(rlat,drlata,drlona,drlatb,drlonb) result(area)! [sarea] use pietc_s, only: u0,u1 @@ -670,12 +797,14 @@ function dtarea_s(rlat,drlata,drlona,drlatb,drlonb) result(area)! [sarea] area=huarea(-sa,sb)+huarea(sc,-sa) end function dtarea_s -!> Compute the area with dtarea_d -!! @param rlat latitude -!! @param drlata latitudes -!! @param drlona longitudes -!! @param drlatb latitudes -!! @param drlonb longitudes +!> Compute the area with dtarea_d. +!! +!! @param[in] rlat latitude +!! @param[in] drlata latitudes +!! @param[in] drlona longitudes +!! @param[in] drlatb latitudes +!! @param[in] drlonb longitudes +!! @param[out] area result !! @author R. J. Purser function dtarea_d(rlat,drlata,drlona,drlatb,drlonb) result(area)! [sarea] use pietc, only: u0,u1 @@ -709,13 +838,15 @@ end function dtarea_d !! estimate even when the quadrilateral is very small, provided the !! diagonal making the B-increment is not disproportionately small compared to !! the characteristic size of the quadrilateral. -!! @param rlat latitude -!! @param drlata latitudes -!! @param drlona longitudes -!! @param drlatb latitudes -!! @param drlonb longitudes -!! @param drlatc latitudes -!! @param drlonc longitudes +!! +!! @param[in] rlat latitude +!! @param[in] drlata latitudes +!! @param[in] drlona longitudes +!! @param[in] drlatb latitudes +!! @param[in] drlonb longitudes +!! @param[in] drlatc latitudes +!! @param[in] drlonc longitudes +!! @param[out] area result !! @author R. J. Purser function dqarea_s &! [sarea] (rlat,drlata,drlona,drlatb,drlonb,drlatc,drlonc) result(area) @@ -726,14 +857,16 @@ function dqarea_s &! [sarea] -sarea(rlat,drlatc,drlonc,drlatb,drlonb) end function dqarea_s -!> Compute the area using dqarea_d -!! @param rlat latitude -!! @param drlata latitudes -!! @param drlona longitudes -!! @param drlatb latitudes -!! @param drlonb longitudes -!! @param drlatc latitudes -!! @param drlonc longitudes +!> Compute the area using dqarea_d. +!! +!! @param[in] rlat latitude +!! @param[in] drlata latitudes +!! @param[in] drlona longitudes +!! @param[in] drlatb latitudes +!! @param[in] drlonb longitudes +!! @param[in] drlatc latitudes +!! @param[in] drlonc longitudes +!! @param[out] area result !! @author R. J. Purser function dqarea_d &! [sarea] (rlat,drlata,drlona,drlatb,drlonb,drlatc,drlonc) result(area) @@ -744,10 +877,12 @@ function dqarea_d &! [sarea] -sarea(rlat,drlatc,drlonc,drlatb,drlonb) end function dqarea_d -!> Calculate dlltoxy_s -!! @param rlat latitude -!! @param drlat latitude -!! @param drlon longitudes +!> Calculate dlltoxy_s. +!! +!! @param[in] rlat latitude +!! @param[in] drlat latitude +!! @param[in] drlon longitudes +!! @param[out] x2 output !! @author R. J. Purser subroutine dlltoxy_s(rlat,drlat,drlon,x2)! [dlltoxy] use pietc_s, only: u2 @@ -759,10 +894,12 @@ subroutine dlltoxy_s(rlat,drlat,drlon,x2)! [dlltoxy] x2=(/clata*sin(drlon),sin(drlat)+u2*sin(rlat)*clata*hav(drlon)/) end subroutine dlltoxy_s -!> Calculate dlltoxy_d -!! @param rlat latitude -!! @param drlat latitude -!! @param drlon longitudes +!> Calculate dlltoxy_d. +!! +!! @param[in] rlat latitude +!! @param[in] drlat latitude +!! @param[in] drlon longitudes +!! @param[out] x2 output !! @author R. J. Purser subroutine dlltoxy_d(rlat,drlat,drlon,x2)! [dlltoxy] use pietc, only: u2 @@ -774,7 +911,10 @@ subroutine dlltoxy_d(rlat,drlat,drlon,x2)! [dlltoxy] x2=(/clata*sin(drlon),sin(drlat)+u2*sin(rlat)*clata*hav(drlon)/) end subroutine dlltoxy_d -!> Haversine function +!> Haversine function. +!! +!! @param[in] t input +!! @param[out] a result !! @author R. J. Purser function hav_s(t) result(a)! [hav] use pietc_s, only: o2 @@ -784,7 +924,10 @@ function hav_s(t) result(a)! [hav] a=(sin(t*o2))**2 end function hav_s -!> hav_d function +!> Doing hav_d function. +!! +!! @param[in] t input +!! @param[out] a result !! @author R. J. Purser function hav_d(t) result(a)! [hav] use pietc, only: o2 @@ -794,8 +937,9 @@ function hav_d(t) result(a)! [hav] a=(sin(t*o2))**2 end function hav_d -!> Normalize the given vector -!! @param v vector +!> Normalize the given vector. +!! +!! @param[inout] v vector !! @author R. J. Purser subroutine normalize_s(v)! [normalize] use pietc_s, only: u0,u1 @@ -805,8 +949,9 @@ subroutine normalize_s(v)! [normalize] s=absv(v); if(s==0)then; v=u0; v(1)=u1; else; v=v/s; endif end subroutine normalize_s -!> normalize_d calculation for given vector -!! @param v vector +!> Doing normalize_d calculation for given vector. +!! +!! @param[inout] v vector !! @author R. J. Purser subroutine normalize_d(v)! [normalize] use pietc, only: u0,u1 @@ -816,7 +961,7 @@ subroutine normalize_d(v)! [normalize] s=absv(v); if(s==u0)then; v=0; v(1)=u1; else; v=v/s; endif end subroutine normalize_d -!> gram_s routine +!> ??? !! @author R. J. Purser subroutine gram_s(as,b,nrank,det)! [gram] use pietc_s, only: u0,u1 @@ -878,7 +1023,7 @@ subroutine gram_s(as,b,nrank,det)! [gram] enddo end subroutine gram_s -!> gram_d routine +!> ??? !! @author R. J. Purser subroutine gram_d(as,b,nrank,det)! [gram] use pietc, only: u0,u1 @@ -945,7 +1090,8 @@ end subroutine gram_d !! matrix is singular, the "sign" of the determinant, detsign, is returned !! as zero (instead of either +1 or -1) and ldet is then just the log of !! the nonzero factors found by the process. -!! @param detsign singular determinant +!! +!! @param[out] detsign singular determinant !! @author R. J. Purser subroutine graml_d(as,b,nrank,detsign,ldet)! [gram] use pietc, only: u0 @@ -1019,6 +1165,9 @@ end subroutine graml_d !> A "plain" (unpivoted) version of Gram-Schmidt, for square matrices only. +!! +!! @param[inout] b matrices +!! @param[out] nrank result !! @author R. J. Purser subroutine plaingram_s(b,nrank)! [gram] use pietc_s, only: u0 @@ -1051,6 +1200,9 @@ subroutine plaingram_s(b,nrank)! [gram] end subroutine plaingram_s !> A "plain" (unpivoted) version of Gram-Schmidt, for square matrices only. +!! +!! @param[inout] b matrices +!! @param[out] nrank result !! @author R. J. Purser subroutine plaingram_d(b,nrank)! [gram] use pietc, only: u0 @@ -1095,10 +1247,11 @@ end subroutine plaingram_d !! (in this module) because the negligibility criterion depends upon an !! "epsilon" value that is fixed (10**(-13)) and assumes elements of a are !! never too different in magnitude from unity, unless they are actually zero. -!! @param a rectangular input matrix -!! @param ipiv pivoting sequence -!! @param tt row-normalization -!! @param b orthonormalized rows +!! +!! @param[in] a rectangular input matrix +!! @param[out] ipiv pivoting sequence +!! @param[out] tt row-normalization +!! @param[out] b orthonormalized rows !! @author R. J. Purser subroutine rowgram(m,n,a,ipiv,tt,b,rank)! [gram] use pietc, only: u0,u1 @@ -1177,8 +1330,9 @@ end subroutine rowgram !> Apply the row-operations, implied by ipiv and tt returned by rowgram, to !! the single column vector, v, to produce the transformed vector vv. -!! @param v vector -!! @param vv vector +!! +!! @param[in] v vector +!! @param[out] vv vector !! @author R. J. Purser subroutine rowops(m,n,ipiv,tt,v,vv)! [rowops] implicit none @@ -1212,10 +1366,10 @@ end subroutine rowops !! Once a solution for D and E is found, return their exponentials, d and e, !! together with the rescaled matrix aa such that a = d.aa.e when d and e are !! interpreted as diagonal matrices. -!! @param d positive diagonals -!! @param e positive diagonals -!! @param a positive diagonals -!! @param f Lagrange multiplier +!! +!! @param[in] a positive diagonals +!! @param[out] d positive diagonals +!! @param[out] e positive diagonals !! @author R. J. Purser subroutine corral(m,n,mask,a,d,aa,e)! [corral] use pietc, only: u0,u1 @@ -1290,8 +1444,9 @@ end subroutine corral !! interpreted as encoding a rotation (as in subroutine axtorot). Note that !! such ax3 are not unique -- adding any multiple of 2*pi times the parallel !! unit vector leads to the same orth33. -!! @param orth33 3*3 proper rotation matrix -!! @param ax3 axial 3-vector +!! +!! @param[in] orth33 3*3 proper rotation matrix +!! @param[out] ax3 axial 3-vector !! @author R. J. Purser subroutine rottoax(orth33,ax3)! [rottoax] implicit none @@ -1319,8 +1474,9 @@ end subroutine rottoax !> Construct the 3*3 orthogonal matrix, orth33, that corresponds to the !! proper rotation encoded by the 3-vector, ax3. The antisymmetric matrix !! ax33 equivalent to the axial vector ax3 is exponentiated to obtain orth33. -!! @param orth33 3*3 orthogonal matrix -!! @param ax3 axial 3-vector +!! +!! @param[in] ax3 axial 3-vector +!! @param[out] orth33 3*3 orthogonal matrix !! @author R. J. Purser subroutine axtorot(ax3,orth33)! [axtorot] implicit none @@ -1332,8 +1488,9 @@ subroutine axtorot(ax3,orth33)! [axtorot] end subroutine axtorot !> Go from the spinor to the quaternion representation. -!! @param cspin spinor representation -!! @param q quaternion representation +!! +!! @param[in] cspin spinor representation +!! @param[out] q quaternion representation !! @author R. J. Purser subroutine spintoq(cspin,q)! [spintoq] implicit none @@ -1343,9 +1500,10 @@ subroutine spintoq(cspin,q)! [spintoq] q(2)=real(cspin(2,1)); q(1)=aimag(cspin(2,1)) end subroutine spintoq -!> Go from the quaternion to the spinor representation -!! @param cspin spinor representation -!! @param q quaternion representation +!> Go from the quaternion to the spinor representation. +!! +!! @param[in] q quaternion representation +!! @param[out] cspin spinor representation !! @author R. J. Purser subroutine qtospin(q,cspin)! [qtospin] implicit none @@ -1357,9 +1515,10 @@ subroutine qtospin(q,cspin)! [qtospin] cspin(2,2)=cmplx( q(0),-q(3)) end subroutine qtospin -!> Go from rotation matrix to quaternion representation -!! @param rot rotation matrix -!! @param q quaternion representation +!> Go from rotation matrix to quaternion representation. +!! +!! @param[in] rot rotation matrix +!! @param[out] q quaternion representation !! @author R. J. Purser subroutine rottoq(rot,q)! [rottoq] use pietc, only: zero=>u0,one=>u1,two=>u2 @@ -1412,9 +1571,10 @@ subroutine rottoq(rot,q)! [rottoq] q(1:3)=t1(3,:)*s end subroutine rottoq -!> Go from quaternion to rotation matrix representations -!! @param q quaternion representation -!! @param rot rotation matrix representations +!> Go from quaternion to rotation matrix representations. +!! +!! @param[in] q quaternion representation +!! @param[out] rot rotation matrix representations !! @author R. J. Purser subroutine qtorot(q,rot)! [qtorot] implicit none @@ -1423,9 +1583,10 @@ subroutine qtorot(q,rot)! [qtorot] call setem(q(0),q(1),q(2),q(3),rot) end subroutine qtorot -!> Go from an axial 3-vector to its equivalent quaternion -!! @param q quaternion -!! @param v axial 3-vector +!> Go from an axial 3-vector to its equivalent quaternion. +!! +!! @param[in] v axial 3-vector +!! @param[out] q quaternion !! @author R. J. Purser subroutine axtoq(v,q)! [axtoq] implicit none @@ -1436,9 +1597,10 @@ subroutine axtoq(v,q)! [axtoq] call rottoq(rot,q) end subroutine axtoq -!> Go from quaternion to axial 3-vector -!! @param q quaternion -!! @param v axial 3-vector +!> Go from quaternion to axial 3-vector. +!! +!! @param[in] q quaternion +!! @param[in] v axial 3-vector !! @author R. J. Purser subroutine qtoax(q,v)! [qtoax] implicit none @@ -1449,7 +1611,7 @@ subroutine qtoax(q,v)! [qtoax] call rottoax(rot,v) end subroutine qtoax -!> setem routine +!> ??? !! @author R. J. Purser subroutine setem(c,d,e,g,r)! [setem] implicit none @@ -1464,9 +1626,11 @@ subroutine setem(c,d,e,g,r)! [setem] r(3,2)=2*(eg+dc); r(1,3)=2*(dg+ec); r(2,1)=2*(de+gc) end subroutine setem -!> Multiply quaternions, a*b, assuming operation performed from right to left -!! @param a real quaternion -!! @param b real quaternion +!> Multiply quaternions, a*b, assuming operation performed from right to left. +!! +!! @param[in] a real quaternion +!! @param[in] b real quaternion +!! @param[out] c result !! @author R. J. Purser function mulqq(a,b)result(c)! [mulqq] implicit none @@ -1482,9 +1646,10 @@ end function mulqq !! Apply the iterated squaring method, m times, to the approximation to !! exp(a/(2**m)) obtained as a Taylor expansion of degree L !! See Fung, T. C., 2004, Int. J. Numer. Meth. Engng, 59, 1273--1286. -!! @param b exponential -!! @param a matrix -!! @param n degree +!! +!! @param[in] n degree +!! @param[in] a matrix +!! @param[out] b exponential matrix of a !! @author R. J. Purser subroutine expmat(n,a,b,detb)! [expmat] use pietc, only: u0,u1,u2,o2 @@ -1516,9 +1681,10 @@ subroutine expmat(n,a,b,detb)! [expmat] end subroutine expmat !> Like expmat, but for the 1st derivatives also. -!! @param b exponential -!! @param a matrix -!! @param n degree +!! +!! @param[in] n degree +!! @param[in] a matrix +!! @param[out] b exponential matrix of a !! @author R. J. Purser subroutine expmatd(n,a,b,bd,detb,detbd)! [expmat] use pietc, only: u0,u1,u2,o2 @@ -1577,9 +1743,10 @@ subroutine expmatd(n,a,b,bd,detb,detbd)! [expmat] end subroutine expmatd !> Like expmat, but for the 1st and 2nd derivatives also. -!! @param b exponential -!! @param a matrix -!! @param n degree +!! +!! @param[in] n degree +!! @param[in] a matrix +!! @param[out] b exponential matrix of a !! @author R. J. Purser subroutine expmatdd(n,a,b,bd,bdd,detb,detbd,detbdd)! [expmat] use pietc, only: u0,u1,u2,o2 @@ -1661,7 +1828,7 @@ subroutine expmatdd(n,a,b,bd,bdd,detb,detbd,detbdd)! [expmat] detbdd=u0; do ki=1,n; do kj=1,n; detbdd(ki,kj)=detb; enddo; enddo end subroutine expmatdd -!> routine zntay +!> ??? !! @author R. J. Purser subroutine zntay(n,z,zn)! [zntay] use pietc, only: u2 @@ -1691,7 +1858,7 @@ subroutine zntay(n,z,zn)! [zntay] print'("In zntay; full complement of iterations used")' end subroutine zntay -!> routine znfun +!> ??? !! @author R. J. Purser subroutine znfun(n,z,zn,znd,zndd,znddd)! [znfun] use pietc, only: u0,u2,u3 @@ -1749,6 +1916,8 @@ end subroutine znfun !! [ cc3, dd3 ] [ cc2, dd2 ] [ cc1, dd1 ] . !! !! Note that the determinant of these matrices is always +1 +!! +!! @param[in] v matric !! @author R. J. Purser subroutine ctoz(v, z,infz)! [ctoz] use pietc, only: u0,u1 @@ -1769,7 +1938,7 @@ subroutine ctoz(v, z,infz)! [ctoz] z=z*zzpi end subroutine ctoz -!> routine ztoc +!> ??? !! @author R. J. Purser subroutine ztoc(z,infz, v)! [ztoc] implicit none @@ -1793,8 +1962,9 @@ end subroutine ztoc !! delta_v = Real(vd*delta_z). !! Thus, by a kind of Cauchy-Riemann relation, Imag(vd)=v CROSS Real(vd). !! THE DERIVATIVE FOR THE IDEAL POINT AT INFINITY HAS NOT BEEN CODED YET!!! -!! @param z delta_z -!! @param v delta_v +!! +!! @param[in] z complex infinitesimal map displacement +!! @param[out] v cartesian vector position !! @author R. J. Purser subroutine ztocd(z,infz, v,vd)! [ztoc] implicit none @@ -1827,10 +1997,14 @@ end subroutine ztocd !! with aa*dd-bb*cc=1, for a standard (north-)polar stereographic transformation !! that takes cartesian point, xc0 to the north pole, xc1 to (lat=0,lon=0), !! xc2 to the south pole (=complex infinity). -!! @param aa Mobius transformation complex coefficients -!! @param bb Mobius transformation complex coefficients -!! @param cc Mobius transformation complex coefficients -!! @param dd Mobius transformation complex coefficients +!! +!! @param[in] xc0 cartesian point +!! @param[in] xc1 cartesian point +!! @param[in] xc2 cartesian point +!! @param[out] aa Mobius transformation complex coefficients +!! @param[out] bb Mobius transformation complex coefficients +!! @param[out] cc Mobius transformation complex coefficients +!! @param[out] dd Mobius transformation complex coefficients !! @author R. J. Purser subroutine setmobius(xc0,xc1,xc2, aa,bb,cc,dd)! [setmobius] implicit none @@ -1903,10 +2077,17 @@ end subroutine setmobius !! the mapping are given in standard complex stereographic form, together !! with the logical codes "infzn" that are TRUE if that point is itself !! the projection pole (i.e., the South Pole for a north polar stereographic). -!! @param aa Mobius transformation complex coefficients -!! @param bb Mobius transformation complex coefficients -!! @param cc Mobius transformation complex coefficients -!! @param dd Mobius transformation complex coefficients +!! +!! @param[in] z0 polar stereographic point +!! @param[in] infz0 point at infinity z0 +!! @param[in] z1 polar stereographic point +!! @param[in] infz1 point at infinity z0 +!! @param[in] z2 polar stereographic point +!! @param[in] infz2 point at infinity z0 +!! @param[out] aa Mobius transformation complex coefficients +!! @param[out] bb Mobius transformation complex coefficients +!! @param[out] cc Mobius transformation complex coefficients +!! @param[out] dd Mobius transformation complex coefficients !! @author R. J. Purser subroutine zsetmobius(z0,infz0, z1,infz1, z2,infz2, aa,bb,cc,dd) implicit none @@ -1969,14 +2150,15 @@ end subroutine zsetmobius !! where the transformation coefficients are the standard aa,bb,cc,dd. !! Infz is .TRUE. only when z is at complex infinity; likewise infw and w. !! For these infinite cases, it is important that numerical z==(0,0). -!! @param z Mobius transformation -!! @param infz Mobius transformation -!! @param w Mobius transformation -!! @param infw Mobius transformation -!! @param aa Mobius transformation complex coefficients -!! @param bb Mobius transformation complex coefficients -!! @param cc Mobius transformation complex coefficients -!! @param dd Mobius transformation complex coefficients +!! +!! @param[in] aa Mobius transformation complex coefficients +!! @param[in] bb Mobius transformation complex coefficients +!! @param[in] cc Mobius transformation complex coefficients +!! @param[in] dd Mobius transformation complex coefficients +!! @param[in] z Mobius transformation complex coefficients +!! @param[in] infz point at infinity z +!! @param[out] w Mobius transformation output +!! @param[out] infw point at infinity w !! @author R. J. Purser subroutine zmobius(aa,bb,cc,dd, z,infz, w,infw)! [mobius] implicit none @@ -2005,12 +2187,13 @@ end subroutine zmobius !> Perform a complex Mobius transformation from cartesian vz to cartesian vw !! where the transformation coefficients are the standard aa,bb,cc,dd. -!! @param aa Mobius transformation coefficients -!! @param bb Mobius transformation coefficients -!! @param cc Mobius transformation coefficients -!! @param dd Mobius transformation coefficients -!! @param vz cartesian -!! @param vw cartesian +!! +!! @param[in] aa Mobius transformation coefficients +!! @param[in] bb Mobius transformation coefficients +!! @param[in] cc Mobius transformation coefficients +!! @param[in] dd Mobius transformation coefficients +!! @param[in] vz Cartesian vaule +!! @param[out] vw Cartesian vaule !! @author R. J. Purser subroutine cmobius(aa,bb,cc,dd, vz,vw)! [mobius] implicit none @@ -2026,10 +2209,12 @@ end subroutine cmobius !> Perform the inverse of the mobius transformation with coefficients, !! {aa,bb,cc,dd}. -!! @param aa Mobius transformation coefficients -!! @param bb Mobius transformation coefficients -!! @param cc Mobius transformation coefficients -!! @param dd Mobius transformation coefficients +!! +!! @param[in] aa Mobius transformation coefficients +!! @param[in] bb Mobius transformation coefficients +!! @param[in] cc Mobius transformation coefficients +!! @param[in] dd Mobius transformation coefficients +!! @param[out] zw Inversed output !! @author R. J. Purser subroutine zmobiusi(aa,bb,cc,dd, zz,infz, zw,infw) ! [mobiusi] implicit none From 68f6eb647c05ecf81c3042150779a7e6f20adfcd Mon Sep 17 00:00:00 2001 From: Edward Hartnett Date: Fri, 5 Mar 2021 06:08:55 -0700 Subject: [PATCH 3/3] doxygen fixes --- .../regional_esg_grid.fd/pmat4.f90 | 180 ++++++++++++------ 1 file changed, 121 insertions(+), 59 deletions(-) diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/pmat4.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/pmat4.f90 index 506d7d99a..3cd11cb43 100644 --- a/sorc/grid_tools.fd/regional_esg_grid.fd/pmat4.f90 +++ b/sorc/grid_tools.fd/regional_esg_grid.fd/pmat4.f90 @@ -1,7 +1,9 @@ !> @file !! @brief Euclidean geometry, geometric (stereographic) projections, -!! related transformations (Mobius). -!! +!! related transformations (Mobius). +!! @author R. J. Purser @date Oct 2005 + +!> Module for handy vector and matrix operations in Euclidean geometry. !! Package for handy vector and matrix operations in Euclidean geometry. !! This package is primarily intended for 3D operations and three of the !! functions (Cross_product, Triple_product and Axial) do not possess simple @@ -40,12 +42,6 @@ !! and some associated mobius transformation utilities, since these complex !! operations have a strong geometrical flavor. !! -!! DIRECT DEPENDENCIES -!! Libraries[their Modules]: pmat[pmat] -!! Additional Modules : pkind, pietc -!! @author R. J. Purser @date Oct 2005 - -!> Module for handy vector and matrix operations in Euclidean geometry. !! @author R. J. Purser module pmat4 !============================================================================ @@ -118,7 +114,7 @@ module pmat4 !> Doing sqrt calculation for absv_s function. !! !! @param[in] a real type input value -!! @param[out] s result +!! @return s result !! @author R. J. Purser function absv_s(a)result(s)! [absv] implicit none @@ -130,7 +126,7 @@ end function absv_s !> Doing sqrt calculation for absv_d function. !! !! @param[in] a real type input value -!! @param[out] s result +!! @return s result !! @author R. J. Purser function absv_d(a)result(s)! [absv] implicit none @@ -142,7 +138,7 @@ end function absv_d !> Doing calculation for normalized_s function. !! !! @param[in] a real type input value -!! @param[out] b result +!! @return b result !! @author R. J. Purser function normalized_s(a)result(b)! [normalized] use pietc_s, only: u0 @@ -156,7 +152,7 @@ end function normalized_s !> Doing calculation for normalized_d function. !! !! @param[in] a real type input value -!! @param[out] b result +!! @return b result !! @author R. J. Purser function normalized_d(a)result(b)! [normalized] use pietc, only: u0 @@ -171,7 +167,7 @@ end function normalized_d !! !! @param[in] u real type input value !! @param[in] a real type input value -!! @param[out] b result +!! @return b result !! @author R. J. Purser function orthogonalized_s(u,a)result(b)! [orthogonalized] implicit none @@ -186,7 +182,7 @@ end function orthogonalized_s !! !! @param[in] u real type input value !! @param[in] a real type input value -!! @param[out] b result +!! @return b result !! @author R. J. Purser function orthogonalized_d(u,a)result(b)! [orthogonalized] implicit none @@ -201,7 +197,7 @@ end function orthogonalized_d !! !! @param[in] a real type input value !! @param[in] b real type input value -!! @param[out] c result +!! @return c result !! @author R. J. Purser function cross_product_s(a,b)result(c)! [cross_product] implicit none @@ -214,7 +210,7 @@ end function cross_product_s !! !! @param[in] a real type input value !! @param[in] b real type input value -!! @param[out] c result +!! @return c result !! @author R. J. Purser function cross_product_d(a,b)result(c)! [cross_product] implicit none @@ -231,7 +227,7 @@ end function cross_product_d !! @param[in] u vector !! @param[in] v vector !! @param[in] w vector -!! @param[out] x triple-cross-product vector +!! @return x triple-cross-product vector !! @author R. J. Purser function triple_cross_product_s(u,v,w)result(x)! [cross_product] implicit none @@ -252,7 +248,7 @@ end function triple_cross_product_s !! @param[in] u vector !! @param[in] v vector !! @param[in] w vector -!! @param[out] x result +!! @return x result !! @author R. J. Purser function triple_cross_product_d(u,v,w)result(x)! [cross_product] implicit none @@ -272,7 +268,7 @@ end function triple_cross_product_d !! !! @param[in] a real type input value !! @param[in] b real type input value -!! @param[out] c result +!! @return c result !! @author R. J. Purser function outer_product_s(a,b)result(c)! [outer_product] implicit none @@ -288,7 +284,7 @@ end function outer_product_s !! !! @param[in] a real type input value !! @param[in] b real type input value -!! @param[out] c result +!! @return c result !! @author R. J. Purser function outer_product_d(a,b)result(c)! [outer_product] implicit none @@ -304,7 +300,7 @@ end function outer_product_d !! !! @param[in] a input value !! @param[in] b input value -!! @param[out] c result +!! @return c result !! @author R. J. Purser function outer_product_i(a,b)result(c)! [outer_product] implicit none @@ -321,7 +317,7 @@ end function outer_product_i !! @param[in] a real type input value !! @param[in] b real type input value !! @param[in] c real type input value -!! @param[out] tripleproduct result +!! @return tripleproduct result !! @author R. J. Purser function triple_product_s(a,b,c)result(tripleproduct)! [triple_product] implicit none @@ -335,7 +331,7 @@ end function triple_product_s !! @param[in] a real type input value !! @param[in] b real type input value !! @param[in] c real type input value -!! @param[out] tripleproduct result +!! @return tripleproduct result !! @author R. J. Purser function triple_product_d(a,b,c)result(tripleproduct)! [triple_product] implicit none @@ -347,7 +343,7 @@ end function triple_product_d !> Calculation for det_s function. !! !! @param[in] a real type input value -!! @param[out] det result +!! @return det result !! @author R. J. Purser function det_s(a)result(det)! [det] use pietc_s, only: u0 @@ -368,7 +364,7 @@ end function det_s !> Calculation for det_d function. !! !! @param[in] a real type input value -!! @param[out] det result +!! @return det result !! @author R. J. Purser function det_d(a)result(det)! [det] use pietc, only: u0 @@ -389,7 +385,7 @@ end function det_d !> Calculation for det_i function. !! !! @param[in] a real type input value -!! @param[out] idet result +!! @return idet result !! @author R. J. Purser function det_i(a)result(idet)! [det] implicit none @@ -403,7 +399,7 @@ end function det_i !> Calculation for det_id function. !! !! @param[in] a real type input value -!! @param[out] idet result +!! @return idet result !! @author R. J. Purser function det_id(a)result(idet)! [det] use pkind, only: dp,dpi @@ -418,7 +414,7 @@ end function det_id !> Calculation for axial3_s function. !! !! @param[in] a real type input value -!! @param[out] b result +!! @return b result !! @author R. J. Purser function axial3_s(a)result(b)! [axial] use pietc_s, only: u0 @@ -431,7 +427,7 @@ end function axial3_s !> Calculation for axial3_d function. !! !! @param[in] a real type input value -!! @param[out] b result +!! @return b result !! @author R. J. Purser function axial3_d(a)result(b)! [axial] use pietc, only: u0 @@ -444,7 +440,7 @@ end function axial3_d !> Calculation for axial33_s function. !! !! @param[in] b real type input value -!! @param[out] a result +!! @return a result !! @author R. J. Purser function axial33_s(b)result(a)! [axial] use pietc_s, only: o2 @@ -457,7 +453,7 @@ end function axial33_s !> Calculation for axial33_d function. !! !! @param[in] b real type input value -!! @param[out] a result +!! @return a result !! @author R. J. Purser function axial33_d(b)result(a)! [axial] use pietc, only: o2 @@ -470,7 +466,7 @@ end function axial33_d !> Calculation for diagn_s function. !! !! @param[in] a real type input value -!! @param[out] b result +!! @return b result !! @author R. J. Purser function diagn_s(a)result(b)! [diag] use pietc, only: u0 @@ -485,7 +481,7 @@ end function diagn_s !> Calculation for diagn_d function. !! !! @param[in] a real type input value -!! @param[out] b result +!! @return b result !! @author R. J. Purser function diagn_d(a)result(b)! [diag] use pietc, only: u0 @@ -500,7 +496,7 @@ end function diagn_d !> Calculation for diagn_i function. !! !! @param[in] a input value -!! @param[out] b result +!! @return b result !! @author R. J. Purser function diagn_i(a)result(b)! [diag] implicit none @@ -514,7 +510,7 @@ end function diagn_i !> Calculation for diagnn_s function. !! !! @param[in] b real type input value -!! @param[out] a result +!! @return a result !! @author R. J. Purser function diagnn_s(b)result(a)! [diag] implicit none @@ -528,7 +524,7 @@ end function diagnn_s !> Calculation for diagnn_d function. !! !! @param[in] b real type input value -!! @param[out] a result +!! @return a result !! @author R. J. Purser function diagnn_d(b)result(a)! [diag] implicit none @@ -542,7 +538,7 @@ end function diagnn_d !> Calculation for diagnn_i function. !! !! @param[in] b integer type input value -!! @param[out] a result +!! @return a result !! @author R. J. Purser function diagnn_i(b)result(a)! [diag] implicit none @@ -556,7 +552,7 @@ end function diagnn_i !> Calculation for trace_s function. !! !! @param[in] b real type input value -!! @param[out] s result +!! @return s result !! @author R. J. Purser function trace_s(b)result(s)! [trace] implicit none @@ -568,7 +564,7 @@ end function trace_s !> Calculation for trace_d function. !! !! @param[in] b real type input value -!! @param[out] s result +!! @return s result !! @author R. J. Purser function trace_d(b)result(s)! [trace] implicit none @@ -580,7 +576,7 @@ end function trace_d !> Calculation for trace_i function. !! !! @param[in] b input value -!! @param[out] s result +!! @return s result !! @author R. J. Purser function trace_i(b)result(s)! [trace] implicit none @@ -592,7 +588,7 @@ end function trace_i !> Calculation for identity_i function. !! !! @param[in] n input value -!! @param[out] a result +!! @return a result !! @author R. J. Purser function identity_i(n)result(a)! [identity] implicit none @@ -604,7 +600,7 @@ end function identity_i !> Calculation for identity3_i function. !! -!! @param[out] a result +!! @return a result !! @author R. J. Purser function identity3_i()result(a)! [identity] implicit none @@ -616,9 +612,9 @@ end function identity3_i !> Spherical area of right-angle triangle whose orthogonal sides have !! orthographic projection dimensions, sa and sb. !! -!! @param[in] sa result -!! @param[in] sb result -!! @param[out] area result +!! @param[in] sa ??? +!! @param[in] sb ??? +!! @return area !! @author R. J. Purser function huarea_s(sa,sb)result(area)! [huarea] implicit none @@ -632,9 +628,9 @@ end function huarea_s !> Calculation for huarea_d function. !! -!! @param[in] sa result -!! @param[in] sb result -!! @param[out] area result +!! @param[in] sa ??? +!! @param[in] sb ??? +!! @return area !! @author R. J. Purser function huarea_d(sa,sb)result(area)! [huarea] implicit none @@ -654,7 +650,7 @@ end function huarea_d !! @param[in] v1 area of the spherical triangle !! @param[in] v2 area of the spherical triangle !! @param[in] v3 area of the spherical triangle -!! @param[out] area result +!! @return area result !! @author R. J. Purser function sarea_s(v1,v2,v3)result(area)! [sarea] use pietc_s, only: zero=>u0 @@ -709,7 +705,7 @@ end function sarea_s !! @param[in] v1 area of the spherical triangle !! @param[in] v2 area of the spherical triangle !! @param[in] v3 area of the spherical triangle -!! @param[out] area result +!! @return area result !! @author R. J. Purser function sarea_d(v1,v2,v3)result(area)! [sarea] use pietc, only: zero=>u0 @@ -771,7 +767,7 @@ end function sarea_d !! @param[in] drlona longitudes !! @param[in] drlatb latitudes !! @param[in] drlonb longitudes -!! @param[out] area result +!! @return area result !! @author R. J. Purser function dtarea_s(rlat,drlata,drlona,drlatb,drlonb) result(area)! [sarea] use pietc_s, only: u0,u1 @@ -804,7 +800,7 @@ end function dtarea_s !! @param[in] drlona longitudes !! @param[in] drlatb latitudes !! @param[in] drlonb longitudes -!! @param[out] area result +!! @return area result !! @author R. J. Purser function dtarea_d(rlat,drlata,drlona,drlatb,drlonb) result(area)! [sarea] use pietc, only: u0,u1 @@ -846,7 +842,7 @@ end function dtarea_d !! @param[in] drlonb longitudes !! @param[in] drlatc latitudes !! @param[in] drlonc longitudes -!! @param[out] area result +!! @return area result !! @author R. J. Purser function dqarea_s &! [sarea] (rlat,drlata,drlona,drlatb,drlonb,drlatc,drlonc) result(area) @@ -866,7 +862,7 @@ end function dqarea_s !! @param[in] drlonb longitudes !! @param[in] drlatc latitudes !! @param[in] drlonc longitudes -!! @param[out] area result +!! @return area !! @author R. J. Purser function dqarea_d &! [sarea] (rlat,drlata,drlona,drlatb,drlonb,drlatc,drlonc) result(area) @@ -914,7 +910,7 @@ end subroutine dlltoxy_d !> Haversine function. !! !! @param[in] t input -!! @param[out] a result +!! @return a result !! @author R. J. Purser function hav_s(t) result(a)! [hav] use pietc_s, only: o2 @@ -927,7 +923,7 @@ end function hav_s !> Doing hav_d function. !! !! @param[in] t input -!! @param[out] a result +!! @return a result !! @author R. J. Purser function hav_d(t) result(a)! [hav] use pietc, only: o2 @@ -962,6 +958,11 @@ subroutine normalize_d(v)! [normalize] end subroutine normalize_d !> ??? +!! +!! @param[in] as ??? +!! @param[out] b ??? +!! @param[out] nrank ??? +!! @param[out] det ??? !! @author R. J. Purser subroutine gram_s(as,b,nrank,det)! [gram] use pietc_s, only: u0,u1 @@ -1024,6 +1025,11 @@ subroutine gram_s(as,b,nrank,det)! [gram] end subroutine gram_s !> ??? +!! +!! @param[in] as ??? +!! @param[out] b ??? +!! @param[out] nrank ??? +!! @param[out] det ??? !! @author R. J. Purser subroutine gram_d(as,b,nrank,det)! [gram] use pietc, only: u0,u1 @@ -1091,7 +1097,11 @@ end subroutine gram_d !! as zero (instead of either +1 or -1) and ldet is then just the log of !! the nonzero factors found by the process. !! +!! @param[in] as ??? +!! @param[out] b ??? +!! @param[out] nrank ??? !! @param[out] detsign singular determinant +!! @param[out] ldet ??? !! @author R. J. Purser subroutine graml_d(as,b,nrank,detsign,ldet)! [gram] use pietc, only: u0 @@ -1248,10 +1258,13 @@ end subroutine plaingram_d !! "epsilon" value that is fixed (10**(-13)) and assumes elements of a are !! never too different in magnitude from unity, unless they are actually zero. !! +!! @param[in] m ??? +!! @param[in] n ??? !! @param[in] a rectangular input matrix !! @param[out] ipiv pivoting sequence !! @param[out] tt row-normalization !! @param[out] b orthonormalized rows +!! @param[in] rank ??? !! @author R. J. Purser subroutine rowgram(m,n,a,ipiv,tt,b,rank)! [gram] use pietc, only: u0,u1 @@ -1331,6 +1344,10 @@ end subroutine rowgram !> Apply the row-operations, implied by ipiv and tt returned by rowgram, to !! the single column vector, v, to produce the transformed vector vv. !! +!! @param[in] m ??? +!! @param[in] n ??? +!! @param[in] ipiv ??? +!! @param[in] tt ??? !! @param[in] v vector !! @param[out] vv vector !! @author R. J. Purser @@ -1367,8 +1384,12 @@ end subroutine rowops !! together with the rescaled matrix aa such that a = d.aa.e when d and e are !! interpreted as diagonal matrices. !! +!! @param[in] m ??? +!! @param[in] n ??? +!! @param[in] mask ??? !! @param[in] a positive diagonals !! @param[out] d positive diagonals +!! @param[in] aa ??? !! @param[out] e positive diagonals !! @author R. J. Purser subroutine corral(m,n,mask,a,d,aa,e)! [corral] @@ -1612,6 +1633,12 @@ subroutine qtoax(q,v)! [qtoax] end subroutine qtoax !> ??? +!! +!! @param[in] c ??? +!! @param[in] d ??? +!! @param[in] e ??? +!! @param[in] g ??? +!! @param[in] r ??? !! @author R. J. Purser subroutine setem(c,d,e,g,r)! [setem] implicit none @@ -1630,8 +1657,8 @@ end subroutine setem !! !! @param[in] a real quaternion !! @param[in] b real quaternion -!! @param[out] c result !! @author R. J. Purser +!! @return c function mulqq(a,b)result(c)! [mulqq] implicit none real(dp),dimension(0:3),intent(IN ):: a,b @@ -1650,6 +1677,7 @@ end function mulqq !! @param[in] n degree !! @param[in] a matrix !! @param[out] b exponential matrix of a +!! @param[out] detb ??? !! @author R. J. Purser subroutine expmat(n,a,b,detb)! [expmat] use pietc, only: u0,u1,u2,o2 @@ -1685,6 +1713,9 @@ end subroutine expmat !! @param[in] n degree !! @param[in] a matrix !! @param[out] b exponential matrix of a +!! @param[out] bd ??? +!! @param[out] detb ??? +!! @param[out] detbd ??? !! @author R. J. Purser subroutine expmatd(n,a,b,bd,detb,detbd)! [expmat] use pietc, only: u0,u1,u2,o2 @@ -1747,6 +1778,11 @@ end subroutine expmatd !! @param[in] n degree !! @param[in] a matrix !! @param[out] b exponential matrix of a +!! @param[out] bd ??? +!! @param[out] bdd ??? +!! @param[out] detb ??? +!! @param[out] detbd ??? +!! @param[out] detbdd ??? !! @author R. J. Purser subroutine expmatdd(n,a,b,bd,bdd,detb,detbd,detbdd)! [expmat] use pietc, only: u0,u1,u2,o2 @@ -1829,6 +1865,10 @@ subroutine expmatdd(n,a,b,bd,bdd,detb,detbd,detbdd)! [expmat] end subroutine expmatdd !> ??? +!! +!! @param[in] n ??? +!! @param[in] z ??? +!! @param[in] zn ??? !! @author R. J. Purser subroutine zntay(n,z,zn)! [zntay] use pietc, only: u2 @@ -1859,6 +1899,13 @@ subroutine zntay(n,z,zn)! [zntay] end subroutine zntay !> ??? +!! +!! @param[in] n ??? +!! @param[in] z ??? +!! @param[out] zn ??? +!! @param[out] znd ??? +!! @param[out] zndd ??? +!! @param[out] znddd ??? !! @author R. J. Purser subroutine znfun(n,z,zn,znd,zndd,znddd)! [znfun] use pietc, only: u0,u2,u3 @@ -1911,13 +1958,17 @@ end subroutine znfun !! aa2 etc to zv, is equivalent to a single mapping zz-->zv by the transformatn !! with coefficients aa3,bb3,cc3,dd3, such that, as 2*2 complex matrices: !! +!! 
 !! [ aa3, bb3 ]   [ aa2, bb2 ]   [ aa1, bb1 ]
 !! [          ] = [          ] * [          ]
 !! [ cc3, dd3 ]   [ cc2, dd2 ]   [ cc1, dd1 ] .
+!!
!! -!! Note that the determinant of these matrices is always +1 +!! Note that the determinant of these matrices is always +1. !! !! @param[in] v matric +!! @param[out] z ??? +!! @param[out] infz ??? !! @author R. J. Purser subroutine ctoz(v, z,infz)! [ctoz] use pietc, only: u0,u1 @@ -1939,6 +1990,10 @@ subroutine ctoz(v, z,infz)! [ctoz] end subroutine ctoz !> ??? +!! +!! @param[in] z ??? +!! @param[in] infz ??? +!! @param[in] v ??? !! @author R. J. Purser subroutine ztoc(z,infz, v)! [ztoc] implicit none @@ -1961,10 +2016,14 @@ end subroutine ztoc !! change of cartesian vector position is delta_v given by: !! delta_v = Real(vd*delta_z). !! Thus, by a kind of Cauchy-Riemann relation, Imag(vd)=v CROSS Real(vd). -!! THE DERIVATIVE FOR THE IDEAL POINT AT INFINITY HAS NOT BEEN CODED YET!!! +!! +!! @note The derivative for the ideal point at infinity has not been +!! coded yet. !! !! @param[in] z complex infinitesimal map displacement +!! @param[in] infz ??? !! @param[out] v cartesian vector position +!! @param[out] vd ??? !! @author R. J. Purser subroutine ztocd(z,infz, v,vd)! [ztoc] implicit none @@ -2214,7 +2273,10 @@ end subroutine cmobius !! @param[in] bb Mobius transformation coefficients !! @param[in] cc Mobius transformation coefficients !! @param[in] dd Mobius transformation coefficients +!! @param[out] zz ??? +!! @param[out] infz ??? !! @param[out] zw Inversed output +!! @param[out] infw ??? !! @author R. J. Purser subroutine zmobiusi(aa,bb,cc,dd, zz,infz, zw,infw) ! [mobiusi] implicit none