spetru.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

                  xm(m)*zwalp*v3(2*m-1,k,latm))*zrcsj
            end do
         end do
#else
         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(ir,k) = d(ir,k) - (zwdalp*v3(2*m-1,k,latm) + &
                  xm(m)*zwalp*u3(2*m  ,k,latp))*zrcsj
               d(ii,k) = d(ii,k) - (zwdalp*v3(2*m  ,k,latm) - &
                  xm(m)*zwalp*u3(2*m-1,k,latp))*zrcsj
               t(ir,k) = t(ir,k) + zwalp*t3(2*m-1,k,latp)
               t(ii,k) = t(ii,k) + zwalp*t3(2*m  ,k,latp)
               vz(ir,k) = vz(ir,k) + (zwdalp*u3(2*m-1,k,latm) - &
                  xm(m)*zwalp*v3(2*m  ,k,latp))*zrcsj
               vz(ii,k) = vz(ii,k) + (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(ir,k) = d(ir,k) - (zwdalp*v3(2*m-1,k,latp) + &
                  xm(m)*zwalp*u3(2*m  ,k,latm))*zrcsj
               d(ii,k) = d(ii,k) - (zwdalp*v3(2*m  ,k,latp) - &
                  xm(m)*zwalp*u3(2*m-1,k,latm))*zrcsj
               t(ir,k) = t(ir,k) + zwalp*t3(2*m-1,k,latm)
               t(ii,k) = t(ii,k) + zwalp*t3(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
150      continue                ! k=1,plev
160      continue                  ! 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)))
#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 360 irow=1,plat/2
            latp = irow
            latm = plat - irow + 1
#if ( defined PVP )
            zcor = ez*alp(nalp(2)+1,irow)
#else
            zcor = ez*alp(2,irow)
#endif
!
! Zero fourier fields
!
            phis(:,latm) = 0.
            phis(:,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.

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

            div(:,:plev,latm) = 0.
            div(:,:plev,latp) = 0.
!
! Compute(phi*,u,v,ln(p*),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
!
                  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
!
                  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
!
                  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
!
                  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)
            end do

            do 270 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)