SUBROUTINE OROSETUP( ngrid, nlayer, ktest, pplev, pplay, pu, pv, pt, zgeom, & pvar,pthe, pgam, & !output in capital IKCRIT, IKCRITH, ICRIT, IKENVH,IKNU,IKNU2, & !ISECT, IKHLIM, not used ZRHO,PRI,BV,ZTAU,ZVPH,ZPSI,ZZDEP,ZNU,ZD1,ZD2,ZDMOD, & PULOW, PVLOW) !--------------------------------------------------------------------------------------------- !!**** *GWSETUP*! Computes low level stresses using subcritical and super critical forms. ! As well as, computes anisotropy coefficient as measure of orographic two-dimensionality ! F.LOTT FOR THE NEW-GWDRAG SCHEME NOVEMBER 1993! !-- ! REFERENCE. ! 1. SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S." ! 2. Lott, F., & Miller, M. J. (1997). A new subgrid‐scale orographic drag parametrization: ! Its formulation and testing.Quarterly Journal of the Royal cMeteorological Society, ! 123(537), 101-127. !-- ! MODIFICATIONS. ! 1.Rewiten by J.liu 03/03/2022 !----------------------------------------------------------------------- use dimradmars_mod, only: ndomainsz use comcstfi_h, only: cpp, g, r, pi use yoegwd_h, only: gfrcrit, grcrit, gsigcr, gssec, gtsec, gvsec use yoegwd_h, only: nktopg implicit none ! 0. DECLARATIONS: ! 0.1 inputs: 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 real, intent(in) :: pplev(ndomainsz,nlayer+1)! Pressure at 1/2 levels(Pa) (has been inversed by DRAG_NORO=inv_pplev) real, intent(in) :: pplay(ndomainsz,nlayer) ! Pressure at full levels(Pa) (has been inversed by DRAG_NORO=inv_pplay) real, intent(in) :: pu(ndomainsz,nlayer) ! Zonal wind at full levels(m/s) (has been inversed by DRAG_NORO, =inv_pu) real, intent(in) :: pv(ndomainsz,nlayer) ! Meridional winds at full levels(m/s)(has been inversed by DRAG_NORO, =inv_pv) real, intent(in) :: pt(ndomainsz,nlayer) ! Temperature at full levels(m/s) (has been inversed by DRAG_NORO=inv_pt) real, intent(in) :: zgeom(ndomainsz,nlayer) ! Geopotetial height real, intent(in) :: pvar(ndomainsz) ! Sub-grid scale standard deviation real, intent(in) :: pthe(ndomainsz) ! Sub-grid scale principal axes angle real, intent(inout) :: pgam(ndomainsz) ! Sub-grid scale anisotropy ! 0.2 outputs: INTEGER,intent(out):: IKCRIT(ndomainsz) ! top of low level flow height INTEGER,intent(out):: IKCRITH(ndomainsz) ! dynamical mixing height for the breaking of gravity waves INTEGER,intent(out):: ICRIT(ndomainsz) ! Critical layer where orographic GW breaks ! INTEGER,intent(out):: ISECT(ndomainsz) ! not used ! INTEGER,intent(out):: IKHLIM(ndomainsz) ! not used INTEGER,intent(out):: IKENVH(ndomainsz) ! Top of the blocked layer INTEGER,intent(out):: IKNU(ndomainsz) ! 4*pvar layer INTEGER,intent(out):: IKNU2(ndomainsz) ! 3*pvar layer REAL, intent(out):: ZRHO(ndomainsz,nlayer+1) ! Density at 1/2 level REAL, intent(out):: PRI(ndomainsz,nlayer+1) ! Mean flow richardson number at 1/2 level REAL, intent(out):: BV(ndomainsz,nlayer+1) ! Brunt–Väisälä frequency at 1/2 level REAL, intent(out):: ZTAU(ndomainsz,nlayer+1) ! Gravity wave stress. Set to 0.0 here and will calculate in GWSTRESS later REAL, intent(out):: ZVPH(ndomainsz,nlayer+1) ! Low level wind speed U_H REAL, intent(out):: ZPSI(ndomainsz,nlayer+1) ! The angle between the incident flow direction and the normal ridge direction pthe REAL, intent(out):: ZZDEP(ndomainsz,nlayer) ! dp by full level REAL, intent(out):: PULOW(ndomainsz) ! Low level zonal wind REAL, intent(out):: PVLOW(ndomainsz) ! Low level meridional wind REAL, intent(out):: ZNU(ndomainsz) ! A critical value see equation 9 REAL, intent(out):: ZD1(ndomainsz) ! Bcos^2(psi)-Csin^2(psi) see equation 17 or 18 REAL, intent(out):: ZD2(ndomainsz) ! (B-C)sin(psi)cos(psi) see equation 17 or 18 REAL, intent(out):: ZDMOD(ndomainsz) ! sqrt(zd1^2+zd2^2) !0.3 Local arrays integer IKNUb(ndomainsz) ! 2*pvar layer integer IKNUl(ndomainsz) ! 1*pvar layer integer kentp(ndomainsz) ! initialized value but never used integer ncount(ndomainsz) ! initialized value but never used REAL ZHCRIT(ndomainsz,nlayer) ! tag for 1*pvar, 2*pvar,3*pvar and 4*pvar, pvar is mu means SD ! REAL ZNCRIT(ndomainsz,nlayer) ! not used REAL ZVPF(ndomainsz,nlayer) ! Flow in plane of low level stress. Seems a unit vector REAL ZDP(ndomainsz,nlayer) ! dp differitial of pressure REAL ZNORM(ndomainsz) ! The norm ridge of a moutain? REAL zb(ndomainsz) ! Parameter B in eqution 17 REAL zc(ndomainsz) ! Parameter C in eqution 17 REAL zulow(ndomainsz) ! initialized value but never used REAL zvlow(ndomainsz) ! initialized value but never used REAL znup(ndomainsz) ! znu in top of 1/2 level REAL znum(ndomainsz) ! znu in bottom of 1/2 level integer jk,jl integer ilevm1 !=nlayer-1 integer ilevm2 !=nlayer-2 integer ilevh !=nalyer/3 INTEGER kidia !=1 INTEGER kfdia !=ngrid real zcons1 !=1/r real zcons2 !=g^2/cpp real zcons3 !=1.5*pi real zphi ! direction of the incident flow real zhgeo ! Height calculated by geopotential/g real zu ! Low level zonal wind (to denfine the dirction of background wind) real zwind,zdwind real zvt1,zvt2,zst,zvar real zdelp !dp differitial of pressure ! variables for bv and density rho real zstabm,zstabp,zrhom,zrhop ! real alpha !=3. but never used real zggeenv,zggeom1,zgvar logical lo LOGICAL LL1(ndomainsz,nlayer+1) !-------------------------------------------------------------------------------- ! 1. INITIALIZATION !-------------------------------------------------------------------------------- ! 100 CONTINUE ! continue tag without source, maybe need delete in future !* 1.1 COMPUTATIONAL CONSTANTS kidia=1 kfdia=ngrid ! 110 CONTINUE ! continue tag without source, maybe need delete in future ILEVM1=nlayer-1 ILEVM2=nlayer-2 ILEVH =nlayer/3 !!!! maybe not enough for Mars, need check later ZCONS1=1./r ZCONS2=g**2/cpp ZCONS3=1.5*PI !------------------------------------------------------------------------------------------------------ ! 2. Compute all the critical levels and the coeffecients of anisotropy !----------------------------------------------------------------------------------------------------- ! 200 CONTINUE ! continue tag without source, maybe need delete in future ! 2.1 Define low level wind, project winds in plane of low level wind, ! determine sector in which to take the variance and set indicator for critical levels. DO JL=kidia,kfdia ! initialize all the height into surface (notice the layers have been inversed by preious rountines) IKNU(JL) =nlayer IKNU2(JL) =nlayer IKNUb(JL) =nlayer IKNUl(JL) =nlayer pgam(JL) =max(pgam(jl),gtsec) ! gtsec is from yoegwd.h which is a common variable ll1(jl,nlayer+1)=.false. end DO ! Define top of low level flow (since pressure, zonal and meridional wind have been inversed ! the process is to find the layer from surface (nlayer) to some levels ) by searching several ! altitude scope ! using 4 times sub-grid scale deviation as the threahold of the critical height DO JK=nlayer,ilevh,-1 ! ilevh=nlayer/3=16 DO JL=kidia,kfdia ! jl=1:ngrid ! To found the layer of the "top of low level flow" LO=(pplev(JL,JK)/pplev(JL,nlayer+1)).GE.GSIGCR IF(LO) THEN IKCRIT(JL)=JK ENDIF ZHCRIT(JL,JK)=4.*pvar(JL) ! ! use geopotetial denfination to get geoheight[in meters] of the layer ZHGEO=zgeom(JL,JK)/g ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK)) IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN IKNU(JL)=JK ENDIF ENDDO end DO ! using 3 times sub-grid scale deviation as the threahold of the critical height DO JK=nlayer,ilevh,-1 DO JL=kidia,kfdia ZHCRIT(JL,JK)=3.*pvar(JL) ZHGEO=zgeom(JL,JK)/g ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK)) IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN IKNU2(JL)=JK ENDIF ENDDO end DO ! using 2 times sub-grid scale deviation as the threahold of the critical height DO JK=nlayer,ilevh,-1 DO JL=kidia,kfdia ZHCRIT(JL,JK)=2.*pvar(JL) ZHGEO=zgeom(JL,JK)/g ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK)) IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN IKNUb(JL)=JK ENDIF ENDDO end DO ! using 1 times sub-grid scale deviation as the threahold of the critical height DO JK=nlayer,ilevh,-1 DO JL=kidia,kfdia ZHCRIT(JL,JK)=pvar(JL) ZHGEO=zgeom(JL,JK)/g ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK)) IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN IKNUl(JL)=JK ENDIF ENDDO end DO ! loop to relocate the critical height to make sure everything is okay if theses ! levels hit the model surface or top. do jl=kidia,kfdia IKNU(jl)=min(IKNU(jl),nktopg) ! nktopg is a common variable from yoegwd.h IKNUb(jl)=min(IKNUb(jl),nktopg) if(IKNUb(jl).eq.nktopg) IKNUl(jl)=nlayer ! Change in here to stop IKNUl=IKNUb if(IKNUl(jl).le.IKNUb(jl)) IKNUl(jl)=nktopg enddo ! 210 CONTINUE ! continue tag without source, maybe need delete in future ! Initialize various arrays for the following computes DO JL=kidia,kfdia ZRHO(JL,nlayer+1) =0.0 BV(JL,nlayer+1) =0.0 BV(JL,1) =0.0 PRI(JL,nlayer+1) =9999.0 ZPSI(JL,nlayer+1) =0.0 PRI(JL,1) =0.0 ZVPH(JL,1) =0.0 PULOW(JL) =0.0 PVLOW(JL) =0.0 zulow(JL) =0.0 zvlow(JL) =0.0 IKCRITH(JL) =nlayer ! surface IKENVH(JL) =nlayer ! surface Kentp(JL) =nlayer ! surface ICRIT(JL) =1 ! topmost layer ncount(JL) =0 ll1(JL,nlayer+1) =.false. ENDDO ! Define low-level flow Brunt–Väisälä frequency N^2, density ZRHO ! The incident flow passes over the mean orography is evaluated by averaging the wind, ! the Brunt–Väisälä frequency, and the fluid density between 1*pvar and 2*pvar over the ! model mean orography DO JK=nlayer,2,-1 ! from surface to topmost-1 layer DO JL=kidia,kfdia IF(ktest(JL).EQ.1) THEN ! if the map of the calling points is true ! calcalate density and BV ZDP(JL,JK)=pplay(JL,JK)-pplay(JL,JK-1) !dp>0 ZRHO(JL,JK)=2.*pplev(JL,JK)*ZCONS1/(pt(JL,JK)+pt(JL,JK-1)) !rho=p/(r*T) ! Brunt–Väisälä frequency N^2. This equation for BV is illness since ! too many variables are used. Use N^2=g/T[1/(cpp*T)+dT/dz] to replace in the future BV(JL,JK)=2.*ZCONS2/(pt(JL,JK)+pt(JL,JK-1))*(1.-cpp*ZRHO(JL,JK)*(pt(JL,JK)-pt(JL,JK-1))/ZDP(JL,JK)) BV(JL,JK)=MAX(BV(JL,JK),GSSEC) ENDIF ENDDO end DO DO JK=nlayer,ilevh,-1 DO JL=kidia,kfdia if(jk.ge.IKNUb(jl).and.jk.le.IKNUl(jl)) then ! IF the layer between 1*pvar and 2*pvar ! calculate the low level wind U_H ! pulow/pvlow at a speicfic location equals to sum of u*dp of all levels ! notice here dp is already a positive number pulow(JL)=pulow(JL)+pu(JL,JK)*(pplev(JL,JK+1)-pplev(JL,JK)) pvlow(JL)=pvlow(JL)+pv(JL,JK)*(pplev(JL,JK+1)-pplev(JL,JK)) end if ENDDO end DO ! averaging the wind DO JL=kidia,kfdia ! by divide dp [p differ between iknul and uknub level] pulow(JL)=pulow(JL)/(pplev(JL,IKNUl(jl)+1)-pplev(JL,IKNUb(jl))) pvlow(JL)=pvlow(JL)/(pplev(JL,IKNUl(jl)+1)-pplev(JL,IKNUb(jl))) ! average U to get background U? ZNORM(JL)=MAX(SQRT(PULOW(JL)**2+PVLOW(JL)**2),GVSEC) ZVPH(JL,nlayer+1)=ZNORM(JL) ! The wind below the surface level (e.g., start of the 1/2 level) ENDDO ! The gravity wave drag caused by the flow passes over an single elliptic mountain can be calculated ! by equation 17 and 18 DO 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 ! Here all physics for equation 17 and 18 ! Direction of the incident flow Zphi=ATAN(PVLOW(JL)/ZU) ! The angle between the incident flow direction and the normal ridge direction pthe ZPSI(jl,nlayer+1)=pthe(jl)*pi/180.-zphi ! equation(17) parameter B and C zb(jl)=1.-0.18*pgam(jl)-0.04*pgam(jl)**2 zc(jl)=0.48*pgam(jl)+0.3*pgam(jl)**2 ! Bcos^2(psi)-Csin^2(psi) ZD1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ZPSI(jl,nlayer+1))**2) ! (B-C)sin(psi)cos(psi) ZD2(jl)=(zb(jl)-zc(jl))*sin(ZPSI(jl,nlayer+1))*cos(ZPSI(jl,nlayer+1)) ! squre root of tao1 and tao2 without the constant see equation 17 or 18 ZDMOD(jl)=sqrt(ZD1(jl)**2+ZD2(jl)**2) ENDDO ! Define blocked flow ! Setup orogrphy axes and define plane of profiles ! Define blocked flow in plane of the low level stress DO JK=1,nlayer DO JL=kidia,kfdia IF(ktest(JL).EQ.1) THEN ZVt1 =PULOW(JL)*pu(JL,JK)+PVLOW(JL)*pv(JL,JK) ZVt2 =-PvLOW(JL)*pu(JL,JK)+PuLOW(JL)*pv(JL,JK) ! zvpf is a normalized variable ZVPF(JL,JK)=(zvt1*ZD1(jl)+zvt2*ZD2(JL))/(znorm(jl)*ZDMOD(jl)) ENDIF ZTAU(JL,JK) =0.0 ZZDEP(JL,JK) =0.0 ZPSI(JL,JK) =0.0 ll1(JL,JK) =.FALSE. ENDDO end DO DO JK=2,nlayer DO JL=kidia,kfdia IF(ktest(JL).EQ.1) THEN ZDP(JL,JK)=pplay(JL,JK)-pplay(JL,JK-1) ! dp ! zvph is the U_H in equation 17 e.g. low level wind speed ZVPH(JL,JK)=((pplev(JL,JK)-pplay(JL,JK-1))*ZVPF(JL,JK)+ & (pplay(JL,JK)-pplev(JL,JK))*ZVPF(JL,JK-1))/ZDP(JL,JK) IF(ZVPH(JL,JK).LT.GVSEC) THEN ZVPH(JL,JK)=GVSEC ICRIT(JL)=JK ENDIF endIF ENDDO end DO ! 2.2 Brunt-vaisala frequency and density at half levels 220 CONTINUE ! continue tag without source, maybe need delete in future DO JK=ilevh,nlayer DO JL=kidia,kfdia IF(ktest(JL).EQ.1) THEN IF(jk.ge.(IKNUb(jl)+1).and.jk.le.IKNUl(jl)) THEN ZST=ZCONS2/pt(JL,JK)*(1.-cpp*ZRHO(JL,JK)* & (pt(JL,JK)-pt(JL,JK-1))/ZDP(JL,JK)) BV(JL,nlayer+1)=BV(JL,nlayer+1)+ZST*ZDP(JL,JK) BV(JL,nlayer+1)=MAX(BV(JL,nlayer+1),GSSEC) ZRHO(JL,nlayer+1)=ZRHO(JL,nlayer+1)+pplev(JL,JK)*2.*ZDP(JL,JK) & *ZCONS1/(pt(JL,JK)+pt(JL,JK-1)) ENDIF endIF ENDDO end DO DO JL=kidia,kfdia !***************************************************************************** ! Okay. There is a possible problem here. If IKNUl=IKNUb then division by zero ! occurs. I have put a fix in here but will ask Francois lott about it in Paris. ! Also if this is the case BV and ZRHO are not defined at nlayer+1 so I have ! added the else. ! by: MAT COLLINS 30.1.96 !***************************************************************************** IF (IKNUL(JL).NE.IKNUB(JL)) THEN BV(JL,nlayer+1)=BV(JL,nlayer+1)/(pplay(JL,IKNUl(jl))-pplay(JL,IKNUb(jl))) ZRHO(JL,nlayer+1)=ZRHO(JL,nlayer+1)/(pplay(JL,IKNUl(jl))-pplay(JL,IKNUb(jl))) ELSE WRITE(*,*) 'OROSETUP: IKNUB=IKNUL= ',IKNUB(JL),' AT JL= ',JL BV(JL,nlayer+1)=BV(JL,nlayer) ZRHO(JL,nlayer+1)=ZRHO(JL,nlayer) ENDIF ZVAR=pvar(JL) ENDDO !JL=kidia,kfdia ! 2.3 Mean flow richardson number and critical height for proude layer ! 230 CONTINUE ! continue tag without source, maybe need delete in future DO JK=2,nlayer DO JL=kidia,kfdia IF(ktest(JL).EQ.1) THEN ! du ZDWIND=MAX(ABS(ZVPF(JL,JK)-ZVPF(JL,JK-1)),GVSEC) ! Mean flow Richardson number Ri=g/rho[drho/dz / (du/dz)^2] ! Here dp maybe dp^2 ? Need ask Francios lott later PRI(JL,JK)=BV(JL,JK)*(ZDP(JL,JK)/(g*ZRHO(JL,JK)*ZDWIND))**2 PRI(JL,JK)=MAX(PRI(JL,JK),GRCRIT) ENDIF ENDDO end DO !* Define top of 'envelope' layer DO JL=kidia,kfdia ZNU (jl)=0.0 znum(jl)=0.0 ENDDO DO JK=2,nlayer-1 DO JL=kidia,kfdia IF(ktest(JL).EQ.1) THEN IF (JK.GE.IKNU2(JL)) THEN ! level lower than 3*par ! all codes here is to calculate equation 9 ZNUM(JL)=ZNU(JL) ZWIND=(pulow(JL)*pu(jl,jk)+pvlow(jl)*pv(jl,jk))/max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec) ZWIND=max(sqrt(zwind**2),gvsec) ZDELP=pplev(JL,JK+1)-pplev(JL,JK) ! dp ZSTABM=SQRT(MAX(BV(JL,JK ),GSSEC)) ZSTABP=SQRT(MAX(BV(JL,JK+1),GSSEC)) ZRHOM=ZRHO(JL,JK ) ZRHOP=ZRHO(JL,JK+1) ! Equation 9. znu is a critical value to find the blocking layer ZNU(JL) = ZNU(JL) + (ZDELP/g)*((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND ! Found the moutain top IF((ZNUM(JL).LE.GFRCRIT).AND.(ZNU(JL).GT.GFRCRIT).AND.(IKENVH(JL).EQ.nlayer)) THEN IKENVH(JL)=JK ENDIF ENDIF ! (JK.GE.IKNU2(JL)) ENDIF !(ktest(JL).EQ.1) ENDDO endDO ! Calculation of a dynamical mixing height for the breaking of gravity waves DO JL=kidia,kfdia znup(jl)=0.0 znum(jl)=0.0 ENDDO DO JK=nlayer-1,2,-1 DO JL=kidia,kfdia IF(ktest(JL).EQ.1) THEN IF (JK.LT.IKENVH(JL)) THEN ZNUM(JL)=ZNUP(JL) ZWIND=(pulow(JL)*pu(jl,jk)+pvlow(jl)*pv(jl,jk))/max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec) ZWIND=max(sqrt(zwind**2),gvsec) ZDELP=pplev(JL,JK+1)-pplev(JL,JK) ZSTABM=SQRT(MAX(BV(JL,JK ),GSSEC)) ZSTABP=SQRT(MAX(BV(JL,JK+1),GSSEC)) ZRHOM=ZRHO(JL,JK ) ZRHOP=ZRHO(JL,JK+1) ZNUP(JL) = ZNUP(JL) + (ZDELP/g)*((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND ! dynamical mixing height for the breaking of gravity waves IF((ZNUM(JL).LE.1.5).AND.(ZNUP(JL).GT.1.5).AND.(IKCRITH(JL).EQ.nlayer)) THEN IKCRITH(JL)=JK ENDIF ENDIF ! (JK.LT.IKENVH(JL)) ENDIF ! (ktest(JL).EQ.1) ENDDO end DO DO JL=KIDIA,KFDIA IKCRITH(JL)=MIN0(IKCRITH(JL),IKNU(JL)) ENDDO ! directional info for flow blocking ************************* DO jk=ilevh,nlayer DO JL=kidia,kfdia IF(jk.ge.IKENVH(jl)) THEN LO=(pu(JL,jk).LT.GVSEC).AND.(pu(JL,jk).GE.-GVSEC) IF(LO) THEN ZU=pu(JL,jk)+2.*GVSEC ELSE ZU=pu(JL,jk) ENDIF Zphi=ATAN(pv(JL,jk)/ZU) ZPSI(jl,jk)=pthe(jl)*pi/180.-zphi end IF ENDDO end DO ! forms the vertical 'leakiness' ************************** DO JK=ilevh,nlayer DO JL=kidia,kfdia IF(jk.ge.IKENVH(jl)) THEN zggeenv=AMAX1(1.,(zgeom(jl,IKENVH(jl))+zgeom(jl,IKENVH(jl)-1))/2.) zggeom1=AMAX1(zgeom(jl,jk),1.) zgvar=amax1(pvar(jl)*g,1.) ZZDEP(jl,jk)=SQRT((zggeenv-zggeom1)/(zggeom1+zgvar)) endIF ENDDO end DO ! 260 CONTINUE ! continue tag without source, maybe need delete in future RETURN END