MODULE gwstress_mod
      
IMPLICIT NONE
      
CONTAINS

      SUBROUTINE GWSTRESS(ngrid,nlayer,ktest,zrho, BV, pvar,psig,zgeom,zdmod,& 
                ! Notice that this 3 variables are actually not used due to illness lines
                 ICRIT,IKNU,ZVPH,                                            &
                ! not used variables
 !                IKCRIT,ISECT,IKHLIM, IKCRITH,IKENVH,PVAR1,pgam,zd1,zd2,znu, &
                ! not defined not used variables
 !                 ZTFR 
                !in(as 0.0)-output:
                ZTAU )
                
      !----------------------------------------------------------------------------------------------
      ! MODULE contains SUBROUTINE gwstress to compute low level stresses using subcritical, super
      ! critical forms.
      ! F. LOTT PUT THE NEW GWD ON IFS      22/11/93
      ! REFERENCE.
      ! SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
      ! Rewirten by J.Liu 03/03/2022
      !----------------------------------------------------------------------------------------------

      use dimradmars_mod, only: ndomainsz   
      use yoegwd_h, only: gkdrag, gtsec, gvcrit

      implicit none
      
      ! 0. DECLARATIONS:
      
      ! 0.1   ARGUMENTS 
      integer,intent(in):: ngrid    ! number of atmospheric columns
      integer,intent(in):: nlayer   ! number of atmospheric layers
      INTEGER,intent(in):: ktest(ndomainsz)   ! map of calling points
!      integer,intent(in):: IKCRIT(ndomainsz) ! not used 
      integer,intent(in):: ICRIT(ndomainsz)   ! actually not used 
 !     integer,intent(in):: IKCRITH(ndomainsz)! not used 
 !     integer,intent(in):: ISECT(ndomainsz)  ! not used
 !     integer,intent(in):: IKHLIM(ndomainsz) ! not used 
 !     integer,intent(in):: IKENVH(ndomainsz) ! The line use this variable has been commented
      integer,intent(in):: IKNU(ndomainsz)    ! actually not used 
      REAL,INTENT(IN):: ZRHO(ndomainsz,nlayer+1) ! Density at 1/2 level
      REAL,INTENT(IN):: BV(ndomainsz,nlayer+1)   ! Brunt–Väisälä frequency at 1/2 level
      REAL,INTENT(IN):: ZVPH(ndomainsz,nlayer+1) ! Low level wind speed U_H
      REAL,INTENT(IN):: zgeom(ndomainsz,nlayer)  ! Geopotetial height
      REAL,INTENT(IN):: pvar(ndomainsz)          ! Sub-grid scale standard deviation
!      REAL,INTENT(IN):: zd1(ndomainsz)       ! not used
!      REAL,INTENT(IN):: zd2(ndomainsz)       ! not used
!      REAL,INTENT(IN):: znu(ndomainsz)       ! not used
      REAL,INTENT(IN):: psig(ndomainsz)       ! SUB-GRID SCALE SLOPE
!      REAL,INTENT(IN):: pgam(ndomainsz)      ! not used
      REAL,INTENT(IN):: zdmod(ndomainsz)      ! Squre root of tao1 and tao2 without the constant, see equation 17 or 18
!      REAL,INTENT(IN):: ZTFR(ndomainsz)      ! not used. It is not even defined in this rountine
      REAL,INTENT(INOUT):: ZTAU(ndomainsz,nlayer+1) !GRAVITY WAVE STRESS.
      
      !0.2   LOCAL ARRAYS
      integer jl
      INTEGER kidia,kfdia
      real zvar    ! Sub-grid scale standard deviation at the calling points
      real zblock,zeff
      logical lo   ! actually not used bucause the if-endif condition that use this
                   ! variable has been commented
                   
!---------------------------------------------------------------------------------------------------
! 1. INITIALIZATION (not important initialization at all may be delete in the future)
!---------------------------------------------------------------------------------------------------
      kidia=1
      kfdia=ngrid
 100   CONTINUE ! continue tag without source, maybe need delete in future
!*         3.1     Gravity wave stress
 300   CONTINUE ! continue tag without source, maybe need delete in future

      DO JL=kidia,kfdia
      	IF(KTEST(JL).EQ.1) THEN
        !Effective mountain height above the blocked flow     
!        IF(IKENVH(JL).EQ.nlayer)THEN
         ZBLOCK=0.0 
!        ELSE
!         ZBLOCK=(zgeom(JL,IKENVH(JL))+zgeom(JL,IKENVH(JL)+1))/2./RG          
!        ENDIF    
        ZVAR=pvar(JL)
        ZEFF=AMAX1(0.,2.*ZVAR-ZBLOCK)
        ! Evaluate equation 17 to get the GW stress
        ZTAU(JL,nlayer+1)=zrho(JL,nlayer+1)*GKDRAG*psig(jl)*ZEFF**2    &
         /4./ZVAR*ZVPH(JL,nlayer+1)*zdmod(jl)*sqrt(BV(jl,nlayer+1))

      !  Too small value of stress or low level flow include critical level 
      !  or low level flow: gravity wave stress nul.                 
        LO=(ZTAU(JL,nlayer+1).LT.GTSEC).OR.(ICRIT(JL).GE.IKNU(JL)).OR. &
        (ZVPH(JL,nlayer+1).LT.GVCRIT)
!       IF(LO) ZTAU(JL,nlayer+1)=0.0      
      	ELSE     
          ZTAU(JL,nlayer+1)=0.0          
      	ENDIF      
      ENDDO

      END SUBROUTINE GWSTRESS
      
END MODULE gwstress_mod