      SUBROUTINE OROSETUP
     *         ( klon   , klev  , KTEST
     *         , KKCRIT, KKCRITH, KCRIT, KSECT , KKHLIM
     *         , kkenvh, kknu  , kknu2
     *         , PAPHM1, PAPM1 , PUM1   , PVM1 , PTM1  , PGEOM1, pvaror
     *         , PRHO  , PRI   , PSTAB  , PTAU , PVPH  ,ppsi, pzdep
     *         , PULOW , PVLOW  
     *         , Ptheta, pgamma,  pnu  ,  pd1  ,  pd2  ,pdmod  )
C
C**** *GWSETUP*
C
C     PURPOSE.
C     --------
C
C**   INTERFACE.
C     ----------
C          FROM *ORODRAG*
C
C        EXPLICIT ARGUMENTS :
C        --------------------
C     ==== INPUTS ===
C     ==== OUTPUTS ===
C
C        IMPLICIT ARGUMENTS :   NONE
C        --------------------
C
C     METHOD.
C     -------
C
C
C     EXTERNALS.
C     ----------
C
C
C     REFERENCE.
C     ----------
C
C        SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
C
C     AUTHOR.
C     -------
C
C     MODIFICATIONS.
C     --------------
C     F.LOTT  FOR THE NEW-GWDRAG SCHEME NOVEMBER 1993
C
C-----------------------------------------------------------------------
      implicit none
C

#include "dimensions.h"
#include "dimphys.h"
#include "dimradmars.h"
      integer klon,klev,kidia,kfdia

#include "comcstfi.h"
#include "yoegwd.h"

C-----------------------------------------------------------------------
C
C*       0.1   ARGUMENTS
C              ---------
C
      INTEGER KKCRIT(NDLO2),KKCRITH(NDLO2),KCRIT(NDLO2),KSECT(NDLO2),
     *        KKHLIM(NDLO2),KTEST(NDLO2),kkenvh(NDLO2)

C
      REAL PAPHM1(NDLO2,KLEV+1),PAPM1(NDLO2,KLEV),PUM1(NDLO2,KLEV),
     *     PVM1(NDLO2,KLEV),PTM1(NDLO2,KLEV),PGEOM1(NDLO2,KLEV),
     *     PRHO(NDLO2,KLEV+1),PRI(NDLO2,KLEV+1),PSTAB(NDLO2,KLEV+1),
     *     PTAU(NDLO2,KLEV+1),PVPH(NDLO2,KLEV+1),ppsi(NDLO2,klev+1),
     *     pzdep(NDLO2,klev)
       REAL PULOW(NDLO2),PVLOW(NDLO2),ptheta(NDLO2),pgamma(NDLO2),
     *     pnu(NDLO2),
     *     pd1(NDLO2),pd2(NDLO2),pdmod(NDLO2)
      real pvaror(NDLO2)
C
C-----------------------------------------------------------------------
C
C*       0.2   LOCAL ARRAYS
C              ------------
C
C
      LOGICAL LL1(NDLO2,nlayermx+1)
      integer kknu(NDLO2),kknu2(NDLO2),kknub(NDLO2),kknul(NDLO2),
     *        kentp(NDLO2),ncount(NDLO2)  
C
      REAL ZHCRIT(NDLO2,nlayermx),ZNCRIT(NDLO2,nlayermx),
     *     ZVPF(NDLO2,nlayermx), ZDP(NDLO2,nlayermx)
      REAL ZNORM(NDLO2),zpsi(NDLO2),zb(NDLO2),zc(NDLO2),
     *      zulow(NDLO2),zvlow(NDLO2),znup(NDLO2),znum(NDLO2)
C
c   declarations pour "implicit none"
      integer jk,jl,ilevm1,ilevm2,ilevh
      real zu,zphi,zcons1,zcons2,zcons3,zwind,zdwind,zhgeo
      real zvt1,zvt2,zst,zvar,zdelp,zstabm,zstabp,zrhom,zrhop
      real alpha,zggeenv,zggeom1,zgvar
      logical lo

C     ------------------------------------------------------------------
C
C*         1.    INITIALIZATION
C                --------------
C
c     print *,' entree gwsetup'
 100  CONTINUE
C
C     ------------------------------------------------------------------
C
C*         1.1   COMPUTATIONAL CONSTANTS
C                -----------------------
C

      kidia=1
	  kfdia=klon

 110  CONTINUE
C
      ILEVM1=KLEV-1
      ILEVM2=KLEV-2
      ILEVH =KLEV/3
C
      ZCONS1=1./r
cold  ZCONS2=G**2/CPD
      ZCONS2=g**2/cpp
cold  ZCONS3=1.5*API
      ZCONS3=1.5*PI
C
C
C     ------------------------------------------------------------------
C
C*         2.
C                --------------
C
 200  CONTINUE
C
C     ------------------------------------------------------------------
C
C*         2.1     DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
C*                 LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
C*                 THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
C
C
C
      DO 2001 JL=kidia,kfdia
      kknu(JL)    =klev
      kknu2(JL)   =klev
      kknub(JL)   =klev
      kknul(JL)   =klev
      pgamma(JL) =max(pgamma(jl),gtsec)
      ll1(jl,klev+1)=.false.
 2001 CONTINUE
C
C*      DEFINE TOP OF LOW LEVEL FLOW
C       ----------------------------
      DO 2002 JK=KLEV,ilevh,-1
      DO 2003 JL=kidia,kfdia
      LO=(PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1)).GE.GSIGCR
      IF(LO) THEN
        KKCRIT(JL)=JK
      ENDIF
      ZHCRIT(JL,JK)=4.*pvaror(JL)
      ZHGEO=PGEOM1(JL,JK)/g
      ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
cccccc g95: XOR does not exist in Fortran  
      IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
        kknu(JL)=JK
      ENDIF
 2003 CONTINUE
 2002 CONTINUE
      DO 2004 JK=KLEV,ilevh,-1
      DO 2005 JL=kidia,kfdia
      ZHCRIT(JL,JK)=3.*pvaror(JL)
      ZHGEO=PGEOM1(JL,JK)/g
      ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
cccccc g95: XOR does not exist in Fortran
      IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
        kknu2(JL)=JK
      ENDIF
 2005 CONTINUE
 2004 CONTINUE
      DO 2006 JK=KLEV,ilevh,-1
      DO 2007 JL=kidia,kfdia
      ZHCRIT(JL,JK)=2.*pvaror(JL)
      ZHGEO=PGEOM1(JL,JK)/g
      ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
cccccc g95: XOR does not exist in Fortran
      IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
        kknub(JL)=JK
      ENDIF
 2007 CONTINUE
 2006 CONTINUE
      DO 2008 JK=KLEV,ilevh,-1
      DO 2009 JL=kidia,kfdia
      ZHCRIT(JL,JK)=pvaror(JL)
      ZHGEO=PGEOM1(JL,JK)/g
      ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
cccccc g95: XOR does not exist in Fortran
      IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
        kknul(JL)=JK
      ENDIF
 2009 CONTINUE
 2008 CONTINUE
C
      do 2010 jl=kidia,kfdia  
      kknu(jl)=min(kknu(jl),nktopg)
      kknub(jl)=min(kknub(jl),nktopg)
      if(kknub(jl).eq.nktopg) kknul(jl)=klev
C
C     CHANGE IN HERE TO STOP KKNUL=KKNUB
C
      if(kknul(jl).le.kknub(jl)) kknul(jl)=nktopg
 2010 continue      
C

 210  CONTINUE
C
C
CC*     INITIALIZE VARIOUS ARRAYS
C
      DO 2107 JL=kidia,kfdia
      PRHO(JL,KLEV+1)  =0.0
      PSTAB(JL,KLEV+1) =0.0
      PSTAB(JL,1)      =0.0
      PRI(JL,KLEV+1)   =9999.0
      ppsi(JL,KLEV+1)  =0.0
      PRI(JL,1)        =0.0
      PVPH(JL,1)       =0.0
      PULOW(JL)        =0.0
      PVLOW(JL)        =0.0
      zulow(JL)        =0.0
      zvlow(JL)        =0.0
      KKCRITH(JL)      =KLEV
      KKenvH(JL)       =KLEV
      Kentp(JL)        =KLEV
      KCRIT(JL)        =1
      ncount(JL)       =0
      ll1(JL,klev+1)   =.false.
 2107 CONTINUE
C
C*     DEFINE LOW-LEVEL FLOW
C      ---------------------
C
      DO 223 JK=KLEV,2,-1
      DO 222 JL=kidia,kfdia
      IF(KTEST(JL).EQ.1) THEN
        ZDP(JL,JK)=PAPM1(JL,JK)-PAPM1(JL,JK-1)
        PRHO(JL,JK)=2.*PAPHM1(JL,JK)*ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
        PSTAB(JL,JK)=2.*ZCONS2/(PTM1(JL,JK)+PTM1(JL,JK-1))*
     *  (1.-cpp*PRHO(JL,JK)*(PTM1(JL,JK)-PTM1(JL,JK-1))/ZDP(JL,JK))
        PSTAB(JL,JK)=MAX(PSTAB(JL,JK),GSSEC)
      ENDIF
  222 CONTINUE
  223 CONTINUE
C
C********************************************************************
C
C*     DEFINE blocked FLOW
C      -------------------
      DO 2115 JK=klev,ilevh,-1
      DO 2116 JL=kidia,kfdia
      if(jk.ge.kknub(jl).and.jk.le.kknul(jl)) then
        pulow(JL)=pulow(JL)+PUM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
        pvlow(JL)=pvlow(JL)+PVM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
      end if
 2116 CONTINUE
 2115 CONTINUE
      DO 2110 JL=kidia,kfdia
      pulow(JL)=pulow(JL)/(PAPHM1(JL,Kknul(jl)+1)-PAPHM1(JL,kknub(jl)))
      pvlow(JL)=pvlow(JL)/(PAPHM1(JL,Kknul(jl)+1)-PAPHM1(JL,kknub(jl)))
      ZNORM(JL)=MAX(SQRT(PULOW(JL)**2+PVLOW(JL)**2),GVSEC)
      PVPH(JL,KLEV+1)=ZNORM(JL)
 2110 CONTINUE
C
C*******  SETUP OROGRAPHY AXES AND DEFINE PLANE OF PROFILES  *******
C
      DO 2112 JL=kidia,kfdia
      LO=(PULOW(JL).LT.GVSEC).AND.(PULOW(JL).GE.-GVSEC)
      IF(LO) THEN
        ZU=PULOW(JL)+2.*GVSEC
      ELSE
        ZU=PULOW(JL)
      ENDIF
      Zphi=ATAN(PVLOW(JL)/ZU)
      ppsi(jl,klev+1)=ptheta(jl)*pi/180.-zphi
      zb(jl)=1.-0.18*pgamma(jl)-0.04*pgamma(jl)**2
      zc(jl)=0.48*pgamma(jl)+0.3*pgamma(jl)**2
      pd1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ppsi(jl,klev+1))**2)
      pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))*cos(ppsi(jl,klev+1))
      pdmod(jl)=sqrt(pd1(jl)**2+pd2(jl)**2)
 2112 CONTINUE
C
C  ************ DEFINE FLOW IN PLANE OF LOWLEVEL STRESS *************
C
      DO 213 JK=1,KLEV
      DO 212 JL=kidia,kfdia
      IF(KTEST(JL).EQ.1)  THEN
        ZVt1       =PULOW(JL)*PUM1(JL,JK)+PVLOW(JL)*PVM1(JL,JK)
        ZVt2       =-PvLOW(JL)*PUM1(JL,JK)+PuLOW(JL)*PVM1(JL,JK)
        ZVPF(JL,JK)=(zvt1*pd1(jl)+zvt2*pd2(JL))/(znorm(jl)*pdmod(jl))
      ENDIF
      PTAU(JL,JK)  =0.0
      Pzdep(JL,JK) =0.0
      Ppsi(JL,JK)  =0.0
      ll1(JL,JK)   =.FALSE.
  212 CONTINUE
  213 CONTINUE
      DO 215 JK=2,KLEV
      DO 214 JL=kidia,kfdia
      IF(KTEST(JL).EQ.1) THEN
        ZDP(JL,JK)=PAPM1(JL,JK)-PAPM1(JL,JK-1)
        PVPH(JL,JK)=((PAPHM1(JL,JK)-PAPM1(JL,JK-1))*ZVPF(JL,JK)+
     *            (PAPM1(JL,JK)-PAPHM1(JL,JK))*ZVPF(JL,JK-1))
     *            /ZDP(JL,JK)
        IF(PVPH(JL,JK).LT.GVSEC) THEN
          PVPH(JL,JK)=GVSEC
          KCRIT(JL)=JK
        ENDIF
      ENDIF
  214 CONTINUE
  215 CONTINUE
C
C
C*         2.2     BRUNT-VAISALA FREQUENCY AND DENSITY AT HALF LEVELS.
C
  220 CONTINUE
C
      DO 2211 JK=ilevh,KLEV
      DO 221 JL=kidia,kfdia
      IF(KTEST(JL).EQ.1) THEN
      IF(jk.ge.(kknub(jl)+1).and.jk.le.kknul(jl)) THEN
           ZST=ZCONS2/PTM1(JL,JK)*(1.-cpp*PRHO(JL,JK)*
     *                   (PTM1(JL,JK)-PTM1(JL,JK-1))/ZDP(JL,JK))
           PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV+1)+ZST*ZDP(JL,JK)
           PSTAB(JL,KLEV+1)=MAX(PSTAB(JL,KLEV+1),GSSEC)
           PRHO(JL,KLEV+1)=PRHO(JL,KLEV+1)+PAPHM1(JL,JK)*2.*ZDP(JL,JK)
     *                   *ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
      ENDIF
      ENDIF
  221 CONTINUE
 2211 CONTINUE
C
      DO 2212 JL=kidia,kfdia
C*****************************************************************************
C
C     O.K. THERE IS A POSSIBLE PROBLEM HERE. IF KKNUL=KKNUB THEN
C     DIVISION BY ZERO OCCURS. I HAVE PUT A FIX IN HERE BUT WILL ASK FRANCOIS 
C     LOTT ABOUT IT IN PARIS.
C
C     MAT COLLINS 30.1.96
C
C     ALSO IF THIS IS THE CASE PSTAB AND PRHO ARE NOT DEFINED AT KLEV+1
C     SO I HAVE ADDED THE ELSE
C
C*****************************************************************************
      IF (KKNUL(JL).NE.KKNUB(JL)) THEN
        PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV+1)/(PAPM1(JL,Kknul(jl))
     *                                          -PAPM1(JL,kknub(jl)))
        PRHO(JL,KLEV+1)=PRHO(JL,KLEV+1)/(PAPM1(JL,Kknul(jl))
     *                                          -PAPM1(JL,kknub(jl)))
      ELSE
      WRITE(*,*) 'KKNUB=KKNUL= ',KKNUB(JL),' AT JL= ',JL
        PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV)
        PRHO(JL,KLEV+1)=PRHO(JL,KLEV)
      ENDIF
        ZVAR=PVARor(JL)
 2212 CONTINUE
C
C*         2.3     MEAN FLOW RICHARDSON NUMBER.
C*                 AND CRITICAL HEIGHT FOR FROUDE LAYER
C
  230 CONTINUE
C
      DO 232 JK=2,KLEV
      DO 231 JL=kidia,kfdia
      IF(KTEST(JL).EQ.1) THEN
        ZDWIND=MAX(ABS(ZVPF(JL,JK)-ZVPF(JL,JK-1)),GVSEC)
        PRI(JL,JK)=PSTAB(JL,JK)*(ZDP(JL,JK)
     *          /(g*PRHO(JL,JK)*ZDWIND))**2
        PRI(JL,JK)=MAX(PRI(JL,JK),GRCRIT)
      ENDIF
  231 CONTINUE
  232 CONTINUE
C
C
C*      DEFINE TOP OF 'envelope' layer
C       ----------------------------

      DO 233 JL=kidia,kfdia
      pnu (jl)=0.0
      znum(jl)=0.0
 233  CONTINUE
      
      DO 234 JK=2,KLEV-1
      DO 234 JL=kidia,kfdia
      
      IF(KTEST(JL).EQ.1) THEN
       
      IF (JK.GE.KKNU2(JL)) THEN
          
            ZNUM(JL)=PNU(JL)
            ZWIND=(pulow(JL)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
     *            max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
            ZWIND=max(sqrt(zwind**2),gvsec)
            ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
            ZSTABM=SQRT(MAX(PSTAB(JL,JK  ),GSSEC))
            ZSTABP=SQRT(MAX(PSTAB(JL,JK+1),GSSEC))
            ZRHOM=PRHO(JL,JK  )
            ZRHOP=PRHO(JL,JK+1)
            PNU(JL) = PNU(JL) + (ZDELP/g)*
     *            ((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND     
            IF((ZNUM(JL).LE.GFRCRIT).AND.(PNU(JL).GT.GFRCRIT)
     *                          .AND.(KKENVH(JL).EQ.KLEV))
     *      KKENVH(JL)=JK
     
      ENDIF    

      ENDIF
      
 234  CONTINUE

C  CALCULATION OF A DYNAMICAL MIXING HEIGHT FOR THE BREAKING
C  OF GRAVITY WAVES:

              
      DO 235 JL=kidia,kfdia
      znup(jl)=0.0
      znum(jl)=0.0
 235  CONTINUE

      DO 236 JK=KLEV-1,2,-1
      DO 236 JL=kidia,kfdia
      
      IF(KTEST(JL).EQ.1) THEN
        
      IF (JK.LT.KKENVH(JL)) THEN

            ZNUM(JL)=ZNUP(JL)
            ZWIND=(pulow(JL)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
     *            max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
            ZWIND=max(sqrt(zwind**2),gvsec)
            ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
            ZSTABM=SQRT(MAX(PSTAB(JL,JK  ),GSSEC))
            ZSTABP=SQRT(MAX(PSTAB(JL,JK+1),GSSEC))
            ZRHOM=PRHO(JL,JK  )
            ZRHOP=PRHO(JL,JK+1)
            ZNUP(JL) = ZNUP(JL) + (ZDELP/g)*
     *            ((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND     
            IF((ZNUM(JL).LE.1.5).AND.(ZNUP(JL).GT.1.5)
     *                          .AND.(KKCRITH(JL).EQ.KLEV))
     *      KKCRITH(JL)=JK
     
      ENDIF
      
      ENDIF
      
 236  CONTINUE
 
      DO 237 JL=KIDIA,KFDIA
      KKCRITH(JL)=MIN0(KKCRITH(JL),KKNU(JL))
 237  CONTINUE
c
c     directional info for flow blocking ************************* 
c
      do 251 jk=ilevh,klev    
      DO 252 JL=kidia,kfdia
      IF(jk.ge.kkenvh(jl)) THEN
      LO=(PUm1(JL,jk).LT.GVSEC).AND.(PUm1(JL,jk).GE.-GVSEC)
      IF(LO) THEN
        ZU=PUm1(JL,jk)+2.*GVSEC
      ELSE
        ZU=PUm1(JL,jk)
      ENDIF
       Zphi=ATAN(PVm1(JL,jk)/ZU)
       ppsi(jl,jk)=ptheta(jl)*pi/180.-zphi
      end if
 252  continue
 251  CONTINUE
c      forms the vertical 'leakiness' **************************

      alpha=3.
      
      DO 254  JK=ilevh,klev
      DO 253  JL=kidia,kfdia
      IF(jk.ge.kkenvh(jl)) THEN
        zggeenv=AMAX1(1.,
     *          (pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)-1))/2.)      
        zggeom1=AMAX1(pgeom1(jl,jk),1.)
        zgvar=amax1(pvaror(jl)*g,1.)     
        pzdep(jl,jk)=SQRT((zggeenv-zggeom1)/(zggeom1+zgvar))      
      end if
 253  CONTINUE
 254  CONTINUE

 260  CONTINUE

      RETURN
      END
