spetru.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

              q(ir,k)  = q(ir,k) + zwalp*q3(2*m-1,k,latm)
              q(ii,k)  = q(ii,k) + zwalp*q3(2*m  ,k,latm)
              vz(ir,k) = vz(ir,k) + (zwdalp*u3(2*m-1,k,latp) - &
                         xm(m)*zwalp*v3(2*m  ,k,latm))*zrcsj
              vz(ii,k) = vz(ii,k) + (zwdalp*u3(2*m  ,k,latp) + &
                         xm(m)*zwalp*v3(2*m-1,k,latm))*zrcsj
           end do
        end do
#endif
     end do
  end do
!
  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)))
#if ( defined PVP )
     do n=1,pmax
        ic = ncoefi(n) - 1
        do m=1,nm(n)
           nspec=m-1+n
           ft = 1. - (float(nspec)/filtlim)**2
           if (float(nspec) .ge. filtlim) ft = 0.
           phi(1,ic+m) = phi(1,ic+m)*ft
           phi(2,ic+m) = phi(2,ic+m)*ft
        end do
     end do
#else    
     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
#endif   
     call hordif1(rearth  ,phi     )
  end if
!
! Compute grid point values of:phi*,u,v,ln(p*),t,q,vz,d and grad(ln(p*)).
!
  do irow=1,plat/2
     latp = irow
     latm = plat - irow + 1
!
! Zero fourier fields
!
     phis(:,latm) = 0.
     phis(:,latp) = 0.

     phisl(:,latm) = 0.
     phisl(:,latp) = 0.

     phism(:,latm) = 0.
     phism(:,latp) = 0.

     ps(:,latm) = 0.
     ps(:,latp) = 0.

     dpsl(:,latm) = 0.
     dpsl(:,latp) = 0.

     dpsm(:,latm) = 0.
     dpsm(:,latp) = 0.

     u3(:,:plev,latm) = 0.
     u3(:,:plev,latp) = 0.

     v3(:,:plev,latm) = 0.
     v3(:,:plev,latp) = 0.

     t3(:,:plev,latm) = 0.
     t3(:,:plev,latp) = 0.

     tl(:,:plev,latm) = 0.
     tl(:,:plev,latp) = 0.
!
     tm(:,:plev,latm) = 0.
     tm(:,:plev,latp) = 0.
!
     ql(:,:plev,latm) = 0.
     ql(:,:plev,latp) = 0.
!
     qm(:,:plev,latm) = 0.
     qm(:,:plev,latp) = 0.

     div(:,:plev,latm) = 0.
     div(:,:plev,latp) = 0.
!
! Compute(phi*,grad(phi*),u,d(u)/d(lamda),v,d(v)/d(lamda),
! ln(p*),t,grad(T),q,vz,d,grad(ln(p*)))m
!
#if ( defined PVP )
     do n=1,pmax,2
        ic = ncoefi(n) - 1
        ialp = nalp(n)
        do m=1,nmreduced(n,irow)
           phialpr = phi(1,ic+m)*alp(ialp+m,irow)
           phialpi = phi(2,ic+m)*alp(ialp+m,irow)
!
           phis(2*m-1,latm) = phis(2*m-1,latm) + phialpr
           phis(2*m  ,latm) = phis(2*m  ,latm) + phialpi
!
           phisl(2*m-1,latm) = phisl(2*m-1,latm) - phialpi*ra
           phisl(2*m  ,latm) = phisl(2*m  ,latm) + phialpr*ra
!
           phdalpr = phi(1,ic+m)*dalp(ialp+m,irow)
           phdalpi = phi(2,ic+m)*dalp(ialp+m,irow)
!
           phism(2*m-1,latp) = phism(2*m-1,latp) + phdalpr*ra
           phism(2*m  ,latp) = phism(2*m  ,latp) + phdalpi*ra
!
           ir = 2*(ic+m) - 1
           ii = ir + 1
           psalpr = alps(ir)*alp(ialp+m,irow)
           psalpi = alps(ii)*alp(ialp+m,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(ir)*dalp(ialp+m,irow)
           psdalpi = alps(ii)*dalp(ialp+m,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 n=2,pmax,2
        ic = ncoefi(n) - 1
        ialp = nalp(n)
        do m=1,nmreduced(n,irow)
           phialpr = phi(1,ic+m)*alp(ialp+m,irow)
           phialpi = phi(2,ic+m)*alp(ialp+m,irow)
!
           phis(2*m-1,latp) = phis(2*m-1,latp) + phialpr
           phis(2*m  ,latp) = phis(2*m  ,latp) + phialpi
           phisl(2*m-1,latp) = phisl(2*m-1,latp) - phialpi*ra
           phisl(2*m  ,latp) = phisl(2*m  ,latp) + phialpr*ra
!
           phdalpr = phi(1,ic+m)*dalp(ialp+m,irow)
           phdalpi = phi(2,ic+m)*dalp(ialp+m,irow)
!
           phism(2*m-1,latm) = phism(2*m-1,latm) + phdalpr*ra
           phism(2*m  ,latm) = phism(2*m  ,latm) + phdalpi*ra
!
           ir = 2*(ic+m) - 1
           ii = ir + 1
           psalpr = alps(ir)*alp(ialp+m,irow)
           psalpi = alps(ii)*alp(ialp+m,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(ir)*dalp(ialp+m,irow)
           psdalpi = alps(ii)*dalp(ialp+m,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 n=1,pmax
        ne = n - 1
        ialp = nalp(n)
        do m=1,nmreduced(n,irow)
           alpn (ialp+m) =  alp(ialp+m,irow)*rsq(ne+m)*xm(m)*ra
           dalpn(ialp+m) = dalp(ialp+m,irow)*rsq(ne+m)      *ra
        end do
     end do
#else
     do m=1,nmmax(irow)
        mr = nstart(m)
        mc = 2*mr
        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
!
           phisl(2*m-1,latm) = phisl(2*m-1,latm) - phialpi*ra
           phisl(2*m  ,latm) = phisl(2*m  ,latm) + phialpr*ra
!
           phdalpr = phi(1,mr+n)*dalp(mr+n,irow)
           phdalpi = phi(2,mr+n)*dalp(mr+n,irow)
!
           phism(2*m-1,latp) = phism(2*m-1,latp) + phdalpr*ra
           phism(2*m  ,latp) = phism(2*m  ,latp) + phdalpi*ra
!
           ir = mc + 2*n - 1
           ii = ir + 1
           psalpr = alps(ir)*alp(mr+n,irow)
           psalpi = alps(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(ir)*dalp(mr+n,irow)
           psdalpi = alps(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
           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
           phisl(2*m-1,latp) = phisl(2*m-1,latp) - phialpi*ra
           phisl(2*m  ,latp) = phisl(2*m  ,latp) + phialpr*ra
!
           phdalpr = phi(1,mr+n)*dalp(mr+n,irow)
           phdalpi = phi(2,mr+n)*dalp(mr+n,irow)
!
           phism(2*m-1,latm) = phism(2*m-1,latm) + phdalpr*ra
           phism(2*m  ,latm) = phism(2*m  ,latm) + phdalpi*ra
!
           ir = mc + 2*n - 1
           ii = ir + 1
           psalpr = alps(ir)*alp(mr+n,irow)
           psalpi = alps(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(ir)*dalp(mr+n,irow)
           psdalpi = alps(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)
        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
#endif
     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)
        phisl(2*m-1,latm) = xm(m)*phisl(2*m-1,latm)
        phisl(2*m  ,latm) = xm(m)*phisl(2*m  ,latm)
        phisl(2*m-1,latp) = xm(m)*phisl(2*m-1,latp)
        phisl(2*m  ,latp) = xm(m)*phisl(2*m  ,latp)
     end do

     do k=1,plev
#if ( defined PVP )
        do n=1,pmax,2
           ic = ncoefi(n) - 1
           ialp = nalp(n)
!DIR$ IVDEP
           do m=1,nmreduced(n,irow)
              ir = 2*(ic+m) - 1
              ii = ir + 1
!
              tmpr = d(ir,k)*alpn(ialp+m)
              tmpi = d(ii,k)*alpn(ialp+m)
              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(ir,k)*dalpn(ialp+m)
              tmpi = d(ii,k)*dalpn(ialp+m)
              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(ir,k)*dalpn(ialp+m)
              tmpi = vz(ii,k)*dalpn(ialp+m)
              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(ir,k)*alpn(ialp+m)
              tmpi = vz(ii,k)*alpn(ialp+m)
              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 = t(ir,k)*alp(ialp+m,irow)
              tmpi = t(ii,k)*alp(ialp+m,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
              tl(2*m-1,k,latm) = tl(2*m-1,k,latm) - tmpi*ra
              tl(2*m  ,k,latm) = tl(2*m  ,k,latm) + tmpr*ra
!
              tmpr = t(ir,k)*dalp(ialp+m,irow)
              tmpi = t(ii,k)*dalp(ialp+m,irow)
              tm(2*m-1,k,latp) = tm(2*m-1,k,latp) + tmpr*ra
              tm(2*m  ,k,latp) = tm(2*m  ,k,latp) + tmpi*ra
!
              tmpr = q(ir,k)*alp(ialp+m,irow)
              tmpi = q(ii,k)*alp(ialp+m,irow)
              ql(2*m-1,k,latm) = ql(2*m-1,k,latm) - tmpi*ra
              ql(2*m  ,k,latm) = ql(2*m  ,k,latm) + tmpr*ra
!
              tmpr = q(ir,k)*dalp(ialp+m,irow)
              tmpi = q(ii,k)*dalp(ialp+m,irow)
              qm(2*m-1,k,latp) = qm(2*m-1,k,latp) + tmpr*ra
              qm(2*m  ,k,latp) = qm(2*m  ,k,latp) + tmpi*ra
!
              tmpr = d(ir,k)*alp(ialp+m,irow)
              tmpi = d(ii,k)*alp(ialp+m,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
           end do
        end do
!
        do n=2,pmax,2