PROGRAM ebu
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !$$$  main program documentation block
  !                .      .    .                                       .
  ! main program: etafcst      early eta model forecast driver
  !   prgmmr: roger            org: np22        date: 97-08-13  
  !
  ! abstract: ebu3 contains the primary runstream for the early eta
  !   forecast model.  after an initial call to subroutine init, calls
  !   are made to subroutines which compute the various dynamical and
  !   physical processes in the model.  the variable 'ntsd' is the
  !   fundamental timestep counter and thus its value determines when
  !   the subroutines are called.  information pertaining to the schemes
  !   used in the model as well as additional references may be found
  !   in "the step-mountain eta coordinate regional model:  a documen-
  !   tation" (black 1988; development division) and "the new nmc meso-
  !   scale eta model:  description and forecast examples (black 1994;
  !   weather and forecasting).
  !
  ! program history log:
  !   87-08-??  janjic,     originator
  !             black
  !   93-05-12  treadon     docblock inserted
  !   93-10-25  black       docblock updated
  !
  !
  ! usage:  main program
  !
  !   input files:  none (see subroutine init)
  !
  !   output files:  none (see subroutine chkout)
  !
  !   subprograms called:
  !     unique:
  !       init     - initialize variables at model start
  !       pfdht    - update pressure gradient force, coroilis,
  !                  divergence, and horizontal part of the
  !                  omega-alpha equation
  !       pdte     - update surface pressure tendency and eta dot
  !       chkomg   - set vertical velocity fields for post processor.
  !                  for the external post this routine generates a
  !                  temporary file containing data needed to compute
  !                  vertical velocity.  for the internal post the
  !                  routine loads data into common blocks passed to
  !                  the post processor.
  !       chkpto   - write history data for omega
  !       vtadv    - vertical advection
  !       hzadv    - horizontal advection
  !       qsmth    - remove negative specific humidity
  !       pdnew    - update surface pressure
  !       radtn    - radiation driver
  !       rdtemp   - apply temperature tendency due to radiation
  !       surface  - update surface temperature, moisture, and other
  !                - ground hydrology
  !       turbl    - perform the vertical turbulent exchange
  !       pblout   - output surface layer diagnostics
  !       gscond   - grid scale precipitation
  !       cucnvc   - convective adjustment for deep or shallow convection
  !       hdiff    - lateral diffusion
  !       bocoh    - update h points on the boundaries
  !       chkout   - post profile data.  for internal post,
  !                  posts model output.  for external post,
  !                  writes temporary file containing all model
  !                  arrays.
  !       ddamp    - apply divergence damping
  !       bocov    - update v points on the boundaries
  !
  !   exit states:
  !     cond =   1 - normal exit
  !
  !     ******************************************************************
  !     *                                                                *
  !     *                       eta coordinate                           *
  !     *                      limited area model                        *
  !     *             with step-mountain (eta) coordinate                *
  !     *                                                                *
  !     *                                                                *
  !     *    noaa / national meteorological center, camp springs, md,    *
  !     *                                                                *
  !     *  geophysical fluid dynamics laboratory / noaa, princeton, nj,  *
  !     *                                                                *
  !     *  university corporation for atmospheric research, boulder, co, *
  !     *                              &                                 *
  !     * department of meteorology, university of belgrade, yugoslavia  *
  !     *                                                                *
  !     ******************************************************************
  !
  !     ******************************************************************
  !     *                                                                *
  !     *                        references                              *
  !     *             for the dynamical part of the model                *
  !     *                                                                *
  !     *  step-mountain eta coordinate:                                 *
  !     *      f. mesinger, 1983, in res. activities in atmos. and       *
  !     *      oceanic modelling, rep. no. 5, wmo, geneva, 4.9-4.10.     *
  !     *                                                                *
  !     *  horizontal advection, continuity equation:                    *
  !     *      z.i. janjic, 1984, mwr, 112, no.6, 1234-1245.             *
  !     *                                                                *
  !     *  internal boundaries, omega-alpha term, coding, performance:   *
  !     *      mesinger et al., 1988, mwr, 116 no.7, 1493-1518.          *
  !     *                                                                *
  !     *  n.b. for more details on these topics see also:               *
  !     *                                                                *
  !     *  1.  mesinger, f., and z.i. janjic, 1985: problems and         *
  !     *        numerical methods of the incorporation of mountains in  *
  !     *        atmospheric models.  lectures in applied mathematics,   *
  !     *        vol 22, amer. math. soc.; also, numerical methods for   *
  !     *        weather prediction, seminar 1983, ecmwf, 103-157;       *
  !     *        also, short- and medium-range weather prediction        *
  !     *        research publ. ser., no. 8, wmo, geneva, 175-233.       *
  !     *                                                                *
  !     *  2.  janjic, z.i., and f. mesinger, 1983: finite-difference    *
  !     *        methods for the shallow water equations on various      *
  !     *        horizontal grids.  numerical methods for weather        *
  !     *        prediction, seminar 1983, ecmwf,29-101.                 *
  !     *                                                                *
  !     *                        references                              *
  !     *             for the dynamical part of the model                *
  !     *                                                                *
  !     *  janjic, z.i., 1990: the step-mountain coordinate:             *
  !     *     physical package.  monthly weather review, vol. 118,       *
  !     *     no. 7, 1429-1443.                                          *
  !     *                                                                *
  !     *  also see references in physical subroutines                   *
  !     *                                                                *
  !     ******************************************************************
  !
  !     ******************************************************************
  !     *                                                                *
  !     *  this version of the program is written in standard ansi       *
  !     *  fortran 77                                                    *
  !     *                                                                *
  !     *  principal programmers:                                        *
  !     *                                                                *
  !     *  z. janjic, university of belgrade,                            *
  !     *  t. black, nmc.                                                *
  !     *                                                                *
  !     ******************************************************************
  !     *                                                                *
  !     *  the model uses the semi-staggered e grid in arakawa notation. *
  !     *  horizontal indexing is one-dimensional. the relation between  *
  !     *  one-dimensional and two-dimensional indices is shown in the   *
  !     *  figure below:                                                 *
  !     *                                                                *
  !     *          ---     ---     ---     ---     ---     ---           *
  !     *      jm / h   v / h   v / h   v / h   v / h   v / h  (v)       *
  !     *         -----   -----   -----   -----   -----   -----          *
  !     *           v*/ h   v / h   v / h   v / h   v / h   v /          *
  !     *          ---     ---     ---     ---     ---     ---           *
  !     *         / h   v / h   v / h   v / h   v / h   v / h  (v)       *
  !     *      .  -----   -----   -----   -----   -----   -----          *
  !     *      .    v*/ h   v / h   v / h   v / h   v / h   v /          *
  !     *      .   ---     ---     ---     ---     ---     ---           *
  !     *         / h   v / h   v / h   v / h   v / h   v / h  (v)       *
  !     *         -----   -----   -----   -----   -----   -----          *
  !     *      2    v*/ h   v / h   v / h   v / h   v / h   v /          *
  !     *          ---     ---     ---     ---     ---     ---           *
  !     *     j=1 / h   v / h   v / h   v / h   v / h   v / h  (v)       *
  !     *         -----   -----   -----   -----   -----   -----          *
  !     *                                                                *
  !     *            i=1      2        ...                  im           *
  !     *                                                                *
  !     *  the one-dimensional index k increases going from left to      *
  !     *  right starting from the left lower corner. after reaching the *
  !     *  right boundary of a horizontal row, it jumps to the left      *
  !     *  boundary of the next horizontal row above.                    *
  !     *                                                                *
  !     *  wind points v* along the left boundary have j values one less *
  !     *  than the line they are in; that is, they are stored in places *
  !     *  denoted by (v).                                               *
  !     *                                                                *
  !     *  the total number of grid points in the horizontal is          *
  !     *       imjm=im*jm-jm/2                                          *
  !     *                                                                *
  !     *  the index increments needed to reach n, s, ne, nw, sw, and se *
  !     *  grid point neighbors are for height points respectively       *
  !     *       khn=im-1; khs=-im;                                       *
  !     *       kne=im; knw=im-1; ksw=-im; kse=-im+1;                    *
  !     *  and for velocity points respectively                          *
  !     *       kvn=im; kvs=-im+1;                                       *
  !     *       kne=im; knw=im-1; ksw=-im; kse=-im+1.                    *
  !     *                                                                *
  !     *  the horizontal loop boundaries are treated differently in     *
  !     *  scalar and vector parts of the code. in scalar parts of the   *
  !     *  code the lower and the upper boundaries for index k are       *
  !     *  defined and stored for each horizontal row (see subroutine    *
  !     *  'econ' in the programe econvk for details). in vector parts   *
  !     *  of the code the upper and the lower boundaries are defined    *
  !     *  for the whole horizontal field. these boundaries are defined  *
  !     *  in 'parameter' statements in each subroutine.                 *
  !     *                                                                *
  !     *  the second letters in loop boundary names indicate whether    *
  !     *  they define the boundaries for height (h) or velocity (v)     *
  !     *  points.  the third letter is l for low or h for high.         *
  !     *  the suffix digit(s) indicate(s) for how many lines the loop   *
  !     *  boundaries are shifted from the boundaries of the integration *
  !     *  domain.                                                       *
  !     *                                                                *
  !     *  boundary masks and topography masks are defined for vector    *
  !     *  processing. the boundary masks hbm2(k) and vbm2(k) are        *
  !     *  equal to one everywhere except at the two outermost rows      *
  !     *  where they are equal to zero. the boundary mask vbm3(k) is    *
  !     *  equal to one everywhere except at the three outermost rows    *
  !     *  where it is equal to zero. the topography masks (htm(k,l)     *
  !     *  and vtm(k,l)) are set to zero underneath the topography and   *
  !     *  to one elswhere. in addition, for treatment of physical       *
  !     *  processes, maximum values of the vertical index are defined   *
  !     *  and stored (lmh(k) and lmv(k), see subroutine econ in the     *
  !     *  program econvk for details).                                  *
  !     *                                                                *
  !     ******************************************************************
  IMPLICIT NONE
  INCLUDE "ctlblk.f90"

  !*** initialization:
  !*** initialize constants and variables

  CALL init

  !*** model history file

  CALL chkout 

  !*** pgf,coriolis,divergence & horizontal part of omega-alpha

  CALL pfdht

  !*** divergence damping and boundary values for restart

  IF (restrt) THEN
     CALL ddamp
     CALL bocov
     restrt =.FALSE.
  ENDIF

  !*** end of initialization 
  !*** (in fact, variable first is not needed for loop control;
  !***  kept for security, since it may be needed elsewhere)

  first = .false. 

  !*** time loop

  DO

     !*** loop control

     IF (ntsd .GE. ntstm) EXIT
     ntsd = ntsd+1

     !*** press. tend.,eta dot & vertical omega-alpha

     CALL pdte

     !*** is it time for a check point on the model history file?

     CALL chkomg

     !*** is it time for advection?

     IF (MOD(ntsd-1,idtad) .EQ. 0) THEN
        CALL vtadv
        CALL hzadv
        CALL hzadv2
     ENDIF

     !*** updating pressure difference

     CALL pdnew

     !*** radiation

     IF (MOD(ntsd-1,nrads) .EQ. 0 .OR. MOD(ntsd-1,nradl) .EQ. 0) THEN
        CALL radtn
     ENDIF

     !*** apply temperature tendency due to radiation

     CALL rdtemp

     !*** turbulent processes and precipitation

     IF (MOD(ntsd-1,nphs ) .EQ. 0) THEN
        CALL turbl
        CALL gscond
     ENDIF
     IF (MOD(ntsd-1,ncnvc) .EQ. 0) THEN
        CALL cucnvc
     ENDIF
     IF (MOD(ntsd-1,nphs ) .EQ. 0) THEN
        CALL precpd
     ENDIF

     !*** lateral diffusion

     CALL hdiff

     !*** updating boundary values at height points

     CALL bocoh

     !*** model history file

     CALL chkout 

     !*** pgf,coriolis,divergence & horizontal part of omega-alpha

     CALL pfdht

     !*** divergence damping

     CALL ddamp

     !*** updating boundary values at velocity points

     CALL bocov

  END DO
END PROGRAM ebu
