spetru.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

                     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
!
                     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
!
                     tmpr = vz(ir,k)*alp(ialp+m,irow)
                     tmpi = vz(ii,k)*alp(ialp+m,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 n=2,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,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
!
                     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
!
                     tmpr = vz(ir,k)*alp(ialp+m,irow)
                     tmpi = vz(ii,k)*alp(ialp+m,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
#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
!
                     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
!
                     tmpr = vz(ir,k)*alp(mr+n,irow)
                     tmpi = vz(ii,k)*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(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
!
                     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
!
                     tmpr = vz(ir,k)*alp(mr+n,irow)
                     tmpi = vz(ii,k)*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
#endif
270            continue                ! k=1,plev
!
! 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 = 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 = 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
               end do
360            continue                  ! irow=1,plat/2
!
! 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 370 lat=1,plat
!     
! Transform Fourier -> grid, obtaining spectrally truncated
! grid point values.
! 1st transform: U,V,T: note contiguity assumptions
! 2nd: ln(PS). 3rd: PHIS. 4th: longitudinal derivative of ln(PS)
! 5th: meridional derivative of ln(PS)
! 6th: vorticity. 7th: 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(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
!
! Convert from ln(ps) to ps
!
                  do i=1,nlon(lat)
                     ps(i,lat) = exp(ps(i,lat))
                  end do
370               continue

                  return
               end subroutine spetru