#include "dims.h"
!
      subroutine lsqdsq(lat,fjp,fj,fjm,fout,lenx)
      use cons_module,only: cs,dphi,dlamda,re
      implicit none
!
! Subroutine to calculate L**2DEL**2(PSI)
! On input:
!   lat is current latitude index
!   fjp,fj,fjm are input arrays at lat+1, lat, and lat-1 respectively
! On output:
!   fout is defined with result
!
#include "params.h"
!
! Args:
      integer,intent(in) :: lat,lenx
      real,intent(in) :: fjp(lenx+1),fj(lenx+1),fjm(lenx+1)
      real,intent(out) :: fout(lenx+1)
!
! Local:
      real :: flenl,flenp,con1,con2,dpsi,dpsip,dpsim
      integer :: i
!
! Longitude and latitude length functions
! flenl = a*abs(cos(lat))*pi/kut(lat)
! flenp = a*dphi
!
      flenl=re*abs(cs(lat))*dlamda
      flenp=re*dphi
      con1=flenl*flenl/(re*cs(lat)*dlamda)**2
      con2=flenp*flenp/((re*dphi)**2*cs(lat))
!
! fout = L*L/(A*CS(J)*DLAMBDA)**2*(PSI(I+1) - 2*PSI(I) + PSI(I-1))
!
      do i=1,lenx
        fout(i) = ((fj(i+1)*(-2.)+fj(i+2))+fj(i))*con1
!
! DPSIP=(PSI(J+1)-PSI(J))*.5*(CS(J+1)+CS(J))
        dpsip = (fjp(i+1)-fj(i+1))*(0.5*(cs(lat+1)+cs(lat)))
!
! DPSIM=(PSI(J)-PSI(J-1))*.5*(CS(J)+CS(J-1))
        dpsim = (fj(i+1)-fjm(i+1))*(0.5*(cs(lat)+cs(lat-1)))
!
! DPSI = (DPSIP-DPSIM)*(LP**2/(A*A*CS(J)*DPHI*2)
        dpsi = (dpsip-dpsim)*con2
!
! S = S + DPSI = (L**2)(D**2)PSI
        fout(i) = fout(i) + dpsi
      enddo
      return
      end