C$$$   SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    HYBINT      VERTICALLY INTERPOLATE HYB-B GRIDS
C
C   PRGMMR:  BENJAMIN, STAN ORG: ERL/PROFS      DATE: 93-01-12
C
C ABSTRACT:  
C ***   This routine reads in a hybrid-b (non-intersecting
C       surfaces) analysis/forecast and interpolates the
C       data to isobaric levels and pressure variables (h,T)
C       and writes out a MAPS pressure analysis (MOA) or forecast 
C       (MDP) file.
C       It also interpolates the hybrid-b data back to hybrid-a
C       coordinates (overlapping theta and sigma levels)
C       and writes that data out as IOA/IDP files.
C       Finally, freezing level and maximum wind level fields
C       are also calculated and written out.
C
C PROGRAM HISTORY LOG:
C    93-01-12      S. BENJAMIN     ORIGINAL VERSION
C    93-10-08      S. BENJAMIN     NMC VERSION WITH FREEZING LEVEL
C                                   AND MAX WIND LEVEL AND
C                                   50-MB ISOBARIC DATA
C    94-02-17      S. BENJAMIN     ADDED BOUNDARY-LAYER FIELDS
C    94-03-04      S. BENJAMIN     CHANGED TO MULTIPLE 30-MB
C                                   MEAN FIELDS
C    94-05-31      S. BENJAMIN     CHANGED TO 25-MB RESOLUTION FOR
C                                   ISOBARIC LEVELS
C    94-05-31      S. BENJAMIN     CHANGED TO WRITE OUT INITIALIZED
C                                   FIELDS 
C    96-01         S. Benjamin     40-km version
C    99-06         S. Benjamin     NCEP/FSL combined version
C
C USAGE:           CALL   HYBINT (ATIME,IHR,G3,MAPSLP,ISTATUS)
C
C   INPUT ARGUMENT LIST:
C     ATIME    - CHAR*9   ANALYSIS TIME IN YYJJJHH00 FORMAT
C     IHR      - CHAR*2   2-DIGIT TIME FOR FCST HOUR
C     G3       - REAL     ARRAY of cycled hybrid-b native grids
C     GSFC     - REAL     ARRAY of cycled 2-d grids
C     GSFC_M   - REAL     ARRAY of extra 2-d grids for MORDS
C     GEX2     - REAL     ARRAY of still more extra 2-d grids for RUC output and for MORDS
C     MAPSLP   - REAL     ARRAY of MAPS sea level pressure
C   OUTPUT ARGUMENT LIST:
C     ISTATUS  - INTEGER  RETURN STATUS
C
C REMARKS: NONE
C
C ATTRIBUTES:
C   MACHINE: DEC - VAX, VMS
C
C$$$
        SUBROUTINE HYBINT (ATIME,IHR,G3,gsfc,gsfc_m,gex2,
     1     isfc,ISTATUS)


C ***   This routine takes  a hybrid-b (non-intersecting
C       surfaces) analysis/forecast and interpolates the
C       data to isobaric levels and pressure variables (h,T)
C       and writes out a MAPS pressure analysis (MOA) or forecast 
C       (MDP) file
C       It also interpolates the hybrid-b data back to hybrid-a
C       coordinates (overlapping theta and sigma levels)
C       and writes that data out as IOA/IDP files.

C       History
C 28 April 1988      S. Benjamin     Original version
C 18 April 1989      S. Benjamin     18-level 80 km version
C 19 July 1989       S. Benjamin    Extrapolated to extra 420K level
C                                    to ensure definition of variables
C                                    at 150mb.  (395K was occurring
C                                    beneath the 150mb level in Canada.)
C 18 Jan 1990        S. Benjamin     Hybrid version
C 21 Mar 1990        S. Benjamin     Hybrid version after sanitizing in
C                                    analysis (hydrostatic consistency,
C                                    theta and pressure both indicating
C                                    which theta level is 1st above
C                                    interface.
C    22 Oct 1990     S. Benjamin     Changed from GRDTRNS to VERINT,
C                                    ESAT to ESW.
C 13 Dec 1990        S. Benjamin     Added thermal wind correction in
C                                    extrapolating to 500 K.
C 17 Dec 1990        S. Benjamin     De-virtualized temps for isobaric 
C                                    levels
C 14 Jan 1991        S. Benjamin     60km version
C 26 Jul 1992        S. Benjamin     For hybrid-b coordinates.
C    Dec 1995        S. Benjamin     40-km version

        INCLUDE 'IMPLICIT'

        INCLUDE 'MAPSCON'
        INCLUDE 'METCON'
        include 'LOCATION'

        INTEGER NZPN_P
C ***   NZPN_P = NUMBER OF 25-MB ISOBARIC LEVELS, 1000-100 MB
        PARAMETER (NZPN_P = NZP25_P)

C ***   NZM_P = number of mean layers
        INTEGER NZMEAN_P
        PARAMETER (NZMEAN_P = 6)

        INTEGER NVARP_P
C ***  NVARP_P = NUMBER OF 3-D VARIABLES for isobaric grids
        PARAMETER (NVARP_P = 6)
        integer nvarp_main_p,nvarp_supp_p
C ***  NVARP_MAIN_P = NO OF 3-D VARS FOR ISOBARIC MAIN FILE AT NCEP
        PARAMETER (NVARP_MAIN_P = 6)

C   extra variables
C      VERTICAL VELOCITY
C      POTENTIAL VORTICITY 
C  --- NX3D_P = number of extra 3-d variables for isobaric grids
        INTEGER NX3D_P
        PARAMETER (NX3D_P = NVARP_P-NVAR_P+1)
C  If NVARP_P is changed, be sure to add to G1 equivalence
C     statement for new variable.

C   --- number of mean-layer variables
        INTEGER NVARM_P
        PARAMETER (NVARM_P=5)
C   --- number of special level variables
        INTEGER NVARS_P
        PARAMETER (NVARS_P=6)
C   --- number of surface variables
        INTEGER NVARSFC_P
        PARAMETER (NVARSFC_P=12)
C   --- number of surface variables - NCEP
        INTEGER NVARSFC_NCEP_P
        PARAMETER (NVARSFC_NCEP_P=25)
C   --- number of total 2-d variables
        INTEGER NVARTOT2D_P
        PARAMETER (NVARTOT2D_P=88)
c       PARAMETER (NVARTOT2D_P=87)
C   --- number of MORDS upper-air variables
        INTEGER NVAR_MORDS_P
        PARAMETER (NVAR_MORDS_P=16)

C     Cloud_def_p       = 0.000001      g/g cloud top threshold for model
        real    Cloud_def_p
        data  Cloud_def_p       / 0.000001/


      
        INTEGER STATUS,ISTATUS,ISFC,nhr

        INTEGER  I,J,K,IVAR,K1

        real rhmax,rhmin
        integer imax,jmax,imin,jmin
        REAL     G(NX_P,NY_P,NVARTOT2D_P)
        REAL     Gx(NX_P,NY_P,NVARm_P+1)
        REAL     GM(NX_P,NY_P,NVAR_MORDS_P)
        REAL     G3(NX_P,NY_P,NZTN_P,NVAR3_P)
     1  ,G3x(NX_P,NY_P,nztn_p)
     1  ,GSFC(NX_P,NY_P,MXVRS_P)
C Gridded 2-d fields (cape/cin, trop fields, MAPS SLP)
     1  ,gsfc_m(nx_p,ny_p,mxvrsgcm_p)

        REAL     G1MEAN(NX_P,NY_P,NZMEAN_P,NVARM_P)
C-- GSFC_EX are extra 2d fields for stations with class=0,1
     1  ,GEX2(nx_p,ny_p,MXVREX_P)

        REAL USTORM(nx_p,ny_p)
     1      ,VSTORM(nx_p,ny_p),HELI(nx_p,ny_p)
     1      ,watericetotal(nztn_p),watericemax
     1      ,zcldbase,zcldtop,zpbltop,windsp,windsp_sfc
     1      ,dz,del_windsp
     1      ,vis(nx_p,ny_p)

C     *** native (hybrid) coordinate VARIABLES
        REAL
     1   TB    (NX_P,NY_P,NZTN_P)
     1  ,RHB   (NX_P,NY_P,NZTN_P)

        REAL
     1   CORIO(NX_P,NY_P)
C Coriolis parameter - gridpts
     1  ,SCALE(NX_P,NY_P)
C Map scale factor
     1  ,TOPO(NX_P,NY_P)
C Mini topography


C *** PRESSURE VARIABLES
        REAL     H     (NX_P,NY_P,NZPN_P)
     1  ,T     (NX_P,NY_P,NZPN_P)
     1  ,RHP   (NX_P,NY_P,NZPN_P)
     1  ,UP    (NX_P,NY_P,NZPN_P)
     1  ,VP    (NX_P,NY_P,NZPN_P)
     1  ,WP    (NX_P,NY_P,NZPN_P,NX3D_P)

C *** FREEZING LEVEL VARIABLES
        REAL
     1   FRZ    (NX_P,NY_P)
     1  ,FRP    (NX_P,NY_P)
     1  ,FRH    (NX_P,NY_P)
        REAL
     1   FRZ2   (NX_P,NY_P)
     1  ,FRP2   (NX_P,NY_P)
     1  ,FRH2   (NX_P,NY_P)

C *** MAXIMUM WIND LEVEL VARIABLES
        REAL
     1   PMX    (NX_P,NY_P)
     1  ,UMX    (NX_P,NY_P)
     1  ,VMX    (NX_P,NY_P)

C *** BOUNDARY LAYER VARIABLES
c       REAL
c    1   RHBL    (NX_P,NY_P)
c    1  ,UBL    (NX_P,NY_P)
c    1  ,VBL    (NX_P,NY_P)
c    1  ,TBL    (NX_P,NY_P)
c    1  ,wBL    (NX_P,NY_P)

C *** SURFACE VARIABLES
        REAL
     1   DUMMY  (NX_P,NY_P)

      INTEGER IPROJ,inum,iq
      CHARACTER*9 P3TIME
      CHARACTER*2 P3PROJ
      REAL Z1,T1,Z5,T5,GAM,TSFC,DEWPT

      real s,t2,dlt,tl,tk,pt,thetae,thetaemax(nx_p,ny_p),
     1   tqk,w,x,xli,pol,
     1   b,ftim,qv,tx,esx,enl,dwpt

      real totprec,snowratio


      real tq(nx_p,ny_p),d(nx_p,ny_p)

      real tmpmax(nx_p,ny_p)

        REAL     G1(NX_P,NY_P,NZPN_P,NVARP_P)

        EQUIVALENCE (Gx(1,1,1),FRZ(1,1))
        EQUIVALENCE (Gx(1,1,2),FRP(1,1))
        EQUIVALENCE (Gx(1,1,3),FRH(1,1))

        EQUIVALENCE (Gx(1,1,4),PMX(1,1))
        EQUIVALENCE (Gx(1,1,5),UMX(1,1))
        EQUIVALENCE (Gx(1,1,6),VMX(1,1))

c       EQUIVALENCE (G (1,1,1),TBL(1,1))
c       EQUIVALENCE (G (1,1,2),RHBL(1,1))
c       EQUIVALENCE (G (1,1,3),UBL(1,1))
c       EQUIVALENCE (G (1,1,4),VBL(1,1))

        EQUIVALENCE (G1(1,1,1,1),H(1,1,1))
        EQUIVALENCE (G1(1,1,1,2),T(1,1,1))
        EQUIVALENCE (G1(1,1,1,3),RHP(1,1,1))
        EQUIVALENCE (G1(1,1,1,4),UP(1,1,1))
        EQUIVALENCE (G1(1,1,1,5),VP(1,1,1))
        EQUIVALENCE (G1(1,1,1,6),WP(1,1,1,1))
c       EQUIVALENCE (G1(1,1,1,7),WP(1,1,1,2))
c       EQUIVALENCE (G1(1,1,1,8),WP(1,1,1,3))
c       EQUIVALENCE (G1(1,1,1,9),WP(1,1,1,4))
c       EQUIVALENCE (G1(1,1,1,10),WP(1,1,1,5))
c       EQUIVALENCE (G1(1,1,1,11),WP(1,1,1,6))
c       EQUIVALENCE (G1(1,1,1,12),WP(1,1,1,7))
c       EQUIVALENCE (G1(1,1,1,13),WP(1,1,1,8))
c       EQUIVALENCE (G1(1,1,1,14),WP(1,1,1,9))

        REAL     EXNER(nx_p,ny_p,nztn_p),Q,TEMP,ES,E,ESW

        REAL     PLEV (NZPN_P)

        DATA     PLEV /
     1          1000.,975.,950.,925.,900.,875.,850.,825.,800.,
     1           775.,750.,725.,700.,675.,650.,625.,600.,575.,550.,
     1           525.,500.,475.,450.,425.,400.,375.,350.,325.,300.,
     1           275.,250.,225.,200.,175.,150.,125.,100./

        CHARACTER ATIME*9 , IHR*2, ADATE*9, APROJ*4
        character bdate*9, bproj*4
        CHARACTER*80 NETDIR
        INTEGER GUNIT,KP1,I1,I2,J1,J2
        REAL  STAB, VORT

c     INTEGER*4 NETVARID2(1)
c     INTEGER*4 NETVARID3(NVAR_P+1)
c     INTEGER*4 NETVARID3(NVAR_P)

        data  b/ 2.6518986/

        character*60 varname2
        character*60 varnamex
        character*60 varnamep
        character*60 varname3p(NVARP_P)
        character*60 varname_mean(NVARM_P)
        character*60 varname_speciallev(NVARS_P)
        character*60 varname_sfc(NVARSFC_P)
        character*60 varname_sfc_ncep(NVARSFC_ncep_P)
        character*60 varname_2d_ncep(NVARtot2d_P)
        character*60 varname_mords_upr(NVAR_MORDS_P)


        data varname_mords_upr /
     1  '700 hPa vertical velocity',
     1  '250 hPa height',
     1  '250 hPa u-wind component',
     1  '250 hPa v-wind component',
     1  '500 hPa height',
     1  '500 hPa absolute vorticity',
     1  '850 hPa height'
     1  ,'850 hPa temperature'
     1  ,'850 hPa specific humidity'
     1  ,'850 hPa u-wind component'
     1  ,'850 hPa v-wind component'
     1  ,'1000 hPa height'
     1  ,'mean temperature in 1st 30hPa layer above ground'
     1  ,'spec hum in 1st 30hPa layer above ground'
     1  ,'u-wind component in 1st 30hPa layer above ground'
     1  ,'v-wind component in 1st 30hPa layer above ground'
     1  /

        data varname3p /
     1  'height - p',
     1  'temperature - p',
     1  'relative humidity - p',
     1  'u-component of wind - p',
     1  'v-component of wind - p',
     1  'vertical velocity - p'
     1  /

        data varname_mean /
     1  'temperature',
     1  'relative humidity',
     1  'u-component of wind',
     1  'v-component of wind',
     1  'pressure vertical velocity'
     1  /

        data varname_speciallev /
     1  'height - frz'
     1 ,'pressure - frz'
     1 ,'relative humidity - frz'
     1 ,'pressure at max wind level'
     1 ,'u-component of wind at max wind level'
     1 ,'v-component of wind at max wind level'
     1  /

        data varname2 /
     1  'MAPS mean sea level pressure'
     1  /

        data varnamep /
     1  'pressure - b'
     1  /

        data varnamex /
     1  '       '
     1  /

        data varname_sfc /
     1  'MAPS mean sea level pressure'
     1 ,'altimeter setting'
     1 ,'3-hr pressure tendency - s'
     1 ,'potential temperature'
     1 ,'dew point temperature'
     1 ,'dew point depression'
     1 ,'u-component of wind - s'
     1 ,'v-component of wind - s'
     1 ,'temperature - s'
     1 ,'specific humidity'
     1 ,'surface pressure'
     1 ,'relative humidity - s'
     1   /

        data varname_sfc_ncep /
     1  'MAPS mean sea level pressure'
     1 ,'altimeter setting'
     1 ,'3-hr pressure tendency'
     1 ,'potential temperature, 2m fixed height'
     1 ,'dew point temperature, 2m fixed height'
     1 ,'dew point depression, 2m fixed height'
     1 ,'u-component of wind, 10m fixed height'
     1 ,'v-component of wind, 10m fixed height'
     1 ,'temperature, 2m fixed height'
     1 ,'specific humidity, 2m fixed height'
     1 ,'equivalent potential temperature'
     1 ,'convective available potential energy'
     1 ,'convective inhibition'
     1 ,'lifted index'
     1 ,'best lifted index'
     1 ,'resolvable scale precipitation - 3h'
     1 ,'sub-grid scale precipitation - 3h'
     1 ,'resolvable scale precipitation - 1h'
     1 ,'sub-grid scale precipitation - 1h'
     1 ,'categorical rain'
     1 ,'categorical freezing rain'
     1 ,'categorical ice pellets'
     1 ,'categorical snow'
     1 ,'snow accumulation'
     1 ,'helicity'
     1   /

       data varname_2d_ncep /
     1  'MAPS mean sea level pressure'
     1 ,'altimeter setting'
     1 ,'3-hr pressure tendency'
     1 ,'potential temperature, 2m fixed height'
     1 ,'dew point temperature, 2m fixed height'
     1 ,'dew point depression, 2m fixed height'
     1 ,'u-component of wind, 10m fixed height'
     1 ,'v-component of wind, 10m fixed height'
     1 ,'temperature, 2m fixed height'
     1 ,'specific humidity, 2m fixed height'
     1 ,'equivalent potential temperature'
     1 ,'convective available potential energy'
     1 ,'convective inhibition'
     1 ,'lifted index'
     1 ,'best lifted index'
     1 ,'resolvable scale precipitation - 3h'
     1 ,'sub-grid scale precipitation - 3h'
     1 ,'resolvable scale precipitation - 1h'
     1 ,'sub-grid scale precipitation - 1h'
     1 ,'categorical rain'
     1 ,'categorical freezing rain'
     1 ,'categorical ice pellets'
     1 ,'categorical snow'
     1 ,'snow accumulation - 3h'
     1 ,'storm-relative helicity'
     1 ,'pressure - frz'
     1 ,'height - frz'
     1 ,'relative humidity - frz'
     1 ,'pressure at tropopause'
     1 ,'potential temperature at tropopause'
     1 ,'u-component of wind at tropopause'
     1 ,'v-component of wind at tropopause'
     1 ,'pressure at max wind level'
     1 ,'u-component of wind at max wind level'
     1 ,'v-component of wind at max wind level'
     1 ,'precipitation rate'
     1 ,'soil moisture availability'
     1 ,'surface runoff'
     1 ,'temperature - m1'
     1 ,'relative humidity - m1'
     1 ,'u-component of wind - m1'
     1 ,'v-component of wind - m1'
     1 ,'vertical velocity - m1'
     1 ,'temperature - m2'
     1 ,'relative humidity - m2'
     1 ,'u-component of wind - m2'
     1 ,'v-component of wind - m2'
     1 ,'vertical velocity - m2'
     1 ,'temperature - m3'
     1 ,'relative humidity - m3'
     1 ,'u-component of wind - m3'
     1 ,'v-component of wind - m3'
     1 ,'vertical velocity - m3'
     1 ,'temperature - m4'
     1 ,'relative humidity - m4'
     1 ,'u-component of wind - m4'
     1 ,'v-component of wind - m4'
     1 ,'vertical velocity - m4'
     1 ,'temperature - m5'
     1 ,'relative humidity - m5'
     1 ,'u-component of wind - m5'
     1 ,'v-component of wind - m5'
     1 ,'vertical velocity - m5'
     1 ,'temperature - m6'
     1 ,'relative humidity - m6'
     1 ,'u-component of wind - m6'
     1 ,'v-component of wind - m6'
     1 ,'vertical velocity - m6'
     1 ,'surface pressure at mini topography'
     1 ,'earth surface height'
     1 ,'pressure at top freezing level'
     1 ,'geopotential height at top freezing level'
     1 ,'relative humidity at top freezing level'
     1 ,'snow depth'
     1 ,'precipitable water'
     1 ,'relative humidity, 2m fixed height'
     1 ,'sub-surface runoff'
     1 ,'u-component of storm motion'
     1 ,'v-component of storm motion'
     1 ,'PBL height'
     1 ,'gust wind speed'
     1 ,'cloud base height'
     1 ,'cloud top height'
     1 ,'visibility'
     1 ,'maximum equivalent potential temperature pressure'
     1 ,'convective cloud top height'
     1 ,'equilibrium level height'
     1 ,'temperature at tropopause'
     1   /


       real const,mslp,gams,gamd,tstd,coninv,gami,govrgam
C    1          const,  ! = g / (R * GAMMA_STD)
C    1          coninv, ! = (R * GAMMA_STD) / g
C    1          mslp,   ! Sea-level pressure for Std. Atmos.
C    1          gams ,  ! = Std. Atmos. lapse rate = 6.5 K/1000 m
C    1          gami ,  ! = isothermal lapse rate = 0.02 K/1000 m
C               gamd    ! = dry adiabatic lapse rate = 10 K/1000 m
C    1          tstd,   ! = Temp at sea-level in Std. Atmos. = 288.15K
      DATA             CONST/5.255303/,
     1          MSLP/1013.25/,
     1          GAMS/0.0065/,
     1          GAMD/0.0100/,
     1          GAMi/0.0002/,
     1          TSTD/288.15/
     1          CONINV /0.190284/


        ATIME = ATIME(1:7)//IHR


C *** READ IN GRDPARM FILE FOR CORIOLIS AND MAP-SCALE FACTOR ARRAYS
      CALL GETLUN(GUNIT,1)
      gunit = 13
        OPEN (UNIT=GUNIT, FILE='GRDPARM',STATUS='OLD',
     1  FORM='UNFORMATTED')
        READ (GUNIT) CORIO,SCALE
        CLOSE (GUNIT)
      CALL GETLUN(GUNIT,2)


         adate = atime(1:7)//'00'
         aproj = '00'//atime(8:9)

        ISTATUS = SUCCESS_P

C ************************************************************
C ************************************************************
C     Calculate RH on native levels
C ************************************************************
C ************************************************************

        DO 20 K = 1,NZTN_P
        DO 20 J = 1,NY_P
        DO 20 I = 1,NX_P

          EXNER(i,j,k) = CPD_P*(g3(I,J,K,p_p)/100000.)**ROVCP_P
          qv = g3(i,j,k,qv_p)
          Q = qv/(1.+qv)
          TEMP = g3(I,J,K,theta_p)*EXNER(i,j,k)/(CPD_P*(1.+0.6078*Q))
          TB(I,J,K) = TEMP
        Tx = temp-273.15
        POL = 0.99999683       + TX*(-0.90826951E-02 +
     1     TX*(0.78736169E-04   + TX*(-0.61117958E-06 +
     2     TX*(0.43884187E-08   + TX*(-0.29883885E-10 +
     3     TX*(0.21874425E-12   + TX*(-0.17892321E-14 +
     4     TX*(0.11112018E-16   + TX*(-0.30994571E-19)))))))))
        esx = 6.1078/POL**8

          ES = esx
          E = g3(I,J,K,p_p)/100.*QV/(0.62197+QV)
          RHB(I,J,K) = 100.*AMIN1(1.,E/ES)
20      CONTINUE


C -- Convert PLEV to Pascals from mb
        do k=1,nzpn_p
          plev(k) = plev(k)*100.
        end do



C -- NX3D is number of extra 3-d files to interpolate to
C     isobaric levels.  
C *** INTERPOLATE TO 25-MB RESOLUTION ISOBARIC LEVELS
        CALL HB2P (g3,RHB,
     1     H,T,RHP,UP,VP,WP,NX3D_P,PLEV,NZPN_P)

C --- Print out isobaric values
        do 30 k=1,nzpn_p
          print *,' isobaric values',k,(g1(15,15,k,ivar),
     1       ivar=1,nvarp_p)
30      continue

c -- Calculate helicity and storm motion vector

        call helic (g3,nx_p,ny_p,nztn_p,nvar_p,heli,
     1     ustorm,vstorm)


C ************************************************************
C ************************************************************
C     Calculate PBL height, wind gust speed
C ************************************************************
C ************************************************************

        do 35 j=1,ny_p
        do 35 i=1,nx_p
C PBL height
          g(i,j,80) = 0.
C gust potential
          g(i,j,81) = sqrt(g3(i,j,1,u_p)**2 + g3(i,j,1,v_p)**2)
          windsp_sfc = g(i,j,81)
C   Is there any mixed layer at all?
          if (g3(i,j,2,theta_p) .lt. g3(i,j,1,theta_p)) then

C   Calculate k1 level as first above PBL top
            do 34 k=3,nztn_p
              k1 = k
              if (g3(i,j,k,theta_p).gt.g3(i,j,1,theta_p)) go to 341
 34         continue
 341        zpbltop = g3(i,j,k1,h_p) + 
     1             (g3(i,j,1,theta_p)-g3(i,j,k1,theta_p))
     1           * (g3(i,j,k1-1,h_p)-g3(i,j,k1,h_p))
     1           / (g3(i,j,k1-1,theta_p) - g3(i,j,k1,theta_p))
            g(i,j,80) = zpbltop - g3(i,j,1,h_p)

C - do some kind of mixing with boundary layer depth.
C     Figure that only 50% of excess with speed can
C     be mixed to surface from 1 km level.

            do 342 k=2,k1-1
              windsp = sqrt(g3(i,j,k,u_p)**2 + g3(i,j,k,v_p)**2)
              del_windsp = windsp - windsp_sfc
              dz = g3(i,j,k,h_p) - g3(i,j,1,h_p)
              del_windsp = del_windsp*(1.0 - 
     1            min(0.5,dz/2000.))
C   2000 for 0.50 (50%) and 1 km depth
              g(i,j,81) = max(g(i,j,81),windsp_sfc+del_windsp)
 342        continue
          end if
 35     continue

c       CALL SMOOTH(g(1,1,80),DUMMY,NX_P,NY_P,0.5)
c       CALL SMOOTH(g(1,1,81),DUMMY,NX_P,NY_P,0.5)

C ************************************************************
C ************************************************************
C     Create cloud top/base and visibility grids
C ************************************************************
C ************************************************************

        do 37 j=1,ny_p
        do 37 i=1,nx_p
          g(i,j,82) = -1.*spval_p
          g(i,j,83) = -1.*spval_p
          watericemax = -99999.
          do k=1,nztn_p
            watericetotal(k) = g3(i,j,k,qc_p) + g3(i,j,k,qi_p)
            watericemax = max(watericemax,watericetotal(k))
          end do
          if (watericemax.lt.cloud_def_p) go to 37
C  Cloud base
C       At surface?
          if (watericetotal(1).gt.cloud_def_p) then
            zcldbase = g3(i,j,1,h_p)
            go to 3788
          end if
C       Aloft?
          do 371 k=2,nztn_p
            k1 = k
            if (watericetotal(k).gt.cloud_def_p) go to 372
 371      continue
          go to 37
 372      zcldbase = g3(i,j,k1,h_p) + (cloud_def_p-watericetotal(k1))
     1           * (g3(i,j,k1-1,h_p)-g3(i,j,k1,h_p))
     1           / (watericetotal(k1-1) - watericetotal(k1))
 3788   continue

C  Cloud top
          if (watericetotal(nztn_p).gt.cloud_def_p) then
            zcldtop = g3(i,j,nztn_p,h_p)
            go to 3799
          end if
          do 373 k=nztn_p,2,-1
            k1 = k
            if (watericetotal(k).gt.cloud_def_p) go to 374
 373      continue
          go to 37
 374      zcldtop = g3(i,j,k1,h_p) + (cloud_def_p-watericetotal(k1))
     1           * (g3(i,j,k1+1,h_p)-g3(i,j,k1,h_p))
     1           / (watericetotal(k1+1) - watericetotal(k1))
 3799     continue
          g(i,j,82) = zcldbase
          g(i,j,83) = zcldtop
 37    continue

C ************************************************************
C ************************************************************
C     Create visibility grid
C ************************************************************
C ************************************************************

       call calvis (g3(1,1,1,qv_p), g3(1,1,1,qc_p),
     1             g3(1,1,1,qr_p), g3(1,1,1,qi_p),
     1             g3(1,1,1,qn_p), g3(1,1,1,qg_p),
     1             g3(1,1,1,theta_p), g3(1,1,1,p_p),
     1             exner(1,1,1),rhb(1,1,2),vis,nx_p,ny_p)

C--- smooth log of visibility
       do j=1,ny_p
       do i=1,nx_p
         g(i,j,84) = alog(vis(i,j))
       end do
       end do

c       CALL SMOOTH(g(1,1,84),DUMMY,NX_P,NY_P,0.5)

C--- change from log(vis) back to vis and convert to meters.
       do j=1,ny_p
       do i=1,nx_p
         g(i,j,84) = 1000.*exp(g(i,j,84))
       end do
       end do

C ************************************************************
C ************************************************************
C     Put extra grids from CAPCIN into g array
C ************************************************************
C ************************************************************
C max theta-e pressure, convective cloud top ht, equilibrium level ht

       do j=1,ny_p
       do i=1,nx_p
         g(i,j,85) = gex2(i,j,16)
         g(i,j,86) = gex2(i,j,17)
         g(i,j,87) = gex2(i,j,18)
       end do
       end do

C ************************************************************
C ************************************************************
C     Create surface grids
C ************************************************************
C ************************************************************

C --- Read in topomini file
        CALL IN2FIL (TOPO,NX_P,NY_P,NX_P,NY_P,'TOPMINI',STATUS)


        print *,' Start loop 120'

        do 120 j=1,ny_p
        do 120 i=1,nx_p
C --- MAPS SLP
          g(i,j,1) = gsfc(i,j,mslp_p)
C      ---   z = surface elevation
          z1 = g3(i,j,1,h_p)
C --- altimeter setting -- account for cockpit height
C        3 m (~0.3 mb) above runway
          g(i,j,2) = (g3(i,j,1,p_p)-30.) * ((tstd + z1*gams)/tstd)
     1          ** const
C --- q = specific humidity
          qv = g3(i,j,1,qv_p)
          Q = QV/(1.+QV)
          T1 = g3(I,J,1,theta_p)*EXNER(i,j,1)/(CPD_P*(1.+0.6078*Q))
          e=g3(i,j,1,p_p)/100.*qv/(0.62197+qv) 
C --- dew-point temperature at original sfc
c         g(i,j,5) = dewpt(e)+273.15
        ENL = ALOG(E)
        DWPT = (243.5*ENL-440.8)/(19.48-ENL)

          g(i,j,5) = dwpt+273.15
C      ---   get values at level 5 for lapse rate calculations
          z5 = g3(i,j,5,h_p)
          Q = g3(I,J,5,qv_p)/(1.+g3(i,j,5,qv_p))
          T5 = g3(I,J,5,theta_p)*EXNER(i,j,5)/(CPD_P*(1.+0.6078*Q))
          GAM = (T1-T5)/(Z5-Z1)
          GAM = MIN(GAMD,MAX(GAM,GAMi))
C      --- temperature at topomini sfc
          tsfc = t1 + (z1-topo(i,j))*gam
C --- pressure at topomini
C --- put topomini pressure as sfc pressure for Isobaric main file
C --- govrgam = g / (R * gam)
          govrgam = g0_p / (rd_p * gam)
          g(i,j,69) = g3(i,j,1,p_p)*(tsfc/t1)**govrgam
C --- sfc pressure at model topo for variable 11 in MAPS_SFC file
C       best for 3h pressure change
          g(i,j,11) = g3(i,j,1,p_p)
C --- temperature
          g(i,j,9) = tsfc
C --- potential temperature at topomini sfc
          g(i,j,4) = tsfc*(100000./g(i,j,11))**rovcp_p
C       Set dewpoint depression to that at original sfc
c --- dewpoint depression
          g(i,j,6) = max(0.,t1 - g(i,j,5) )
C --- dew-pt at topomini
          g(i,j,5) = tsfc - g(i,j,6)
C --- surface winds
          g(i,j,7) = g3(i,j,1,u_p)
          g(i,j,8) = g3(i,j,1,v_p)
C --- specific humidity
          g(i,j,10) = Q
          tmpmax(i,j) = -99999.
C --- sfc RH
          g(i,j,12) = rhb(i,j,1)
120     continue

C -- Calculate max temp up to level 20
        do 31 k=2,20
        do 31 j=1,ny_p
        do 31 i=1,nx_p
          tmpmax(i,j) = max(tmpmax(i,j),tb(i,j,k))
31      continue

C *******
C  Read in sfc pressure field from 3 h ago for 3-h pressure change

C --- Get integer for projection time
       read (atime(8:9),'(i2)') iproj
       print *,' iproj',iproj
       if (iproj.lt.3) then
         call chngtm (atime(1:7)//'00',-10800 + iproj*3600,
     +      p3time,istatus)
       else
         write (p3proj,'(i2.2)') iproj-3
         p3time = atime(1:7)//p3proj
       end if
         print *,' p3time',p3time
         bdate = p3time(1:7)//'00'
         bproj = '00'//p3time(8:9)
         netdir = 'BACK_SFC '
         call readg (netdir,bdate,bproj,
c    +     nx_p,ny_p,nztn_p,1,nvarsfc_p,1,
     +     nx_p,ny_p,1,1,nvarsfc_p,1,
     +     varname_sfc,
     +     varnamex,gm,g3x,istatus)
          if (istatus.ne.0 .or. g3(2,2,1,1).eq.0.) then
            print *,'istatus',istatus
            print *,' file not found', bdate,bproj
            do 110 j=1,ny_p
            do 110 i=1,nx_p
              g(i,j,3) = spval_p
110         continue
          else
            do 115 j=1,ny_p
            do 115 i=1,nx_p
C --- 3-h pressure change in Pa/second
c             g(i,j,3) = (g3(i,j,1,p_p) - g3x(i,j,1))/10800.
              g(i,j,3) = (g(i,j,11) - gm(i,j,11))/10800.
115         continue

        print *,' old/new',g(60,60,11),gm(60,60,11)


        CALL SMOOTH(g(1,1,3),DUMMY,NX_P,NY_P,0.5)
        CALL SMOOTH(g(1,1,3),DUMMY,NX_P,NY_P,0.5)

          end if

        CALL SMOOTH(g(1,1,2),DUMMY,NX_P,NY_P,0.5)
        CALL SMOOTH(g(1,1,2),DUMMY,NX_P,NY_P,0.5)

        do 130 ivar=1,nvarsfc_p
          print *,' sfc field, ivar',ivar,g(50,50,ivar)
130     continue

        if (facility_p.eq.fsl_p)
     1   call writeg ('MAPS_SFC ',adate,aproj,9,nx_p,ny_p,
     1                  1,1,nvarsfc_p,1,
     1                  varname_sfc,varnamex,g,g1)

        call writegg ('MAPSG_SFC ',adate,aproj,9,nx_p,ny_p,
     1                  1,1,nvarsfc_p,1,
     1                  varname_sfc,varnamex,g,g1)

        do 132 j=1,ny_p
        do 132 i=1,nx_p
          s = g(i,j,6)
          t2 = g(i,j,9)-273.15

C --- Calculate equivalent potential temperature
C ----------------------------------------------
C - condensation temperature.
C   THE APPROXIMATION BELOW, A THIRD ORDER POLYNOMIAL IN S AND T,
C   IS DUE TO HERMAN WOBUS. THE SOURCE OF DATA FOR FITTING THE  
C   POLYNOMIAL IS UNKNOWN.  (Source - Tom Schlatter - FSL)     
                                                              
          DLT = S*(1.2185+1.278E-03*T2+                       
     1        S*(-2.19E-03+1.173E-05*S-5.2E-06*T2))          
          tl = t2-dlt+273.15
          tk = g(i,j,9)
          w = g3(i,j,1,qv_p)*1000. 
          PT = TK*(100000./g3(i,j,1,p_p))**(0.2854*(1.-0.00028*W))        
          g(i,j,11)= PT*EXP((3.376/TL-0.00254)*W*(1.+0.00081*W))             
C -- CAPE
          g(i,j,12)= gsfc(i,j,14)
C -- CIN
          g(i,j,13)= gsfc(i,j,15)
C  variables for next sets of calculations
          thetaemax(i,j) = 0.
          tq(i,j) = 253.16         
          d(i,j)  = 120.
132     continue

C -- LI/best LI
c       do k=1,nztn_p
C  Assume 20 levels from ground covers up to at least 300 mb
C         above ground, which is the depth that Storm Prediction
C         Center wants for calculating LIs

C  - 1st, calculate theta-e max in lowest 20 levels
        do 1325 k=1,20
        do 1325 j=1,ny_p
        do 1325 i=1,nx_p
          t2 = tb(i,j,k)-273.15
          DLT = S*(1.2185+1.278E-03*T2+
     1      S*(-2.19E-03+1.173E-05*S-5.2E-06*T2))
          tl = t2-dlt+273.15
          tk = tb(i,j,k)
          w = g3(i,j,k,qv_p)*1000.
          PT = TK*(100000./g3(i,j,k,p_p))**(0.2854*(1.-0.00028*W))
          thetae= PT*EXP((3.376/TL-0.00254)*W*(1.+0.00081*W))
          if ((g3(i,j,1,p_p)-g3(i,j,k,p_p)).le.30000.) 
     1        thetaemax(i,j) = max(thetaemax(i,j),thetae)
1325    continue


C -- code to get temp at 500mb for thetaemax
C       first guess 

          do 122 iq = 1,10
            do 1215 j=1,ny_p
            do 1215 i=1,nx_p
c           tqk = tq-273.15
            thetae = thetaemax(i,j)
            d(i,j)   = d(i,j)/2.
        Tx = tq(i,j)-273.15
        POL = 0.99999683       + TX*(-0.90826951E-02 +
     1     TX*(0.78736169E-04   + TX*(-0.61117958E-06 +
     2     TX*(0.43884187E-08   + TX*(-0.29883885E-10 +
     3     TX*(0.21874425E-12   + TX*(-0.17892321E-14 +
     4     TX*(0.11112018E-16   + TX*(-0.30994571E-19)))))))))
        esx = 6.1078/POL**8

c           e   = esw(tqk)
            e   = esx
            w   = 0.62197*e/(500.-e)
C        1.219013 = (1000/500)**(2/7)
            x   = thetae * exp(-b*w*1000./tq(i,j) )
     1            - tq(i,j)*1.219013
c           if (abs(x).lt.1.e-7) go to 222
c           tq(i,j) = tq(i,j) + sign(d,x)
            tq(i,j)=tq(i,j) + sign(d(i,j),x)
1215      continue
122       continue

        do 123 j=1,ny_p
        do 123 i=1,nx_p
C - level 21 should be 500 mb level
          g(i,j,15) = t(i,j,21)-tq(i,j)
          tq(i,j) = 253.16         
          d(i,j)  = 120.
123     continue
C   best LI


C - now get sfc LI
          do 722 iq = 1,10
            do 7215 j=1,ny_p
            do 7215 i=1,nx_p
c           tqk = tq-273.15
            thetae = g(i,j,11)
            d(i,j)   = d(i,j)/2.
        Tx = tq(i,j)-273.15
        POL = 0.99999683       + TX*(-0.90826951E-02 +
     1     TX*(0.78736169E-04   + TX*(-0.61117958E-06 +
     2     TX*(0.43884187E-08   + TX*(-0.29883885E-10 +
     3     TX*(0.21874425E-12   + TX*(-0.17892321E-14 +
     4     TX*(0.11112018E-16   + TX*(-0.30994571E-19)))))))))
        esx = 6.1078/POL**8

c           e   = esw(tqk)
            e   = esx
            w   = 0.62197*e/(500.-e)
C        1.219013 = (1000/500)**(2/7)
            x   = thetae * exp(-b*w*1000./tq(i,j) )
     1            - tq(i,j)*1.219013
c           if (abs(x).lt.1.e-7) go to 222
            tq(i,j) = tq(i,j) + sign(d(i,j),x)
7215      continue
722       continue

C - level 21 should be 500 mb level
        do 723 j=1,ny_p
        do 723 i=1,nx_p
          g(i,j,14) = t(i,j,21)-tq(i,j)
C   sfc LI

C ************************************************************
C ************************************************************
C     Calculate precipitation type
C ************************************************************
C ************************************************************

C -- resolvable scale precip - 3h
          g(i,j,16)= gsfc(i,j,7)
C -- sub-grid scale precip - 3h
          g(i,j,17)= gsfc(i,j,8)
C -- resolvable scale precip - 1h 
          g(i,j,18)= gex2(i,j,13)
C -- sub-grid scale precip - 1h
          g(i,j,19)=gex2(i,j,14)
C -- precip types
C       20 - rain
C       21 - freezing rain
C       22 - ice pellets
C       23 - snow

          g(i,j,20) = 0.
          g(i,j,21) = 0.
          g(i,j,22) = 0.
          g(i,j,23) = 0.
723    continue

        do 133 j=1,ny_p
        do 133 i=1,nx_p
C   accum rain
          totprec = (gsfc(i,j,7)+gsfc(i,j,8))
          snowratio = 0.0
          if (totprec.gt.0.01) snowratio =
     1       gsfc(i,j,33)/totprec

C --   snow
c         if (gex2(i,j,7).gt.0.2e-9) g(i,j,23) = 1.
          if (gex2(i,j,7).gt.0.2e-9
     1        .or. (totprec.gt.0.1.and.snowratio.ge.0.25)) g(i,j,23) = 1.

C --   rain/freezing rain
c         if (gex2(i,j,6).gt.1.e-2) then
          if (gex2(i,j,6).gt.1.e-2
     1        .or. (totprec.gt.0.1 .and.snowratio.lt.0.75)) then

            if (g(i,j,9).ge.273.15) then
C              rain
              g(i,j,20) = 1.
            else if (tmpmax(i,j).gt.273.15) then
C              freezing rain
C              Only allow frz rain if some level above
C               ground is above freezing.
              g(i,j,21) = 1.
            end if
          end if

C -- graupel/ice pellets/snow
          if (gex2(i,j,8).gt.1.e-9) then
            if (g(i,j,9).le.273.15) then
C              check for max rain mixing ratio
C              if it's > 0.05 g/kg, => ice pellets
              if (gex2(i,j,15).gt.0.00005) then
                if (gex2(i,j,8).gt.gex2(i,j,7)) then
C              ice pellets
                  g(i,j,22) = 1.
                else
C              snow
                  g(i,j,23) = 1.
                end if
              else
C              snow
                g(i,j,23) = 1.
              end if
            else
              if (g(i,j,9).ge.275.15) then
C              rain
                g(i,j,20) = 1.
              else
C              snow
                g(i,j,23) = 1.
              end if
            end if
          end if
              
C -- snow accumulation
          g(i,j,24) = gsfc(i,j,33)
C -- helicity
          g(i,j,25) = heli(i,j)            
        
C ---- Fill in fields for BUFR soundings
C    snow y/n
        gex2(i,j,6) = g(i,j,23)
C    ice pellets y/n
        gex2(i,j,7) = g(i,j,22)
C    freezing rain y/n
        gex2(i,j,8) = g(i,j,21)
C    rain y/n
        gex2(i,j,9) = g(i,j,20)
        gex2(i,j,10) = ustorm(i,j)
        gex2(i,j,11) = vstorm(i,j)
        gex2(i,j,12) = heli(i,j)

133     continue

c - set grids without good values at analysis times
C     equal to missing at analysis time
        if (ihr.eq.'00') then
        do 1331 j=1,ny_p
        do 1331 i=1,nx_p
          g(i,j,18) = spval_p
          g(i,j,19) = spval_p
          g(i,j,20) = spval_p
          g(i,j,21) = spval_p
          g(i,j,22) = spval_p
          g(i,j,23) = spval_p
          g(i,j,24) = spval_p
1331    continue
        end if
          
c       call writegg ('MAPSG_SFC ',adate,aproj,9,nx_p,ny_p,
c    1                  1,1,nvarsfc_ncep_p,1,
c    1                  varname_sfc_ncep,varnamex,g,g1)


C ***
C   Write out 25-mb isobaric data
c       if (facility_p.eq.fsl_p)
c    1    call writeg ('MAPS_PRES ',adate,aproj,2,nx_p,ny_p,
c    1                  nzpn_p,1,1,nvarp_p,
c    1                  varnamex,varname3p,g,g1)

c       call writegg ('MAPSG_PRES ',adate,aproj,2,nx_p,ny_p,
c    1                  nzpn_p,1,1,nvarp_p,
c    1                  varnamex,varname3p,g,g1)

C   Write out RUC sfc file w/ ~25 2-d variables
        call writegg ('RUC_SFC ',adate,aproj,2,nx_p,ny_p,
     1                  nzpn_p,1,nvarsfc_ncep_p,1,
     1                  varname_sfc_ncep,varnamex,
     1                  g,g1)


        do 735 j=1,ny_p
        do 735 i=1,nx_p
          g(i,j,70) = g3(i,j,1,h_p)
c snow depth
          g(i,j,74) = gsfc(i,j,34)
C precipitable water
          g(i,j,75) = gsfc(i,j,pw_p)
C --- sfc RH
        g(i,j,76) = rhb(i,j,1)
C sub-surface runoff
        g(i,j,77) = gsfc(i,j,36)
C storm motion components
        g(i,j,78) = ustorm(i,j)
        g(i,j,79) = vstorm(i,j)
C 700mb vert vel
          gm(i,j,1) = g1(i,j,13,6)
C 250mb height/u/v
          gm(i,j,2) = g1(i,j,31,1)
          gm(i,j,3) = g1(i,j,31,4)
          gm(i,j,4) = g1(i,j,31,5)
C 500mb height/abs vort
          gm(i,j,5) = g1(i,j,21,1)
          J1 = MAX(1,J-1)
          J2 = MIN(NY_P,J+1)
          I1 = MAX(1,I-1)
          I2 = MIN(NX_P,I+1)
          k=21
          gm(i,j,6) = 0.5/SCALE(I,J)* (VP(I2,J,K)-VP(I1,J,K)
     1       -UP(I,J2,K)+UP(I,J1,K)) + corio(i,j)
C 850mb height / temp / spec hum /u/v
          gm(i,j,7) = g1(i,j,7,1)
          gm(i,j,8) = g1(i,j,7,2)
        Tx = g1(i,j,7,2)-273.15
        POL = 0.99999683       + TX*(-0.90826951E-02 +
     1     TX*(0.78736169E-04   + TX*(-0.61117958E-06 +
     2     TX*(0.43884187E-08   + TX*(-0.29883885E-10 +
     3     TX*(0.21874425E-12   + TX*(-0.17892321E-14 +
     4     TX*(0.11112018E-16   + TX*(-0.30994571E-19)))))))))
        esx = 6.1078/POL**8
c         e = g1(i,j,9,3)*esw(g1(i,j,9,2)-273.15)
          e = g1(i,j,7,3)*esx/100.

          gm(i,j,9) = 0.62197*e/(850.-0.37803*e)
          gm(i,j,10) = g1(i,j,7,4)
          gm(i,j,11) = g1(i,j,7,5)
C 1000mb height
          gm(i,j,12) = g1(i,j,1,1)
735    continue



C *** INTERPOLATE TO FREEZING LEVEL
        CALL HB2FRZ (g3(1,1,1,p_p),g3(1,1,1,h_p),TB,
     1              g3(1,1,1,theta_p),RHB,NX_P,NY_P,NZTN_P,
     1              FRZ,FRP,FRH,frz2,frp2,frh2)

        rhmax=0.
        rhmin=999.
        do i=1,nx_p
        do j=1,ny_p
          if(frh(i,j).ge.rhmax)then
             rhmax=frh(i,j)
             imax=i
             jmax=j
          endif
          if(frh(i,j).le.rhmin)then
            rhmin=frh(i,j)
            imin=i
            jmin=j
          endif
        enddo
        enddo
        print *,'RHMAX/RHMIN =',rhmax,imax,jmax,rhmin,imin,jmin
C *** INTERPOLATE TO MAXIMUM WIND LEVEL
        CALL HB2MW   (g3(1,1,1,p_p),g3(1,1,1,u_p),
     1              g3(1,1,1,v_p),NX_P,NY_P,NZTN_P,
     1              PMX,UMX,VMX)

c       if (facility_p.eq.fsl_p)
c    1   call writeg ('MAPS_SPEC ',adate,aproj,7,nx_p,ny_p,
c    1                  1,1,nvars_p,1,
c    1                  varname_speciallev,varnamex,gx,g1)

c       call writegg ('MAPSG_SPEC ',adate,aproj,7,nx_p,ny_p,
c    1                  1,1,nvars_p,1,
c    1                  varname_speciallev,varnamex,gx,g1)

c       open (17,file='frz2',form='unformatted',type='unknown')
c       write (17) frp
c       write (17) frp2
c       close (17)

        do j=1,ny_p
        do i=1,nx_p
          g(i,j,26) = frp(i,j)
          g(i,j,27) = frz(i,j)
          g(i,j,28) = frh(i,j)
          g(i,j,71) = frp2(i,j)
          g(i,j,72) = frz2(i,j)
          g(i,j,73) = frh2(i,j)
          g(i,j,29) = gsfc(i,j,tropp_p)
          g(i,j,30) = gsfc(i,j,tropth_p)
          g(i,j,88) = gsfc(i,j,tropth_p)
     1        *(gsfc(i,j,tropp_p)/100000.)**rovcp_p
          g(i,j,31) = gsfc(i,j,tropu_p)
          g(i,j,32) = gsfc(i,j,tropv_p)
          g(i,j,33) = pmx(i,j)
          g(i,j,34) = umx(i,j)
          g(i,j,35) = vmx(i,j)
C -- precipitation rate
          g(i,j,36) = gsfc(i,j,precin_p)
C -- soil moisture availability
          g(i,j,37) = gex2(i,j,5)
C -- surface runoff
          g(i,j,38) = gsfc(i,j,35)
        end do
        end do

C *** INTERPOLATE TO 30-MB LAYER MEAN FIELDS
C *** 1) 30-MB BOUNDARY LAYER
        CALL HB230MB (g3(1,1,1,p_p),TB,RHB,g3(1,1,1,u_p),
     1             g3(1,1,1,v_p),g3(1,1,1,vv_p),
     1             NX_P,NY_P,NZTN_P,1,nvarm_p,gx)

        DO 60 J=1,NY_P
        DO 60 I=1,NX_P
          gm(i,j,13) = gx(i,j,1)
        Tx = gx(i,j,1)-273.15
        POL = 0.99999683       + TX*(-0.90826951E-02 +
     1     TX*(0.78736169E-04   + TX*(-0.61117958E-06 +
     2     TX*(0.43884187E-08   + TX*(-0.29883885E-10 +
     3     TX*(0.21874425E-12   + TX*(-0.17892321E-14 +
     4     TX*(0.11112018E-16   + TX*(-0.30994571E-19)))))))))
        esx = 6.1078/POL**8

c         e = gx(i,j,2)*esw(gx(i,j,1)-273.15)
          e = gx(i,j,2)*esx/100.
C  --- 30-mb layer centered at 15mb above surface
          gm(i,j,14) = 0.62197*e/(g3(i,j,1,p_p)/100.-15.-0.37803*e)
          gm(i,j,15) = gx(i,j,3)
          gm(i,j,16) = gx(i,j,4)
60      CONTINUE
        DO 61 IVAR=1,nvarm_p
        DO 61 J=1,NY_P
        DO 61 I=1,NX_P
          G1MEAN(I,J,1,IVAR) = Gx(I,J,IVAR)
          g (i,j,38+ivar) = gx(i,j,ivar)
61      continue

        do ivar = 1,nvar_mords_p
          print *, ' MORDS field',ivar,' ',gm(60,60,ivar)
        end do

c       call writeg ('MAPS_MORDS ',adate,'98'//aproj(3:4),13,nx_p,
       if (facility_p.eq.fsl_p)
     1 call writegg ('MAPS_MORDS ',adate,'98'//aproj(3:4),13,nx_p,
     1     ny_p,1,1,nvar_mords_p,
     1     1,varname_mords_upr,varnamex,gm,g1)

C =============================
C   Calculate mean layer fields for levels 2-6
C =============================
        do 779 k=2,NZMEAN_P
C ***  30-60-MB ,60-90, 90-120, 120-150, 150-180 mb LAYERs
        CALL HB230MB (g3(1,1,1,p_p),TB,RHB,g3(1,1,1,u_p),
     1             g3(1,1,1,v_p),g3(1,1,1,vv_p),
     1             NX_P,NY_P,NZTN_P,k,nvarm_p,gx)
        DO 65 IVAR=1,nvarm_p
        DO 65 J=1,NY_P
        DO 65 I=1,NX_P
          G1MEAN(I,J,k,IVAR) = Gx(I,J,IVAR)
          g (i,j,38+(k-1)*nvarm_p+ivar) = gx(i,j,ivar)
65      CONTINUE
779     continue


c       if (facility_p.eq.fsl_p)
c    1   call writeg ('MAPS_MEAN ',adate,aproj,8,nx_p,ny_p,
c    1                  NZMEAN_P,1,1,nvarm_p,
c    1                  varnamex,varname_mean,g,g1mean)

c        call writegg ('MAPSG_MEAN ',adate,aproj,8,nx_p,ny_p,
c     1                  NZMEAN_P,1,1,nvarm_p,
c     1                  varnamex,varname_mean,g,g1mean)

C   Write out RUC isobaric main file - ~6 3-d vars and ~79 2-d vars

       if (facility_p.eq.fsl_p)
     1   call writeg ('RUC_PRESM ',adate,aproj,4,nx_p,ny_p,
     1                  nzpn_p,1,nvartot2d_p,nvarp_main_p,
     1                  varname_2d_ncep,varname3p,g,g1)

        call writegg ('RUC_PRESM ',adate,aproj,4,nx_p,ny_p,
     1                  nzp25_p,1,nvartot2d_p,nvarp_main_p,
     1                  varname_2d_ncep,varname3p,g,g1)


C   Write out isobaric supplemental file
c       call writegg ('RUC_PRESS ',adate,aproj,2,nx_p,ny_p,
c    1                  nzpn_p,1,1,nvarp_supp_p,
c    1                  varnamex,varname3p(nvarp_main_p+1),
c    1                  g,g1(1,1,1,nvarp_main_p+1) )


C--- PROFOUT is a routine for providing output profiles and
C--- surface fields in BUFR for National Weather Service and GCIP

C get ftim for BUFR
        read (aproj,'(I4.0)') NHR
        ftim=nhr*3600
        print *,' ftim= ', ftim

        CALL PROFOUT (G3,gex2,GSFC,GSFC_M,FTIM,adate,aproj)

        ISTATUS=SUCCESS_P
        RETURN
        END
