#include "dims.h"
      SUBROUTINE GLP(ZL,ZP)
      use cons_module,only: len1,len2,len3,kmax,kmaxp1,imax,imaxp2,t0,
     |  freq_semidi,dt,freq_3m3,dzgrav,grav,re_inv
      use input_module,only: ncep,nmc
      use init_module,only: iter
      use bndry_module,only: zb,zb2,zba,bnd,bnd2,bnda,ci
      implicit none
#include "params.h"
#include "fcom.h"
#include "cons3m3.h"
#include "buff.h"
#include "vscr.h"
#include "index.h"
#include "phys.h"
#include "ncep.h"
#include "bndz.h"
!
! Args:
      real,intent(out) :: ZL(ZIMXP,ZKMXP),ZP(ZIMXP,ZKMXP)
!
! Local:
      COMPLEX EXPU,EXPT
      real :: wt
      integer :: ntk,ntnmk,ntnpk,i,nmsk,k,nzk
!
      WT=0.225
C     ****     S6=TN+WT*(TNM-2.*TN+TNP)=TBAR
      NTK=NJ+NT
      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)
      enddo
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
      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
C     ****     STORE NZ TEMPORARILY
      NZK=NJ+NZ
      DO 2 I=1,LEN3
      S6(I,1)=F(I,NZK)
    2 CONTINUE

C     ****     Z(1)=ZB
!
! Do not use ncep boundary:
      if (ncep <= 0 .and. nmc <= 0) then
        EXPT=CEXP(CI*freq_semidi*dt*ITER)
        DO 3 I=1,LEN1
	  F(I,NZK)=ZBND(I,J)+REAL(ZB(J)*BND(I)*EXPT)
C         F(I,NZK)=REAL(ZB(J)*BND(I)*EXPT)
    3   CONTINUE
C     ****     Z(1) = ZB+ZB2
        EXPT = CEXP(CI*.5*freq_semidi*dt*ITER)
        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.
        DO 8 I = 1,LEN1
          F(I,NZK)  = F(I,NZK)+REAL(ZBA(J)*BNDA(I)*EXPT)
    8   CONTINUE
!
! Use ncep Z boundary:
      else
        EXPU = CEXP(CI*freq_semidi*dt*ITER)
        EXPT = CEXP(CI*.5*freq_semidi*dt*ITER)
        if (ncep > 0) then
          DO I = 1,LEN1
            F(I,NZK)  = zncep(I,J,2)+REAL(ZB2(J)*BND2(I)*EXPT) ! new ncep
     +                  +REAL(ZB(J)*BND(I)*EXPU)
          enddo
        elseif (nmc > 0) then
          DO I = 1,LEN1
            F(I,NZK)  = znmc(I,J,2)+REAL(ZB2(J)*BND2(I)*EXPT) ! old nmc
     +                  +REAL(ZB(J)*BND(I)*EXPU)
          enddo
        endif 
      endif
C     ****
C     ****     Add effect of (3,-3) 2-day wave
C     ****
      EXPT = CEXP(CI*freq_3m3*dt*ITER)
      DO I = 1,LEN1
	F(I,NZK) = F(I,NZK)+ REAL(ZB3M3(J)*BND3M3(I)*EXPT)
      ENDDO
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)
        enddo
      NZK=NZK+1
    4 CONTINUE
C     ****     INSERT PERIODIC POINTS
      NZK=NJ+NZ-1
      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 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