#include "dims.h"
      SUBROUTINE GLP(ZL,ZP)
      use bndry_module,only: zb,zb2,zba,bnd,bnd2,bnda,ci
      use init_module,only: iter
      use cons_module,only: len1,len2,len3,kmax,kmaxp1,imax,imaxp2,
     |  t0,dt,freq_semidi,re_inv,grav,dzgrav
      implicit none
#include "params.h"
#include "fcom.h"
#include "buff.h"
#include "vscr.h"
#include "index.h"
#include "phys.h"
!     include "strt.h" ! for iter
!
! Args:
      real,intent(out) :: ZL(ZIMXP,ZKMXP),ZP(ZIMXP,ZKMXP)
!
! Local:
      real :: wt
      real :: ztmp(zimxp,zkmxp),barmi(zimxp,zkmxp),tbar(zimxp,zkmxp),
     |  tn(zimxp,zkmxp),tn_nm(zimxp,zkmxp),tn_upd(zimxp,zkmxp),
     |  zlbc0(zimxp),zlbc1(zimxp)
      integer :: ntk,ntnmk,ntnpk,i,nzk,k,nmsk,lat
      COMPLEX EXPT,expt1,expt2,expta
!
      DATA WT/0.225/
!
      NTK=NJ+NT
C     ****     S6=TN+WT*(TNM-2.*TN+TNP)=TBAR
      NTNMK=NJ+NTNM
      NTNPK=NJNP+NT
C     ****     S6=-2.*TN+TNM
      do i=1,len2
        s6(i,1) = f(i,ntk)*(-2.)+f(i,ntnmk)
      enddo
C     ****     S6=(S6+TNP)*WT
      do i=1,len2
        s6(i,1) = (s6(i,1)+f(i,ntnpk))*wt
      enddo
C     ****     S6=S6+TN
      do i=1,len2
        s6(i,1) = s6(i,1)+f(i,ntk)

        tn    (i,1) = f(i,ntk)
        tn_nm (i,1) = f(i,ntnmk)
        tn_upd(i,1) = f(i,ntnpk)
        tbar  (i,1) = s6(i,1)
      enddo

!     call addfsech('TBAR',' ',' ',s6,zimxp,zkmxp,zkmx,j)

C     ****     S7(K)=(MS(K)+MS(K+1))*.5=MBAR
      NMSK=NJ+NMS
      do i=1,len2
        s7(i,1) = (f(i,nmsk)+f(i,nmsk+1))*0.5
        barmi(i,1) = s7(i,1)
      enddo
C     ****     S6=TBAR+T0
      DO 1 K=1,KMAX
        do i=1,len1
          s6(i,k) = s6(i,k)+(0.5*(t0(k)+t0(k+1)))
        enddo
    1 CONTINUE
C     ****     S7=S6/S7=(TBAR+T0)/M
      do i=1,len2
        s7(i,1) = s6(i,1) / s7(i,1)
      enddo
C     ****     S7=(DS*R/G)*S7
      do i=1,len2
        s7(i,1) = (dz/dzgrav)*s7(i,1)
      enddo

!     call addfsech('DZTBAR',' ',' ',s7,zimxp,zkmxp,zkmx,j)

C     ****     STORE NZ TEMPORARILY
      NZK=NJ+NZ
      DO 2 I=1,LEN3
      S6(I,1)=F(I,NZK)
    2 CONTINUE
      do i=1,len1
        zlbc0(i) = s6(i,1) ! save lower boundary
      enddo
C     ****     Z(1)=ZB
      EXPT=CEXP(CI*freq_semidi*dt*ITER)
      expt1 = expt
      DO 3 I=1,LEN1
      F(I,NZK)=REAL(ZB(J)*BND(I)*EXPT)
    3 CONTINUE
C     ****     Z(1) = ZB+ZB2
      EXPT = CEXP(CI*.5*freq_semidi*dt*ITER)
      expt2 = expt
      DO 7 I = 1,LEN1
        F(I,NZK)  = F(I,NZK)+REAL(ZB2(J)*BND2(I)*EXPT)
    7 CONTINUE
C     *****     Z(1) = ZB+ZB2+ZBA
      EXPT = 1.
      expta = expt
      DO 8 I = 1,LEN1
        F(I,NZK)  = F(I,NZK)+REAL(ZBA(J)*BNDA(I)*EXPT)
        ztmp(i,1) = f(i,nzk)
        zlbc1(i) = f(i,nzk) ! save lower boundary
    8 CONTINUE

! Compare lower boundaries (should be the same):
!     write(6,"('glp: j=',i3,' zlbc0-zlbc1=',/,(6e12.4))")
!    |  j,zlbc0-zlbc1 

C     ****     Z(K+1)=S7(K)+Z(K)
      DO 4 K=1,KMAX
        do i=1,len1
          f(i,nzk+1) = s7(i,k)+f(i,nzk)
          ztmp(i,k+1) = f(i,nzk+1)
        enddo
        NZK=NZK+1
    4 CONTINUE
C     ****     INSERT PERIODIC POINTS
      NZK=NJ+NZ-1
      DO 6 I=1,2
      DO 6 K=1,KMAXP1
      F(I,NZK+K)=F(I+IMAX,NZK+K)
      F(I+IMAXP2,NZK+K)=F(I+2,NZK+K)
    6 CONTINUE

!     call addfsech('ZTMP',' ',' ',f(1,nj+nz),zimxp,zkmxp,zkmx,j)

      CALL DLDP(NZ,ZL,ZP,LEN3-4)
C     ****     RESTORE NZ
      NZK=NJ+NZ
      DO 5 I=1,LEN3
      F(I,NZK)=S6(I,1)
    5 CONTINUE
      do i=1,len2-4
        zl(i+2,1) = (zl(i+2,1)+zl(i+2,2))*(.5*grav*racs)
      enddo
      do i=1,len2-4
        zp(i+2,1) = (zp(i+2,1)+zp(i+2,2))*(.5*grav*re_inv)
      enddo
C     ****
C     ****     PERIODIC POINTS FOR ZL AND ZP
C     ****
      DO I = 1,2
        DO K = 1,KMAX
          ZL(I,K)= ZL(I+IMAX,K)
          ZL(I+IMAXP2,K)= ZL(I+2,K)
          ZP(I,K)= ZP(I+IMAX,K)
          ZP(I+IMAXP2,K)= ZP(I+2,K)
        ENDDO
      ENDDO
      RETURN
      END
C