#include <misc.h>
#include <params.h>
module spetru 4
!-----------------------------------------------------------------------
!
! Purpose: Spectrally truncate initial data fields.
!
! Method: Truncate one or a few fields at a time, to minimize the
! memory requirements
!
! Original version: J. Rosinski
! Standardized: J. Rosinski, June 1992
! Reviewed: B. Boville, J. Hack, August 1992
! Modified to implement processing of subsets of fields: P. Worley, May 2003
!
!-----------------------------------------------------------------------
contains
subroutine spetru_phis(phis ,phis_hires) 2,18
!-----------------------------------------------------------------------
!
! Purpose:
!
! Method:
! Spectrally truncate PHIS input field.
!
! Author:
! Original version: J. Rosinski
! Standardized: J. Rosinski, June 1992
! Reviewed: B. Boville, J. Hack, August 1992
! Modified: P. Worley, May 2003
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid, only: plon, plond, plev, plat
use pspect
use comspe
use rgrid, only: nlon, nmmax
use commap, only: w, xm
use dynconst, only: rearth
implicit none
#include <comctl.h>
#include <comfft.h>
!
! Input/Output arguments
!
real(r8), intent(inout) :: phis(plond,plat) ! Fourier -> spec. coeffs. for sfc geo.
logical, intent(in) :: phis_hires ! true => PHIS came from hi res topo file
!
!---------------------------Local workspace-----------------------------
!
real(r8) phi(2,psp/2) ! used in spectral truncation of phis
real(r8) tmp1 ! vector temporary
real(r8) tmp2 ! vector temporary
real(r8) phialpr,phialpi ! phi*alp (real and imaginary)
real(r8) zwalp ! zw*alp
real(r8) zw ! w**2
real(r8) filtlim ! filter function
real(r8) ft ! filter multiplier for spectral coefficients
#if ( ! defined USEFFTLIB )
real(r8) work((plon+1)*plev) ! Workspace for fft
#else
real(r8) work((plon+1)*pcray) ! Workspace for fft
#endif
integer i ! longitude index
integer irow ! latitude pair index
integer latm,latp ! symmetric latitude indices
integer lat
integer m ! longitudinal wavenumber index
integer n ! latitudinal wavenumber index
integer nspec
integer mr ! spectral indices
!
!-----------------------------------------------------------------------
!
! Zero spectral array
!
phi(:,:) = 0.
!
! Transform grid -> fourier
!
do lat=1,plat
irow = lat
if (lat.gt.plat/2) irow = plat - lat + 1
call fft991(phis(1,lat),work,trig(1,irow),ifax(1,irow),1,plond, &
nlon(lat),1,-1)
end do ! lat=1,plat
!
! Loop over latitude pairs
!
do irow=1,plat/2
latp = irow
latm = plat - irow + 1
zw = w(irow)*2.
do i=1,2*nmmax(irow)
!
! Compute symmetric and antisymmetric components
!
tmp1 = 0.5*(phis(i,latm) - phis(i,latp))
tmp2 = 0.5*(phis(i,latm) + phis(i,latp))
phis(i,latm) = tmp1
phis(i,latp) = tmp2
end do
!
! Compute phi*mn
!
do m=1,nmmax(irow)
mr = nstart(m)
do n=1,nlen(m),2
zwalp = zw*alp(mr+n,irow)
phi(1,mr+n) = phi(1,mr+n) + zwalp*phis(2*m-1,latp)
phi(2,mr+n) = phi(2,mr+n) + zwalp*phis(2*m ,latp)
end do
do n=2,nlen(m),2
zwalp = zw*alp(mr+n,irow)
phi(1,mr+n) = phi(1,mr+n) + zwalp*phis(2*m-1,latm)
phi(2,mr+n) = phi(2,mr+n) + zwalp*phis(2*m ,latm)
end do
end do
enddo ! irow=1,plat/2
!
if (phis_hires) then
!
! Apply spectral filter to phis
! filter is a function of n
! if n < filter limit then
! spectral_coeff = spectral_coeff * (1. - (float(n)/filtlim)**2)
! else
! spectral_coeff = 0.
! endif
! where filter limit = 1.4*PTRN
!
filtlim = float(int(1.4*float(ptrn)))
do m=1,pmmax
mr = nstart(m)
do n=1,nlen(m)
nspec=m-1+n
ft = 1. - (float(nspec)/filtlim)**2
if (float(nspec) .ge. filtlim) ft = 0.
phi(1,mr+n) = phi(1,mr+n)*ft
phi(2,mr+n) = phi(2,mr+n)*ft
end do
end do
call hordif1(rearth,phi)
end if
!
! Compute grid point values of phi*.
!
do irow=1,plat/2
latp = irow
latm = plat - irow + 1
!
! Zero fourier fields
!
phis(:,latm) = 0.
phis(:,latp) = 0.
!
! Compute(phi*)m
!
do m=1,nmmax(irow)
mr = nstart(m)
do n=1,nlen(m),2
phialpr = phi(1,mr+n)*alp(mr+n,irow)
phialpi = phi(2,mr+n)*alp(mr+n,irow)
!
phis(2*m-1,latm) = phis(2*m-1,latm) + phialpr
phis(2*m ,latm) = phis(2*m ,latm) + phialpi
end do
end do
do m=1,nmmax(irow)
mr = nstart(m)
do n=2,nlen(m),2
phialpr = phi(1,mr+n)*alp(mr+n,irow)
phialpi = phi(2,mr+n)*alp(mr+n,irow)
!
phis(2*m-1,latp) = phis(2*m-1,latp) + phialpr
phis(2*m ,latp) = phis(2*m ,latp) + phialpi
!
end do
end do
!
! Recompute real fields from symmetric and antisymmetric parts
!
do i=1,nlon(latm)+2
tmp1 = phis(i,latm) + phis(i,latp)
tmp2 = phis(i,latm) - phis(i,latp)
phis(i,latm) = tmp1
phis(i,latp) = tmp2
end do
enddo ! irow=1,plat/2
!
do lat=1,plat
!
! Transform Fourier -> grid, obtaining spectrally truncated
! grid point values.
!
irow = lat
if (lat.gt.plat/2) irow = plat - lat + 1
call fft991(phis(1,lat),work,trig(1,irow),ifax(1,irow),1,plond, &
nlon(lat),1,+1)
enddo
return
end subroutine spetru_phis
!************************************************************************
subroutine spetru_ps(ps ,dpsl ,dpsm) 2,7
!-----------------------------------------------------------------------
!
! Purpose:
!
! Method:
! Spectrally truncate PS input field.
!
! Author:
! Original version: J. Rosinski
! Standardized: J. Rosinski, June 1992
! Reviewed: B. Boville, J. Hack, August 1992
! Modified: P. Worley, May 2003
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid, only: plon, plond, plev, plat
use pspect
use comspe
use rgrid, only: nlon, nmmax
use commap, only: w, xm
use dynconst, only: ra
implicit none
#include <comctl.h>
#include <comfft.h>
!
! Input/Output arguments
!
real(r8), intent(inout) :: ps(plond,plat) ! Fourier -> spec. coeffs. for ln(ps)
!
! Output arguments
!
real(r8), intent(out) :: dpsl(plond,plat) ! Spectrally trunc d(ln(ps))/d(longitude)
real(r8), intent(out) :: dpsm(plond,plat) ! Spectrally trunc d(ln(ps))/d(latitude)
!
!---------------------------Local workspace-----------------------------
!
real(r8) alps_tmp(psp) ! used in spectral truncation of phis
real(r8) tmp1 ! vector temporary
real(r8) tmp2 ! vector temporary
real(r8) zwalp ! zw*alp
real(r8) psdalpr,psdalpi ! alps (real and imaginary)*dalp
real(r8) psalpr,psalpi ! alps (real and imaginary)*alp
real(r8) zw ! w**2
#if ( ! defined USEFFTLIB )
real(r8) work((plon+1)*plev) ! Workspace for fft
#else
real(r8) work((plon+1)*pcray) ! Workspace for fft
#endif
integer ir,ii ! indices complex coeffs. of spec. arrs.
integer i,k ! longitude, level indices
integer irow ! latitude pair index
integer latm,latp ! symmetric latitude indices
integer lat
integer m ! longitudinal wavenumber index
integer n ! latitudinal wavenumber index
integer nspec
integer mr,mc ! spectral indices
!
!-----------------------------------------------------------------------
!
! Zero spectral array
!
alps_tmp(:) = 0.
!
! Compute the 2D quantities which are transformed to spectral space:
! ps= ln(ps).
!
do lat=1,plat
irow = lat
if (lat.gt.plat/2) irow = plat - lat + 1
do i=1,nlon(lat)
ps(i,lat) = log(ps(i,lat))
end do
!
! Transform grid -> fourier
!
call fft991(ps(1,lat),work,trig(1,irow),ifax(1,irow),1,plond, &
nlon(lat),1,-1)
end do ! lat=1,plat
!
! Loop over latitude pairs
!
do irow=1,plat/2
latp = irow
latm = plat - irow + 1
zw = w(irow)*2.
do i=1,2*nmmax(irow)
!
! Compute symmetric and antisymmetric components
!
tmp1 = 0.5*(ps(i,latm) - ps(i,latp))
tmp2 = 0.5*(ps(i,latm) + ps(i,latp))
ps(i,latm) = tmp1
ps(i,latp) = tmp2
end do
!
! Compute ln(p*)mn
!
do m=1,nmmax(irow)
mr = nstart(m)
mc = 2*mr
do n=1,nlen(m),2
zwalp = zw*alp(mr+n,irow)
ir = mc + 2*n - 1
ii = ir + 1
alps_tmp(ir) = alps_tmp(ir) + zwalp*ps(2*m-1,latp)
alps_tmp(ii) = alps_tmp(ii) + zwalp*ps(2*m ,latp)
end do
do n=2,nlen(m),2
zwalp = zw*alp(mr+n,irow)
ir = mc + 2*n - 1
ii = ir + 1
alps_tmp(ir) = alps_tmp(ir) + zwalp*ps(2*m-1,latm)
alps_tmp(ii) = alps_tmp(ii) + zwalp*ps(2*m ,latm)
end do
end do
enddo ! irow=1,plat/2
!
! Compute grid point values of:ln(p*) and grad(ln(p*)).
!
do irow=1,plat/2
latp = irow
latm = plat - irow + 1
!
! Zero fourier fields
!
ps(:,latm) = 0.
ps(:,latp) = 0.
dpsl(:,latm) = 0.
dpsl(:,latp) = 0.
dpsm(:,latm) = 0.
dpsm(:,latp) = 0.
!
! Compute(ln(p*),grad(ln(p*)))m
!
do m=1,nmmax(irow)
mr = nstart(m)
mc = 2*mr
do n=1,nlen(m),2
ir = mc + 2*n - 1
ii = ir + 1
psalpr = alps_tmp(ir)*alp(mr+n,irow)
psalpi = alps_tmp(ii)*alp(mr+n,irow)
!
ps(2*m-1,latm) = ps(2*m-1,latm) + psalpr
ps(2*m ,latm) = ps(2*m ,latm) + psalpi
dpsl(2*m-1,latm) = dpsl(2*m-1,latm) - psalpi*ra
dpsl(2*m ,latm) = dpsl(2*m ,latm) + psalpr*ra
!
psdalpr = alps_tmp(ir)*dalp(mr+n,irow)
psdalpi = alps_tmp(ii)*dalp(mr+n,irow)
!
dpsm(2*m-1,latp) = dpsm(2*m-1,latp) + psdalpr*ra
dpsm(2*m ,latp) = dpsm(2*m ,latp) + psdalpi*ra
end do
end do
do m=1,nmmax(irow)
mr = nstart(m)
mc = 2*mr
do n=2,nlen(m),2
ir = mc + 2*n - 1
ii = ir + 1
psalpr = alps_tmp(ir)*alp(mr+n,irow)
psalpi = alps_tmp(ii)*alp(mr+n,irow)
!
ps(2*m-1,latp) = ps(2*m-1,latp) + psalpr
ps(2*m ,latp) = ps(2*m ,latp) + psalpi
dpsl(2*m-1,latp) = dpsl(2*m-1,latp) - psalpi*ra
dpsl(2*m ,latp) = dpsl(2*m ,latp) + psalpr*ra
!
psdalpr = alps_tmp(ir)*dalp(mr+n,irow)
psdalpi = alps_tmp(ii)*dalp(mr+n,irow)
!
dpsm(2*m-1,latm) = dpsm(2*m-1,latm) + psdalpr*ra
dpsm(2*m ,latm) = dpsm(2*m ,latm) + psdalpi*ra
end do
end do
do m=1,nmmax(irow)
dpsl(2*m-1,latm) = xm(m)*dpsl(2*m-1,latm)
dpsl(2*m ,latm) = xm(m)*dpsl(2*m ,latm)
dpsl(2*m-1,latp) = xm(m)*dpsl(2*m-1,latp)
dpsl(2*m ,latp) = xm(m)*dpsl(2*m ,latp)
end do
!
! Recompute real fields from symmetric and antisymmetric parts
!
do i=1,nlon(latm)+2
!
tmp1 = ps(i,latm) + ps(i,latp)
tmp2 = ps(i,latm) - ps(i,latp)
ps(i,latm) = tmp1
ps(i,latp) = tmp2
!
tmp1 = dpsl(i,latm) + dpsl(i,latp)
tmp2 = dpsl(i,latm) - dpsl(i,latp)
dpsl(i,latm) = tmp1
dpsl(i,latp) = tmp2
!
tmp1 = dpsm(i,latm) + dpsm(i,latp)
tmp2 = dpsm(i,latm) - dpsm(i,latp)
dpsm(i,latm) = tmp1
dpsm(i,latp) = tmp2
end do
!
enddo ! irow=1,plat/2
!
do lat=1,plat
!
! Transform Fourier -> grid, obtaining spectrally truncated
! grid point values.
!
irow = lat
if (lat.gt.plat/2) irow = plat - lat + 1
call fft991(ps(1,lat),work,trig(1,irow),ifax(1,irow),1,plond, &
nlon(lat),1,+1)
call fft991(dpsl(1,lat),work,trig(1,irow),ifax(1,irow),1,plond, &
nlon(lat),1,+1)
call fft991(dpsm(1,lat),work,trig(1,irow),ifax(1,irow),1,plond, &
nlon(lat),1,+1)
!
! Convert from ln(ps) to ps
!
do i=1,nlon(lat)
ps(i,lat) = exp(ps(i,lat))
end do
!
enddo
return
end subroutine spetru_ps
!************************************************************************
subroutine spetru_t(t3) 2,6
!-----------------------------------------------------------------------
!
! Purpose:
!
! Method:
! Spectrally truncate T input field.
!
! Author:
! Original version: J. Rosinski
! Standardized: J. Rosinski, June 1992
! Reviewed: B. Boville, J. Hack, August 1992
! Modified: P. Worley, May 2003
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid, only: plon, plond, plev, plat
use pspect
use comspe
use rgrid, only: nlon, nmmax
use commap, only: w, xm
implicit none
#include <comctl.h>
#include <comfft.h>
!
! Input/Output arguments
!
real(r8), intent(inout) :: t3(plond,plev,plat) ! Fourier -> spec. coeffs. for temperature
!
!---------------------------Local workspace-----------------------------
!
real(r8) t_tmp(psp) ! used in spectral truncation of t
real(r8) tmp1 ! vector temporary
real(r8) tmp2 ! vector temporary
real(r8) tmpr ! vector temporary (real)
real(r8) tmpi ! vector temporary (imaginary)
real(r8) zwalp ! zw*alp
real(r8) zw ! w**2
#if ( ! defined USEFFTLIB )
real(r8) work((plon+1)*plev) ! Workspace for fft
#else
real(r8) work((plon+1)*pcray) ! Workspace for fft
#endif
integer ir,ii ! indices complex coeffs. of spec. arrs.
integer i,k ! longitude, level indices
integer irow ! latitude pair index
integer latm,latp ! symmetric latitude indices
integer lat
integer m ! longitudinal wavenumber index
integer n ! latitudinal wavenumber index
integer nspec
integer mr,mc ! spectral indices
!
!-----------------------------------------------------------------------
!
! Transform grid -> fourier
!
do lat=1,plat
irow = lat
if (lat.gt.plat/2) irow = plat - lat + 1
call fft991(t3(1,1,lat),work,trig(1,irow),ifax(1,irow),1,plond, &
nlon(lat),plev,-1)
end do ! lat=1,plat
!
! Loop over vertical levels
!
do k=1,plev
!
! Zero spectral array
!
t_tmp(:) = 0.
!
! Loop over latitude pairs
!
do irow=1,plat/2
latp = irow
latm = plat - irow + 1
zw = w(irow)*2.
!
! Multi-level field: T
!
do i=1,2*nmmax(irow)
tmp1 = 0.5*(t3(i,k,latm) - t3(i,k,latp))
tmp2 = 0.5*(t3(i,k,latm) + t3(i,k,latp))
t3(i,k,latm) = tmp1
t3(i,k,latp) = tmp2
end do
!
! Compute tmn
!
do m=1,nmmax(irow)
mr = nstart(m)
mc = 2*mr
do n=1,nlen(m),2
zwalp = zw*alp (mr+n,irow)
ir = mc + 2*n - 1
ii = ir + 1
t_tmp(ir) = t_tmp(ir) + zwalp*t3(2*m-1,k,latp)
t_tmp(ii) = t_tmp(ii) + zwalp*t3(2*m ,k,latp)
end do
end do
do m=1,nmmax(irow)
mr = nstart(m)
mc = 2*mr
do n=2,nlen(m),2
zwalp = zw*alp (mr+n,irow)
ir = mc + 2*n - 1
ii = ir + 1
t_tmp(ir) = t_tmp(ir) + zwalp*t3(2*m-1,k,latm)
t_tmp(ii) = t_tmp(ii) + zwalp*t3(2*m ,k,latm)
end do
end do
enddo ! irow=1,plat/2
!
! Compute grid point values of:t.
!
do irow=1,plat/2
latp = irow
latm = plat - irow + 1
!
! Zero fourier fields
!
t3(:,k,latm) = 0.
t3(:,k,latp) = 0.
do m=1,nmmax(irow)
mr = nstart(m)
mc = 2*mr
do n=1,nlen(m),2
ir = mc + 2*n - 1
ii = ir + 1
!
tmpr = t_tmp(ir)*alp(mr+n,irow)
tmpi = t_tmp(ii)*alp(mr+n,irow)
t3(2*m-1,k,latm) = t3(2*m-1,k,latm) + tmpr
t3(2*m ,k,latm) = t3(2*m ,k,latm) + tmpi
end do
end do
do m=1,nmmax(irow)
mr = nstart(m)
mc = 2*mr
do n=2,nlen(m),2
ir = mc + 2*n - 1
ii = ir + 1
!
tmpr = t_tmp(ir)*alp(mr+n,irow)
tmpi = t_tmp(ii)*alp(mr+n,irow)
t3(2*m-1,k,latp) = t3(2*m-1,k,latp) + tmpr
t3(2*m ,k,latp) = t3(2*m ,k,latp) + tmpi
end do
end do
!
! Recompute real fields from symmetric and antisymmetric parts
!
do i=1,nlon(latm)+2
tmp1 = t3(i,k,latm) + t3(i,k,latp)
tmp2 = t3(i,k,latm) - t3(i,k,latp)
t3(i,k,latm) = tmp1
t3(i,k,latp) = tmp2
end do
enddo ! irow=1,plat/2
enddo ! k=1,plev
!
do lat=1,plat
!
! Transform Fourier -> grid, obtaining spectrally truncated
! grid point values.
irow = lat
if (lat.gt.plat/2) irow = plat - lat + 1
call fft991(t3(1,1,lat),work,trig(1,irow),ifax(1,irow),1,plond, &
nlon(lat),plev,+1)
enddo
return
end subroutine spetru_t
!***********************************************************************
subroutine spetru_uv(u3 ,v3 ,vort ,div) 2,8
!-----------------------------------------------------------------------
!
! Purpose:
!
! Method:
! Spectrally truncate U, V input fields.
!
! Author:
! Original version: J. Rosinski
! Standardized: J. Rosinski, June 1992
! Reviewed: B. Boville, J. Hack, August 1992
! Modified: P. Worley, May 2003
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid, only: plon, plond, plev, plat
use pspect
use comspe
use rgrid, only: nlon, nmmax
use commap, only: w, xm, rsq, cs
use dynconst, only: ez, ra
implicit none
#include <comctl.h>
#include <comfft.h>
!
! Input/Output arguments
!
real(r8), intent(inout) :: u3(plond,plev,plat) ! Fourier -> spec. coeffs. for u-wind
real(r8), intent(inout) :: v3(plond,plev,plat) ! Fourier -> spec. coeffs. for v-wind
!
! Output arguments
!
real(r8), intent(out) :: vort(plond,plev,plat) ! Spectrally truncated vorticity
real(r8), intent(out) :: div(plond,plev,plat) ! Spectrally truncated divergence
!
!---------------------------Local workspace-----------------------------
!
real(r8) d_tmp(psp) ! used in spectral truncation of div
real(r8) vz_tmp(psp) ! used in spectral truncation of vort
real(r8) alpn(pspt) ! alp*rsq*xm*ra
real(r8) dalpn(pspt) ! dalp*rsq*ra
real(r8) tmp1 ! vector temporary
real(r8) tmp2 ! vector temporary
real(r8) tmpr ! vector temporary (real)
real(r8) tmpi ! vector temporary (imaginary)
real(r8) zcor ! correction for absolute vorticity
real(r8) zwalp ! zw*alp
real(r8) zwdalp ! zw*dalp
real(r8) zrcsj ! ra/(cos**2 latitude)
real(r8) zw ! w**2
#if ( ! defined USEFFTLIB )
real(r8) work((plon+1)*plev) ! Workspace for fft
#else
real(r8) work((plon+1)*pcray) ! Workspace for fft
#endif
real(r8) zsqcs
integer ir,ii ! indices complex coeffs. of spec. arrs.
integer i,k ! longitude, level indices
integer irow ! latitude pair index
integer latm,latp ! symmetric latitude indices
integer lat
integer m ! longitudinal wavenumber index
integer n ! latitudinal wavenumber index
integer nspec
integer mr,mc ! spectral indices
call print_memusage ('spetru_uv')
!
!-----------------------------------------------------------------------
!
! Compute the quantities which are transformed to spectral space:
! 1. u = u*sqrt(1-mu*mu), u * cos(phi)
! 2. v = v*sqrt(1-mu*mu), v * cos(phi)
!
do lat=1,plat
irow = lat
if (lat.gt.plat/2) irow = plat - lat + 1
zsqcs = sqrt(cs(irow))
do k=1,plev
do i=1,nlon(lat)
u3(i,k,lat) = u3(i,k,lat)*zsqcs
v3(i,k,lat) = v3(i,k,lat)*zsqcs
end do
end do
!
! Transform grid -> fourier
! 1st transform: U,V,T: note contiguity assumptions
! 2nd transform: LN(PS). 3rd transform: surface geopotential
!
call fft991(u3(1,1,lat),work,trig(1,irow),ifax(1,irow),1,plond, &
nlon(lat),plev,-1)
call fft991(v3(1,1,lat),work,trig(1,irow),ifax(1,irow),1,plond, &
nlon(lat),plev,-1)
end do ! lat=1,plat
!
! Multi-level fields: U, V
!
do k=1,plev
!
! Zero spectral arrays
!
vz_tmp(:) = 0.
d_tmp(:) = 0.
!
! Loop over latitude pairs
!
do irow=1,plat/2
latp = irow
latm = plat - irow + 1
zrcsj = ra/cs(irow)
zw = w(irow)*2.
do i=1,2*nmmax(irow)
tmp1 = 0.5*(u3(i,k,latm) - u3(i,k,latp))
tmp2 = 0.5*(u3(i,k,latm) + u3(i,k,latp))
u3(i,k,latm) = tmp1
u3(i,k,latp) = tmp2
tmp1 = 0.5*(v3(i,k,latm) - v3(i,k,latp))
tmp2 = 0.5*(v3(i,k,latm) + v3(i,k,latp))
v3(i,k,latm) = tmp1
v3(i,k,latp) = tmp2
end do
!
! Compute vzmn and dmn
!
do m=1,nmmax(irow)
mr = nstart(m)
mc = 2*mr
do n=1,nlen(m),2
zwdalp = zw*dalp(mr+n,irow)
zwalp = zw*alp (mr+n,irow)
ir = mc + 2*n - 1
ii = ir + 1
d_tmp(ir) = d_tmp(ir) - (zwdalp*v3(2*m-1,k,latm) + &
xm(m)*zwalp*u3(2*m ,k,latp))*zrcsj
d_tmp(ii) = d_tmp(ii) - (zwdalp*v3(2*m ,k,latm) - &
xm(m)*zwalp*u3(2*m-1,k,latp))*zrcsj
vz_tmp(ir) = vz_tmp(ir) + (zwdalp*u3(2*m-1,k,latm) - &
xm(m)*zwalp*v3(2*m ,k,latp))*zrcsj
vz_tmp(ii) = vz_tmp(ii) + (zwdalp*u3(2*m ,k,latm) + &
xm(m)*zwalp*v3(2*m-1,k,latp))*zrcsj
end do
end do
do m=1,nmmax(irow)
mr = nstart(m)
mc = 2*mr
do n=2,nlen(m),2
zwdalp = zw*dalp(mr+n,irow)
zwalp = zw*alp (mr+n,irow)
ir = mc + 2*n - 1
ii = ir + 1
d_tmp(ir) = d_tmp(ir) - (zwdalp*v3(2*m-1,k,latp) + &
xm(m)*zwalp*u3(2*m ,k,latm))*zrcsj
d_tmp(ii) = d_tmp(ii) - (zwdalp*v3(2*m ,k,latp) - &
xm(m)*zwalp*u3(2*m-1,k,latm))*zrcsj
vz_tmp(ir) = vz_tmp(ir) + (zwdalp*u3(2*m-1,k,latp) - &
xm(m)*zwalp*v3(2*m ,k,latm))*zrcsj
vz_tmp(ii) = vz_tmp(ii) + (zwdalp*u3(2*m ,k,latp) + &
xm(m)*zwalp*v3(2*m-1,k,latm))*zrcsj
end do
end do
enddo ! irow=1,plat/2
!
! Compute grid point values of:u,v,vz, and d.
!
do irow=1,plat/2
latp = irow
latm = plat - irow + 1
zcor = ez*alp(2,irow)
!
! Compute(u,v,vz,d)m
!
do m=1,nmmax(irow)
mr = nstart(m)
do n=1,nlen(m)
!
! These statements will likely not be bfb since xm*ra is now a scalar
!
alpn (mr+n) = alp(mr+n,irow)*rsq(n+m-1)*xm(m)*ra
dalpn(mr+n) = dalp(mr+n,irow)*rsq(n+m-1) *ra
end do
end do
!
! Zero fourier fields
!
u3(:,k,latm) = 0.
u3(:,k,latp) = 0.
v3(:,k,latm) = 0.
v3(:,k,latp) = 0.
vort(:,k,latm) = 0.
vort(:,k,latp) = 0.
div(:,k,latm) = 0.
div(:,k,latp) = 0.
do m=1,nmmax(irow)
mr = nstart(m)
mc = 2*mr
do n=1,nlen(m),2
ir = mc + 2*n - 1
ii = ir + 1
!
tmpr = d_tmp(ir)*alpn(mr+n)
tmpi = d_tmp(ii)*alpn(mr+n)
u3(2*m-1,k,latm) = u3(2*m-1,k,latm) + tmpi
u3(2*m ,k,latm) = u3(2*m ,k,latm) - tmpr
!
tmpr = d_tmp(ir)*dalpn(mr+n)
tmpi = d_tmp(ii)*dalpn(mr+n)
v3(2*m-1,k,latp) = v3(2*m-1,k,latp) - tmpr
v3(2*m ,k,latp) = v3(2*m ,k,latp) - tmpi
!
tmpr = vz_tmp(ir)*dalpn(mr+n)
tmpi = vz_tmp(ii)*dalpn(mr+n)
u3(2*m-1,k,latp) = u3(2*m-1,k,latp) + tmpr
u3(2*m ,k,latp) = u3(2*m ,k,latp) + tmpi
!
tmpr = vz_tmp(ir)*alpn(mr+n)
tmpi = vz_tmp(ii)*alpn(mr+n)
v3(2*m-1,k,latm) = v3(2*m-1,k,latm) + tmpi
v3(2*m ,k,latm) = v3(2*m ,k,latm) - tmpr
!
tmpr = d_tmp(ir)*alp(mr+n,irow)
tmpi = d_tmp(ii)*alp(mr+n,irow)
div(2*m-1,k,latm) = div(2*m-1,k,latm) + tmpr
div(2*m ,k,latm) = div(2*m ,k,latm) + tmpi
!
tmpr = vz_tmp(ir)*alp(mr+n,irow)
tmpi = vz_tmp(ii)*alp(mr+n,irow)
vort(2*m-1,k,latm) = vort(2*m-1,k,latm) + tmpr
vort(2*m ,k,latm) = vort(2*m ,k,latm) + tmpi
end do
end do
do m=1,nmmax(irow)
mr = nstart(m)
mc = 2*mr
do n=2,nlen(m),2
ir = mc + 2*n - 1
ii = ir + 1
!
tmpr = d_tmp(ir)*alpn(mr+n)
tmpi = d_tmp(ii)*alpn(mr+n)
u3(2*m-1,k,latp) = u3(2*m-1,k,latp) + tmpi
u3(2*m ,k,latp) = u3(2*m ,k,latp) - tmpr
!
tmpr = d_tmp(ir)*dalpn(mr+n)
tmpi = d_tmp(ii)*dalpn(mr+n)
v3(2*m-1,k,latm) = v3(2*m-1,k,latm) - tmpr
v3(2*m ,k,latm) = v3(2*m ,k,latm) - tmpi
!
tmpr = vz_tmp(ir)*dalpn(mr+n)
tmpi = vz_tmp(ii)*dalpn(mr+n)
u3(2*m-1,k,latm) = u3(2*m-1,k,latm) + tmpr
u3(2*m ,k,latm) = u3(2*m ,k,latm) + tmpi
!
tmpr = vz_tmp(ir)*alpn(mr+n)
tmpi = vz_tmp(ii)*alpn(mr+n)
v3(2*m-1,k,latp) = v3(2*m-1,k,latp) + tmpi
v3(2*m ,k,latp) = v3(2*m ,k,latp) - tmpr
!
tmpr = d_tmp(ir)*alp(mr+n,irow)
tmpi = d_tmp(ii)*alp(mr+n,irow)
div(2*m-1,k,latp) = div(2*m-1,k,latp) + tmpr
div(2*m ,k,latp) = div(2*m ,k,latp) + tmpi
!
tmpr = vz_tmp(ir)*alp(mr+n,irow)
tmpi = vz_tmp(ii)*alp(mr+n,irow)
vort(2*m-1,k,latp) = vort(2*m-1,k,latp) + tmpr
vort(2*m ,k,latp) = vort(2*m ,k,latp) + tmpi
end do
end do
!
! Correction to get the absolute vorticity.
!
vort(1,k,latp) = vort(1,k,latp) + zcor
!
! Recompute real fields from symmetric and antisymmetric parts
!
do i=1,nlon(latm)+2
tmp1 = u3(i,k,latm) + u3(i,k,latp)
tmp2 = u3(i,k,latm) - u3(i,k,latp)
u3(i,k,latm) = tmp1
u3(i,k,latp) = tmp2
!
tmp1 = v3(i,k,latm) + v3(i,k,latp)
tmp2 = v3(i,k,latm) - v3(i,k,latp)
v3(i,k,latm) = tmp1
v3(i,k,latp) = tmp2
!
tmp1 = vort(i,k,latm) + vort(i,k,latp)
tmp2 = vort(i,k,latm) - vort(i,k,latp)
vort(i,k,latm) = tmp1
vort(i,k,latp) = tmp2
!
tmp1 = div(i,k,latm) + div(i,k,latp)
tmp2 = div(i,k,latm) - div(i,k,latp)
div(i,k,latm) = tmp1
div(i,k,latp) = tmp2
end do
enddo ! irow=1,plat/2
enddo ! k=1,plev
!
do lat=1,plat
!
! Transform Fourier -> grid, obtaining spectrally truncated
! grid point values.
!
irow = lat
if (lat.gt.plat/2) irow = plat - lat + 1
call fft991(u3(1,1,lat),work,trig(1,irow),ifax(1,irow),1,plond, &
nlon(lat),plev,+1)
call fft991(v3(1,1,lat),work,trig(1,irow),ifax(1,irow),1,plond, &
nlon(lat),plev,+1)
call fft991(vort(1,1,lat),work,trig(1,irow),ifax(1,irow),1, &
plond,nlon(lat),plev,+1)
call fft991(div(1,1,lat),work,trig(1,irow),ifax(1,irow),1,plond, &
nlon(lat),plev,+1)
!
! Convert U,V to u,v
!
zsqcs = sqrt(cs(irow))
do k=1,plev
do i=1,nlon(lat)
u3(i,k,lat) = u3(i,k,lat)/zsqcs
v3(i,k,lat) = v3(i,k,lat)/zsqcs
end do
end do
enddo
return
end subroutine spetru_uv
end module spetru