subroutine INTLZ (IM,JM,DX,DY,DXJ,PHI)

!     ******************************************************************
!     *                                                                *
!     *  ROUTINE FOR INITIALIZATION OF CONSTANTS AND VARIABLES         *
!     *                                                                *
!     *               (ETA MODEL)                                      *
!     ******************************************************************
!----------------------------------------------------------------------

 DIMENSION DXJ(JM)
 DIMENSION PHI(IM,JM)

!--------------HORIZONTAL GRID CONSTANTS--------------------------------

 INCLUDE 'griddef'

!--------------UNIVERSAL CONSTANTS--------------------------------------
                             DATA &
     A/6376000./,G/9.800/,TWOM/.00014584/,R/287.04/,CP/1004.6/ &
    ,PI/3.141592654/

!--------------DERIVED GEOMETRICAL CONSTANTS----------------------------

      DTR=PI/180.
      DPR=180./PI
      TPH0=TPH0D*DTR
      WB=WBD*DTR
      SB=SBD*DTR
      DLM=DLMD*DTR
      DPH=DPHD*DTR
      TDLM=DLM+DLM
      TDPH=DPH+DPH
      RDLM=1./DLM
      RDPH=1./DPH

      WBI=WB+TDLM
      SBI=SB+TDPH
      EBI=WB+IM2*TDLM
      ANBI=SB+JM3*DPH

      STPH0=SIN(TPH0)
      CTPH0=COS(TPH0)
!--------------DERIVED HORIZONTAL GRID CONSTANTS------------------------
      DY=A*DPH

      TPH=SB-DPH
 DO J=1,JM
      TPH=TPH+DPH
      DXP=A*DLM*COS(TPH)
      DXJ(J)=DXP
 END DO

! NOTE: for E-grid we switch to 2*dx for (i+1)-(i) distance 

      DX=2.*DXJ(JM/2+1)
      DXMAX=2.*maxval(DXJ)
      write(*,*) 'DXMAX=',DXMAX,'  CHOSEN DX=',DX

!-  get the latitude of the original E-grid
!   (this should be for the analysis grid)
!   (for now, use only an approximation - for use in getuv)

       TPH=SB-DPH
              do J=1,JM
       TPH=TPH+DPH
       TLMV = WB - TDLM + MOD(J,2) * DLM
               do I=1,IM
       TLMV = TLMV + TDLM
       CALL invTLLv(TLMV,TPH,TLM0D,DTR,CTPH0,STPH0,ALM,APH,1)
       PHI(i,j)=APH
               enddo
              enddo

      return
      end subroutine INTLZ

                          SUBROUTINE TLLv &
     (ALMD,APHD,TLM0D,DTR,CTPH0,STPH0,TLM,TPH,n)
!---------
         dimension almd(n),aphd(n),tlm(n),tph(n)
!---------------------------------------------------------------------
         do i=1,n
      RELM=(ALMD(i)-TLM0D)*DTR
      SRLM=SIN(RELM)
      CRLM=COS(RELM)
      APH=APHD(i)*DTR
      SPH=SIN(APH)
      CPH=COS(APH)
      CC=CPH*CRLM
      ANUM=CPH*SRLM
      DENOM=CTPH0*CC+STPH0*SPH
      TLM(i)=ATAN2(ANUM,DENOM)/DTR
      TPH(i)=ASIN(CTPH0*SPH-STPH0*CC)/DTR
         end do
      RETURN
      END subroutine TLLv

                          SUBROUTINE invTLLv &
     (ALMD,APHD,TLM0D,DTR,CTPH0,STPH0,TLM,TPH,n)
!---------
!--------- inverse of TLLv:  input ALMD,APHD is grid lon,lat
!---------                   output is earth lon,lat, TLM,TPH
!---------
         dimension almd(n),aphd(n),tlm(n),tph(n)
!---------------------------------------------------------------------
         do i=1,n
      RELM=ALMD(i)*DTR
      SRLM=SIN(RELM)
      CRLM=COS(RELM)
      APH=APHD(i)*DTR
      SPH=SIN(APH)
      CPH=COS(APH)
      CC=CPH*CRLM
      ANUM=CPH*SRLM
      DENOM=CTPH0*CC-STPH0*SPH
      TLM(i)=tlm0d+ATAN2(ANUM,DENOM)/DTR
      TPH(i)=ASIN(CTPH0*SPH+STPH0*CC)/DTR
         end do
      RETURN
      END SUBROUTINE invTLLv