tools_sp.F90

Go to the documentation of this file.

# 1 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/external/tools_sp.fypp"
# 1 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/external/../instrumentation.fypp" 1
# 1 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/external/../subr_list.fypp" 1
!----------------------------------------------------------------------
! Header: subr_list
!> Subroutines/functions list
! Author: Benjamin Menetrier
! Licensing: this code is distributed under the CeCILL-C license
! Copyright 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
!----------------------------------------------------------------------
 
# 926 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/external/../subr_list.fypp"
# 2 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/external/../instrumentation.fypp" 2
!----------------------------------------------------------------------
! Header: instrumentation
!> Instrumentation functions
! Author: Benjamin Menetrier
! Licensing: this code is distributed under the CeCILL-C license
! Copyright 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
!----------------------------------------------------------------------
 
# 112 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/external/../instrumentation.fypp"
# 2 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/external/tools_sp.fypp" 2
!----------------------------------------------------------------------
! (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_const, only: zero,quarter,half,one,two,ten,pi
use tools_kinds, only: kind_real
 
 
implicit none
 
interface splat
   module procedure sp_splat
end interface
interface lubksb
   module procedure sp_lubksb
end interface
interface ludcmp
   module procedure sp_ludcmp
end interface
 
private
public :: splat
 
contains
 
!-----------------------------------------------------------------
! Subroutine: sp_splat
!> Compute latitude functions
!-----------------------------------------------------------------
subroutine sp_splat(IDRT,JMAX,SLAT,WLAT)
 
implicit none
 
! Passed variables
integer,intent(in) :: IDRT                !< Grid identifier
integer,intent(in) :: JMAX                !< Number of latitudes
real(kind_real),intent(out) :: SLAT(JMAX) !< Sines of latitude
real(kind_real),intent(out) :: WLAT(JMAX) !< Gaussian weights
 
! Local variables
real(kind_real):: PK(JMAX/2),PKM1(JMAX/2),PKM2(JMAX/2)
real(kind_real):: SLATD(JMAX/2),SP,SPMAX,EPS
integer,parameter:: JZ=50
real(kind_real):: BZ(JZ)
real(kind_real):: DLT,D1,d,R
real(kind_real):: AWORK((JMAX+1)/2,((JMAX+1)/2))
real(kind_real):: BWORK(((JMAX+1)/2))
integer:: JH,JHE,JHO,J0
integer:: IPVT((JMAX+1)/2)
real(kind_real),parameter:: C=(one-(two/pi)**2)*quarter
integer:: J,JS,N
 
! Set name
 
 
! Probe in
 
 
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
 
! Probe out
 
 
end subroutine sp_splat
 
!-----------------------------------------------------------------
! Subroutine: sp_lubksb
!> Solves a system of linear equations, follows call to LUDCMP
!-----------------------------------------------------------------
subroutine sp_lubksb(A,N,NP,INDX,B)
 
implicit none
 
! Passed varaibles
integer,intent(in):: NP               !< ?
integer,intent(in):: N                !< ?
real(kind_real),intent(in):: A(NP,NP) !< ?
real(kind_real),intent(inout):: B(N)  !< ?
integer,intent(in):: INDX(N)          !< ?
 
! Local variables
real(kind_real):: SUM
integer:: I,II,J,LL
 
! Set name
 
 
! Probe in
 
 
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
 
! Probe out
 
 
end subroutine sp_lubksb
 
!-----------------------------------------------------------------
! Subroutine: sp_ludcmp
!> Replaces an NxN matrix A with the LU decomposition
!-----------------------------------------------------------------
subroutine sp_ludcmp(A,N,NP,INDX,D)
 
implicit none
 
! Passed variables
integer,intent(in):: N                   !< ?
integer,intent(in):: NP                  !< ?
real(kind_real),intent(inout):: A(NP,NP) !< ?
integer,intent(out):: INDX(N)            !< ?
real(kind_real),intent(out),optional:: D !< ?
 
! Local variables
real(kind_real),parameter:: TINY=1.0e-20_kind_real
real(kind_real):: AAMAX,SUM,DUM
real(kind_real):: VV(N)
integer:: IMAX
integer:: I,J,K
 
! Set name
 
 
! Probe in
 
 
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
 
! Probe out
 
 
end subroutine sp_ludcmp
 
end module tools_sp