diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90
index 3c5c1e4fb..3316b2460 100644
--- a/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90
+++ b/sorc/grid_tools.fd/regional_esg_grid.fd/pesg.f90
@@ -431,7 +431,7 @@ end subroutine xmtoxc_vak1
 !! @param a ESG mapping parameterization
 !! @param k ESG mapping parameterization
 !! @param xm map-space vector
-!! @param xm derivative
+!! @param xc derivative
 !! @param xcd jacobian matrix
 !! @param ff error flag
 !! @author R. J. Purser
@@ -451,7 +451,7 @@ subroutine xmtoxc_ak(a,k,xm,xc,xcd,ff)!                            [xmtoxc_ak]
 xcd=matmul(xcd,xsd)
 end subroutine xmtoxc_ak
 
-!! Inverse mapping of xmtoxc_ak. That is, go from given cartesian unit
+!> Inverse mapping of xmtoxc_ak. That is, go from given cartesian unit
 !! 3-vector, xc, to map coordinate 2-vector xm (or return a raised
 !! failure flag, FF, if the attempt fails).
 !!
@@ -1598,7 +1598,7 @@ subroutine hgrid_ak_rc(lx,ly,nx,ny,A,K,plat,plon,pazi, & !       [hgrid_ak_rc]
 enddo
 end subroutine hgrid_ak_rc
 
-!! Use a and k as the parameters of an Extended Schmidt-transformed
+!> Use a and k as the parameters of an Extended Schmidt-transformed
 !! Gnomonic (ESG) mapping centered at (pdlat,pdlon) and twisted about
 !! this center by an azimuth angle of pdazi counterclockwise (these
 !! angles in degrees).
@@ -1644,7 +1644,7 @@ subroutine hgrid_ak_dd(lx,ly,nx,ny,a,k,pdlat,pdlon,pdazi, & !    [hgrid_ak_dd]
 gdlon=gdlon*rtod !
 end subroutine hgrid_ak_dd
 
-!! Like hgrid_ak_rr_c, except all the angle arguments (but not
+!> Like hgrid_ak_rr_c, except all the angle arguments (but not
 !! delx,dely) are in degrees instead of radians.
 !!
 !! @param[in] lx x grid index for left-lower corner of the grid at center
@@ -1973,6 +1973,8 @@ end subroutine gtoxm_ak_dd_m
 !! @param[in] pdazi ???
 !! @param[in] delx central x-spacing grid point
 !! @param[in] dely central y-spacing grid point
+!! @param[in] dlat ???
+!! @param[in] dlon ???
 !! @param[out] xm ???
 !! @param[out] ff failure flag
 !! @author R. J. Purser
diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/pfun.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/pfun.f90
index e5415ef5c..dbaad10fc 100644
--- a/sorc/grid_tools.fd/regional_esg_grid.fd/pfun.f90
+++ b/sorc/grid_tools.fd/regional_esg_grid.fd/pfun.f90
@@ -1,8 +1,10 @@
 !> @file
+!! ???
 !! @author R. J. Purser
-!! Direct dependencies:
-!! Modules: pkind, pietc_s, pietc
+
+!> This module is for ???
 !!
+!! @author R. J. Purser
 module pfun
 !=============================================================================
 use pkind, only: sp,dp
@@ -27,55 +29,74 @@ module pfun
 
 contains
 
-!=============================================================================
+!> ???
+!!
+!! @param[in] x ???
+!! @return y ???
+!! @author R. J. Purser  
 function gd_s(x) result(y)!                                               [gd]
-!=============================================================================
 ! Gudermannian function
 implicit none
 real(sp),intent(in ):: x
 real(sp)            :: y
 y=atan(sinh(x))
 end function gd_s
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] x ???
+!! @return y ???
+!! @author R. J. Purser  
 function gd_d(x) result(y)!                                               [gd]
-!=============================================================================
 implicit none
 real(dp),intent(in ):: x
 real(dp)            :: y
 y=atan(sinh(x))
 end function gd_d
 
-!=============================================================================
+!> Inverse Gudermannian function for single precision real.
+!!
+!! @param[in] y ???
+!! @return x ???
+!! @author R. J. Purser  
 function gdi_s(y) result(x)!                                             [gdi]
-!=============================================================================
-! Inverse Gudermannian function
 implicit none
 real(sp),intent(in ):: y
 real(sp)            :: x
 x=atanh(sin(y))
 end function gdi_s
-!=============================================================================
+
+!> Inverse Gudermannian function for double precision real.
+!!
+!! @param[in] y ???
+!! @return x ???
+!! @author R. J. Purser  
 function gdi_d(y) result(x)!                                             [gdi]
-!=============================================================================
 implicit none
 real(dp),intent(in ):: y
 real(dp)            :: x
 x=atanh(sin(y))
 end function gdi_d
 
-!=============================================================================
+!> Haversine function for single precision real.
+!!
+!! @param[in] t ???
+!! @return a ???
+!! @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
-!=============================================================================
+
+!>  Haversine function for double precision real.
+!!
+!! @param[in] t ???
+!! @return a ???
+!! @author R. J. Purser  
 function hav_d(t) result(a)!                                             [hav]
-!=============================================================================
 use pietc, only: o2
 implicit none
 real(dp),intent(in ):: t
@@ -83,19 +104,27 @@ function hav_d(t) result(a)!                                             [hav]
 a=(sin(t*o2))**2
 end function hav_d
 
-!=============================================================================
+!> Hyperbolic-haversine for single precision real.
+!!
+!! @note The minus sign in the hyperbolic-haversine definition.
+!!
+!! @param[in] t ???
+!! @return a ???
+!! @author R. J. Purser  
 function havh_s(t) result(a)!                                           [havh]
-!=============================================================================
-! Note the minus sign in the hyperbolic-haversine definition
 use pietc_s, only: o2
 implicit none
 real(sp),intent(in ):: t
 real(sp)            :: a
 a=-(sinh(t*o2))**2
 end function havh_s
-!=============================================================================
+
+!> Hyperbolic-haversine for double precision real.
+!!
+!! @param[in] t ???
+!! @return a ???
+!! @author R. J. Purser  
 function havh_d(t) result(a)!                                           [havh]
-!=============================================================================
 use pietc, only: o2
 implicit none
 real(dp),intent(in ):: t
@@ -103,19 +132,25 @@ function havh_d(t) result(a)!                                           [havh]
 a=-(sinh(t*o2))**2
 end function havh_d
 
-!=============================================================================
+!> Arc-haversine function for single precision real.
+!!
+!! @param[in] a ???
+!! @return t ???
+!! @author R. J. Purser  
 function ahav_s(a) result(t)!                                           [ahav]
-!=============================================================================
 use pietc_s, only: u2
-! Arc-haversine function
 implicit none
 real(sp),intent(in ):: a
 real(sp)            :: t
 t=u2*asin(sqrt(a))
 end function ahav_s
-!=============================================================================
+
+!> Arc-haversine function for double precision real.
+!!
+!! @param[in] a ???
+!! @return t ???
+!! @author R. J. Purser  
 function ahav_d(a) result(t)!                                           [ahav]
-!=============================================================================
 use pietc, only: u2
 implicit none
 real(dp),intent(in ):: a
@@ -123,19 +158,27 @@ function ahav_d(a) result(t)!                                           [ahav]
 t=u2*asin(sqrt(a))
 end function ahav_d
 
-!=============================================================================
+!> Hyperbolic arc-haversine for single precision real.
+!!
+!! @note The minus sign in the hyperbolic arc-haversine definition.
+!!
+!! @param[in] a ???
+!! @return t ???
+!! @author R. J. Purser  
 function ahavh_s(a) result(t)!                                         [ahavh]
-!=============================================================================
 use pietc_s, only: u2
-! Note the minus sign in the hyperbolic arc-haversine definition
 implicit none
 real(sp),intent(in ):: a
 real(sp)            :: t
 t=u2*asinh(sqrt(-a))
 end function ahavh_s
-!=============================================================================
+
+!> Hyperbolic arc-haversine for double precision real.
+!!
+!! @param[in] a ???
+!! @return t ???
+!! @author R. J. Purser  
 function ahavh_d(a) result(t)!                                         [ahavh]
-!=============================================================================
 use pietc, only: u2
 implicit none
 real(dp),intent(in ):: a
@@ -143,9 +186,12 @@ function ahavh_d(a) result(t)!                                         [ahavh]
 t=u2*asinh(sqrt(-a))
 end function ahavh_d
 
-!=============================================================================
+!> ???
+!!
+!! @param[in] t ???
+!! @return a ???
+!! @author R. J. Purser  
 function atanh_s(t) result(a)!                                         [atanh]
-!=============================================================================
 use pietc_s, only: u1,o2,o3,o5
 implicit none
 real(sp),intent(IN ):: t
@@ -157,9 +203,13 @@ function atanh_s(t) result(a)!                                         [atanh]
 else; tt=t*t;            a=t*(u1+tt*(o3+tt*(o5+tt*(o7+tt*o9))))
 endif
 end function atanh_s
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] t ???
+!! @return a ???
+!! @author R. J. Purser  
 function atanh_d(t) result(a)!                                         [atanh]
-!=============================================================================
 use pietc, only: u1,o2,o3,o5
 implicit none
 real(dp),intent(IN ):: t
@@ -172,9 +222,12 @@ function atanh_d(t) result(a)!                                         [atanh]
 endif
 end function atanh_d
 
-!=============================================================================
+!> ???
+!!
+!! @param[in] x ???
+!! @return t ???
+!! @author R. J. Purser  
 function sech_s(x)result(r)!                                            [sech]
-!=============================================================================
 ! This indirect way of computing 1/cosh(x) avoids overflows at large x
 use pietc_s, only: u1,u2
 implicit none
@@ -185,9 +238,13 @@ function sech_s(x)result(r)!                                            [sech]
 e=exp(-ax)
 r=e*u2/(u1+e*e)
 end function sech_s
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] x ???
+!! @return r ???
+!! @author R. J. Purser  
 function sech_d(x)result(r)!                                            [sech]
-!=============================================================================
 use pietc, only: u1,u2
 implicit none
 real(dp),intent(in ):: x
@@ -198,27 +255,36 @@ function sech_d(x)result(r)!                                            [sech]
 r=e*u2/(u1+e*e)
 end function sech_d
 
-!=============================================================================
+!> ???
+!!
+!! @param[in] x ???
+!! @return r ???
+!! @author R. J. Purser  
 function sechs_s(x)result(r)!                                          [sechs]
-!=============================================================================
 implicit none
 real(sp),intent(in ):: x
 real(sp)            :: r
 r=sech(x)**2
 end function sechs_s
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] x ???
+!! @return r ???
+!! @author R. J. Purser  
 function sechs_d(x)result(r)!                                          [sechs]
-!=============================================================================
 implicit none
 real(dp),intent(in ):: x
 real(dp)            :: r
 r=sech(x)**2
 end function sechs_d
 
-!=============================================================================
+!> Evaluate the symmetric real function sin(x)/x-1.
+!!
+!! @param[in] x ???
+!! @return r ???
+!! @author R. J. Purser  
 function sinoxm_d(x) result(r)!                                       [sinoxm]
-!=============================================================================
-! Evaluate the symmetric real function sin(x)/x-1
 use pietc, only: u1
 implicit none
 real(dp),intent(in ):: x
@@ -233,10 +299,12 @@ function sinoxm_d(x) result(r)!                                       [sinoxm]
 endif
 end function sinoxm_d
 
-!=============================================================================
+!> Evaluate the symmetric real function sin(x)/x.
+!!
+!! @param[in] x ???
+!! @return r ???
+!! @author R. J. Purser  
 function sinox_d(x) result(r)!                                         [sinox]
-!=============================================================================
-! Evaluate the symmetric real function sin(x)/x
 use pietc, only: u1
 implicit none
 real(dp),intent(in ):: x
@@ -245,10 +313,12 @@ function sinox_d(x) result(r)!                                         [sinox]
 r=sinoxm(x)+u1
 end function sinox_d
 
-!=============================================================================
+!> Evaluate the symmetric real function sinh(x)/x-1.
+!!
+!! @param[in] x ???
+!! @return r ???
+!! @author R. J. Purser  
 function sinhoxm_d(x) result(r)!                                     [sinhoxm]
-!=============================================================================
-! Evaluate the symmetric real function sinh(x)/x-1
 use pietc, only: u1
 implicit none
 real(dp),intent(in ):: x
@@ -263,10 +333,12 @@ function sinhoxm_d(x) result(r)!                                     [sinhoxm]
 endif
 end function sinhoxm_d
 
-!=============================================================================
+!> Evaluate the symmetric real function sinh(x)/x.
+!!
+!! @param[in] x ???
+!! @return r ???
+!! @author R. J. Purser  
 function sinhox_d(x) result(r)!                                       [sinhox]
-!=============================================================================
-! Evaluate the symmetric real function sinh(x)/x
 use pietc, only: u1
 implicit none
 real(dp),intent(in ):: x
diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/pkind.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/pkind.f90
index 054201fd6..e5879c2d5 100644
--- a/sorc/grid_tools.fd/regional_esg_grid.fd/pkind.f90
+++ b/sorc/grid_tools.fd/regional_esg_grid.fd/pkind.f90
@@ -1,10 +1,11 @@
 !> @file
 module pkind
-integer,parameter:: spi=selected_int_kind(6),&
-                    dpi=selected_int_kind(12),&
-                    sp =selected_real_kind(6,30),&
-                    dp =selected_real_kind(15,300),&
-                    spc=sp,dpc=dp
+integer,parameter:: spi=selected_int_kind(6) !< Single precision integer kind.
+integer,parameter:: dpi=selected_int_kind(12) !< Double precision integer kind.
+integer,parameter:: sp =selected_real_kind(6,30) !< Single precision real kind.
+integer,parameter:: dp =selected_real_kind(15,300) !< Double precision real kind.
+integer,parameter:: spc=sp  !< Single precision real kind.
+integer,parameter:: dpc=dp  !< Double precision real kind.
 !private:: one_dpi; integer(8),parameter:: one_dpi=1
 !integer,parameter:: dpi=kind(one_dpi)
 !integer,parameter:: sp=kind(1.0)
diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/pmat2.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/pmat2.f90
index ccb7c0ecb..34ec2b2ca 100644
--- a/sorc/grid_tools.fd/regional_esg_grid.fd/pmat2.f90
+++ b/sorc/grid_tools.fd/regional_esg_grid.fd/pmat2.f90
@@ -1,13 +1,14 @@
 !> @file
+!! @brief Routines dealing with the operations of banded matrices.
 !! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1994/1999
-!!
-!! Routines dealing with the operations of banded matrices
+
+!> Routines dealing with the operations of banded matrices.
 !! The three special routines allow the construction of compact or
 !! conventional interpolation and differencing stencils to a general
 !! order of accuracy. These are:
-!! - AVCO:  Averaging, or interpolating;
-!! - DFCO:  Differentiating (once);
-!! - DFCO2: Differentiating (twice).
+!! - avco()  Averaging, or interpolating;
+!! - dfco()  Differentiating (once);
+!! - dfco2() Differentiating (twice).
 !!
 !! Other routines provide the tools for applying compact schemes, and for
 !! the construction and application of recursive filters.
@@ -16,20 +17,41 @@
 !!  added nonredundant ldltb and ltdlbv routines for symmetric matrices,
 !!  and remove obsolescent routines.
 !!
-!! DIRECT DEPENDENCIES
-!! Libraries[their modules]: pmat[pmat]
-!! Additional Modules      : pkind
-!!
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1994/1999
 module pmat2
 !============================================================================
 use    pkind, only: spi,sp,dp,dpc
 implicit none
 private
-public:: avco,dfco,dfco2, clipb,cad1b,csb1b,cad2b,csb2b,                    &
-         ldub,ldltb,udlb,l1ubb,l1ueb,ltdlbv,                                &
-         udlbv,udlbx,udlby,udlvb,udlxb,udlyb,u1lbv,u1lbx,u1lby,u1lvb,u1lxb, &
-         u1lyb,linbv,wrtb
-real(dp),parameter:: zero=0
+public:: avco !< ???
+public:: dfco !< ???
+public:: dfco2 !< ???
+public:: clipb !< ???
+public:: cad1b !< ???
+public:: csb1b !< ???
+public:: cad2b !< ???
+public:: csb2b !< ???
+public:: ldub !< ???
+public:: ldltb !< ???
+public:: udlb !< ???
+public:: l1ubb !< ???
+public:: l1ueb !< ???
+public:: ltdlbv !< ???
+public:: udlbv !< ???
+public:: udlbx !< ???
+public:: udlby !< ???
+public:: udlvb !< ???
+public:: udlxb !< ???
+public:: udlyb !< ???
+public:: u1lbv !< ???
+public:: u1lbx !< ???
+public:: u1lby !< ???
+public:: u1lvb !< ???
+public:: u1lxb !< ???
+public:: u1lyb !< ???
+public:: linbv !< ???
+public:: wrtb !< ???
+real(dp),parameter:: zero=0  !< ???
 
 interface AVCO;   module procedure AVCO,  DAVCO,  TAVCO;   end interface
 interface DFCO;   module procedure DFCO,  DDFCO,  TDFCO;   end interface
@@ -61,26 +83,20 @@ module pmat2
 interface WRTB;   module procedure WRTB;                   end interface
 contains
 
-!=============================================================================
+!>  Compute one row of the coefficients for the compact mid-interval
+!!  interpolation scheme characterized by matrix equation of the form,
+!!			 A.t = B.s			       (*)
+!!  Where s is the vector of "source" values, t the staggered "target" values.
+!!
+!! @param[in] na number of t-points operated on by this row of the A of (*)
+!! @param[in] nb number of s-points operated on by this row of the B of (*)
+!! @param[in] za coordinates of t-points used in this row of (*)
+!! @param[in] zb coordinates of s-points used in this row of (*)
+!! @param[in] z0 nominal point of application of this row of (*)
+!! @param[out] a the NA coefficients A for this scheme
+!! @param[out] b the NB coefficients B for this scheme
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine AVCO(na,nb,za,zb,z0,a,b) !                                   [AVCO]
-!=============================================================================
-!		    SUBROUTINE AVCO
-!   R.J.Purser, National Centers for Environmental Prediction, Washington D.C.
-!   jim.purser@noaa.gov					      1999
-!
-!  Compute one row of the coefficients for the compact mid-interval
-!  interpolation scheme characterized by matrix equation of the form,
-!			 A.t = B.s			       (*)
-!  Where s is the vector of "source" values, t the staggered "target" values.
-!
-! --> NA:   number of t-points operated on by this row of the A of (*)
-! --> NB:   number of s-points operated on by this row of the B of (*)
-! --> ZA:   coordinates of t-points used in this row of (*)
-! --> ZB:   coordinates of s-points used in this row of (*)
-! --> Z0:   nominal point of application of this row of (*)
-! <-- A:    the NA coefficients A for this scheme
-! <-- B:    the NB coefficients B for this scheme
-!=============================================================================
 use pietc, only: u0,u1
 use pmat,  only: inv
 implicit none
@@ -103,9 +119,18 @@ subroutine AVCO(na,nb,za,zb,z0,a,b) !                                   [AVCO]
 call INV(w,ab)
 a=ab(1:na); b=ab(na1:nab)
 end subroutine AVCO 
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] na ???
+!! @param[in] nb ???
+!! @param[in] za ???
+!! @param[in] zb ???
+!! @param[in] z0 ???
+!! @param[in] a ???
+!! @param[in] b ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine DAVCO(na,nb,za,zb,z0,a,b) !                                  [AVCO]
-!=============================================================================
 use pietc, only: u0,u1
 use pmat, only: inv
 implicit none
@@ -128,9 +153,15 @@ subroutine DAVCO(na,nb,za,zb,z0,a,b) !                                  [AVCO]
 call INV(w,ab)
 a=ab(1:na); b=ab(na1:nab)
 end subroutine DAVCO
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] xa ???
+!! @param[in] xb ???
+!! @param[in] a ???
+!! @param[in] b ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine TAVCO(xa,xb,a,b)!                                            [AVCO]
-!=============================================================================
 implicit none
 real(dp),dimension(:),intent(IN ):: xa,xb
 real(dp),dimension(:),intent(OUT):: a,b
@@ -142,26 +173,20 @@ subroutine TAVCO(xa,xb,a,b)!                                            [AVCO]
 call DAVCO(na,nb,xa,xb,zero,a,b)
 end subroutine TAVCO
 
-!=============================================================================
+!>  Compute one row of the coefficients for either the compact differencing or
+!!  quadrature scheme characterized by matrix equation of the form,
+!!			 A.d = B.c			       (*)
+!!  In either case, d is the derivative of c.
+!!
+!! @param[in] na number of d-points operated on by this row of the A of (*)
+!! @param[in] nb number of c-points operated on by this row of the B of (*)
+!! @param[in] za coordinates of d-points used in this row of (*)
+!! @param[in] zb coordinates of c-points used in this row of (*)
+!! @param[in] z0 nominal point of application of this row of (*)
+!! @param[in] A the A-coefficients for this scheme
+!! @param[in] B the B-coefficients for this scheme
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine DFCO(na,nb,za,zb,z0,a,b)!                                    [DFCO]
-!=============================================================================
-!   R.J.Purser, National Centers for Environmental Prediction, Washington D.C.
-!   jim.purser@noaa.gov					      1999
-!		    SUBROUTINE DFCO
-!
-!  Compute one row of the coefficients for either the compact differencing or
-!  quadrature scheme characterized by matrix equation of the form,
-!			 A.d = B.c			       (*)
-!  In either case, d is the derivative of c.
-!
-! --> NA:   number of d-points operated on by this row of the A of (*)
-! --> NB:   number of c-points operated on by this row of the B of (*)
-! --> ZA:   coordinates of d-points used in this row of (*)
-! --> ZB:   coordinates of c-points used in this row of (*)
-! --> Z0:   nominal point of application of this row of (*)
-! <-- A:    the A-coefficients for this scheme
-! <-- B:    the B-coefficients for this scheme
-!=============================================================================
 use pietc_s, only: u0,u1
 use pmat,    only: inv
 implicit none
@@ -185,9 +210,18 @@ subroutine DFCO(na,nb,za,zb,z0,a,b)!                                    [DFCO]
 call INV(w,ab)
 a=ab(1:na); b=ab(na1:nab)
 end subroutine DFCO 
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] na ???
+!! @param[in] nb ???
+!! @param[in] za ???
+!! @param[in] zb ???
+!! @param[in] z0 ???
+!! @param[in] a ???
+!! @param[in] b ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine DDFCO(na,nb,za,zb,z0,a,b) ! Real(dp) version of              [DFCO]
-!=============================================================================
 use pietc, only: u0,u1
 use pmat,  only: inv
 implicit none
@@ -211,9 +245,15 @@ subroutine DDFCO(na,nb,za,zb,z0,a,b) ! Real(dp) version of              [DFCO]
 call INV(w,ab)
 a=ab(1:na); b=ab(na1:nab)
 end subroutine DDFCO 
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] xa ???
+!! @param[in] xb ???
+!! @param[in] a ???
+!! @param[in] b ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine TDFCO(xa,xb,a,b)!                                            [DFCO]
-!=============================================================================
 implicit none
 real(dp),dimension(:),intent(IN ):: xa,xb
 real(dp),dimension(:),intent(OUT):: a,b
@@ -225,26 +265,20 @@ subroutine TDFCO(xa,xb,a,b)!                                            [DFCO]
 call DDFCO(na,nb,xa,xb,zero,a,b)
 end subroutine TDFCO
 
-!=============================================================================
+!>  Compute one row of the coefficients for either the compact second-
+!!  differencing scheme characterized by matrix equation of the form,
+!!			 A.d = B.c			       (*)
+!!  Where d is the second-derivative of c.
+!!
+!! @param[in] na number of d-points operated on by this row of the A of (*)
+!! @param[in] nb number of c-points operated on by this row of the B of (*)
+!! @param[in] za coordinates of d-points used in this row of (*)
+!! @param[in] zb coordinates of c-points used in this row of (*)
+!! @param[in] z0 nominal point of application of this row of (*)
+!! @param[in] a the NA coefficients A for this scheme
+!! @param[in] b the NB coefficients B for this scheme
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine DFCO2(na,nb,za,zb,z0,a,b)!                                  [DFCO2] 
-!=============================================================================
-!		    SUBROUTINE DFCO2
-!   R.J.Purser, National Centers for Environmental Prediction, Washington D.C.
-!   jim.purser@noaa.gov					      1999
-!
-!  Compute one row of the coefficients for either the compact second-
-!  differencing scheme characterized by matrix equation of the form,
-!			 A.d = B.c			       (*)
-!  Where d is the second-derivative of c.
-!
-! --> NA:   number of d-points operated on by this row of the A of (*)
-! --> NB:   number of c-points operated on by this row of the B of (*)
-! --> ZA:   coordinates of d-points used in this row of (*)
-! --> ZB:   coordinates of c-points used in this row of (*)
-! --> Z0:   nominal point of application of this row of (*)
-! <-- A:    the NA coefficients A for this scheme
-! <-- B:    the NB coefficients B for this scheme
-!=============================================================================
 use pietc_s, only: u0,u1
 use pmat, only: inv
 implicit none
@@ -268,9 +302,18 @@ subroutine DFCO2(na,nb,za,zb,z0,a,b)!                                  [DFCO2]
 call INV(w,ab)
 a=ab(1:na); b=ab(na1:nab)
 end subroutine DFCO2 
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] na ???
+!! @param[in] nb ???
+!! @param[in] za ???
+!! @param[in] zb ???
+!! @param[in] z0 ???
+!! @param[in] a ???
+!! @param[in] b ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine DDFCO2(na,nb,za,zb,z0,a,b) ! Real(dp) version of            [DFCO2]
-!=============================================================================
 use pietc, only: u0,u1
 use pmat, only: inv
 implicit none
@@ -294,7 +337,14 @@ subroutine DDFCO2(na,nb,za,zb,z0,a,b) ! Real(dp) version of            [DFCO2]
 call INV(w,ab)
 a=ab(1:na); b=ab(na1:nab)
 end subroutine ddfco2 
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] xa ???
+!! @param[in] xb ???
+!! @param[in] a ???
+!! @param[in] b ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine TDFCO2(xa,xb,a,b)!                                          [DFCO2]
 !=============================================================================
 real(dp),dimension(:),intent(IN ):: xa,xb
@@ -307,10 +357,15 @@ subroutine TDFCO2(xa,xb,a,b)!                                          [DFCO2]
 call DDFCO2(na,nb,xa,xb,zero,a,b)
 end subroutine TDFCO2
 
-
-!=============================================================================
+!> ???
+!!
+!! @param[in] m1 ???
+!! @param[in] m2 ???
+!! @param[in] mah1 ???
+!! @param[in] mah2 ???
+!! @param[in] a ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 pure subroutine CLIB(m1,m2,mah1,mah2,a)!                               [CLIPB]
-!=============================================================================
 use pietc_s, only: u0
 implicit none
 integer(spi), intent(IN   ) :: m1, m2, mah1, mah2
@@ -319,9 +374,16 @@ pure subroutine CLIB(m1,m2,mah1,mah2,a)!                               [CLIPB]
 do j=1,mah1; a(1:min(m1,j),-j)=u0; enddo
 do j=m2-m1+1,mah2; a(max(1,m2-j+1):m1,j)=u0; enddo
 end subroutine CLIB
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] m1 ???
+!! @param[in] m2 ???
+!! @param[in] mah1 ???
+!! @param[in] mah2 ???
+!! @param[in] a ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 pure subroutine clib_d(m1,m2,mah1,mah2,a)!                             [CLIPB]
-!=============================================================================
 use pietc, only: u0
 implicit none
 integer(spi),intent(IN   ) :: m1, m2, mah1, mah2
@@ -330,9 +392,16 @@ pure subroutine clib_d(m1,m2,mah1,mah2,a)!                             [CLIPB]
 do j=1,mah1; a(1:min(m1,j),-j)=u0; enddo
 do j=m2-m1+1,mah2; a(max(1,m2-j+1):m1,j)=u0; enddo
 end subroutine clib_d
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] m1 ???
+!! @param[in] m2 ???
+!! @param[in] mah1 ???
+!! @param[in] mah2 ???
+!! @param[in] a ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 pure subroutine clib_c(m1,m2,mah1,mah2,a)!                              [CLIPB]
-!=============================================================================
 use pietc, only: c0
 implicit none
 integer(spi), intent(IN   ) :: m1, m2, mah1, mah2
@@ -342,18 +411,17 @@ pure subroutine clib_c(m1,m2,mah1,mah2,a)!                              [CLIPB]
 do j=m2-m1+1,mah2; a(max(1,m2-j+1):m1,j)=c0; enddo
 end subroutine clib_c
 
-!=============================================================================
+!> Incorporate operand symmetry near end-1 of a band matrix operator.
+!!
+!! @param[in] m1 Size of implied full matrix
+!! @param[in] mah1 Left semi-bandwidths of A.
+!! @param[in] mah2 Right semi-bandwidths of A.
+!! @param[in] mirror2 2*location of symmetry axis relative to end-1 operand element.
+!! @param[inout] a Input as unclipped operator, output as symmetrized and clipped.
+!!      Note: although m2 is not used here, it IS used in companion routines
+!!            cad2b and csb2b; it is retained in the interests of uniformity.
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine CAD1B(m1,mah1,mah2,mirror2,a)!                              [CAD1B]
-!=============================================================================
-! Incorporate operand symmetry near end-1 of a band matrix operator
-!
-! <-> A:      Input as unclipped operator, output as symmetrized and clipped.
-! m1, m2:     Sizes of implied full matrix
-! mah1, mah2: Left and right semi-bandwidths of A.
-! mirror2:    2*location of symmetry axis relative to end-1 operand element.
-!      Note: although m2 is not used here, it IS used in companion routines
-!            cad2b and csb2b; it is retained in the interests of uniformity.
-!=============================================================================
 use pietc_s, only: u0
 implicit none
 integer(spi),intent(IN   ):: m1,mah1,mah2,mirror2
@@ -370,11 +438,15 @@ subroutine CAD1B(m1,mah1,mah2,mirror2,a)!                              [CAD1B]
 enddo
 end subroutine CAD1B
 
-!=============================================================================
+!> Like cad1b, but for antisymmetric operand.
+!!
+!! @param[in] m1 ???
+!! @param[in] mah1 ???
+!! @param[in] mah2 ???
+!! @param[in] mirror2 ???
+!! @param[in] a ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine CSB1B(m1,mah1,mah2,mirror2,a)!                              [CSB1B]
-!=============================================================================
-! Like cad1b, but for antisymmetric operand
-!=============================================================================
 use pietc_s, only: u0
 implicit none
 integer(spi),intent(IN   ):: m1,mah1,mah2,mirror2
@@ -391,16 +463,16 @@ subroutine CSB1B(m1,mah1,mah2,mirror2,a)!                              [CSB1B]
 enddo
 end subroutine CSB1B
 
-!=============================================================================
+!> Incorporate operand symmetry near end-2 of a band matrix operator
+!!
+!! @param[in] m1 Sizes of implied full matrix
+!! @param[in] m2 Sizes of implied full matrix
+!! @param[in] mah1 Left semi-bandwidths of A.
+!! @param[in] mah2 Right semi-bandwidths of A.
+!! @param[in] mirror2 2*location of symmetry axis relative to end-2 operand element.
+!! @param[inout] a Input as unclipped operator, output as symmetrized and clipped.
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine CAD2B(m1,m2,mah1,mah2,mirror2,a)!                           [CAD2B]
-!=============================================================================
-! Incorporate operand symmetry near end-2 of a band matrix operator
-!
-! <-> A:      Input as unclipped operator, output as symmetrized and clipped.
-! m1, m2:     Sizes of implied full matrix
-! mah1, mah2: Left and right semi-bandwidths of A.
-! mirror2:    2*location of symmetry axis relative to end-2 operand element.
-!=============================================================================
 use pietc_s, only: u0
 implicit none
 integer(spi),intent(IN   ):: m1,m2,mah1,mah2,mirror2
@@ -418,9 +490,16 @@ subroutine CAD2B(m1,m2,mah1,mah2,mirror2,a)!                           [CAD2B]
 enddo
 end subroutine CAD2B
 
-!=============================================================================
+!> ???
+!!
+!! @param[in] m1 ???
+!! @param[in] m2 ???
+!! @param[in] mah1 ???
+!! @param[in] mah2 ???
+!! @param[in] mirror2 ???
+!! @param[in] a ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine CSB2B(m1,m2,mah1,mah2,mirror2,a)!                           [CSB2B]
-!=============================================================================
 use pietc_s, only: u0
 implicit none
 integer(spi),intent(IN   ):: m1,m2,mah1,mah2,mirror2
@@ -438,22 +517,18 @@ subroutine CSB2B(m1,m2,mah1,mah2,mirror2,a)!                           [CSB2B]
 enddo
 end subroutine CSB2B
 
-!=============================================================================
+!> Compute [L]*[D**-1]*[U] decomposition of asymmetric band-matrix.
+!!
+!! @param[inout] a input as the asymmetric band matrix. On output, it contains
+!! the [L]*[D**-1]*[U] factorization of the input matrix, where
+!!     - [L] is lower triangular with unit main diagonal
+!!     - [D] is a diagonal matrix
+!!     - [U] is upper triangular with unit main diagonal
+!! @param[in] m The number of rows of array A
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[in] mah2 the right half-bandwidth of fortran array A
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1994
 subroutine LDUB(m,mah1,mah2,a)!                                         [LDUB]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
-!		    SUBROUTINE LDUB
-!  Compute [L]*[D**-1]*[U] decomposition of asymmetric band-matrix
-!
-! <-> A: input as the asymmetric band matrix. On output, it contains
-!     the [L]*[D**-1]*[U] factorization of the input matrix, where
-!     [L] is lower triangular with unit main diagonal
-!     [D] is a diagonal matrix
-!     [U] is upper triangular with unit main diagonal
-! --> M:    The number of rows of array A
-! --> MAH1: the left half-bandwidth of fortran array A
-! --> MAH2: the right half-bandwidth of fortran array A
-!=============================================================================
 use pietc_s, only: u0,u1
 implicit none
 integer(spi),intent(IN   ):: m,mah1, mah2 
@@ -481,9 +556,15 @@ subroutine LDUB(m,mah1,mah2,a)!                                         [LDUB]
    a(j,1:jmost-j)=ajji*a(j,1:jmost-j)
 enddo
 end subroutine LDUB
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] m ???
+!! @param[in] mah1 ???
+!! @param[in] mah2 ???
+!! @param[in] a ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine DLDUB(m,mah1,mah2,a) ! Real(dp) version of                   [LDUB]
-!=============================================================================
 use pietc, only: u0,u1
 implicit none
 integer(spi),intent(IN   ):: m,mah1, mah2 
@@ -512,9 +593,13 @@ subroutine DLDUB(m,mah1,mah2,a) ! Real(dp) version of                   [LDUB]
 enddo
 end subroutine DLDUB
 
-!=============================================================================
+!> ???
+!!
+!! @param[in] m ???
+!! @param[in] mah1 ???
+!! @param[in] a ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine LDLTB(m,mah1,a) ! Real(sp) version of                       [LDLTB]
-!=============================================================================
 use pietc_s, only: u0,u1
 integer(spi),intent(IN   ):: m,mah1
 real(sp),    intent(INOUT):: a(m,-mah1:0) 
@@ -541,9 +626,14 @@ subroutine LDLTB(m,mah1,a) ! Real(sp) version of                       [LDLTB]
   enddo
 enddo
 end subroutine LDLTB
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] m ???
+!! @param[in] mah1 ???
+!! @param[in] a ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine DLDLTB(m,mah1,a) ! Real(dp) version of                   [LDLTB]
-!=============================================================================
 use pietc, only: u0,u1
 integer(spi),intent(IN   ) :: m,mah1
 real(dp),    intent(INOUT) :: a(m,-mah1:0) 
@@ -571,9 +661,14 @@ subroutine DLDLTB(m,mah1,a) ! Real(dp) version of                   [LDLTB]
 enddo
 end subroutine DLDLTB
 
-!=============================================================================
+!> ???
+!!
+!! @param[in] m ???
+!! @param[in] mah1 ???
+!! @param[in] mah2 ???
+!! @param[in] a ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine UDLB(m,mah1,mah2,a) ! Reversed-index version of ldub        [UDLB]
-!=============================================================================
 implicit none
 integer(spi),                    intent(IN   ) :: m,mah1,mah2
 real(sp),dimension(m,-mah1:mah2),intent(INOUT) :: a(m,-mah1:mah2)
@@ -583,9 +678,15 @@ subroutine UDLB(m,mah1,mah2,a) ! Reversed-index version of ldub        [UDLB]
 at=a(m:1:-1,mah2:-mah1:-1); call LDUB(m,mah2,mah1,at)
 a=at(m:1:-1,mah1:-mah2:-1)
 end subroutine UDLB 
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] m ???
+!! @param[in] mah1 ???
+!! @param[in] mah2 ???
+!! @param[in] a ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine DUDLB(m,mah1,mah2,a) ! real(dp) version of udlb              [UDLB]
-!=============================================================================
 implicit none
 integer(spi),                    intent(IN   ) :: m,mah1,mah2
 real(dp),dimension(m,-mah1:mah2),intent(INOUT) :: a(m,-mah1:mah2)
@@ -596,25 +697,21 @@ subroutine DUDLB(m,mah1,mah2,a) ! real(dp) version of udlb              [UDLB]
 a=at(m:1:-1,mah1:-mah2:-1)
 end subroutine DUDLB 
 
-!=============================================================================
+!>  Form the [L]*[D]*[U] decomposition of asymmetric band-matrix  [A] replace
+!!  lower triangular elements of [A] by [D**-1]*[L]*[D], the upper by [U],
+!!  replace matrix [B] by [D**-1]*[B].
+!!
+!! @param[inout] a input as band matrix, output as lower and upper triangulars with 1s
+!! implicitly assumed to lie on the main diagonal. The product of these
+!! triangular matrices is [D**-1]*[A], where [D] is a diagonal matrix.
+!! @param[inout] b in as band matrix, out as same but premultiplied by diagonal [D**-1]
+!! @param[in] m Number of rows of A and B
+!! @param[in] mah1 left half-width of fortran array A
+!! @param[in] mah2 right half-width of fortran array A
+!! @param[in] mbh1 left half-width of fortran array B
+!! @param[in] mbh2 right half-width of fortran array B
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1996
 subroutine L1UBB(m,mah1,mah2,mbh1,mbh2,a,b)!                           [L1UBB] 
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1996
-!		    SUBROUTINE L1UBB
-!  Form the [L]*[D]*[U] decomposition of asymmetric band-matrix  [A] replace
-!  lower triangular elements of [A] by [D**-1]*[L]*[D], the upper by [U],
-!  replace matrix [B] by [D**-1]*[B].
-!
-! <-> A input as band matrix, output as lower and upper triangulars with 1s
-!     implicitly assumed to lie on the main diagonal. The product of these
-!     triangular matrices is [D**-1]*[A], where [D] is a diagonal matrix.
-! <-> B in as band matrix, out as same but premultiplied by diagonal [D**-1]
-! --> M    Number of rows of A and B
-! --> MAH1 left half-width of fortran array A
-! --> MAH2 right half-width of fortran array A
-! --> MBH1 left half-width of fortran array B
-! --> MBH2 right half-width of fortran array B
-!=============================================================================
 use pietc_s, only: u0,u1
 implicit none
 integer(spi), intent(IN   ) ::  m,mah1, mah2, mbh1, mbh2 
@@ -640,9 +737,18 @@ subroutine L1UBB(m,mah1,mah2,mbh1,mbh2,a,b)!                           [L1UBB]
    b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
 enddo
 end subroutine L1UBB
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] m ???
+!! @param[in] mah1 ???
+!! @param[in] mah2 ???
+!! @param[in] mbh1 ???
+!! @param[in] mbh2 ???
+!! @param[in] a ???
+!! @param[in] b ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine DL1UBB(m,mah1,mah2,mbh1,mbh2,a,b) ! Real(dp) version of     [L1UBB]
-!=============================================================================
 use pietc, only: u0,u1
 implicit none
 integer(spi),intent(IN   ) :: m,mah1, mah2, mbh1, mbh2 
@@ -669,29 +775,27 @@ subroutine DL1UBB(m,mah1,mah2,mbh1,mbh2,a,b) ! Real(dp) version of     [L1UBB]
 enddo
 end subroutine DL1UBB
 
-!=============================================================================
+!>  Form the [L]*[D]*[U] decomposition of asymmetric band-matrix [A]
+!!  replace all but row zero of the lower triangular elements of [A]
+!!  by [D**-1]*[L]*[D], the upper by [U], replace matrix [B] by
+!!  [D**-1]*[B].
+!!
+!!  This is a special adaptation of L1UBB used to process quadarature
+!!  weights for QEDBV etc in which the initial quadrature value is
+!!  provided as input instead of being implicitly assumed zero (which
+!!  is the case for QZDBV etc).
+!!
+!! @param[in] m number of rows of B, one less than the rows of A (which has "row 0")
+!! @param[in] mah1 left half-width of fortran array A
+!! @param[in] mah2 right half-width of fortran array A
+!! @param[in] mbh1 left half-width of fortran array B
+!! @param[in] mbh2 right half-width of fortran array B
+!! @param[inout] a input as band matrix, output as lower and upper triangulars with 1s
+!! implicitly assumed to lie on the main diagonal. The product of these
+!! triangular matrices is [D**-1]*[A], where [D] is a diagonal matrix.
+!! @param[inout] b in as band matrix, out as same but premultiplied by diagonal [D**-1]
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1998
 subroutine L1UEB(m,mah1,mah2,mbh1,mbh2,a,b)!                           [L1UEB]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1998
-!		    SUBROUTINE L1UEB
-!  Form the [L]*[D]*[U] decomposition of asymmetric band-matrix  [A] replace
-!  all but row zero of the
-!  lower triangular elements of [A] by [D**-1]*[L]*[D], the upper by [U],
-!  replace matrix [B] by [D**-1]*[B].
-!  This is a special adaptation of L1UBB used to process quadarature weights
-!  for QEDBV etc in which the initial quadrature value is provided as input
-!  instead of being implicitly assumed zero (which is the case for QZDBV etc).
-!
-! <-> A input as band matrix, output as lower and upper triangulars with 1s
-!     implicitly assumed to lie on the main diagonal. The product of these
-!     triangular matrices is [D**-1]*[A], where [D] is a diagonal matrix.
-! <-> B in as band matrix, out as same but premultiplied by diagonal [D**-1]
-! --> M    number of rows of B, one less than the rows of A (which has "row 0")
-! --> MAH1 left half-width of fortran array A
-! --> MAH2 right half-width of fortran array A
-! --> MBH1 left half-width of fortran array B
-! --> MBH2 right half-width of fortran array B
-!=============================================================================
 use pietc_s, only: u0,u1
 implicit none
 integer(spi),intent(IN   ) :: m,mah1, mah2, mbh1, mbh2 
@@ -717,9 +821,18 @@ subroutine L1UEB(m,mah1,mah2,mbh1,mbh2,a,b)!                           [L1UEB]
    b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
 enddo
 end subroutine L1UEB
-!=============================================================================
+
+!> ???
+!!
+!! @param[in] m ???
+!! @param[in] mah1 ???
+!! @param[in] mah2 ???
+!! @param[in] mbh1 ???
+!! @param[in] mbh2 ???
+!! @param[in] a ???
+!! @param[in] b ???
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine DL1UEB(m,mah1,mah2,mbh1,mbh2,a,b) ! Real(dp) version of     [L1UEB]
-!=============================================================================
 use pietc, only: u0,u1
 implicit none
 integer(spi),intent(IN   ):: m,mah1, mah2, mbh1, mbh2 
@@ -746,21 +859,17 @@ subroutine DL1UEB(m,mah1,mah2,mbh1,mbh2,a,b) ! Real(dp) version of     [L1UEB]
 enddo
 end subroutine DL1UEB
 
-!=============================================================================
+!>  Back-substitution step of linear inversion involving banded matrix
+!!  and Vector.
+!!
+!! @param[in] m the number of rows assumed for A and for V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[in] mah2 the right half-bandwidth of fortran array A
+!! @param[inout] a encodes the (L)*(D**-1)*(U) factorization of the linear-system
+!! matrix, as supplied by subroutine LDUB
+!! @param[inout] v input as right-hand-side vector, output as solution vector
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1994
 subroutine UDLBV(m,mah1,mah2,a,v)!                                    [UDLBV]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
-!		    SUBROUTINE UDLBV
-!  BACk-substitution step of linear inversion involving
-!  Banded matrix and Vector.
-!
-! --> A encodes the (L)*(D**-1)*(U) factorization of the linear-system
-!     matrix, as supplied by subroutine LDUB
-! <-> V input as right-hand-side vector, output as solution vector
-! --> M the number of rows assumed for A and for V
-! --> MAH1 the left half-bandwidth of fortran array A
-! --> MAH2 the right half-bandwidth of fortran array A
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ):: m, mah1, mah2
 real(sp),    intent(IN   ):: a(m,-mah1:mah2)
@@ -778,9 +887,18 @@ subroutine UDLBV(m,mah1,mah2,a,v)!                                    [UDLBV]
    do i=max(1,j-mah2),j-1; v(i)=v(i)-a(i,j-i)*vj; enddo
 enddo
 end subroutine UDLBV
-!=============================================================================
+
+!>  Back-substitution step of linear inversion involving banded matrix
+!!  and Vector.
+!!
+!! @param[in] m the number of rows assumed for A and for V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[in] mah2 the right half-bandwidth of fortran array A
+!! @param[inout] a encodes the (L)*(D**-1)*(U) factorization of the linear-system
+!! matrix, as supplied by subroutine LDUB
+!! @param[inout] v input as right-hand-side vector, output as solution vector
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1994
 subroutine dudlbv(m,mah1,mah2,a,v)!                                    [udlbv]
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ) :: m, mah1, mah2
 real(dp),    intent(IN   ) :: a(m,-mah1:mah2)
@@ -799,12 +917,17 @@ subroutine dudlbv(m,mah1,mah2,a,v)!                                    [udlbv]
 enddo
 end subroutine dudlbv
 
-!=============================================================================
+!> Like udlbv, except assuming a is the ltdl decomposition of a
+!! SYMMETRIC banded matrix, with only the non-upper part provided (to
+!! avoid redundancy)
+!!
+!! @param[in] m the number of rows assumed for A and for V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[inout] a encodes the (L)*(D**-1)*(U) factorization of the linear-system
+!! matrix, as supplied by subroutine LDUB
+!! @param[inout] v input as right-hand-side vector, output as solution vector
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine ltdlbv(m,mah1,a,v)!                                        [ltdlbv]
-!=============================================================================
-! Like udlbv, except assuming a is the ltdl decomposition of a SYMMETRIC
-! banded matrix, with only the non-upper part provided (to avoid redundancy)
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ) :: m, mah1
 real(sp),    intent(IN   ) :: a(m,-mah1:0)
@@ -822,12 +945,19 @@ subroutine ltdlbv(m,mah1,a,v)!                                        [ltdlbv]
    do i=max(1,j-mah1),j-1; v(i)=v(i)-a(j,i-j)*vj; enddo
 enddo
 end subroutine ltdlbv
-!=============================================================================
+
+
+!> Like udlbv, except assuming a is the ltdl decomposition of a
+!! SYMMETRIC banded matrix, with only the non-upper part provided (to
+!! avoid redundancy)
+!!
+!! @param[in] m the number of rows assumed for A and for V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[inout] a encodes the (L)*(D**-1)*(U) factorization of the linear-system
+!! matrix, as supplied by subroutine LDUB
+!! @param[inout] v input as right-hand-side vector, output as solution vector
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine dltdlbv(m,mah1,a,v)!                                       [ltdlbv]
-!=============================================================================
-! Like udlbv, except assuming a is the ltdl decomposition of a SYMMETRIC
-! banded matrix, with only the non-upper part provided (to avoid redundancy)
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ) :: m, mah1
 real(dp),    intent(IN   ) :: a(m,-mah1:0)
@@ -846,23 +976,19 @@ subroutine dltdlbv(m,mah1,a,v)!                                       [ltdlbv]
 enddo
 end subroutine dltdlbv
 
-!=============================================================================
+!>  Back-substitution step of parallel linear inversion involving
+!!  Banded matrix and X-Vectors.
+!!
+!! @param[in] mx the number of rows assumed for A and length of
+!!     X-vectors stored in V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[in] mah2 the right half-bandwidth of fortran array A
+!! @param[in] my number of parallel X-vectors inverted
+!! @param[in] a encodes the (L)*(D**-1)*(U) factorization of the linear-system
+!!     matrix, as supplied by subroutine LDUB or, if N=NA, by LDUB
+!! @param[inout] v input as right-hand-side vectors, output as solution vectors
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1994
 subroutine UDLBX(mx,mah1,mah2,my,a,v)!                                [UDLBX]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
-!		    SUBROUTINE UDLBX
-!  BACk-substitution step of parallel linear inversion involving
-!  Banded matrix and X-Vectors.
-!
-! --> A encodes the (L)*(D**-1)*(U) factorization of the linear-system
-!     matrix, as supplied by subroutine LDUB or, if N=NA, by LDUB
-! <-> V input as right-hand-side vectors, output as solution vectors
-! --> MX the number of rows assumed for A and length of
-!     X-vectors stored in V
-! --> MAH1 the left half-bandwidth of fortran array A
-! --> MAH2 the right half-bandwidth of fortran array A
-! --> MY number of parallel X-vectors inverted
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ) :: mx, mah1, mah2, my
 real(sp),    intent(IN   ) :: a(mx,-mah1:mah2)
@@ -879,23 +1005,19 @@ subroutine UDLBX(mx,mah1,mah2,my,a,v)!                                [UDLBX]
 enddo
 end subroutine UDLBX
 
-!=============================================================================
+!>  Back-substitution step of parallel linear inversion involving
+!!  Banded matrix and Y-Vectors.
+!!
+!! @param[in] my the number of rows assumed for A and length of
+!!     Y-vectors stored in V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[in] mah2 the right half-bandwidth of fortran array A
+!! @param[in] mx number of parallel Y-vectors inverted
+!! @param[in] a encodes the (L)*(D**-1)*(U) factorization of the linear-system
+!!     matrix, as supplied by subroutine LDUB or, if N=NA, by LDUB
+!! @param[inout] v input as right-hand-side vectors, output as solution vectors
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1994
 subroutine UDLBY(my,mah1,mah2,mx,a,v)!                                [UDLBY]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
-!		    SUBROUTINE UDLBY
-!  BACk-substitution step of parallel linear inversion involving
-!  Banded matrix and Y-Vectors.
-!
-! --> A encodes the (L)*(D**-1)*(U) factorization of the linear-system
-!     matrix, as supplied by subroutine LDUB or, if N=NA, by LDUB
-! <-> V input as right-hand-side vectors, output as solution vectors
-! --> MY the number of rows assumed for A and length of
-!     Y-vectors stored in V
-! --> MAH1 the left half-bandwidth of fortran array A
-! --> MAH2 the right half-bandwidth of fortran array A
-! --> MX number of parallel Y-vectors inverted
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ) :: my, mah1, mah2, mx
 real(sp),    intent(IN   ) :: a(my,-mah1:mah2)
@@ -912,21 +1034,17 @@ subroutine UDLBY(my,mah1,mah2,mx,a,v)!                                [UDLBY]
 enddo
 end subroutine UDLBY
 
-!=============================================================================
+!>  Back-substitution step of linear inversion involving
+!!  row-Vector and Banded matrix.
+!!
+!! @param[in] m the number of rows assumed for A and columns for V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[in] mah2 the right half-bandwidth of fortran array A
+!! @param[inout] v input as right-hand-side row-vector, output as solution vector
+!! @param[in] a encodes the (L)*(D**-1)*(U) factorization of the linear-system
+!!     matrix, as supplied by subroutine LDUB
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1994
 subroutine UDLVB(m,mah1,mah2,v,a)!                                    [UDLVB]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
-!		    SUBROUTINE UDLVB
-!  BACk-substitution step of linear inversion involving
-!  row-Vector and Banded matrix.
-!
-! <-> V input as right-hand-side row-vector, output as solution vector
-! --> A encodes the (L)*(D**-1)*(U) factorization of the linear-system
-!     matrix, as supplied by subroutine LDUB
-! --> M the number of rows assumed for A and columns for V
-! --> MAH1 the left half-bandwidth of fortran array A
-! --> MAH2 the right half-bandwidth of fortran array A
-!=============================================================================
 implicit none
 integer(spi), intent(IN   ) :: m, mah1, mah2
 real(sp),     intent(IN   ) :: a(m,-mah1:mah2)
@@ -946,23 +1064,19 @@ subroutine UDLVB(m,mah1,mah2,v,a)!                                    [UDLVB]
 enddo
 end subroutine UDLVB
 
-!=============================================================================
+!>  Back-substitution step of parallel linear inversion involving
+!!  Banded matrix and row-X-Vectors.
+!!
+!! @param[in] mx the number of rows assumed for A and length of
+!!     X-vectors stored in V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[in] mah2 the right half-bandwidth of fortran array A
+!! @param[in] my number of parallel X-vectors inverted
+!! @param[inout] v input as right-hand-side vectors, output as solution vectors
+!! @param[in] a encodes the (L)*(D**-1)*(U) factorization of the linear-system
+!!     matrix, as supplied by subroutine LDUB
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1994
 subroutine UDLXB(mx,mah1,mah2,my,v,a)!                                [UDLXB]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
-!		    SUBROUTINE UDLXB
-!  BACk-substitution step of parallel linear inversion involving
-!  Banded matrix and row-X-Vectors.
-!
-! <-> V input as right-hand-side vectors, output as solution vectors
-! --> A encodes the (L)*(D**-1)*(U) factorization of the linear-system
-!     matrix, as supplied by subroutine LDUB
-! --> MX the number of rows assumed for A and length of
-!     X-vectors stored in V
-! --> MAH1 the left half-bandwidth of fortran array A
-! --> MAH2 the right half-bandwidth of fortran array A
-! --> MY number of parallel X-vectors inverted
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ) :: mx, mah1, mah2, my
 real(sp),    intent(IN   ) :: a(mx,-mah1:mah2)
@@ -979,23 +1093,19 @@ subroutine UDLXB(mx,mah1,mah2,my,v,a)!                                [UDLXB]
 enddo
 end subroutine UDLXB
 
-!=============================================================================
+!>  BACk-substitution step of parallel linear inversion involving
+!!  Banded matrix and row-Y-Vectors.
+!!
+!! @param[in] my the number of rows assumed for A and length of
+!!     Y-vectors stored in V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[in] mah2 the right half-bandwidth of fortran array A
+!! @param[in] mx number of parallel Y-vectors inverted
+!! @param[inout] v input as right-hand-side vectors, output as solution vectors
+!! @param[in] a encodes the (L)*(D**-1)*(U) factorization of the linear-system
+!!     matrix, as supplied by subroutine LDUB
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1994
 subroutine UDLYB(my,mah1,mah2,mx,v,a)!                                [UDLYB]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
-!		    SUBROUTINE UDLYB
-!  BACk-substitution step of parallel linear inversion involving
-!  Banded matrix and row-Y-Vectors.
-!
-! <-> V input as right-hand-side vectors, output as solution vectors
-! --> A encodes the (L)*(D**-1)*(U) factorization of the linear-system
-!     matrix, as supplied by subroutine LDUB
-! --> MY the number of rows assumed for A and length of
-!     Y-vectors stored in V
-! --> MAH1 the left half-bandwidth of fortran array A
-! --> MAH2 the right half-bandwidth of fortran array A
-! --> MX number of parallel Y-vectors inverted
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ) :: my, mah1, mah2, mx
 real(sp),    intent(IN   ) :: a(my,-mah1:mah2)
@@ -1012,21 +1122,17 @@ subroutine UDLYB(my,mah1,mah2,mx,v,a)!                                [UDLYB]
 enddo
 end subroutine UDLYB
 
-!=============================================================================
+!>  BACk-substitution step ((U**-1)*(L**-1)) of linear inversion involving
+!!  special Banded matrix and right-Vector.
+!!
+!! @param[in] m the number of rows assumed for A and for V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[in] mah2 the right half-bandwidth of fortran array A
+!! @param[in] a encodes the [L]*[U] factorization of the linear-system
+!!     matrix, as supplied by subroutine L1UBB
+!! @param[inout] v input as right-hand-side vector, output as solution vector
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1996
 subroutine U1LBV(m,mah1,mah2,a,v)!                                    [U1LBV]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1996
-!		    SUBROUTINE U1LBV
-!  BACk-substitution step ((U**-1)*(L**-1)) of linear inversion involving
-!  special Banded matrix and right-Vector.
-!
-! --> A encodes the [L]*[U] factorization of the linear-system
-!     matrix, as supplied by subroutine L1UBB
-! <-> V input as right-hand-side vector, output as solution vector
-! --> M the number of rows assumed for A and for V
-! --> MAH1 the left half-bandwidth of fortran array A
-! --> MAH2 the right half-bandwidth of fortran array A
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ) :: m, mah1, mah2
 real(sp),    intent(IN   ) :: a(m,-mah1:mah2)
@@ -1045,23 +1151,19 @@ subroutine U1LBV(m,mah1,mah2,a,v)!                                    [U1LBV]
 enddo
 end subroutine U1LBV
 
-!=============================================================================
+!>  Special BaCk-substitution step of parallel linear inversion involving
+!!  Banded matrix and X-right-Vectors.
+!!
+!! @param[in] mx the number of rows assumed for A and length of
+!!     X-vectors stored in V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[in] mah2 the right half-bandwidth of fortran array A
+!! @param[in] my number of parallel X-vectors inverted
+!! @param[in] a encodes the [L]*[U] factorization of the linear-system
+!!     matrix, as supplied by subroutine L1UBB
+!! @param[inout] v input as right-hand-side vectors, output as solution vectors
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1996
 subroutine U1LBX(mx,mah1,mah2,my,a,v)!                                [U1LBX]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1996
-!		    SUBROUTINE U1LBX
-!  Special BaCk-substitution step of parallel linear inversion involving
-!  Banded matrix and X-right-Vectors.
-!
-! --> A encodes the [L]*[U] factorization of the linear-system
-!     matrix, as supplied by subroutine L1UBB
-! <-> V input as right-hand-side vectors, output as solution vectors
-! --> MX the number of rows assumed for A and length of
-!     X-vectors stored in V
-! --> MAH1 the left half-bandwidth of fortran array A
-! --> MAH2 the right half-bandwidth of fortran array A
-! --> MY number of parallel X-vectors inverted
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ) :: mx, mah1, mah2, my
 real(sp),    intent(IN   ) :: a(mx,-mah1:mah2)
@@ -1077,23 +1179,19 @@ subroutine U1LBX(mx,mah1,mah2,my,a,v)!                                [U1LBX]
 enddo
 end subroutine U1LBX
 
-!=============================================================================
+!>  Special BaCk-substitution step of parallel linear inversion involving
+!!  Banded matrix and Y-right-Vectors.
+!!
+!! @param[in] my the number of rows assumed for A and length of
+!!     Y-vectors stored in V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[in] mah2 the right half-bandwidth of fortran array A
+!! @param[in] mx number of parallel Y-vectors inverted
+!! @param[in] a encodes the [L]*[U] factorization of the linear-system
+!!     matrix, as supplied by subroutine L1UBB
+!! @param[inout] v input as right-hand-side vectors, output as solution vectors
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1996
 subroutine U1LBY(my,mah1,mah2,mx,a,v)!                                [U1LBY]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1996
-!		    SUBROUTINE U1LBY
-!  Special BaCk-substitution step of parallel linear inversion involving
-!  Banded matrix and Y-right-Vectors.
-!
-! --> A encodes the [L]*[U] factorization of the linear-system
-!     matrix, as supplied by subroutine L1UBB
-! <-> V input as right-hand-side vectors, output as solution vectors
-! --> MY the number of rows assumed for A and length of
-!     Y-vectors stored in V
-! --> MAH1 the left half-bandwidth of fortran array A
-! --> MAH2 the right half-bandwidth of fortran array A
-! --> MX number of parallel Y-vectors inverted
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ) :: my, mah1, mah2, mx
 real(sp),    intent(IN   ) :: a(my,-mah1:mah2)
@@ -1109,21 +1207,17 @@ subroutine U1LBY(my,mah1,mah2,mx,a,v)!                                [U1LBY]
 enddo
 end subroutine U1LBY
 
-!=============================================================================
+!>  Special BaCk-substitution step of linear inversion involving
+!!  left-Vector and Banded matrix.
+!!
+!! @param[in] m the number of rows assumed for A and columns for V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[in] mah2 the right half-bandwidth of fortran array A
+!! @param[inout] v input as right-hand-side row-vector, output as solution vector
+!! @param[in] a encodes the special [L]*[U] factorization of the linear-system
+!!     matrix, as supplied by subroutine L1UBB
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1996
 subroutine U1LVB(m,mah1,mah2,v,a)!                                    [U1LVB]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1996
-!		    SUBROUTINE U1LVB
-!  Special BaCk-substitution step of linear inversion involving
-!  left-Vector and Banded matrix.
-!
-! <-> V input as right-hand-side row-vector, output as solution vector
-! --> A encodes the special [L]*[U] factorization of the linear-system
-!     matrix, as supplied by subroutine L1UBB
-! --> M the number of rows assumed for A and columns for V
-! --> MAH1 the left half-bandwidth of fortran array A
-! --> MAH2 the right half-bandwidth of fortran array A
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ) :: m, mah1, mah2
 real(sp),    intent(IN   ) :: a(m,-mah1:mah2)
@@ -1142,23 +1236,19 @@ subroutine U1LVB(m,mah1,mah2,v,a)!                                    [U1LVB]
 enddo
 end subroutine U1LVB
 
-!=============================================================================
+!>  Special BaCk-substitution step of parallel linear inversion involving
+!!  Banded matrix and X-left-Vectors.
+!!
+!! @param[in] mx the number of rows assumed for A and length of
+!!     X-vectors stored in V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[in] mah2 the right half-bandwidth of fortran array A
+!! @param[in] my number of parallel X-vectors inverted
+!! @param[in] v input as right-hand-side vectors, output as solution vectors
+!! @param[in] a encodes the special [L]*[U] factorization of the linear-system
+!!     matrix, as supplied by subroutine L1UBB
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1996
 subroutine U1LXB(mx,mah1,mah2,my,v,a)!                                [U1LXB]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1996
-!		    SUBROUTINE U1LXB
-!  Special BaCk-substitution step of parallel linear inversion involving
-!  Banded matrix and X-left-Vectors.
-!
-! <-> V input as right-hand-side vectors, output as solution vectors
-! --> A encodes the special [L]*[U] factorization of the linear-system
-!     matrix, as supplied by subroutine L1UBB
-! --> MX the number of rows assumed for A and length of
-!     X-vectors stored in V
-! --> MAH1 the left half-bandwidth of fortran array A
-! --> MAH2 the right half-bandwidth of fortran array A
-! --> MY number of parallel X-vectors inverted
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ) :: mx, mah1, mah2, my
 real(sp),    intent(IN   ) :: a(mx,-mah1:mah2)
@@ -1174,23 +1264,19 @@ subroutine U1LXB(mx,mah1,mah2,my,v,a)!                                [U1LXB]
 enddo
 end subroutine U1LXB
 
-!=============================================================================
+!>  Special BaCk-substitution step of parallel linear inversion involving
+!!  special Banded matrix and Y-left-Vectors.
+!!
+!! @param[in] my the number of rows assumed for A and length of
+!!     Y-vectors stored in V
+!! @param[in] mah1 the left half-bandwidth of fortran array A
+!! @param[in] mah2 the right half-bandwidth of fortran array A
+!! @param[in] mx number of parallel Y-vectors inverted
+!! @param[inout] v input as right-hand-side vectors, output as solution vectors
+!! @param[in] a encodes the [L]*[U] factorization of the linear-system
+!!     matrix, as supplied by subroutine L1UBB
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1996
 subroutine U1LYB(my,mah1,mah2,mx,v,a)!                                [U1LYB]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1996
-!		    SUBROUTINE U1LYB
-!  Special BaCk-substitution step of parallel linear inversion involving
-!  special Banded matrix and Y-left-Vectors.
-!
-! <-> V input as right-hand-side vectors, output as solution vectors
-! --> A encodes the [L]*[U] factorization of the linear-system
-!     matrix, as supplied by subroutine L1UBB
-! --> MY the number of rows assumed for A and length of
-!     Y-vectors stored in V
-! --> MAH1 the left half-bandwidth of fortran array A
-! --> MAH2 the right half-bandwidth of fortran array A
-! --> MX number of parallel Y-vectors inverted
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ) :: my, mah1, mah2, mx
 real(sp),    intent(IN   ) :: a(my,-mah1:mah2)
@@ -1206,19 +1292,15 @@ subroutine U1LYB(my,mah1,mah2,mx,v,a)!                                [U1LYB]
 enddo
 end subroutine U1LYB
 
-!=============================================================================
+!> Solve LINear system with square Banded-matrix and vector V.
+!!
+!! @param[in] m order of matrix A
+!! @param[in] mah1 left half-bandwidth of A
+!! @param[in] mah2 right half-bandwidth of A
+!! @param[inout] a system matrix on input, its [L]*[D**-1]*[U] factorization on exit
+!! @param[inout] v vector of right-hand-sides on input, solution vector on exit
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1994
 subroutine LINBV(m,mah1,mah2,a,v)!                                     [LINBV]
-!=============================================================================
-!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
-!		    SUBROUTINE LINBV
-!   Solve LINear system with square Banded-matrix and vector V
-!
-! <-> A system matrix on input, its [L]*[D**-1]*[U] factorization on exit
-! <-> V vector of right-hand-sides on input, solution vector on exit
-! --> M order of matrix A
-! --> MAH1 left half-bandwidth of A
-! --> MAH2 right half-bandwidth of A
-!=============================================================================
 implicit none
 integer(spi),intent(IN   ) :: m, mah1, mah2
 real(sp),    intent(INOUT) :: a(m,-mah1:mah2), v(m)
@@ -1227,9 +1309,15 @@ subroutine LINBV(m,mah1,mah2,a,v)!                                     [LINBV]
 call UDLBV(m,mah1,mah2,a,v)
 end subroutine LINBV
 
-!=============================================================================
+!> ???
+!!
+!! @param[in] m1
+!! @param[in] m2
+!! @param[in] mah1
+!! @param[in] mah2
+!! @param[in] a
+!! @author R. J. Purser, Tsukasa Fujita (JMA) @date 1999
 subroutine WRTB(m1,m2,mah1,mah2,a)!                                     [WRTB]
-!=============================================================================
 implicit none
 integer(spi),intent(IN) :: m1, m2, mah1, mah2
 real(sp),    intent(IN) :: a(m1,-mah1:mah2)