#include "dims.h"
!
      subroutine addiag(lat,ixt,iprint)
      use bndry_module,only: zb,zb2,zba,bnd,bnd2,bnda,ci
      use init_module,only: iter
      use cons_module,only: kmax,kmaxp1,t0,len1,len2,len3,cs,rmassinv,
     |  dzgrav,dt,freq_semidi
      implicit none
!
! Calculate z, m, and cos(phi)*v at latitude lat. Store
!   results in fg. This is same as addiag, except it reads
!   and defines fg rather than f (and lat rather than lat+2)
! This is called at all latitudes from a multi-tasked loop
!   in advnce. 
!
#include "params.h"
#include "fgcom.h"
#include "vscr.h"
#include "index.h"
!     include "strt.h" ! for iter
!
! Args:
      integer,intent(in) :: lat,ixt,iprint
!
! Local:
      complex :: expt
      integer :: i,k,nmsk,ntk,nzk,npsk,nps2k
      real :: fmin,fmax
!
      if (iprint > 0) write(6,"('enter addiag: lat=',i2)") lat
!
! nvc = cos(phi)*v
!
      do i=1,len2
        fg(i,nvc+1,lat,ixt) = cs(lat)*fg(i,nv+1,lat,ixt)
      enddo
!
! mbar = mean molecular weight (k+1/2):
!
      do i=1,len3
        fg(i,nms+1,lat,ixt) = 1./(fg(i,nps+1,lat,ixt)*rmassinv(1) +
     |    fg(i,nps2+1,lat,ixt)*rmassinv(2)+(1.-fg(i,nps+1,lat,ixt)-
     |    fg(i,nps2+1,lat,ixt))*rmassinv(3))
      enddo
!
! t1 = mbar(k=0) (linear extrapolation)
!
      do i=1,len1
        t1(i) = 1.5*fg(i,nms+1,lat,ixt)-0.5*fg(i,nms+2,lat,ixt)
      enddo
!
! mbar(k) = 0.5*(mbar(k+1/2)+mbar(k-1/2)), k = kmaxp1,2,1
!
      nmsk = nms+kmaxp1+1
      do k=1,kmax
        nmsk = nmsk-1
        do i=1,len1
          fg(i,nmsk,lat,ixt) = 0.5*(fg(i,nmsk,  lat,ixt)+
     |                              fg(i,nmsk-1,lat,ixt))
        enddo
      enddo
!
! mbar(1) = t1
!
      do i=1,len1
        fg(i,nms+1,lat,ixt) = t1(i)
      enddo
!     call addfsech('FNMS',' ',' ',fg(1,nms+1,lat,ixt),
!    |  zimxp,zkmxp,zkmx,lat)
!
! Calculate Z:
!
      do i=1,len2
        s1(i,1) = (fg(i,nms+1,lat,ixt)+fg(i,nms+2,lat,ixt))*0.5
      enddo
!     call addfsech('W1a',' ',' ',s1,zimxp,zkmxp,zkmx,lat)
!
! s2 = t+t0
!
      ntk = nt
      do k=1,kmax
        ntk = ntk+1
        do i=1,len1
          s2(i,k) = fg(i,ntk,lat,ixt)+(.5*(t0(k)+t0(k+1)))
        enddo
      enddo
!
! s1=s2/s1=(t+t0)/m
!
      do i=1,len2
        s1(i,1) = s2(i,1) / s1(i,1)
      enddo
!     call addfsech('TNa',' ',' ',fg(1,nt+1,lat,ixt),zimxp,zkmxp,zkmx,
!    |  lat)
!     call addfsech('W1b',' ',' ',s1,zimxp,zkmxp,zkmx,lat)
!
! s1=(ds*r/g)*s1
!
      do i=1,len2
        s1(i,1) = (dz/dzgrav) * s1(i,1)
      enddo
!     call addfsech('W1c',' ',' ',s1,zimxp,zkmxp,zkmx,lat)
!
! z(1) = zb
      expt=cexp(ci*freq_semidi*dt*iter)
      do i=1,len1
        fg(i,nz+1,lat,ixt) = real(zb(lat)*bnd(i)*expt)
      enddo
!
! Add in effect of (1,1) tidal component to lbc
!
      expt = cexp(ci*.5*freq_semidi*dt*iter)
      do i=1,len1
        fg(i,nz+1,lat,ixt) = fg(i,nz+1,lat,ixt)+
     |                       real(zb2(lat)*bnd2(i)*expt)
      enddo
!
! Add in effect of annual tidal component to lbc
!
      expt = 1.
      do i=1,len1
        fg(i,nz+1,lat,ixt) = fg(i,nz+1,lat,ixt)+
     |                       real(zba(lat)*bnda(i)*expt)
      enddo
!
! z(k+1)=s1(k)+z(k)
!
      nzk = nz
      do k=1,kmax
        nzk = nzk+1
        do i=1,len1
          fg(i,nzk+1,lat,ixt) = s1(i,k)+fg(i,nzk,lat,ixt)
        enddo
      enddo

!     call addfsech('ADIAG_Z',' ',' ',fg(1,nz+1,lat,ixt),
!    |  zimxp,zkmxp,zkmx,lat)

      end subroutine addiag