C*******************************************************************
C
C SUBROUTINE DRAGON
C
C PURPOSE:  TO CALCULATE DRAG BASED ON 1ST GUESS SEA LEVEL WINDS, 1ST
C           GUESS GROUND TEMPERATURE, AND ANALYZED SEA LEVEL TEMP.
C
C INPUT/OUTPUT:  PHISA  - MODEL 'PHIS'   (TOPOGRAPHY*GRAVITY)
C                ZLWFLAG- LAND/WATER FLAG (0=LAND, 1=WATER)
C                USLMF  - 1ST GUESS SEA LEVEL U
C                VSLMF  - 1ST GUESS SEA LEVEL V
C                TSURF  - ANALYZED SEA LEVEL TEMPERATURE
C                GT     - 1ST GUESS GROUND TEMPERATURE
C
C ON OUTPUT :    DRAG :  CALCULATED DRAGS
C              IWATER :  INTEGER ARRAY WHICH SIGNALS IF GRID POINT IS
C                        OVER LAND OR OCEAN (0=LAND, 1=WATER)
C
C LOCAL VARIABLES:  DELT - DIFFERENCE BETWEEN GROUND AND SEA LEVEL TEMP
C                   WIND - MAGNITUDE OF SEA LEVEL WIND VELOCITY
C                   HOLD - WORK ARRAY
C                   HOLDBIT- WORK BIT ARRAY
C                   CD   - DRAG COEFFICIENT
C                   DRG  - INTERACTION COEFFICIENT
C
C NOTES:   1) THE DRAG IS CALCULATED IN THE SAME MANNER AS THAT OF
C             THE MODEL. SEE MODEL DOCUMENTATION FOR DETAILS.
C
C         *2) DRAGS ARE CALCULATED OVER OCEANS ONLY FOR THIS ANALYSIS
C             SINCE DRAGS COMPUTED OVER LAND ARE UNRELIABLE. NOTE THAT
C             THIS IMPLIES THAT 'DELT' IS A DIFFERENCE BETWEEN GROUND
C             AND SEA LEVEL TEMPERATURES WHEN OVER OCEANS.
C
C********************************************************************
C
      SUBROUTINE dragon(phisa,zlwflag,uslmf,vslmf,tsurf,tg)
       include 'params.h'
       include 'processors.h'
      DIMENSION TG(IM,JNP),TSURF(IM,JNP),VSLMF(IM,JNP),USLMF(IM, 
     &JNP), ZLWFLAG(IM,JNP),PHISA(IM,JNP)
!hpf$ distribute phisa(*,block) onto PROC_ARRAY_1D
!hpf$ distribute tg(*,block) onto PROC_ARRAY_1D
!hpf$ distribute tsurf(*,block) onto PROC_ARRAY_1D
!hpf$ distribute uslmf(*,block) onto PROC_ARRAY_1D
!hpf$ distribute vslmf(*,block) onto PROC_ARRAY_1D
!hpf$ distribute zlwflag(*,block) onto PROC_ARRAY_1D
      integer iw(imr,jmr+1)
      COMMON / sfcdrag / drag (imr,jmr+1), iwater (13104)
!hpf$ distribute drag(*,block) onto PROC_ARRAY_1D
!hpf$ distribute iwater(block) onto PROC_ARRAY_1D

      COMMON / CONTRL / IM, JNP, NB, NDALT, NHMS, NLAY, NLAYM1, NLAYP1, 
     &NYMD, STRAT, IDIMEN, SLP, FDCOS, FEXP, SLT, SIGDAT, MAXREP
      LOGICAL STRAT, SLP, FDCOS, FEXP, SLT, SIGDAT
      COMMON / PHYSIC / BETA, CP, GRAV, OMEGA2, RADE, RGAS, ROCP, CPOG, 
     &GAMMA, ROCP1, GOCP
      REAL  DELT(im,jnp)
      REAL  CD(im,jnp)
      REAL  WIND(im,jnp)
      REAL  DRG(im,jnp)
!hpf$ align with phisa:: iw,DELT,CD,WIND,DRG
      DATA ZERO / 0.0 /
      DATA ONE / 1.0 /
      DATA TWO / 2.0 /
      DATA FIVE / 5.0 /
      DATA SEVEN / 7.0 /
      DATA EIGHTEEN / 18.0 /
      DATA HAFF / 0.5 /
      DATA TENM3 / .001 /
      DATA PT00125 / .00125 /
      DATA PT0025 / .0025 /
      DATA PT07 / .07 /
C----------------------------------------------------------------------
C----------------------------------------------------------------------
      write(6,*) 'jnp= ',jnp
c      iwater=RESHAPE(zlwflag + 0.001,(/im*jnp/))
      do j=1,jnp
        forall(i=1:im) iwater(i+(j-1)*im)=zlwflag(i,j)+0.001
      end do
      iw=zlwflag + 0.001
      PRINT 6020
6020  FORMAT (' IWATER FLAG IN DRAGON:')
      k = im * jnp
6021  FORMAT (' ',72I1)
c      DO j = 1 , jnp
c        PRINT 6021, (iwater(i),i = -(j - 1) * im + k - im + 1,-(j - 1) 
c     &      * im + k)
c        k = k - im
c      ENDDO
      delt = tg - tsurf
      wind = sqrt(uslmf* uslmf + vslmf * vslmf)
      WHERE (wind .LT. tenm3)
        wind = zero
      ENDWHERE
      cd = (one + pt07 * wind) * tenm3

C  TEMPERATURE DIFFERENCE BETWEEN GROUND AND SURFACE
C  -------------------------------------------------
C  WIND SPEED CALCULATIONS
C  ----------------------
C  CALCULATE DRAG COEFFICIENT CD
C  -----------------------------
C  CALCULATE STABILITY-DEPENDENT AIR-SURFACE INTERACTION COEFFICIENT
C  -----------------------------------------------------------------
C FOR SEA LEVEL PRESSURE ANALYSIS, CALCULATE DRAG OVER OCEANS ONLY
C ----------------------------------------------------------
      where(cd.lt.pt0025)
        cd=cd*haff
      elsewhere
        cd=pt00125
      endwhere
      where(delt.ge.zero)
        drg = cd * (wind + sqrt(delt))
      elsewhere
        drg = cd * wind*wind*wind/(wind* wind - seven * delt)
      endwhere
      where(iw.eq.1)
       drag=eighteen*grav*drg/(rgas*tsurf)
      elsewhere
       drag=zero
      endwhere

      RETURN
      END
