spetru.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

           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,latp) = u3(2*m-1,k,latp) + tmpi
              u3(2*m  ,k,latp) = u3(2*m  ,k,latp) - tmpr
!     
              tmpr = d(ir,k)*dalpn(ialp+m)
              tmpi = d(ii,k)*dalpn(ialp+m)
              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(ir,k)*dalpn(ialp+m)
              tmpi = vz(ii,k)*dalpn(ialp+m)
              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(ir,k)*alpn(ialp+m)
              tmpi = vz(ii,k)*alpn(ialp+m)
              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 = t(ir,k)*alp(ialp+m,irow)
              tmpi = t(ii,k)*alp(ialp+m,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
              tl(2*m-1,k,latp) = tl(2*m-1,k,latp) - tmpi*ra
              tl(2*m  ,k,latp) = tl(2*m  ,k,latp) + tmpr*ra
!
              tmpr = t(ir,k)*dalp(ialp+m,irow)
              tmpi = t(ii,k)*dalp(ialp+m,irow)
              tm(2*m-1,k,latm) = tm(2*m-1,k,latm) + tmpr*ra
              tm(2*m  ,k,latm) = tm(2*m  ,k,latm) + tmpi*ra
!
              tmpr = q(ir,k)*alp(ialp+m,irow)
              tmpi = q(ii,k)*alp(ialp+m,irow)
              ql(2*m-1,k,latp) = ql(2*m-1,k,latp) - tmpi*ra
              ql(2*m  ,k,latp) = ql(2*m  ,k,latp) + tmpr*ra
!
              tmpr = q(ir,k)*dalp(ialp+m,irow)
              tmpi = q(ii,k)*dalp(ialp+m,irow)
              qm(2*m-1,k,latm) = qm(2*m-1,k,latm) + tmpr*ra
              qm(2*m  ,k,latm) = qm(2*m  ,k,latm) + tmpi*ra
!
              tmpr = d(ir,k)*alp(ialp+m,irow)
              tmpi = d(ii,k)*alp(ialp+m,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
           end do
        end do
#else
        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(ir,k)*alpn(mr+n)
              tmpi = d(ii,k)*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(ir,k)*dalpn(mr+n)
              tmpi = d(ii,k)*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(ir,k)*dalpn(mr+n)
              tmpi = vz(ii,k)*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(ir,k)*alpn(mr+n)
              tmpi = vz(ii,k)*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 = t(ir,k)*alp(mr+n,irow)
              tmpi = t(ii,k)*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
              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(mr+n,irow)
              tmpi = t(ii,k)*dalp(mr+n,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(mr+n,irow)
              tmpi = q(ii,k)*alp(mr+n,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(mr+n,irow)
              tmpi = q(ii,k)*dalp(mr+n,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(mr+n,irow)
              tmpi = d(ii,k)*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
           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(ir,k)*alpn(mr+n)
              tmpi = d(ii,k)*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(ir,k)*dalpn(mr+n)
              tmpi = d(ii,k)*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(ir,k)*dalpn(mr+n)
              tmpi = vz(ii,k)*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(ir,k)*alpn(mr+n)
              tmpi = vz(ii,k)*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 = t(ir,k)*alp(mr+n,irow)
              tmpi = t(ii,k)*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
              tl(2*m-1,k,latp) = tl(2*m-1,k,latp) - tmpi*ra
              tl(2*m  ,k,latp) = tl(2*m  ,k,latp) + tmpr*ra
!
              tmpr = t(ir,k)*dalp(mr+n,irow)
              tmpi = t(ii,k)*dalp(mr+n,irow)
              tm(2*m-1,k,latm) = tm(2*m-1,k,latm) + tmpr*ra
              tm(2*m  ,k,latm) = tm(2*m  ,k,latm) + tmpi*ra
!
              tmpr = q(ir,k)*alp(mr+n,irow)
              tmpi = q(ii,k)*alp(mr+n,irow)
              ql(2*m-1,k,latp) = ql(2*m-1,k,latp) - tmpi*ra
              ql(2*m  ,k,latp) = ql(2*m  ,k,latp) + tmpr*ra
!
              tmpr = q(ir,k)*dalp(mr+n,irow)
              tmpi = q(ii,k)*dalp(mr+n,irow)
              qm(2*m-1,k,latm) = qm(2*m-1,k,latm) + tmpr*ra
              qm(2*m  ,k,latm) = qm(2*m  ,k,latm) + tmpi*ra
!
              tmpr = d(ir,k)*alp(mr+n,irow)
              tmpi = d(ii,k)*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
           end do
        end do
#endif
!
! d(T)/d(lamda)
! d(U)/d(lamda)
! d(V)/d(lamda)
!
!DIR$ IVDEP
        do m=1,nmmax(irow)
           tl(2*m-1,k,latm) = xm(m)*tl(2*m-1,k,latm)
           tl(2*m  ,k,latm) = xm(m)*tl(2*m  ,k,latm)
           tl(2*m-1,k,latp) = xm(m)*tl(2*m-1,k,latp)
           tl(2*m  ,k,latp) = xm(m)*tl(2*m  ,k,latp)
           ql(2*m-1,k,latm) = xm(m)*ql(2*m-1,k,latm)
           ql(2*m  ,k,latm) = xm(m)*ql(2*m  ,k,latm)
           ql(2*m-1,k,latp) = xm(m)*ql(2*m-1,k,latp)
           ql(2*m  ,k,latp) = xm(m)*ql(2*m  ,k,latp)
        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
!
        tmp1 = phisl(i,latm) + phisl(i,latp)
        tmp2 = phisl(i,latm) - phisl(i,latp)
        phisl(i,latm) = tmp1
        phisl(i,latp) = tmp2
!
        tmp1 = phism(i,latm) + phism(i,latp)
        tmp2 = phism(i,latm) - phism(i,latp)
        phism(i,latm) = tmp1
        phism(i,latp) = tmp2
!
        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
!
     do k=1,plev
        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 = 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
!
           tmp1 = tl(i,k,latm) + tl(i,k,latp)
           tmp2 = tl(i,k,latm) - tl(i,k,latp)
           tl(i,k,latm) = tmp1
           tl(i,k,latp) = tmp2
!
           tmp1 = tm(i,k,latm) + tm(i,k,latp)
           tmp2 = tm(i,k,latm) - tm(i,k,latp)
           tm(i,k,latm) = tmp1
           tm(i,k,latp) = tmp2
!
           tmp1 = ql(i,k,latm) + ql(i,k,latp)
           tmp2 = ql(i,k,latm) - ql(i,k,latp)
           ql(i,k,latm) = tmp1
           ql(i,k,latp) = tmp2
!
           tmp1 = qm(i,k,latm) + qm(i,k,latp)
           tmp2 = qm(i,k,latm) - qm(i,k,latp)
           qm(i,k,latm) = tmp1
           qm(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
     end do
  end do
!
! 2nd pass through initial data to obtain and merge all untruncated
! fields:read in and store the data which do not need to be spectrally
! truncated, skipping over header records first.  Also complete
! initialization of prognostics.F
!     
!
  do lat=1,plat
!     
! Transform Fourier -> grid, obtaining spectrally truncated
! grid point values.
! 1st transform: U,V,T
! 2nd: ln(PS). 3rd: PHIS. 4th: longitudinal derivative of ln(PS)
! 5th: meridional derivative of ln(PS)
! 6th: divergence
!
     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 (t3   (1,1,lat) ,work     ,trig(1,irow),ifax(1,irow),1       , &
                     plond       ,nlon(lat),plev        ,+1      )
     call fft991 (ps   (1,lat)   ,work     ,trig(1,irow),ifax(1,irow),1       , &
                     plond       ,nlon(lat),1           ,+1      )
     call fft991 (phis (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      )
     call fft991 (div  (1,1,lat) ,work     ,trig(1,irow),ifax(1,irow),1       , &
                     plond       ,nlon(lat),plev        ,+1      )
!
! Still more fft's
!
! 1st: zonal t derivative
! 2nd: meridional t derivative
! 3rd: zonal phis derivative
! 4th: meridional phis derivative
!
     call fft991 (tl   (1,1,lat) ,work     ,trig(1,irow),ifax(1,irow),1       , &
                     plond       ,nlon(lat),plev        ,+1      )
     call fft991 (tm   (1,1,lat) ,work     ,trig(1,irow),ifax(1,irow),1       , &
                     plond       ,nlon(lat),plev        ,+1      )
     call fft991 (ql   (1,1,lat) ,work     ,trig(1,irow),ifax(1,irow),1       , &
                     plond       ,nlon(lat),plev        ,+1      )
     call fft991 (qm   (1,1,lat) ,work     ,trig(1,irow),ifax(1,irow),1       , &
                     plond       ,nlon(lat),plev        ,+1      )
     call fft991 (phisl(1,lat)   ,work     ,trig(1,irow),ifax(1,irow),1       , &
                     plond       ,nlon(lat),1           ,+1      )
     call fft991 (phism(1,lat)   ,work     ,trig(1,irow),ifax(1,irow),1       , &
                     plond       ,nlon(lat),1           ,+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
!
! Convert from ln(ps) to ps
!
     do i=1,nlon(lat)
        ps(i,lat) = exp(ps(i,lat))
     end do
  end do

  return
end subroutine spetru