              PROGRAM EBU
C$$$  MAIN PROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C MAIN PROGRAM: ETAFCST      EARLY ETA MODEL FORECAST DRIVER
C   PRGMMR: ROGER            ORG: NP22        DATE: 97-08-13  
C
C ABSTRACT: EBU3 CONTAINS THE PRIMARY RUNSTREAM FOR THE EARLY ETA
C   FORECAST MODEL.  AFTER AN INITIAL CALL TO SUBROUTINE INIT, CALLS
C   ARE MADE TO SUBROUTINES WHICH COMPUTE THE VARIOUS DYNAMICAL AND
C   PHYSICAL PROCESSES IN THE MODEL.  THE VARIABLE 'NTSD' IS THE
C   FUNDAMENTAL TIMESTEP COUNTER AND THUS ITS VALUE DETERMINES WHEN
C   THE SUBROUTINES ARE CALLED.  INFORMATION PERTAINING TO THE SCHEMES
C   USED IN THE MODEL AS WELL AS ADDITIONAL REFERENCES MAY BE FOUND
C   IN "THE STEP-MOUNTAIN ETA COORDINATE REGIONAL MODEL:  A DOCUMEN-
C   TATION" (BLACK 1988; DEVELOPMENT DIVISION) AND "THE NEW NMC MESO-
C   SCALE ETA MODEL:  DESCRIPTION AND FORECAST EXAMPLES (BLACK 1994;
C   WEATHER AND FORECASTING).
C
C PROGRAM HISTORY LOG:
C   87-08-??  JANJIC,     ORIGINATOR
C             BLACK
C   93-05-12  TREADON     DOCBLOCK INSERTED
C   93-10-25  BLACK       DOCBLOCK UPDATED
C
C
C USAGE:  MAIN PROGRAM
C
C   INPUT FILES:  NONE (SEE SUBROUTINE INIT)
C
C   OUTPUT FILES:  NONE (SEE SUBROUTINE CHKOUT)
C
C   SUBPROGRAMS CALLED:
C     UNIQUE:
C       INIT     - INITIALIZE VARIABLES AT MODEL START
C       PFDHT    - UPDATE PRESSURE GRADIENT FORCE, COROILIS,
C                  DIVERGENCE, AND HORIZONTAL PART OF THE
C                  OMEGA-ALPHA EQUATION
C       PDTE     - UPDATE SURFACE PRESSURE TENDENCY AND ETA DOT
C       CHKOMG   - SET VERTICAL VELOCITY FIELDS FOR POST PROCESSOR.
C                  FOR THE EXTERNAL POST THIS ROUTINE GENERATES A
C                  TEMPORARY FILE CONTAINING DATA NEEDED TO COMPUTE
C                  VERTICAL VELOCITY.  FOR THE INTERNAL POST THE
C                  ROUTINE LOADS DATA INTO COMMON BLOCKS PASSED TO
C                  THE POST PROCESSOR.
C       CHKPTO   - WRITE HISTORY DATA FOR OMEGA
C       VTADV    - VERTICAL ADVECTION
C       HZADV    - HORIZONTAL ADVECTION
C       QSMTH    - REMOVE NEGATIVE SPECIFIC HUMIDITY
C       PDNEW    - UPDATE SURFACE PRESSURE
C       RADTN    - RADIATION DRIVER
C       RDTEMP   - APPLY TEMPERATURE TENDENCY DUE TO RADIATION
C       SURFACE  - UPDATE SURFACE TEMPERATURE, MOISTURE, AND OTHER
C                - GROUND HYDROLOGY
C       TURBL    - PERFORM THE VERTICAL TURBULENT EXCHANGE
C       PBLOUT   - OUTPUT SURFACE LAYER DIAGNOSTICS
C       GSCOND   - GRID SCALE PRECIPITATION
C       CUCNVC   - CONVECTIVE ADJUSTMENT FOR DEEP OR SHALLOW CONVECTION
C       HDIFF    - LATERAL DIFFUSION
C       BOCOH    - UPDATE H POINTS ON THE BOUNDARIES
C       CHKOUT   - POST PROFILE DATA.  FOR INTERNAL POST,
C                  POSTS MODEL OUTPUT.  FOR EXTERNAL POST,
C                  WRITES TEMPORARY FILE CONTAINING ALL MODEL
C                  ARRAYS.
C       DDAMP    - APPLY DIVERGENCE DAMPING
C       BOCOV    - UPDATE V POINTS ON THE BOUNDARIES
C
C   EXIT STATES:
C     COND =   1 - NORMAL EXIT
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN
C   MACHINE : CRAY Y-MP 8/32
C
C$$$
C     ******************************************************************
C     *                                                                *
C     *                       ETA COORDINATE                           *
C     *                      LIMITED AREA MODEL                        *
C     *             WITH STEP-MOUNTAIN (ETA) COORDINATE                *
C     *                                                                *
C     *                                                                *
C     *    NOAA / NATIONAL METEOROLOGICAL CENTER, CAMP SPRINGS, MD,    *
C     *                                                                *
C     *  GEOPHYSICAL FLUID DYNAMICS LABORATORY / NOAA, PRINCETON, NJ,  *
C     *                                                                *
C     *  UNIVERSITY CORPORATION FOR ATMOSPHERIC RESEARCH, BOULDER, CO, *
C     *                              &                                 *
C     * DEPARTMENT OF METEOROLOGY, UNIVERSITY OF BELGRADE, YUGOSLAVIA  *
C     *                                                                *
C     ******************************************************************
C
C     ******************************************************************
C     *                                                                *
C     *                        REFERENCES                              *
C     *             FOR THE DYNAMICAL PART OF THE MODEL                *
C     *                                                                *
C     *  STEP-MOUNTAIN ETA COORDINATE:                                 *
C     *      F. MESINGER, 1983, IN RES. ACTIVITIES IN ATMOS. AND       *
C     *      OCEANIC MODELLING, REP. NO. 5, WMO, GENEVA, 4.9-4.10.     *
C     *                                                                *
C     *  HORIZONTAL ADVECTION, CONTINUITY EQUATION:                    *
C     *      Z.I. JANJIC, 1984, MWR, 112, NO.6, 1234-1245.             *
C     *                                                                *
C     *  INTERNAL BOUNDARIES, OMEGA-ALPHA TERM, CODING, PERFORMANCE:   *
C     *      MESINGER ET AL., 1988, MWR, 116 NO.7, 1493-1518.          *
C     *                                                                *
C     *  N.B. FOR MORE DETAILS ON THESE TOPICS SEE ALSO:               *
C     *                                                                *
C     *  1.  MESINGER, F., AND Z.I. JANJIC, 1985: PROBLEMS AND         *
C     *        NUMERICAL METHODS OF THE INCORPORATION OF MOUNTAINS IN  *
C     *        ATMOSPHERIC MODELS.  LECTURES IN APPLIED MATHEMATICS,   *
C     *        VOL 22, AMER. MATH. SOC.; ALSO, NUMERICAL METHODS FOR   *
C     *        WEATHER PREDICTION, SEMINAR 1983, ECMWF, 103-157;       *
C     *        ALSO, SHORT- AND MEDIUM-RANGE WEATHER PREDICTION        *
C     *        RESEARCH PUBL. SER., NO. 8, WMO, GENEVA, 175-233.       *
C     *                                                                *
C     *  2.  JANJIC, Z.I., AND F. MESINGER, 1983: FINITE-DIFFERENCE    *
C     *        METHODS FOR THE SHALLOW WATER EQUATIONS ON VARIOUS      *
C     *        HORIZONTAL GRIDS.  NUMERICAL METHODS FOR WEATHER        *
C     *        PREDICTION, SEMINAR 1983, ECMWF,29-101.                 *
C     *                                                                *
C     *                        REFERENCES                              *
C     *             FOR THE DYNAMICAL PART OF THE MODEL                *
C     *                                                                *
C     *  JANJIC, Z.I., 1990: THE STEP-MOUNTAIN COORDINATE:             *
C     *     PHYSICAL PACKAGE.  MONTHLY WEATHER REVIEW, VOL. 118,       *
C     *     NO. 7, 1429-1443.                                          *
C     *                                                                *
C     *  ALSO SEE REFERENCES IN PHYSICAL SUBROUTINES                   *
C     *                                                                *
C     ******************************************************************
C
C     ******************************************************************
C     *                                                                *
C     *  THIS VERSION OF THE PROGRAM IS WRITTEN IN STANDARD ANSI       *
C     *  FORTRAN 77                                                    *
C     *                                                                *
C     *  PRINCIPAL PROGRAMMERS:                                        *
C     *                                                                *
C     *  Z. JANJIC, UNIVERSITY OF BELGRADE,                            *
C     *  T. BLACK, NMC.                                                *
C     *                                                                *
C     ******************************************************************
C     *                                                                *
C     *  THE MODEL USES THE SEMI-STAGGERED E GRID IN ARAKAWA NOTATION. *
C     *  HORIZONTAL INDEXING IS ONE-DIMENSIONAL. THE RELATION BETWEEN  *
C     *  ONE-DIMENSIONAL AND TWO-DIMENSIONAL INDICES IS SHOWN IN THE   *
C     *  FIGURE BELOW:                                                 *
C     *                                                                *
C     *          ---     ---     ---     ---     ---     ---           *
C     *      JM / H   V / H   V / H   V / H   V / H   V / H  (V)       *
C     *         -----   -----   -----   -----   -----   -----          *
C     *           V*/ H   V / H   V / H   V / H   V / H   V /          *
C     *          ---     ---     ---     ---     ---     ---           *
C     *         / H   V / H   V / H   V / H   V / H   V / H  (V)       *
C     *      .  -----   -----   -----   -----   -----   -----          *
C     *      .    V*/ H   V / H   V / H   V / H   V / H   V /          *
C     *      .   ---     ---     ---     ---     ---     ---           *
C     *         / H   V / H   V / H   V / H   V / H   V / H  (V)       *
C     *         -----   -----   -----   -----   -----   -----          *
C     *      2    V*/ H   V / H   V / H   V / H   V / H   V /          *
C     *          ---     ---     ---     ---     ---     ---           *
C     *     J=1 / H   V / H   V / H   V / H   V / H   V / H  (V)       *
C     *         -----   -----   -----   -----   -----   -----          *
C     *                                                                *
C     *            I=1      2        ...                  IM           *
C     *                                                                *
C     *  THE ONE-DIMENSIONAL INDEX K INCREASES GOING FROM LEFT TO      *
C     *  RIGHT STARTING FROM THE LEFT LOWER CORNER. AFTER REACHING THE *
C     *  RIGHT BOUNDARY OF A HORIZONTAL ROW, IT JUMPS TO THE LEFT      *
C     *  BOUNDARY OF THE NEXT HORIZONTAL ROW ABOVE.                    *
C     *                                                                *
C     *  WIND POINTS V* ALONG THE LEFT BOUNDARY HAVE J VALUES ONE LESS *
C     *  THAN THE LINE THEY ARE IN; THAT IS, THEY ARE STORED IN PLACES *
C     *  DENOTED BY (V).                                               *
C     *                                                                *
C     *  THE TOTAL NUMBER OF GRID POINTS IN THE HORIZONTAL IS          *
C     *       IMJM=IM*JM-JM/2                                          *
C     *                                                                *
C     *  THE INDEX INCREMENTS NEEDED TO REACH N, S, NE, NW, SW, AND SE *
C     *  GRID POINT NEIGHBORS ARE FOR HEIGHT POINTS RESPECTIVELY       *
C     *       KHN=IM-1; KHS=-IM;                                       *
C     *       KNE=IM; KNW=IM-1; KSW=-IM; KSE=-IM+1;                    *
C     *  AND FOR VELOCITY POINTS RESPECTIVELY                          *
C     *       KVN=IM; KVS=-IM+1;                                       *
C     *       KNE=IM; KNW=IM-1; KSW=-IM; KSE=-IM+1.                    *
C     *                                                                *
C     *  THE HORIZONTAL LOOP BOUNDARIES ARE TREATED DIFFERENTLY IN     *
C     *  SCALAR AND VECTOR PARTS OF THE CODE. IN SCALAR PARTS OF THE   *
C     *  CODE THE LOWER AND THE UPPER BOUNDARIES FOR INDEX K ARE       *
C     *  DEFINED AND STORED FOR EACH HORIZONTAL ROW (SEE SUBROUTINE    *
C     *  'ECON' IN THE PROGRAME ECONVK FOR DETAILS). IN VECTOR PARTS   *
C     *  OF THE CODE THE UPPER AND THE LOWER BOUNDARIES ARE DEFINED    *
C     *  FOR THE WHOLE HORIZONTAL FIELD. THESE BOUNDARIES ARE DEFINED  *
C     *  IN 'PARAMETER' STATEMENTS IN EACH SUBROUTINE.                 *
C     *                                                                *
C     *  THE SECOND LETTERS IN LOOP BOUNDARY NAMES INDICATE WHETHER    *
C     *  THEY DEFINE THE BOUNDARIES FOR HEIGHT (H) OR VELOCITY (V)     *
C     *  POINTS.  THE THIRD LETTER IS L FOR LOW OR H FOR HIGH.         *
C     *  THE SUFFIX DIGIT(S) INDICATE(S) FOR HOW MANY LINES THE LOOP   *
C     *  BOUNDARIES ARE SHIFTED FROM THE BOUNDARIES OF THE INTEGRATION *
C     *  DOMAIN.                                                       *
C     *                                                                *
C     *  BOUNDARY MASKS AND TOPOGRAPHY MASKS ARE DEFINED FOR VECTOR    *
C     *  PROCESSING. THE BOUNDARY MASKS HBM2(K) AND VBM2(K) ARE        *
C     *  EQUAL TO ONE EVERYWHERE EXCEPT AT THE TWO OUTERMOST ROWS      *
C     *  WHERE THEY ARE EQUAL TO ZERO. THE BOUNDARY MASK VBM3(K) IS    *
C     *  EQUAL TO ONE EVERYWHERE EXCEPT AT THE THREE OUTERMOST ROWS    *
C     *  WHERE IT IS EQUAL TO ZERO. THE TOPOGRAPHY MASKS (HTM(K,L)     *
C     *  AND VTM(K,L)) ARE SET TO ZERO UNDERNEATH THE TOPOGRAPHY AND   *
C     *  TO ONE ELSWHERE. IN ADDITION, FOR TREATMENT OF PHYSICAL       *
C     *  PROCESSES, MAXIMUM VALUES OF THE VERTICAL INDEX ARE DEFINED   *
C     *  AND STORED (LMH(K) AND LMV(K), SEE SUBROUTINE ECON IN THE     *
C     *  PROGRAM ECONVK FOR DETAILS).                                  *
C     *                                                                *
C     ******************************************************************
                             L O G I C A L
     & RUN,FIRST,RESTRT,SIGMA
C----------------------------------------------------------------------
      INCLUDE "CTLBLK"
C
C
C      CALL W3LOG('$S97225.65','ETAFCST ')
Chou      CALL W3TAGB('ETAFCST ',0097,0225,0065,'NP22   ')                   






C
C
C--------INITIALIZE CONSTANTS AND VARIABLES-----------------------------
                    C A L L    I N I T
C--------JUMP INTO THE TIME LOOP FOR PRELIMINARY CALCULATIONS----------
                   G O   T O  1200
C********ENTRY INTO THE TIME LOOP***************************************
 2000 CONTINUE
      NTSD=NTSD+1
      WRITE(6,2001) NTSD,(NTSD-1)*DT
 2001 FORMAT('EBU:  TIMESTEP NTSD=',I5,'  FCST TIME=',F8.1)
C***********************************************************************
C--------PRESS. TEND.,ETA DOT & VERTICAL OMEGA-ALPHA--------------------
                   C A L L    P D T E
C--------IS IT TIME FOR A CHECK POINT ON THE MODEL HISTORY FILE?--------
                   C A L L    C H K O M G
C--------IS IT TIME FOR ADVECTION?--------------------------------------
      IF(MOD(NTSD-1,IDTAD).EQ.0)THEN
                   C A L L    V T A D V
                   C A L L    H Z A D V
                   C A L L    H Z A D V 2
      ENDIF
C--------UPDATING PRESSURE DIFFERENCE-----------------------------------
                   C A L L    P D N E W
C--------RADIATION------------------------------------------------------
      IF(MOD(NTSD-1,NRADS).EQ.0.OR.MOD(NTSD-1,NRADL).EQ.0)THEN
            C A L L    R A D T N
      ENDIF
C--------APPLY TEMPERATURE TENDENCY DUE TO RADIATION--------------------
                   C A L L    R D T E M P
C--------TURBULENT PROCESSES AND PRECIPITATION--------------------------
      IF(MOD(NTSD-1,NPHS).EQ.0)        T H E N
                   C A L L    T U R B L
      ENDIF
      IF(MOD(NTSD-1,NPHS).EQ.0)       T H E N
                   C A L L    G S C O N D
      ENDIF
      IF(MOD(NTSD-1,NCNVC).EQ.0)        T H E N
                   C A L L    C U C N V C
      ENDIF
      IF(MOD(NTSD-1,NPHS).EQ.0)       T H E N
                   C A L L    P R E C P D
      ENDIF
C--------LATERAL DIFFUSION---------------------------------------------
                   C A L L    H D I F F
C--------UPDATING BOUNDARY VALUES AT HEIGHT POINTS----------------------
                   C A L L    B O C O H
C--------IS IT TIME FOR A CHECK POINT ON THE MODEL HISTORY FILE?-------
 1200 CONTINUE
                   C A L L    C H K O U T 
C--------PGF,CORIOLIS,DIVERGENCE & HORIZONTAL PART OF OMEGA-ALPHA-------
                   C A L L    P F D H T
C--------CLEAN UP AFTER 1ST STEP PRELIMINARY CALCULATIONS---------------
      IF(FIRST.AND..NOT.RESTRT)        T H E N
      FIRST=.FALSE.
                   GO TO 2000
      ENDIF
C--------DIVERGENCE DAMPING---------------------------------------------
                   C A L L    D D A M P
C--------UPDATING BOUNDARY VALUES AT VELOCITY POINTS--------------------
                   C A L L    B O C O V
C--------CLEAN UP AFTER RESTART-----------------------------------------
      IF(RESTRT)                       T H E N
                   FIRST=.FALSE.
                   RESTRT=.FALSE.
      ENDIF
C***********************************************************************
      IF(NTSD.LT.NTSTM)                GO TO 2000
C********EXIT FROM THE TIME LOOP****************************************
C
C      CALL W3LOG('$E')
Chou      CALL W3TAGE('ETAFCST ') 
C
      STOP1
      END
