MODULE gwprofil_mod IMPLICIT NONE CONTAINS SUBROUTINE GWPROFIL( ngrid, nlayer, kgwd ,kdx, ktest, & IKCRITH, ICRIT, IKENVH, IKNU,IKNU2,pplev, & ZRHO, BV, ZVPH, PRI,zdmod, psig , pvar, & !not used variables !IKCRIT,ZTAUf,ZTFR,znu,pgam, ! in-output ZTAU ) !--------------------------------------------------------------------------- ! THe stress profile for gravity waves is computed as follows: ! It is constant (no GWD) at the levels between the ground and ! the top of the blocked layer (IKENVH). It decreases linearly ! with heights from the top of the blocked layer to 3*varor(IKNU) ! to simulates lee waves or nonlinear gravity wave breaking. ! Above it is constant, except when the wave encounters a critical ! level(ICRIT) or when it breaks.! ! PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93) ! REFERENCE. ! SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S." ! MODIFICATIONS. ! Rewrited by J.Liu 03/03/2022 !--------------------------------------------------------------------------- use dimradmars_mod, only: ndomainsz use yoegwd_h, only: gkdrag, grcrit, gssec, gtsec 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) :: kgwd ! Number of points at which the scheme is called ! INTEGER, INTENT(IN) :: IKCRIT(ndomainsz ) ! Top of low level flow height (not used) INTEGER, INTENT(IN) :: IKCRITH(ndomainsz ) ! dynamical mixing height for the breaking of gravity waves INTEGER, INTENT(IN) :: ICRIT(ndomainsz ) ! Critical layer where orographic GW breaks INTEGER, INTENT(IN) :: kdx(ndomainsz ) ! Points at which to call the scheme INTEGER, INTENT(IN) :: ktest(ndomainsz ) ! Map of calling points INTEGER, INTENT(IN) :: IKENVH(ndomainsz ) ! Top of the blocked layer INTEGER, INTENT(IN) :: IKNU(ndomainsz ) ! 4*pvar layer INTEGER, INTENT(IN) :: IKNU2(ndomainsz ) ! 3*pvar layer REAL, INTENT(IN):: pplev(ndomainsz ,nlayer+1) ! Pressure at 1/2 levels(Pa) REAL, INTENT(IN):: BV(ndomainsz ,nlayer+1) ! Brunt–Väisälä frequency at 1/2 level (output from OROSETUP input for GWSTRESS and GWPROFIL) REAL, INTENT(IN):: ZRHO (ndomainsz ,nlayer+1) ! Atmospheric density at 1/2 level REAL, INTENT(IN):: ZVPH (ndomainsz ,nlayer+1) ! SUB-GRID SCALE STANDARD DEVIATION at 1/2 level REAL, INTENT(IN):: PRI (ndomainsz ,nlayer+1) ! Mean flow richardson number at 1/2 level ! REAL, INTENT(IN):: ZTFR (ndomainsz ) ! not used ! REAL, INTENT(IN):: ZTAUf (ndomainsz ,nlayer+1) ! not used ! REAL, INTENT(IN):: zdmod (ndomainsz ) ! REAL, INTENT(IN):: znu (ndomainsz ) ! not used REAL, INTENT(IN):: psig(ndomainsz ) ! Sub-grid scale slope ! REAL, INTENT(IN):: pgam(ndomainsz ) ! Sub-grid scale anisotropy(not used) REAL, INTENT(IN):: pvar(ndomainsz ) ! Sub-grid scale standard deviation REAL, INTENT(inout):: ZTAU(ndomainsz ,nlayer+1) ! Gravity waves induced stress ! 0.2 LOCAL ARRAYS real zsqr,zalfa,zriw,zalpha,zb,zdel,zdz2n,zdelp,zdelpt integer ji,jk,jl,ilevh REAL ZDZ2 (ndomainsz ,nlayer) , ZNORM(ndomainsz ) , zoro(ndomainsz ) REAL Z_TAU (ndomainsz ,nlayer+1) !----------------------------------------------------------------------- !* 1. Initialization !----------------------------------------------------------------------- ! kidia=1 ! kfdia=ngrid 100 CONTINUE ! 1.1 Computional constants . ilevh=nlayer/3 DO ji=1,kgwd jl=kdx(ji) Zoro(JL)=psig(JL)*zdmod(JL)/4./max(pvar(jl),1.0) Z_TAU(JL,IKNU(JL)+1)=ZTAU(JL,IKNU(JL)+1) Z_TAU(JL,nlayer+1)=ZTAU(JL,nlayer+1) ENDDO ! all start here with this long loop DO JK=nlayer,2,-1 ! 4.1 Constant wave stress until top of the blocking layer 410 CONTINUE DO ji=1,kgwd jl=kdx(ji) IF(JK.GE.IKNU2(JL)) THEN ZTAU(JL,JK)=Z_TAU(JL,nlayer+1) ENDIF ENDDO !* 4.15 Constant shear stress until the top of the low level flow layer 415 CONTINUE ! 4.2 Wave displacement at next level. 420 CONTINUE DO ji=1,kgwd jl=kdx(ji) IF(JK.LT.IKNU2(JL)) THEN ZNORM(JL)=gkdrag*ZRHO(JL,JK)*SQRT(BV(JL,JK))*ZVPH(JL,JK)*zoro(jl) ZDZ2(JL,JK)=ZTAU(JL,JK+1)/max(ZNORM(JL),gssec) ENDIF ENDDO ! 4.3 Wave Richardson number, new wave displacement and stress: ! Breaking evaluation and critical level C DO ji=1,kgwd jl=kdx(ji) IF(JK.LT.IKNU2(JL)) THEN IF((ZTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.ICRIT(JL))) THEN ZTAU(JL,JK)=0.0 ELSE ZSQR=SQRT(PRI(JL,JK)) ZALFA=SQRT(BV(JL,JK)*ZDZ2(JL,JK))/ZVPH(JL,JK) ZRIW=PRI(JL,JK)*(1.-ZALFA)/(1+ZALFA*ZSQR)**2 IF(ZRIW.LT.GRCRIT) THEN ZDEL=4./ZSQR/GRCRIT+1./GRCRIT**2+4./GRCRIT ZB=1./GRCRIT+2./ZSQR ZALPHA=0.5*(-ZB+SQRT(ZDEL)) ZDZ2N=(ZVPH(JL,JK)*ZALPHA)**2/BV(JL,JK) ZTAU(JL,JK)=ZNORM(JL)*ZDZ2N ELSE ZTAU(JL,JK)=ZNORM(JL)*ZDZ2(JL,JK) ENDIF ZTAU(JL,JK)=MIN(ZTAU(JL,JK),ZTAU(JL,JK+1)) ENDIF ENDIF ENDDO end DO !JK=nlayer,2,-1 !------------------------------------------------------------------------------------------- 440 CONTINUE ! write(*,*) 'ZTAU' ! write(*,99) ((ji,ilevh,ZTAU(ji,ilevh),ji=1,ndomainsz ), ! . ilevh=1,nlayer+1) 99 FORMAT(i3,i3,f15.5) ! not used ! Proganisation of the stress profile if breaking occurs at low level DO ji=1,kgwd jl=kdx(ji) Z_TAU(JL,IKENVH(JL))=ZTAU(JL,IKENVH(JL)) Z_TAU(JL,IKCRITH(JL))=ZTAU(JL,IKCRITH(JL)) ENDDO DO JK=1,nlayer DO ji=1,kgwd jl=kdx(ji) IF(JK.GT.IKCRITH(JL).AND.JK.LT.IKENVH(JL))THEN ZDELP=pplev(JL,JK)-pplev(JL,IKENVH(JL)) ZDELPT=pplev(JL,IKCRITH(JL))-pplev(JL,IKENVH(JL)) ZTAU(JL,JK)=Z_TAU(JL,IKENVH(JL))+(Z_TAU(JL,IKCRITH(JL))-Z_TAU(JL,IKENVH(JL)))*ZDELP/ZDELPT ENDIF ENDDO end DO END SUBROUTINE GWPROFIL END MODULE gwprofil_mod