      subroutine eta2eta(hin,t,u,v,q,psfc,zsfc,wet,sst,si,erefse,idate,
     &           keyref,ref,etop)
C
C       SUBPROGRAM DOCUMENTAION BLOCK
C
C  SUBPROGRAM:  ETA2ETA               SIGMA-TO-ETA INTERPOLATION PROGRAM
C    AUTHOR: FEDOR MESINGER           ORG: W/NMC22        DATE: ??-??-??
C
C
C  ABSTRACT:  SLIGHTLY ADJUSTS ETA MODEL GUESS HISTORY VARIABLES NAD
C             RECOMPUTES HEIGHT, TEMPERATURES, AND WINDS ON ETA LEVELS
C             TP BE CONSISTENT WITH SURFACE PRESSURE
C
C
C  PROGRAM HISTORY LOG:
C    ??-??-??  F MESINGER WROTE ORIGINAL PTETA CODE
C    ??-??-??  T BLACK MADE NUMEROUS REVISONS
C    90-03-14  E ROGERS REWROTE CODE TO MAKE INDEXING CONSISTENT WITH
C              ETA MODEL FORECAST CODES
C    ??-??-??  D DEAVEN AND E ROGERS REVISED OLD PTETA CODE 
C
C---------------------------------------------------------------------
      INCLUDE "parmeta.h"
C---------------------------------------------------------------------
                             D I M E N S I O N
     1  H(IMJM,LDMZ), T(IMJM,LDM), U(IMJM,LDM), V(IMJM,LDM)
     2, Q(IMJM,LDM), HETA(IMJM,LM), QETA(IMJM,LM), QETAIJ(LM)
     3, PD(IMJM), UETA(IMJM,LM), VETA(IMJM,LM), TETA(IMJM,LM)
     4, HETAIJ(LM), UETAIJ(LM), VETAIJ(LM), PMVP(LDM), HIN(IMJM,LDMZ)
     5, ALPR(LDMZ), HREF(LDM), WF(IMJM), SM(IMJM),hsave(imjm,ldmz)
     6, PSFC(IMJM), WET(IMJM), SST(IMJM), SI(IMJM),psave(imjm)
     7, PINT(IMJM,LDMZ), PMID(IMJM,LDM), EXNL(IMJM),lrof(imjm)
     8, ZSFC(IMJM), ALP(LDMZ), ALPSQ(LDMZ), LDT1(LM)
     9, DLT(LM), PDB(KB,2), TB(KB,LM,2), UB(KB,LM,2)
     O, VB(KB,LM,2), QB(KB,LM,2), FIX(IMJM), IMATCH(IMJM)
     1, HGT(IMJM), PHIS(IMJM), SIGMID(LDM), SIGINT(LDMZ)
     2, ETA(LMP1), ETAL(LM), DETA(LM), ZETA(LMP1), REF(IMJM), IDAT(3)
     3, HFP(2,2), HFPA(IMJM), HRSD(IMJM), DZ(IMJM), LHGT(IMJM)
     4, PETAS(LM),PERIC(IMJM),IDATE(4),erefse(imjm),tsave(imjm,ldm)

C
    
    
                             L O G I C A L
     &  PRINT, TEST, RUN, IDTBLS, FPMNTS, SILUET, SIGMA, MRMSXL
     &, NOMNTS
C
C                          D A T A
C    1 SIGMID/.22329E-1,.827399E-1,.1463738,.2144695,.2859137,.3595226
C    2, .4341393,.5086423,.581962,.6530743,.7210065,.7848307,.8436723
C    3, .89671,.9431677,.9823040/
C
C                          D A T A
C    1 SIGINT/0.,.0538132,.1143270,.1801915,.2501068,.3228228,.3971393
C    2, .4719059,.5460227,.6184392,.6881547,.7542188,.8157315,.8718414
C    3, .9217484,.9647014,1./
C
                           D A T A
     1 LIN /11/, LOROG/14/, LDATE/13/, LOUT/51/
C
c                          E Q U I V A L E N C E
c    & (HETA(1,1),H(1,1)),(UETA(1,1),U(1,1))
c    &, (VETA(1,1),V(1,1)),(h(1,1),hin(1,1))
c    &, (TETA(1,1),T(1,1)),(QETA(1,1),Q(1,1))
C
    4 FORMAT('  HEIGHT OF ISOBARIC SURFACE')
    5 FORMAT('  ZONAL COMPONENT OF VELOCITY')
    6 FORMAT('  MERIDIONAL COMPONENT OF VELOCITY')
    8 FORMAT(' ',//)
    9 FORMAT('1')
C
C--------------UNIVERSAL CONSTANTS--------------------------------------
C
      PI = 3.141592654
      PIHF = 0.5 * PI
      DTR = PI / 180.
C
C--------------HORIZONTAL GRID CONSTANTS--------------------------------
C
C
         do 3349 j = 1,ldmz
         do 3349 i = 1,imjm
         hsave(i,j) = hin(i,j)
         if(j.eq.ldmz) go to 3349
         tsave(i,j) = t(i,j)
3349     continue
C
C --SET UP SWITCHES FOR PRINTING OF OUTPUT, REMOVING MOUNTAINS, ETC.----
C
C        HOUR00 PRESSURE DATA TABLES TO BE PRINTED/NOT PRINTED
C
      print = .true.
      IDTBLS=.FALSE.
C
C        MOUNTAINS TO BE FOUR-POINT ('REGULAR' OR 'SILHOUETTE')
C        AVERAGED, SO AS TO HAVE THE SAME GROUND TOPOGRAPHY AT GROUPS
C        OF FOUR NEIGHBORING HEIGHT POINTS, OR NOT FOUR-POINT AVERAGED
C
      FPMNTS=.false.
      SILUET=.TRUE.
C
C        HOUR00 FIELDS TO BE CREATED FOR A SIGMA MODE/NOT A SIGMA MODE
C        EXPERIMENT
C
      SIGMA=.FALSE.
C
C        IF A SIGMA MODE EXPERIMENT, MOUNTAINS TO BE REMOVED 6 GRIDS
C        LINES ALONG THE BOUNDARY/NOT TO BE REMOVED
C
      MRMSXL=.FALSE.
C
C        NO MOUNTAINS EXPERIMENT/NOT A NO MOUNTAINS EXPERIMENT
C
      NOMNTS=.FALSE.
C
C--------HEIGHTS OF REFERENCE CONTOURS, AND 200 MB CONTOUR INTERVAL------
C
CCC   HREF( 1) = 29920.    CHECK
CCC   HREF( 2) = 17840.
      HREF( 1) = 16160.
      HREF( 2) = 13200.
      HREF( 3) = 11680.
      HREF( 4) = 10160.
      HREF( 5) =  8960.
      HREF( 6) =  7080.
      HREF( 7) =  5440.
      HREF( 8) =  2880.
      HREF( 9) =  1320.
      HREF(10) =     0.
      DIS = 80.
C
c     REWIND LIN
C
C------------- READ SIGMA DATA AND DATE RECORD---------------------------
C
c     READ(LDATE) IDATE
      IDAT(1) = IDATE(2)
      IDAT(2) = IDATE(3)
      IDAT(3) = 1900 + IDATE(4)
      IHRST  =  IDATE(1)
      WRITE(6,9909) (IDAT(I), I = 1, 3), IHRST
 9909 FORMAT(2X,4I6)
c     READ(LIN) HIN, T, U, V, Q, PSFC, ZSFC, WET, SST, SI, erefse
      DO 55 L = 1, LM
      DO K = 1, IMJM
        UETA(K,L) = U(K,L)
        VETA(K,L) = V(K,L)
        QETA(K,L) = Q(K,L)
        TETA(K,L) = T(K,L)
      ENDDO
   55 CONTINUE
C
C
      DO 78845 K = 1, IMJM
        PSFC(K) = PSFC(K) * 100.0
78845 CONTINUE
C
C------- DEFINITION OF CONSTANTS NEEDED FOR VERTICAL INTERPOLATION ---
C
      G = 9.80
      R = 287.04
      CP = 1004.6
      P00 = 1.E5
      GAMMA = 0.0065
      PRF0 = 101325.0
      T0 = 288.0
      GORG = G / (R * GAMMA)
      RGOG = 1.0 / GORG
      ROG = R / G
      CPOG = CP / G
      CAPA = R / CP
      PT= etop * 100.0
      PETAT = PT
      DETA2 = 0.0
      ALPT = ALOG(PT)
C
C------- DEPTH OF ETA LAYERS ----------------------------------------
C
      INCLUDE "deta"
C
C------- COMPUTE ETA AT THE INTERFACES OF THE MODEL LAYERS ------------------
C
      ZMAX = -100000.0
      ZMIN = 100000.0
      DO 7589 K = 1, IMJM
c       IF(ZSFC(K) .LE. 0.004) THEN
c         FIX(K) = 1.0
c       ELSE IF(ZSFC(K) .GT. 0.004) THEN
c         FIX(K) = 0.0
c       END IF
        fix(k) = 0.0
        PD(K) = PSFC(K) - PT
        psave(k) = pd(k)
        PERIC(K) = PSFC(K)
 7589 CONTINUE
      ETA(1) = 0.
      SIGINT(1) = 0.
C
      DO 272 L = 1, LMM1
         ETA(L+1) = ETA(L) + DETA(L)
         SIGINT(L+1) = ETA(L+1)
 272  CONTINUE
C
      ETA(LMP1) = 1.0
      SIGINT(LMP1) = 1.0
      DETA(LM) = 1.0 - ETA(LM)
C
      DO 901 L = 1, LM
         ETAL(L) = 0.5 * (ETA(L) + ETA(L+1))
         SIGMID(L) = ETAL(L)
 901  CONTINUE
C
      DO 904 L = 1, LMP1
         ZETA(L) = T0 * (1.-((PT + ETA(L)*(PRF0 - PT)) / PRF0) ** RGOG)
     1            /GAMMA
 904  CONTINUE



C
C------------- OUTPUT OF READ PRESSURE DATA ---------------------------
C
      IF (.NOT. IDTBLS) GO TO 909
    1 FORMAT(' ',21F6.0)
    2 FORMAT(' ',3X,20F6.0)
    3 FORMAT('  LEVEL (FROM ABOVE) ',I2)
C     PRINT 9
C     PRINT 8
C
C     DO 195 L = 1, LDM, 2
         DO 196 K = 1, IMJM
            WF(K) = H(K,L)
  196    CONTINUE
         Z0 = HREF(L)
         IF(L. GE. 8)  DIS = 40.
         IF(L. GE. 11) DIS = 20.
  195 CONTINUE
C
C------------- VERTICAL INTERPOLATION TO CONSTANT ETA LAYERS ---------
C
C------- DEFINITION OF PRESSURE DATA AND ETA LAYERS CONSTANTS --------
C
c  read in the pre-computed adjusted steps, sea mask and ref
c  file is stored in full cray words
c
  909 rewind lorog
      read (lorog) hgt,sm
      read (lorog) ref
C     PRINT 9
C
      IF(SIGMA) THEN
          DO 929 K = 1, IMJM
          REF(K) = 1.0
  929     CONTINUE
      END IF
      if(keyref.ne.0) return








C
C     FIND HEIGHTS ON SIGMA LAYER INTERFACES
C
          kksetem = 0
      DO 914 K = 1, IMJM
         PRFS = PRF0 * ((T0 - GAMMA * HGT(K)) / T0) ** GORG
         ETAS = (PRFS - PT) / (PRF0 - PT)
         DO 9145 L = 1, LDM
            PMID(K,L) = (ETAL(L) / erefse(k)) * (psfc(k) - PT) + PT
 9145    CONTINUE
*ODIR SKIP
         DO 9146 L = 1, LDMZ
            h(k,l) = 9.e99
            PINT(K,L) = (ETA(L) / erefse(k)) * (psfc(k) - PT) + PT
            if(eta(l)/erefse(k).ge.0.99999) then
                h(k,l) = zsfc(k)
                kksetem = kksetem + 1
            end if
            if(eta(l)/erefse(k).ge.1.0001) then
                kksetem = kksetem - 1
            pint(k,l) = -10.0
            h(k,l) = 9.e99
            end if 
9146    CONTINUE
c        H(K,LDMZ) = 0.004
C        H(K,LDMZ) = ZSFC(K)
  914 CONTINUE
c         print *,"kksetem points is ",kksetem
C
      DO 7001 LL = 1, LDM
         L = LDMZ - LL
         DO 7771 K=1,IMJM
C---------------------------------------------------------------------
            IF(ZSFC(K) .EQ. h(k,l+1)) THEN
               EXNL(K) = (PSFC(K) / P00) ** CAPA
               else if(h(k,l+1).gt.1.e50) then
               go to 7771
            END IF
C---------------------------------------------------------------------
            EXNT = (PINT(K,L) / P00) ** CAPA
c           H(K,L) = H(K,L+1) - CPOG * T(K,L) * (EXNT - EXNL(K)) *
c    1              (P00 / PMID(K,L)) ** CAPA
            h(k,l) = h(k,l+1) + t(k,l)*r*(log(pint(k,l+1)/pint(k,l)))/g
C
            IF(K .EQ. 14 .OR. K .EQ. 100 .OR. K .EQ. 300) THEN
c             WRITE(6,60440) K,L,ETA(L),REF(K),PSFC(K),H(K,L),
c    1           EXNL(K), EXNT, PINT(K,LDMZ)
60440         FORMAT(2X,'K,L= ',2I5,' ETA= ',E11.4,' REF= ',E11.4,
     1         ' PSFC= ',E11.4,' H= ',E11.4,
     2         ' EXNL= ',E11.4,' EXNT= ',E11.4,' PINT(LDMZ)= ',E11.4)
            END IF
C
            EXNL(K) = EXNT
7771       continue
 7001 CONTINUE
C     DO 7009 L = 1, LM
C       WRITE(6,87310) L,ETA(L),ETAL(L),SIGINT(L),SIGMID(L),ZETA(L)
C7310   FORMAT(2X,I3,1X,5(E12.5,1X))
C7009   CONTINUE
c     DO 88844 K = 1, IMJM
c       IF(K .eq. -9 ) then
c         WRITE(6,76574) K, PSFC(K), ZSFC(K), HGT(K), SM(K), FIX(K)
76574     FORMAT(2X,'K(INPUT) =',I5,5(1X,E12.5))
c         DO 88849 L = 1, LM
c          WRITE(6,88845) L,T(K,L),Q(K,L),U(K,L),V(K,L),H(K,L)
c    1       ,PMID(K,L),PINT(K,L)
88845      FORMAT(2X,I2,1X,7(E11.4,1X))
88849     CONTINUE
c       END IF
88844 CONTINUE
C
C
C------------- VERTICAL INTERPOLATION: QUADRATIC INTERPOLATION OF
C------------- HEIGHTS (EQUIVALENT TO TEMPERATURE BEING LINEAR IN LN P)
C------------- AND LINEAR INTERPOLATION OF WINDS
C
C        COMPUTATION OF THE 'SEA LEVEL' PRESSURE DIFFERENCE
C
      DO 514 K = 1, IMJM
c        SPLM = PINT(K,LDMZ)
c        RSPL = SPLM / PINT(K,LDM)
cwas     DO 515 LD = 2, LDM
         DO 515 LD = 1, LDM
            if(pint(k,ld+1).gt.0.0.and.pint(k,ld).gt.0.0) then
            ALPR(LD) = ALOG(PINT(K,LD+1) / PINT(K,LD))
            else
            alpr(ld) = 9.e99
            end if
  515    CONTINUE
cwas     DO 530 LD = 2, LDMZ
         DO 530 LD = 1, LDMZ
            if(abs(pint(k,ld)-psfc(k)).lt.0.5)  then
               splm = psfc(k)
               rspl = splm / pint(k,ld-1)
            end if
            if(pint(k,ld).gt.0.0) then
            ALP(LD) = ALOG(PINT(K,LD))
            ALPSQ(LD) = ALP(LD)**2
            else
            alp(ld) = 9.e99
            alpsq(ld) = 9.e99
            end if
  530    CONTINUE
c        IF(FIX(K) .EQ. 1.0) GO TO 514
c         WRITE(6,75698) K
c75698     FORMAT(2X,'FIXED POINT = ',I5)
C
C
         LDMZ1 = LDMZ - 1
cwas         DO 5301 LS = 2, LDMZ1
c        DO 5301 LS = 2, LDMZ
c           IF(K .EQ. 5214) THEN
c             WRITE(6,22976) K, LS, H(K,LS), HGT(K)
22976         FORMAT(2X,2(I5,1X),2(E12.5,1X))
c           END IF
cwas            IF(H(K,LS) .LT. HGT(K)) THEN
c           IF(H(K,LS) .LE. HGT(K)) THEN
cwas              L1 = LS - 1
c              l1 = ls - 2
c              GO TO 5302
c           END IF
5301    CONTINUE
         L1 = LDMZ-2
5302     continue
          if(pint(k,l1+2).le.0.0) then
                 l1 = l1 - 1
                 go to 5302
          end if
         DLTPD = ALPR(L1+1)*ALOG(PINT(K,L1+2) / PINT(K,L1))
     1          * (-ALPR(L1))
C
         L2 = L1 + 1
         L3 = L1 + 2
         H1a = H(K,L1)
         H2 = H(K,L2)
         H3 = H(K,L3)
         H3M2 = H3 - H2
         H3M1 = H3 - H1a
         B = (H3M2*(ALPSQ(L3)-ALPSQ(L1))-H3M1*(ALPSQ(L3)-ALPSQ(L2)))/
     1       DLTPD
         C = (H3M1*(ALP(L3)-ALP(L2))-H3M2*(ALP(L3)-ALP(L1)))/DLTPD
C
         IF(K .EQ. 5214) THEN
c           WRITE(6,33788) K, L, B, C, DLT(L)
c           WRITE(6,33789) ALPSQ(L1), ALPSQ(L2), H3M2, H3M1
c           WRITE(6,33790) H1a, H2, H3, L1
33788       FORMAT(' K=',I4,' L=',I2,' B=',E12.5,' C=',E12.5,
     1            ' DLT=',E12.5)
33789       FORMAT(' ALPSQ(L1)=',E12.5,' ALPSQ(L2)=',E12.5,' H3M2=',
     1            E12.5,' H3M1=',E12.5)
33790       FORMAT(' H1=',E12.5,' H2=',E12.5,' H3=',E12.5,' L1=',I2)
         END IF
C
         IF (ABS(C) .LT. 0.00001) GO TO 532
         ALPGT = (-B-SQRT(B**2-4.*C*(-C*ALPSQ(L3)-B*ALP(L3)+H3-
     1           HGT(K))))/(2.*C)
         PDIJ = (EXP(ALPGT)-PT)/REF(K)
C
         IF(K .EQ. 5214) THEN
c           WRITE(6,88411) K, L, L1, L2, PDIJ
c           WRITE(6,88412) ALPSQ(L3), ALP(L3), HGT(K)
c           WRITE(6,88413) ALPGT, REF(K)
88411       FORMAT(' 1ST  K=',I4,' L=',I2,' L1=',I2,' L2=',I2,
     1       ' PDIJ=',E12.5)
88412       FORMAT(' ALPSQ(L3)=',E12.5,' ALP(L3)=',E12.5,' HGT=',E12.5)
88413       FORMAT(' ALPGT=',E12.5,' REF=',E12.5)
         END IF
C
         GO TO 533
  532    HLDM = H(K,LDM)
         write(6,*)"oh oh sfc trouble"
CWAS     PDIJ = (SPLM*RSPL**((H(K,LDMM)-HGT(K))/
CWAS 1       (H(K,LDM-2)-H(K,LDMM)))-PT)/REF(K)
         PDIJ = (SPLM*RSPL**((H(K,L2)-HGT(K))/
     1       (H(K,L3)-H(K,L2)))-PT)/REF(K)
C
c        IF(K .EQ. 5214) THEN
            WRITE(6,94765) K, L, L1, L2, L3, H1a, H2, H3, PDIJ
            WRITE(6,44765) K, PDIJ, LDMM, SPLM, HLDM
            WRITE(6,44766) RSPL, HGT(K), H(K,LDMM), PT
44765       FORMAT(' 2ND  K=',I4,' PDIJ=',E12.5,' LDMM=',I2,
     1             ' SPLM=',E12.5,' HLDM=',E12.5)
44766       FORMAT(' RSPL=',E12.5,' HGT=',E12.5,' H(LDMM)=',E12.5,
     1            ' PT=',E12.5)
94765       FORMAT(2X,5(I5,1X),4(E12.5,1X))
c        END IF
C
C        PDIJ HAS BEEN CHANGED TO MAKE IT EQUAL TO THE DIFFERENCE
C        BETWEEN P*, PRESSURE EXTRAPOLATED TO THE ETA=1 SURFACE, AND PT
C
C------- COMPUTATION OF HEIGHTS AT UPPER BOUNDARIES OF ETA LAYERS -----
C------- AND TEMPERATURE AND SPEC HUMIDITY AT MIDDLE OF ETA LAYERS -----
C
  533    CONTINUE
         IF(K .EQ. 5214) THEN
c           WRITE(6,94766) K, L, L1, L2, L3, H1a, H2, H3, PDIJ
c           WRITE(6,44767) K, PDIJ, LDMM, SPLM, HLDM
c           WRITE(6,44768) RSPL, HGT(K), H(K,LDMM), PT
44767       FORMAT(' 2ND  K=',I4,' PDIJ=',E12.5,' LDMM=',I2,
     1             ' SPLM=',E12.5,' HLDM=',E12.5)
44768       FORMAT(' RSPL=',E12.5,' HGT=',E12.5,' H(LDMM)=',E12.5,
     1            ' PT=',E12.5)
94766       FORMAT(2X,5(I5,1X),4(E12.5,1X))
         END IF
         DETA1 = 0.0
         KNTSIG = 1
         PETA = PT + 0.0001
         DO 6315 L = 1, LM
          IF(L .GT. 1) DETA1 = DETA(L-1)
          PETA = PETA + DETA1 * PDIJ
C         PETAS(L) = PETA
            kntsig = 1
6313     CONTINUE
         IF (PETA .GT. PINT(K,KNTSIG) .AND. PETA.LE.PINT(K,KNTSIG+1)
     1         .OR.KNTSIG.GE.LDMM) THEN
              LDT1(L) = KNTSIG
              GO TO 63188
            END IF
            KNTSIG = KNTSIG + 1
            GO TO 6313
63188         IF(K .EQ. 5214 .or.k .eq.500) THEN
c         WRITE(6,63155) K,L,PETA,PDIJ,PINT(K,KNTSIG),PINT(K,KNTSIG+1),
c    1     KNTSIG
63155     FORMAT(2X,'K,L= ',2I5,' PETA= ',E11.4,' PDIJ=' E11.4,
     1     ' PINT(KNTSIG)= ',E11.4,' PINT(KNTSIG+1)= ',E11.4,
     2     ' KNTSIG= ',I3)
         END IF
63156    continue
            if(pint(k,kntsig+2).le.0.0) then
                   kntsig = kntsig - 1
                   ldt1(l) = kntsig
                   go to 63156
            end if
6315    CONTINUE
         do 531 l = 1,lm
         l1 = ldt1(l)
         dlt(l) = alpr(l1+1)*log(pint(k,l1+2) / pint(k,l1))
     $            * (-alpr(l1))
531      continue
C
         L = LM
  517    L1 = LDT1(L)
         L2 = L1 + 1
         L3 = L1 + 2
         H1a = H(K,L1)
         H2 = H(K,L2)
         H3 = H(K,l3)
         H3M2 = H3 - H2
         H3M1 = H3 - H1a
         B = (h3m2*(alpsq(l3)-alpsq(l1))-h3m1*(alpsq(l3)-alpsq(l2)))/
     $       dlt(l)
         c = (h3m1*(alp(l3)-alp(l2))-h3m2*(alp(l3)-alp(l1)))/dlt(l)
518      arglog = pt + pdij * eta(l)
         alpeta = log(arglog)
         arglg2 = pt + pdij * (eta(l) + 0.5 * deta(l))
         alpm = log(arglg2)
         hetaij(l) = h2 + b*(alpeta-alp(l2)) + c*(alpeta**2-alpsq(l2))
c        alnpik = log(pint(k,l1))
c        alnpim = log(pint(k,l2))
c        alnp3 = 0.5 * (alnpik+alnpim)
c        dlnp = alnpik - alnpim
c        rlnp = 1.0 / dlnp
c        aaa = (h1a - h2) * rlnp
c        bbb = rog * (t(k,l1) - t(k,l2)) * rlnp
c        h3 = 0.5*(h1a+h2) + 0.125 * bbb * dlnp*dlnp
c        peta = eta(l) * pdij + pt
c        dlnp = (log(peta) - alnp3)
c        HETAIJ(L) = h3 + (aaa -0.5 * bbb*dlnp)*dlnp
C
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
         IF(L .EQ. 1) GO TO 516
         L = L - 1
         GO TO 517
C
  516    PD(K) = PDIJ
         DO 5140  L = 1, LM
c           if(abs(hetaij(l)-hgt(k)).lt.9.0) hetaij(l) = hgt(k)
            HETA(K,L) = HETAIJ(L)
 5140    CONTINUE
  514 CONTINUE


C     CALL CMAS2D(PD,'PD  ',7448,2,.TRUE.)
C
      TEST = .FALSE.
      IF (.NOT.PRINT) GO TO 525
  529 DO 916 L = 1, LM
C        PRINT 9
C        PRINT 3, L
C        PRINT 526
 526     FORMAT('  HEIGHT AT THE UPPER BOUNDARY')
C        PRINT 8
  916 CONTINUE
  525 CONTINUE
      IF(TEST) GO TO 630
C
C        COMPUTATION OF WINDS WITHIN ETA LAYERS
C
C    KSL IS THE NORTHWEST CORNER POINT
C
      DO 5190 K = 1, IMJM
         KMIM = K + KSW
         KPIM = K + KNE
         KMIMM1 = K + KSE
         IF(K .GT. IM) KT = MOD(KMIM, IMT)
         KT1 = MOD(K, IMT)
         DO 5215 L = 1, LDM
            IF(K .EQ. IMJM) GO TO 5216
C
C      PMVP ALONG EDGES
C
            IF(K .LT. IM) THEN
               PMVP(L) = 0.5 * (PMID(K,L) + PMID(K+1,L))
               IF(L.EQ.1) PDVP = 0.5 * (PD(K) + PD(K+1))
C
            ELSE IF(K .GE. KSL .AND. K .LT. IMJM)THEN
               PMVP(L) = 0.5 * (PMID(K,L) + PMID(K+1,L))
               IF(L.EQ.1) PDVP = 0.5 * (PD(K) + PD(K+1))
C
            ELSE IF(K.EQ.IM .OR. KT.EQ.0 .OR. KT1.EQ.0) THEN
               PMVP(L) = 0.5 * (PMID(KMIMM1,L) + PMID(KPIM,L))
               IF(L.EQ.1) PDVP = 0.5 * (PD(KMIMM1) + PD(KPIM))
C
            ELSE IF(K .LT. KSL) THEN
C
C      PMVP AT ALL INTERIOR POINTS
C
               PMVP(L) = 0.25 * (PMID(K,L)+PMID(K+KNE,L)+PMID(K+1,L)+
     1                   PMID(K+KSE,L))
               IF(L.EQ.1)
     1            PDVP = 0.25 * (PD(K)+PD(K+KNE)+PD(K+1)+PD(K+KSE))
             END IF
         GO TO 5215
C
C     PMVP AT NORTHEAST CORNER POINT
C
 5216    PMVP(L) = 0.5 * (PMID(K,L) + PMID(K-1,L))
         IF(L.EQ.1) PDVP = 0.5 * (PD(K) + PD(K-1))
 5215    CONTINUE
C
         DO 815 LD = 1, LDMM
            ALPR(LD) = ALOG(PMVP(LD + 1) / PMVP(LD))
  815    CONTINUE
         DO 830 LD = 1, LDM
            ALP(LD) = ALOG(PMVP(LD))
  830    CONTINUE
C
C    COMPUTE WINDS
C
         IF(FIX(K) .EQ. 1.0) GO TO 5190
         DO 523 L = 1 , LM
            IF(L.EQ.1) ALPLB = ALPT
            ALPUB = ALPLB
            ALPLB = ALOG(PT + PDVP * ETA(L+1))
            ALPETA = SQRT(0.5 * (ALPUB**2 + ALPLB**2))
C
            DO 524 LD=2,LDM
C              IF(K .EQ. 1) THEN
C               WRITE(6,2976) L,LD,ALPUB,ALPLB,ALPETA,ALP(LD)
C2976           FORMAT(1X,'L,LD= ',2I3,' ALPUB= ',E11.4,' ALPLB= ',
C    1           E11.4,' ALPETA= ',E11.4,' ALP(LD)= ',E11.4)
C              END IF
               IF (ALPETA.GT.ALP(LD) .AND. LD.LT.LDM) GO TO 524
               ULD = U(K,LD)
               VLD = V(K,LD)
               LDM1 = LD-1
               CF = (ALP(LD) - ALPETA) / ALPR(LDM1)
               UETAIJ(L) = ULD + (U(K,LDM1) - ULD) * CF
               VETAIJ(L) = VLD + (V(K,LDM1) - VLD) * CF
C              IF(K .EQ. 1) THEN
C               WRITE(6,2977) L,LDM1,ULD,VLD,UETAIJ(L),VETAIJ(L),
C    1                        ALPR(LDM1),CF,U(K,LDM1),V(K,LDM1)
C2977           FORMAT(1X,'L,LDM= ',2I3,' ULD,VLD= ',2(E11.4,1X),
C    1           ' UETA,VETA= ',2(E11.4,1X),' ALPR(LDM1)= ',E11.4,
C    2           ' CF= 'E11.4,' ULDM1,VLDM1= ',2(E11.4,1X))
C              END IF
               GO TO 523
  524       CONTINUE
C
  523    CONTINUE
C
         DO 519 L = 1, LM
            UETA(K,L) = UETAIJ(L)
            VETA(K,L) = VETAIJ(L)
  519    CONTINUE
C
 5190 CONTINUE
C
C------------- GROUND SURFACE HEIGHTS CONVERTED TO GEOPOTENTIALS, AND --
C------------- GEOPOTENTIAL TO TEMPERATURE INVERSION -------------------
C
      DO 501 K = 1, IMJM
         PHIS(K) = G * HGT(K)
C
         IF(FIX(K) .EQ. 1.0) GO TO 501
         PHUB=0.
         IF(SIGMA) PHUB = PHIS(K)
         DO 502 IVI = 1, LM
            L = LMP1 - IVI
            PHLB = PHUB
            PHUB = G * HETA(K,L)
c           TETA(K,L)=-(PHLB-PHUB)*(PT+ETAL(L)*PD(K))/(R*DETA(L)*PD(K))
c           delpp = log(eta(l+1)*pd(k)+pt) - log(eta(l)*pd(k)+pt)
c           teta(k,l)=(phub-phlb)/(r*delpp)
            if(k.eq.-9)write(6,*)"teta,phlb,phub,pt,etal,l,pd,deta  ",
     $         teta(k,l),phlb,phub,pt,etal(l),l,pd(k),deta(l)
            IF(TETA(K,L) .LT. 150.0) TETA(K,L) = 150.0
C
C           IF(K .EQ. 1) THEN
C             WRITE(6,38441) L, TETA(K,L), PHLB, HETA(K,L), PHUB
C             WRITE(6,38442) ETA(L), PD(K), DETA(L), R
C8441         FORMAT(' L=',I2,' TETA=',E12.5,' PHLB=',E12.5,
C    1               ' HETA=',E12.5,' PHUB=',E12.5)
C8442         FORMAT(' ETA=',E12.5,' PD=',E12.5,' DETA=',E12.5,
C    1               ' R=',E12.5)
C           END IF
C
  502    CONTINUE
C
         DO 8150 LD = 1, LDMM
            ALPR(LD) = ALOG(SIGMID(LD+1) / SIGMID(LD))
 8150    CONTINUE
         DO 8300 LD = 1, LDM
cwas        ALP(LD) = ALOG(SIGMID(LD) * PSFC(K))
            ALP(LD) = ALOG(pmid(k,ld))
 8300    CONTINUE
C
C     INTERPOLATE Q LINEARLY IN LN(P) FROM SIGMA TO ETA.
C
         ALPLB = ALPT
         DO 823 L = 1, LM
         IF(FIX(K) .EQ. 1.0) GO TO 823
            ALPUB = ALPLB
            ALPLB = ALOG(PT + PD(K) * ETA(L+1))
            ALPETA = SQRT(0.5 * (ALPUB**2 + ALPLB**2))
            qetaij(l) = q(k,l)
            DO 824 LD = 2, LDM
C           IF(K .EQ. 1) THEN
C            WRITE(6,82439) L,LD,ALPUB,ALPLB,ALPETA,ALP(LD)
C2439        FORMAT(1X,'L,LD=',2I4,' ALPUB= ',E11.4,' ALPLB= ',E11.4,
C    1        ' ALPETA= ',E11.4,' ALP(LD)= ',E11.4)
C           END IF
               IF(ALPETA.GT.ALP(LD).AND.LD.LT.LDM) GO TO 824
               IF(ALPETA.LE.ALP(LD)) THEN
c                if(pint(k,ld).lt.0.0) then
c                ldtmp = ld
c24397            ldtmp = ldtmp - 1
c                   if(pint(k,ldtmp).lt.0.0) go to 24397
c                   qetaij(l) = q(k,ldtmp)
c                else
                    QLD = Q(K,LD)
                    LDM1 = LD - 1
                    CF = (ALP(LD) - ALPETA) / ALPR(LDM1)
c                   QETAIJ(L) = QLD + (Q(K,LDM1) - QLD) * CF
c                   QETAIJ(L) = Q(k,ld+1)
c                end if
C                IF(K .EQ. 1) THEN
C                 WRITE(6,82430) K,L,LD,QLD,Q(K,LDM1),ALP(LD),ALPETA,
C    1             ALPR(LDM1),CF,QETAIJ(L)
C2430             FORMAT(1X,'K,L,LD= ',3I5,' Q(LD)= ',E11.4,
C    1          ' Q(LD-1)= ',E11.4,' ALP(LD)= ',E11.4,' ALPETA= ',E11.4,
C    2          ' ALPR= ',E11.4,' CF= ',E11.4,' QETAIJ= ',E11.4)
C                END IF
                 CLOGES = -CM1/TETA(K,L)-CM2*ALOG10(TETA(K,L))+CM3
                 ESE = 10.0 ** (CLOGES+2.)
                 QS = EPS*ESE/(ETAL(L)*PD(K)+PT-ESE*(1.-EPS))
                 IF(K .EQ. 1) THEN
c                 WRITE(6,82431) K,L,LD,CLOGES,ESE,QS
82431              FORMAT(1X,'K,L,LD= ',3I5,' CLOGES= ',E11.4,' ESE= ',
     1              E11.4,' QS= ',E11.4)
                 END IF
                 QSMX = 0.95 * QS
                 QETAIJ(L) = AMIN1(QETAIJ(L),QSMX)
                 GO TO 823
C
C     IF THE ETA SURFACE IS BELOW THE LOWEST SIGMA SURFACE,
C     ASSUME THE RH (GIVEN TETA) IS EQUAL TO THAT ON THE
C     LOWEST SIGMA SURFACE. NOT APPLICABLE FOR ETA-ETA CONVERSION
C
               ELSE
                 ldtmp = ld
c                if(pint(k,ld).lt.0.0) then
c24398            ldtmp = ldtmp - 1
c                   if(pint(k,ldtmp).lt.0.0) go to 24398
c                end if
c                   qetaij(l) = q(k,ldtmp+1)
                 CLOGES = -CM1/TETA(K,L)-CM2*ALOG10(TETA(K,L))+CM3
                 ESE = 10.0 ** (CLOGES+2.)
                 QS = EPS*ESE/(ETAL(L)*PD(K)+PT-ESE*(1.-EPS))
                 QSMX = 0.95 * QS
                 qsnow = amin1(qetaij(l),qsmx)
              if(qsnow.ne.qetaij(l).and.ref(k).ge.eta(l+1)) then
c                write(6,*)"q was ,is ",qetaij(l),qsnow,l,ldtmp,ld,k
c                write(6,*)"ref,eta ",ref(k),eta(l+1)
              end if
                 QETAIJ(L) = qsnow
                 go to 823
               END IF
  824       CONTINUE
  823    CONTINUE
C
         DO 819 L = 1, LM
            QETA(K,L) = QETAIJ(L)
            IF(QETA(K,L) .LT. 0.0) QETA(K,L)=0.
            IF(QETA(K,L) .GT. 25.E-3) QETA(K,L) = 25.E-3
  819    CONTINUE
  501 CONTINUE
         do 4350 i = 2705,2967,262
c           write(6,*)"k,erefse,ref  ",i,erefse(i),ref(i)
         do 4350 j = 1,ldm
         tmpdif = teta(i,j) - tsave(i,j)
         tsadif = tsadif + tmpdif
c           write(6,*)"t,tsave,tmpdif,i,j ",teta(i,j),tsave(i,j),
c    $                   tmpdif,i,j
4350     continue
          do 4357 j = 1,ldm
          hin(1,j) = 0.0
          hin(2,j) = 0.0
          do 4357 i = 1,imjm
          hin(1,j) = hin(1,j) + tsave(i,j)
          hin(2,j) = hin(2,j) + teta(i,j)
4357      continue
          do 4356 j = 1,ldm
          hin(1,j) = hin(1,j) / imjm
          hin(2,j) = hin(2,j) / imjm
c         write(6,*)"tsave bar,teta bar,level",hin(1,j),hin(2,j),j
4356      continue
C
C     CONVERT FROM VIRTUAL TO TRUE TEMPERATURE.
C
      DO 6025 K = 1, IMJM
      DO 6025 L = 4, LM
         IF(TETA(K,L) .GT. 315.0) TETA(K,L)=315.0
         IF(TETA(K,L) .LT. 150.0) TETA(K,L)=150.0
         TV = TETA(K,L)
         DO 602 ITER = 1, 3
            TETA(K,L) = TV / (1. + 0.608 * QETA(K,L))
            CLOGES = -CM1/TETA(K,L)-CM2*ALOG10(TETA(K,L))+CM3
            ESE = 10.0 ** (CLOGES+2.)
            DENOM2 = ETAL(L)*PD(K)+PT-ESE*(1.-EPS)
            QS = EPS*ESE/(ETAL(L)*PD(K)+PT-ESE*(1.-EPS))
            QSMX = 0.95 * QS
            IF(QETA(K,L) .LT. QSMX) GO TO 6025
            QETA(K,L) = QSMX
  602    CONTINUE
 6025 CONTINUE
          do 4355 j = 1,ldm
          hin(1,j) = 0.0
          hin(2,j) = 0.0
          do 4355 i = 1,imjm
          hin(1,j) = hin(1,j) + q(i,j)
          hin(2,j) = hin(2,j) + qeta(i,j)
4355      continue
          do 4358 j = 1,ldm
          hin(1,j) = hin(1,j) / imjm
          hin(2,j) = hin(2,j) / imjm
c         write(6,*)"qsave bar,qeta bar,level",hin(1,j),hin(2,j),j
4358      continue
C
C     NOW REDEFINE PD TO HAVE IT EQUAL TO PSFC-PT.
C
  603 DO 605 K = 1, IMJM
         IF(FIX(K) .NE. 1.0) THEN
          PD(K) = REF(K) * PD(K)
         ELSE IF(FIX(K) .EQ. 1.0) THEN
          PD(K) = PSFC(K) - PT
         END IF
         PERIC(K) = PD(K) + PT
  605 CONTINUE
C
C     PRINT 9
C     DO 917 L = 1, LM
C       PRINT 3, L
C       PRINT 510
C 510   FORMAT ('  TEMPERATURE')
C       KSTO = 1
C       KSTE = JM + 1
C       KENDO = (JM * 20) + 1
C       KENDE = KENDO + IM
C       DO 918 I = 1, IM
C        IF(MOD(I,2).EQ.1) THEN
C           WRITE(6,1)  (TETA(K,L),K = KSTO, KENDO, JM)
C           KSTO = KSTO + 1
C           KENDO = KENDO + 1
C         ELSE
C           WRITE(6,1)  (TETA(K,L),K = KSTE, KENDE, JM)
C           KSTE = KSTE + 1
C           KENDE = KENDE + 1
C         END IF
C 918   CONTINUE
C 917 CONTINUE
C
      DO 65432 L=1,LM
      DO 65432 K=1,IMJM
      IF(TETA(K,L) .GT. 350. OR.TETA(K,L).LT.150.E0)THEN
c       WRITE(6,87722)K,L,TETA(K,L)
87722   FORMAT(' PTETA K=',I4,' L=',I2,' T=',E12.5)
      ENDIF
      IF(QETA(K,L).GT.1.E0)THEN
c       WRITE(6,87723)K,L,QETA(K,L)
87723   FORMAT(' PTETA K=',I4,' L=',I2,' Q=',E12.5)
      ENDIF
      IF(UETA(K,L).GT.150.E0.OR.UETA(K,L).LT.-150.E0)THEN
c       WRITE(6,87724)K,L,UETA(K,L)
87724   FORMAT(' PTETA K=',I4,' L=',I2,' U=',E12.5)
      ENDIF
      IF(VETA(K,L).GT.150.E0.OR.VETA(K,L).LT.-150.E0)THEN
c       WRITE(6,87725)K,L,VETA(K,L)
87725   FORMAT(' PTETA K=',I4,' L=',I2,' V=',E12.5)
      ENDIF
65432 CONTINUE
C       WRITE(6,7725)
C7725   FORMAT(' BEFORE TOPOG = 0')
C     DO 65433 L=1,LM
      DO 65433 K=1,8580
      IF(PD(K) .LT. 8000.) THEN
c       WRITE(6,87729) K,PD(K)
87729   FORMAT(2X,'K = ',I5,' PD(K) = ',E12.5)
      END IF
65433 CONTINUE
C
C------- SETTING U,V AND T EQUAL TO ZERO AT POINTS BELOW THE GROUND ---
C------- SURFACE ------------------------------------------------------
C
      DO 905 K = 1, IMJM
      DO 905 L = 1 ,LM
        if(eta(l+1)-ref(k).gt.1.e-6) then
c       IF(ETA(L+1) .GT. REF(K)) THEN
           KKSW = K + KSW
           KKNW = K + KNW
c          IF(K .Le. -9) THEN
c           WRITE(6,90550) K,L,ETA(L+1),REF(K),KKNW,KKSW
90550       FORMAT(2X,'K= ',I4,' L= ',I2,' ETA(L+1)= ',E11.4,
     1       ' REF(K) = ',E11.4,' KKNW,KKSW= ',2I4)
c          END IF
           TETA(K,L) = 0.
           QETA(K,L) = 0.
           UETA(K,L) = 0.
           VETA(K,L) = 0.
           UETA(K-1,L) = 0.
           VETA(K-1,L) = 0.
           UETA(K+KNW,L) = 0.
           VETA(K+KNW,L) = 0.
           UETA(K+KSW,L) = 0.
           VETA(K+KSW,L) = 0.
        END IF
 905  CONTINUE
C       WRITE(6,6725)
C6725   FORMAT(' AFTER TOPOG = 0')
C     DO 35433 L=1,LM
C     DO 35433 K=1,8580
C     IF(K .GT. 2950 .AND. K .LT. 3000 .AND.  L .GE. 10) THEN
C       WRITE(6,37722) K,L,TETA(K,L),QETA(K,L),UETA(K,L),VETA(K,L)
C7722   FORMAT(2X,2(I5,1X),4(E12.5,1X))
C     END IF
C5433 CONTINUE
C
      DO 90500 L = 1, LM
C        WRITE(6,49940) L, ETA(L+1), REF(7448)
C        WRITE(6,49941) L, TETA(1,L), QETA(1,L)
C        WRITE(6,49942) L, TETA(7448,L), QETA(7448,L)
C        WRITE(6,49943) L, TETA(8581,L), QETA(8581,L)
C9940    FORMAT(' PTETA  L=',I2,' ETA(+1)=',E12.5,' REF=',E12.5)
C9941    FORMAT(' L=',I2,' TETA(1)=',E12.5,' QETA(1)=',E12.5)
C9942    FORMAT(' L=',I2,' TETA(7748)=',E12.5,' QETA(7748)=',E12.5)
C9943    FORMAT(' L=',I2,' TETA(8581)=',E12.5,' QETA(8581)=',E12.5)
c        WRITE(6,98849) L
98849    FORMAT(2X,' OUT OF PTETA  L = ',I2)
C        DO 90501 K = 1, IMJM
C          IF(MOD(K,100) .EQ. 0 .OR. K .EQ. 1) THEN
C          IF(L .EQ. 8) THEN
C            WRITE(6,90502) K,L,TETA(K,L),QETA(K,L),UETA(K,L),
C    1                      VETA(K,L),HETA(K,L),HGT(K)
C0502        FORMAT(2X,2(I5,1X),6(E12.5,1X))
C          END IF
C0501    CONTINUE
90500 CONTINUE
         DO 90501 K = 1, IMJM
        IF(K .GE. 5214 .AND. K .LE. 5217 .OR. K .GE. 5279 .AND.
     1      K .LE. 5283) THEN
C          IF(MOD(K,100) .EQ. 0 .OR. K .EQ. 1 .OR. K .GT. 8500) THEN
c           WRITE(6,98610) K,PERIC(K),REF(K),HGT(K),SM(K),FIX(K)
98610       FORMAT(2X,'K(OUTPUT)= ',I5,5(1X,E12.5))
            DO 90504 L = 1, LM
c            WRITE(6,90502) L,TETA(K,L),QETA(K,L),UETA(K,L),
c    1                      VETA(K,L),HETA(K,L),ETA(L)
90502        FORMAT(2X,I5,1X,6(E12.5,1X))
90504       CONTINUE
           END IF
90501    CONTINUE
C
C     CALL CMAS3D(TETA,'TETA',8,3427,16,.TRUE.)
C     CALL CMAS3D(UETA,'UETA',8,3427,16,.TRUE.)
C     CALL CMAS3D(VETA,'VETA',8,3427,16,.TRUE.)
C
C
C        SETTING U,V AND T EQUAL TO ZERO AT POINTS ALONG THE
C        NORTHERNMOST (OUTSIDE) LINE, IN ORDER TO ENABLE VECTORIZATION
C        OF SOME OF THE FORECAST CODE
C
C        NOT NEEDED IN CRAY VERSION
C
C     DO 920 L = 1, LM
C       DO 919 K = KSL, IMJM
C         UETA(K,L) = 0.
C         VETA(K,L) = 0.
C         TETA(K,L) = 0.
C         QETA(K,L) = 0.
C919     CONTINUE
C920  CONTINUE
C
C------------- THE SEPARATION OF THE BOUNDARY VALUES -------------------
C
      KHSTWB = 2 * KNE
      KHENWB = KSL + KSW + KSE
      KHICWB = KNW + KNE
      KHSTEB = (2 * KNE) + KNW
      KHENEB = IMJM + KSW + KSE
      KHICEB = KNE + KNW
c     WRITE(6, 84589) KHSTWB, KHENWB, KHICWB
c     WRITE(6, 84589) KHSTEB, KHENEB, KHICEB
84589 FORMAT(2X, 3I7)
      KVSTWB = IM
      KVENWB = KSL + KSW
      KVICWB = KNE + KNW
      KVSTEB = KNE + KNW
      KVENEB = KSLM1
      KVICEB = KNE + KNW
c     WRITE(6, 84589) KVSTWB, KVENWB, KVICWB
c     WRITE(6, 84589) KVSTEB, KVENEB, KVICEB
C


























































C ---------------------------------------------------------------------


C       SOUTH BOUNDARY - H POINTS
C
      N = 1
      DO 541 K = 1, IM
         PDB(N,1) = PD(K)
         PDB(N,2) = 0.
         DO 540 L = 1, LM
            TB(N,L,1) = TETA(K,L)
            TB(N,L,2) = 0.
            QB(N,L,1) = QETA(K,L)
            QB(N,L,2) = 0.
 540     CONTINUE
         N = N + 1
 541  CONTINUE
C
C       NORTH BOUNDARY - H POINTS
C
      DO 543 K = KSL, IMJM
         PDB(N,1) = PD(K)
         PDB(N,2) = 0.
         DO 542 L = 1, LM
            TB(N,L,1) = TETA(K,L)
            TB(N,L,2) = 0.
            QB(N,L,1) = QETA(K,L)
            QB(N,L,2) = 0.
 542     CONTINUE
         N = N + 1
 543  CONTINUE
C
C       WEST BOUNDARY - H POINTS
C
      DO 545 K = KHSTWB, KHENWB, KHICWB
         PDB(N,1) = PD(K)
         PDB(N,2) = 0.
         DO 544 L = 1, LM
            TB(N,L,1) = TETA(K,L)
            TB(N,L,2) = 0.
            QB(N,L,1) = QETA(K,L)
            QB(N,L,2) = 0.
 544     CONTINUE
      N = N + 1
 545  CONTINUE
C
C       EAST BOUNDARY - H POINTS
C
      DO 547 K = KHSTEB, KHENEB, KHICEB
         PDB(N,1) = PD(K)
         PDB(N,2) = 0.
         DO 546 L = 1, LM
            TB(N,L,1) = TETA(K,L)
            TB(N,L,2) = 0.
            QB(N,L,1) = QETA(K,L)
            QB(N,L,2) = 0.
 546     CONTINUE
          N = N + 1
 547  CONTINUE
C
C       SOUTH BOUNDARY - V POINTS
C
      N = 1
      DO 549 K = 1, IMM1
         DO 548 L = 1, LM
            UB(N,L,1) = UETA(K,L)
            UB(N,L,2) = 0.
            VB(N,L,1) = VETA(K,L)
            VB(N,L,2) = 0.
 548     CONTINUE
         N = N + 1
 549  CONTINUE
C
C       NORTH BOUNDARY - V POINTS
C
      DO 551 K = KSL, IMJMM1
         DO 550 L = 1, LM
            UB(N,L,1) = UETA(K,L)
            UB(N,L,2) = 0.
            VB(N,L,1) = VETA(K,L)
            VB(N,L,2) = 0.
 550     CONTINUE
         N = N + 1
 551  CONTINUE
C
C       WEST BOUNDARY - V POINTS
C
      DO 553 K = KVSTWB, KVENWB, KVICWB
         DO 552 L = 1, LM
            UB(N,L,1) = UETA(K,L)
            UB(N,L,2) = 0.
            VB(N,L,1) = VETA(K,L)
            VB(N,L,2) = 0.
 552     CONTINUE
         N = N + 1
 553  CONTINUE
C
C       EAST BOUNDARY - V POINTS
C
      DO 555 K = KVSTEB, KVENEB, KVICEB
         DO 554 L = 1, LM
            UB(N,L,1) = UETA(K,L)
            UB(N,L,2) = 0.
            VB(N,L,1) = VETA(K,L)
            VB(N,L,2) = 0.
 554     CONTINUE
         N = N + 1
 555  CONTINUE
C
      REWIND LOUT
      RUN = .TRUE.
      NTSD = 0
C704     WRITE (LOUT) TETA,PD,PHIS,REF,ETA,PT,DETA,ETAL,ZETA,PDB,TB,UB, VB
         WRITE (LOUT) RUN,IDAT,IHRST,NTSD,UETA,VETA,WET,SST,SI
         WRITE (LOUT) TETA,QETA,PD,PHIS,SM,REF,ETA,PT,DETA,ETAL,ZETA
     1                ,PDB,TB,QB,UB,VB
 704  CONTINUE
C
      TEST = .FALSE.
      IF (TEST) GO TO 529
  630 CONTINUE
         tsadif = 0.0
         kprdif = 0
         do 4349 i = 2705,2967,262
         do 4349 j = 1,ldm
         tmpdif = teta(i,j) - tsave(i,j)
         tsadif = tsadif + tmpdif
c           write(6,*)"t,tsave,tmpdif,i,j ",teta(i,j),tsave(i,j),
c    $                   tmpdif,i,j
         kprdif = kprdif + 1
4349     continue
         do 4348 i = 2705,2967,262
         pdiffff = psave(i) - pd(i)
c        write(6,*)"psave,pd,pdif ",psave(i),pd(i),pdiffff
         do 4348 j = 1,ldm
         zmpdif = hsave(i,j) - heta(i,j)
         if(abs(zmpdif).gt.1.0) then
c           write(6,*)"h,hsave,zmpdif,i,j ",heta(i,j),hsave(i,j),
c    $                   zmpdif,i,j
         end if
4348     continue
            tsadif = tsadif/(imjm*ldm)
c           write(6,*)"tsadif,kprdif  ",tsadif,kprdif
C
         do 4351 i = 2705,2967,262
         do 4351 j = 1,ldm
         qmpdif = qeta(i,j) - q(i,j)
c           write(6,*)"q,qeta,qmpdif,i,j ",q(i,j),qeta(i,j),
c    $                   qmpdif,i,j
4351     continue
      return
      END
