tools_sp.F90

Go to the documentation of this file.

! (C) Copyright 2018-2019 UCAR
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!----------------------------------------------------------------------
! Module: tools_sp
!> Spectral transforms
!----------------------------------------------------------------------
module tools_sp
 
use tools_kinds, only: kind_real
use tools_const, only: pi
 
implicit none
 
private
real(kind=kind_real),parameter:: zero    = 0.0_kind_real
real(kind=kind_real),parameter:: one     = 1.0_kind_real
real(kind=kind_real),parameter:: two     = 2.0_kind_real
real(kind=kind_real),parameter:: ten     = 10.0_kind_real
real(kind=kind_real),parameter:: half    = 0.5_kind_real
real(kind=kind_real),parameter:: quarter = 0.25_kind_real
 
public :: splat
 
contains
 
!-----------------------------------------------------------------
! Subroutine: splat
!> Compute latitude functions
!-----------------------------------------------------------------
subroutine splat(IDRT,JMAX,SLAT,WLAT)
 
implicit none
 
! Passed variables
integer,intent(in ) :: idrt                    !< Grid identifier
integer,intent(in ) :: jmax                    !< Number of latitudes
real(kind=kind_real),intent(out) :: slat(jmax) !< Sines of latitude 
real(kind=kind_real),intent(out) :: wlat(jmax) !< Gaussian weights-
 
! Local variables
real(kind=kind_real):: pk(jmax/2),pkm1(jmax/2),pkm2(jmax/2)
real(kind=kind_real):: slatd(jmax/2),sp,spmax,eps
integer,parameter:: jz=50
real(kind=kind_real):: bz(jz)
real(kind=kind_real):: dlt,d1,d,r
real(kind=kind_real):: awork((jmax+1)/2,((jmax+1)/2))
real(kind=kind_real):: bwork(((jmax+1)/2))
integer:: jh,jhe,jho,j0
integer:: ipvt((jmax+1)/2)
real(kind=kind_real),parameter:: c=(one-(two/pi)**2)*quarter
integer:: j,js,n
 
data bz / 2.4048255577_kind_real,  5.5200781103_kind_real, &
 &        8.6537279129_kind_real, 11.7915344391_kind_real, &
 &       14.9309177086_kind_real, 18.0710639679_kind_real, &
 &       21.2116366299_kind_real, 24.3524715308_kind_real, &
 &       27.4934791320_kind_real, 30.6346064684_kind_real, &
 &       33.7758202136_kind_real, 36.9170983537_kind_real, &
 &       40.0584257646_kind_real, 43.1997917132_kind_real, &
 &       46.3411883717_kind_real, 49.4826098974_kind_real, &
 &       52.6240518411_kind_real, 55.7655107550_kind_real, &
 &       58.9069839261_kind_real, 62.0484691902_kind_real, &
 &       65.1899648002_kind_real, 68.3314693299_kind_real, &
 &       71.4729816036_kind_real, 74.6145006437_kind_real, &
 &       77.7560256304_kind_real, 80.8975558711_kind_real, &
 &       84.0390907769_kind_real, 87.1806298436_kind_real, &
 &       90.3221726372_kind_real, 93.4637187819_kind_real, &
 &       96.6052679510_kind_real, 99.7468198587_kind_real, &
 &       102.888374254_kind_real, 106.029930916_kind_real, &
 &       109.171489649_kind_real, 112.313050280_kind_real, &
 &       115.454612653_kind_real, 118.596176630_kind_real, &
 &       121.737742088_kind_real, 124.879308913_kind_real, &
 &       128.020877005_kind_real, 131.162446275_kind_real, &
 &       134.304016638_kind_real, 137.445588020_kind_real, &
 &       140.587160352_kind_real, 143.728733573_kind_real, &
 &       146.870307625_kind_real, 150.011882457_kind_real, &
 &       153.153458019_kind_real, 156.295034268_kind_real /
 
! Initialization
eps=ten*epsilon(sp)
d1=one
j0=0
 
if(idrt==4) then
   ! Gaussian latitudes
   jh=jmax/2
   jhe=(jmax+1)/2
   r=one/sqrt((real(jmax,kind_real)+half)**2+c)
   do j=1,min(jh,jz)
      slatd(j)=cos(bz(j)*r)
   end do
   do j=jz+1,jh
      slatd(j)=cos((bz(jz)+(j-jz)*pi)*r)
   end do
   spmax=one
   do while(spmax>eps)
      spmax=zero
      do j=1,jh
         pkm1(j)=one
         pk(j)=slatd(j)
      end do
      do n=2,jmax
         do j=1,jh
            pkm2(j)=pkm1(j)
            pkm1(j)=pk(j)
            pk(j)=real((2*n-1)*slatd(j)*pkm1(j)-(n-1)*pkm2(j),kind_real)/real(n,kind_real)
         end do
      end do
      do j=1,jh
         sp=pk(j)*(one-slatd(j)**2)/(jmax*(pkm1(j)-slatd(j)*pk(j)))
         slatd(j)=slatd(j)-sp
         spmax=max(spmax,abs(sp))
      end do
   end do
   do j=1,jh
      slat(j)=slatd(j)
      wlat(j)=(two*(one-slatd(j)**2))/(jmax*pkm1(j))**2
      slat(jmax+1-j)=-slat(j)
      wlat(jmax+1-j)=wlat(j)
   end do
   if(jhe>jh) then
      slat(jhe)=zero
      wlat(jhe)=two/jmax**2
      do n=2,jmax,2
         wlat(jhe)=wlat(jhe)*n**2/(n-1)**2
      end do
   end if
else if(idrt==0) then
   ! Equally-spaced latitudes including poles
   jh=jmax/2
   jhe=(jmax+1)/2
   jho=jhe-1
   dlt=pi/(jmax-1)
   slat(1)=one
   do j=2,jh
      slat(j)=cos(real((j-1),kind_real)*dlt)
   end do
   do js=1,jho
      do j=1,jho
         awork(js,j)=cos(real(2*(js-1)*j,kind_real)*dlt)
      end do
   end do
   do js=1,jho
      bwork(js)=-d1/real((4*(js-1)**2-1),kind_real)
   end do
   call ludcmp(awork,jho,jhe,ipvt)
   call lubksb(awork,jho,jhe,ipvt,bwork)
   wlat(1)=zero
   do j=1,jho
      wlat(j+1)=bwork(j)
   end do
   do j=1,jh
      slat(jmax+1-j)=-slat(j)
      wlat(jmax+1-j)=wlat(j)
   end do
   if(jhe>jh) then
      slat(jhe)=zero
      wlat(jhe)=two*wlat(jhe)
   end if
else if(idrt==256) then
   ! Equally-spaced latitudes excluding poles
   jh=jmax/2
   jhe=(jmax+1)/2
   jho=jhe
   dlt=pi/jmax
   slat(1)=one
   do j=1,jh
      slat(j)=cos((real(j,kind_real)-half)*dlt)
   end do
   do js=1,jho
      do j=1,jho
         awork(js,j)=cos(real(2*(js-1),kind_real)*(real(j,kind_real)-half)*dlt)
      end do
   end do
   do js=1,jho
      bwork(js)=-d1/real(4*(js-1)**2-1,kind_real)
   end do
   call ludcmp(awork,jho,jhe,ipvt,d)
   call lubksb(awork,jho,jhe,ipvt,bwork)
   wlat(1)=zero
   do j=1,jho
      wlat(j)=bwork(j)
   end do
   do j=1,jh
      slat(jmax+1-j)=-slat(j)
      wlat(jmax+1-j)=wlat(j)
   end do
   if(jhe>jh) then
      slat(jhe)=zero
      wlat(jhe)=two*wlat(jhe)
   end if
end if
 
end subroutine splat
 
!-----------------------------------------------------------------
! Subroutine: lubksb
!> Solves a system of linear equations, follows call to LUDCMP
!-----------------------------------------------------------------
subroutine lubksb(A,N,NP,INDX,B)
 
implicit none
 
! Passed varaibles
integer,intent(in):: NP                    !< ?
integer,intent(in):: N                     !< ?
real(kind=kind_real),intent(in):: a(np,np) !< ?
real(kind=kind_real),intent(inout):: b(n)  !< ?
integer,intent(in):: INDX(N)               !< ?
 
! Local variables
real(kind=kind_real):: sum
integer:: I,II,J,LL
 
ii=0
do i=1,n
   ll=indx(i)
   sum=b(ll)
   b(ll)=b(i)
   if(ii/=0) then
      do j=ii,i-1
         sum=sum-a(i,j)*b(j)
      end do
   else if(abs(sum)>zero) then
      ii=i
   end if
   b(i)=sum
end do
do i=n,1,-1
   sum=b(i)
   if(i<n) then
      do j=i+1,n
         sum=sum-a(i,j)*b(j)
      end do
   end if
   b(i)=sum/a(i,i)
end do
 
end subroutine lubksb
 
!-----------------------------------------------------------------
! Subroutine: ludcmp
!> Replaces an NxN matrix A with the LU decomposition
!-----------------------------------------------------------------
subroutine ludcmp(A,N,NP,INDX,D)
 
implicit none
 
! Passed variables
integer,intent(in):: N                        !< ?
integer,intent(in):: NP                       !< ?
real(kind=kind_real),intent(inout):: a(np,np) !< ?
integer,intent(out):: INDX(N)                 !< ? 
real(kind=kind_real),intent(out),optional:: d !< ?
 
! Local variables
real(kind=kind_real),parameter:: tiny=1.0e-20
real(kind=kind_real):: aamax,sum,dum
real(kind=kind_real):: vv(n)
integer:: IMAX
integer:: I,J,K
 
d=one
do i=1,n
   aamax=zero
   do j=1,n
      if(abs(a(i,j))>aamax) aamax=abs(a(i,j))
   end do
   if(.NOT.(abs(aamax)>zero)) print *, 'SINGULAR MATRIX.'
   vv(i)=one/aamax
end do
do j=1,n
   if(j>1) then
      do i=1,j-1
         sum=a(i,j)
         if(i>1) then
            do k=1,i-1
               sum=sum-a(i,k)*a(k,j)
            end do
            a(i,j)=sum
         end if
      end do
   end if
   aamax=zero
   do i=j,n
      sum=a(i,j)
      if(j>1) then
         do k=1,j-1
           sum=sum-a(i,k)*a(k,j)
         end do
         a(i,j)=sum
      end if
      dum=vv(i)*abs(sum)
      if(dum>=aamax) then
         imax=i
         aamax=dum
      end if
   end do
   if(j/=imax) then
      do k=1,n
         dum=a(imax,k)
         a(imax,k)=a(j,k)
         a(j,k)=dum
      end do
      d=-d
      vv(imax)=vv(j)
   end if
   indx(j)=imax
   if(j/=n) then
      if(.NOT.(abs(a(j,j))>zero)) a(j,j)=tiny
      dum=one/a(j,j)
      do i=j+1,n
         a(i,j)=a(i,j)*dum
      end do
   end if
end do
if(.NOT.(abs(a(n,n))>zero)) a(n,n)=tiny
 
end subroutine ludcmp
 
end module tools_sp