Index: /trunk/LMDZ.MARS/libf/phymars/sugwd.F90 =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/sugwd.F90 (revision 3753) +++ /trunk/LMDZ.MARS/libf/phymars/sugwd.F90 (revision 3754) @@ -1,6 +1,6 @@ SUBROUTINE SUGWD(nlayer,sigtest) -! ============================================================================== -! Initialize common variables in yoegwd.h to control the orographic -! graivty wave drag parameterization. That means, all the tunable parameters +! ============================================================================== +! Initialize common variables in yoegwd.h to control the orographic +! gravity wave drag parameterization. That means, all the tunable parameters ! for oro-GW scheme are in this subroutine. ! MARTIN MILLER *ECMWF* ORIGINAL : 90-01-01 Index: /trunk/LMDZ.PLUTO/libf/phypluto/calldrag_noro_mod.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/calldrag_noro_mod.F90 (revision 3754) +++ /trunk/LMDZ.PLUTO/libf/phypluto/calldrag_noro_mod.F90 (revision 3754) @@ -0,0 +1,186 @@ +MODULE calldrag_noro_mod + +IMPLICIT NONE + +CONTAINS + + SUBROUTINE calldrag_noro(ngrid,nlayer,ptimestep,pplay,pplev,pt,pu,pv, & + !output + pdtgw,pdugw,pdvgw) + +!===================================================================================================== +! MODULE designed to call SUBROUTINE drag_noro_mod Interface for sub-grid scale orographic scheme +! The purpose of this subroutine is: +! 1) Make some initial calculation at first call +! 2) Split the calculation in several sub-grid ("sub-domain") to save memory and be able run on +! a workstation at high resolution.The sub-grid size is defined in dimradmars_mod. +! 3) Transfer the increment output of drag_noro_mod into tendencies +! Christophe Hourdin + Francois Forget +! +! VERSION Update: This version uses the variable's splitting, which can be usefull when performing +! very high resolution simulation like LES. ! +! J.-B. Madeleine 10W12 +! +! UPDATE Rewirten into .F90 J.Liu 3/3/2022 Translate the code into .F90. The name of the +! variables are uniformed. Comments are made. +! Unused variables are deleted. +!======================================================================================================= + + use surfdat_h, only: zstd, zsig, zgam, zthe !sub-grid scale standard devitation,slope,anisotropy, and priciple axes angle + ! which will be coded as pvar, psig,pgam,pthe in the called subroutines. + use dimphy, only: ndomainsz + use drag_noro_mod, only: drag_noro + USE ioipsl_getin_p_mod, ONLY : getin_p + + IMPLICIT NONE + + ! 0.1 Declarations : + + ! 0.1.0 Input + INTEGER, intent(in):: ngrid ! Number of atmospheric columns + INTEGER, intent(in):: nlayer ! Number of atmospheric layers + REAL, intent(in):: ptimestep ! Time step of the Physics(s) + REAL, intent(in):: pplev(ngrid,nlayer+1) ! Pressure at 1/2 levels(Pa) + REAL, intent(in):: pplay(ngrid,nlayer) ! Pressure at full levels(Pa) + REAL, intent(in):: pt(ngrid,nlayer) ! Temp at full levels(K) + REAL, intent(in):: pu(ngrid,nlayer) ! Zonal wind at full levels(m/s) + REAL, intent(in):: pv(ngrid,nlayer) ! Meridional wind at full levels(m/s) + + ! 0.1.1 Output + REAL, intent(out):: pdtgw(ngrid,nlayer) ! Tendency on temperature (K/s) due to orographic gravity waves + REAL, intent(out):: pdugw(ngrid,nlayer) ! Tendency on zonal wind (m/s/s) due to orographic gravity waves + REAL, intent(out):: pdvgw(ngrid,nlayer) ! Tendency on meridional wind (m/s/s) due to orographic gravity waves + + !0.1.2 Local variables : + REAL sigtest(nlayer+1) ! sigma coordinate at 1/2 level + INTEGER kgwd ! Number of points at which the scheme is called + INTEGER kgwdim ! kgwdIM=MAX(1,kgwd) + INTEGER ktest(ngrid) ! Map of calling points + INTEGER kdx(ndomainsz) ! Points at which to call the scheme + INTEGER ndomain ! Parameter (ndomain = (ngrid-1) / ndomainsz + 1) + INTEGER l,ig + INTEGER jd,ig0,nd + REAL,save:: zstdthread ! Detecting points concerned by the scheme by height +!$OMP THREADPRIVATE(zstdthread) + REAL zulow(ngrid) ! Low level zonal wind + REAL zvlow(ngrid) ! Low level Meridional wind + REAL zustr(ngrid) ! Low level zonal stress + REAL zvstr(ngrid) ! Low level meridional stress + REAL zplev(ndomainsz,nlayer+1) ! pplev in sub domain + REAL zplay(ndomainsz,nlayer) ! pplay in sub domain + REAL zt(ndomainsz,nlayer) ! pt in sub domain + REAL zu(ndomainsz,nlayer) ! pu in sub domain + REAL zv(ndomainsz,nlayer) ! pv in sub domain + REAL d_t(ndomainsz,nlayer) ! Temperature increment(K) from DRAG_NORO + REAL d_u(ndomainsz,nlayer) ! Zonal wind increment(K=m/s) from DRAG_NORO + REAL d_v(ndomainsz,nlayer) ! Meridional wind increment(K=m/s) from DRAG_NORO + logical,save:: ll=.false. +!$OMP THREADPRIVATE(ll) + LOGICAL,save:: firstcall=.true. +!$OMP THREADPRIVATE(firstcall) + +!----------------------------------------------------------------------------------------------------------------------- +! 1. INITIALISATIONS +!----------------------------------------------------------------------------------------------------------------------- + IF (firstcall) THEN + write(*,*) "calldrag_noro: orographic GW scheme is active!" + zstdthread = 50.0 ! Mars' value !! + call getin_p("oro_gwd_zstdthread", zstdthread) + write(*,*) "oro_gwd_zstdthread: zstdthread=", zstdthread + + do l=1,nlayer+1 + sigtest(l)=pplev(1,l)/pplev(1,1) + enddo + + call sugwd(nlayer,sigtest) + + if (ngrid .EQ. 1) then + IF (ndomainsz .NE. 1) THEN + print* + print*,'ATTENTION !!!' + print*,'pour tourner en 1D, meme pour drag_noro ' + print*,'fixer ndomainsz=1 dans phypluto/dimphy' + print* + call exit(1) + ENDIF + endif + firstcall=.false. + END IF !firstcall + + ! Initialization the tendency of this scheme into zero +! pdtgw(ngrid,nlayer)=0.0 ! Tendency on temperature (K/s) due to orographic gravity waves +! pdugw(ngrid,nlayer)=0.0 ! Tendency on zonal wind (m/s/s) due to orographic gravity waves +! pdvgw(ngrid,nlayer)=0.0 ! Tendency on meridional wind (m/s/s) due to orographic gravity waves + + ! AS: moved out of firstcall to allow nesting+evoluting horiz domain + ndomain = (ngrid-1) / ndomainsz + 1 +!----------------------------------------------------------------------------------------------------------------------- +! 2. Starting loop on sub-domain +!----------------------------------------------------------------------------------------------------------------------- + ! start the loop from the first sub-domain to the ndomain, in which the orographic gravity waves induced + ! wind/temperature tendencies are calculated + DO jd=1,ndomain + ig0=(jd-1)*ndomainsz + if (jd.eq.ndomain) then + nd=ngrid-ig0 + else + nd=ndomainsz + endif + + ! 2.1 Detecting points concerned by the scheme + kgwd=0 + DO ig=ig0+1,ig0+nd + ktest(ig)=0 + ! only the points have a zstd greater than the thread value(default is 50) are concerned. + ll=zstd(ig).gt.zstdthread + IF(ll) then + ! The map of calling is set 1 (which means the orographic map in this points is called?) + ktest(ig)=1 + ! The numbers of the calling points added + kgwd=kgwd+1 + ! The position of the calling points are picked up + kdx(kgwd)=ig - ig0 + ENDIF + ENDDO + kgwdIM=MAX(1,kgwd) + + ! 2.2 Spliting input variable in sub-domain input variables in each level + DO l=1,nlayer+1 + do ig = 1,nd + zplev(ig,l) = pplev(ig0+ig,l) + enddo + end DO + DO l=1,nlayer + do ig = 1,nd + zplay(ig,l) = pplay(ig0+ig,l) + zt(ig,l) = pt(ig0+ig,l) + zu(ig,l) = pu(ig0+ig,l) + zv(ig,l) = pv(ig0+ig,l) + enddo + end DO + + ! 2.3 Calling gravity wave and subgrid scale topo parameterization(all physics in this subroutine) to get + ! zonal wind, meridional wind, and temperature increament + call drag_noro (nd,nlayer,ptimestep,zplay,zplev,zstd(ig0+1),zsig(ig0+1),zgam(ig0+1),zthe(ig0+1),& + kgwd,kgwdim,kdx,ktest(ig0+1),zt, zu, zv, & + ! output + zulow(ig0+1),zvlow(ig0+1),zustr(ig0+1),zvstr(ig0+1),d_t,d_u,d_v) + + ! 2.4 Un-spliting output variable from sub-domain input variables (and devide by ptimestep -> true tendencies) + DO l=1,nlayer + do ig = 1,nd + pdtgw(ig0+ig,l) = d_t(ig,l)/ptimestep + pdugw(ig0+ig,l) = d_u(ig,l)/ptimestep + pdvgw(ig0+ig,l) = d_v(ig,l)/ptimestep + enddo + end DO + + end DO ! (jd=1, ndomain) +! 2. ending loop on sub-domain + +!****************************************************************************** + + END SUBROUTINE calldrag_noro ! end of the subroutine + +END MODULE calldrag_noro_mod ! end of the module + Index: /trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 (revision 3753) +++ /trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 (revision 3754) @@ -4,6 +4,6 @@ logical,save :: callrad,corrk,calldifv,UseTurbDiff !$OMP THREADPRIVATE(callrad,corrk,calldifv,UseTurbDiff) - logical,save :: calladj,calltherm,n2cond,callsoil -!$OMP THREADPRIVATE(calladj,calltherm,n2cond,callsoil) + logical,save :: calladj,calltherm,n2cond,callsoil,calllott +!$OMP THREADPRIVATE(calladj,calltherm,n2cond,callsoil,calllott) logical,save :: callconduct,callmolvis,callmoldiff !$OMP THREADPRIVATE(callconduct,callmolvis,callmoldiff) @@ -199,5 +199,5 @@ real,save :: deltap ! width of transition to alpha_top (Pa) !$OMP THREADPRIVATE(alpha_top,pref,deltap) - + !! Microphysics-specific variables logical,save :: callmufi, call_haze_prod_pCH4 Index: /trunk/LMDZ.PLUTO/libf/phypluto/dimphy.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/dimphy.F90 (revision 3753) +++ /trunk/LMDZ.PLUTO/libf/phypluto/dimphy.F90 (revision 3754) @@ -1,4 +1,4 @@ MODULE dimphy - + INTEGER,SAVE :: klon ! number of atmospheric columns (for this OpenMP subgrid) INTEGER,SAVE :: klev ! number of atmospheric layers, read by master @@ -6,27 +6,28 @@ INTEGER,SAVE :: klevm1 ! number of atmospheric layers-1, read by master ! INTEGER,SAVE :: kflev + integer,save :: ndomainsz !=(ngrid-1)/20 + 1 !$OMP THREADPRIVATE(klon) CONTAINS - + SUBROUTINE Init_dimphy(klon0,klev0) IMPLICIT NONE - + INTEGER, INTENT(in) :: klon0 INTEGER, INTENT(in) :: klev0 - + klon=klon0 - -!$OMP MASTER + +!$OMP MASTER klev=klev0 klevp1=klev+1 klevm1=klev-1 ! kflev=klev -!$OMP END MASTER +!$OMP END MASTER !$OMP BARRIER - + END SUBROUTINE Init_dimphy - + END MODULE dimphy Index: /trunk/LMDZ.PLUTO/libf/phypluto/drag_noro_mod.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/drag_noro_mod.F90 (revision 3754) +++ /trunk/LMDZ.PLUTO/libf/phypluto/drag_noro_mod.F90 (revision 3754) @@ -0,0 +1,149 @@ +MODULE drag_noro_mod + +IMPLICIT NONE + +CONTAINS + + SUBROUTINE drag_noro (ngrid,nlayer,ptimestep,pplay,pplev,pvar, psig, pgam, pthe, & + kgwd,kgwdim,kdx,ktest,t, u, v, & + ! output + pulow, pvlow, pustr, pvstr, d_t, d_u, d_v) + +!-------------------------------------------------------------------------------- +! MODULE contains DRAG_NORO SUBROUNTINE for sub-grid scale orographic sheme. +! 1) Initialization the variables +! 2) Overturn the levels and computes geopotential height +! 3) Call the and ORODRAG subroutine to updates the tendencies +! 4) Overturn back to the normal level of the tendencies and transfer it into +! increments and sum the oro-gw stress. +! the scheme has been called. +! Z.X. Li + F.Lott (LMD/CNRS/ENS) 1995-02-01 +! +! UPDATE: J.Liu 03/03/2022 Translate the code into .F90. The name of the +! variables are uniformed. Comments are made. +! Unused variables are deleted. +!------------------------------------------------------------------------------- + + use dimphy, only: ndomainsz + USE orodrag_mod, ONLY: orodrag + USE comcstfi_mod, ONLY: g, r + + IMPLICIT none + + ! 0. DECLARATIONS: + + ! 0.1 Input: + INTEGER, intent(in):: ngrid ! Number of atmospheric columns [only nd] + INTEGER, intent(in):: nlayer ! Number of atmospheric layers + REAL, intent(in):: ptimestep ! Time step of the Physics(s) + real, intent(in):: pplay(ndomainsz,nlayer) ! Pressure at full levels(Pa) + real, intent(in):: pplev(ndomainsz,nlayer+1) ! Pressure at 1/2 levels(Pa) + REAL, intent(in):: pvar(ndomainsz) ! sub-grid scale standard deviation + REAL, intent(in):: psig(ndomainsz) ! sub-grid scale standard slope + REAL, intent(inout):: pgam(ndomainsz) ! sub-grid scale standard anisotropy + REAL, intent(in):: pthe(ndomainsz) ! sub-grid scale principal axes angle + REAL, intent(in):: u(ndomainsz,nlayer) ! Zonal wind at full levels(m/s) + REAL, intent(in):: v(ndomainsz,nlayer) ! Meridional wind at full levels(m/s) + REAL, intent(in):: t(ndomainsz,nlayer) ! Temperature at full levels(m/s) + INTEGER, intent(in):: kgwd ! Number of points at which the scheme is called + INTEGER, intent(in):: kgwdim ! kgwdim=MAX(1,kgwd) + INTEGER, intent(in):: kdx(ndomainsz) ! Points at which to call the scheme + INTEGER, intent(in):: ktest(ndomainsz) ! Map of calling points + + ! 0.2 Output: + REAL, intent(out):: pulow(ndomainsz) ! Low level zonal wind + REAL, intent(out):: pvlow(ndomainsz) ! Low level meridional wind + REAL, intent(out):: pustr(ndomainsz) ! Low level zonal stress + REAL, intent(out):: pvstr(ndomainsz) ! Low level meridional stress + REAL, intent(out):: d_t(ndomainsz,nlayer) ! Temperature increment(K) due to orographic gravity waves + REAL, intent(out):: d_u(ndomainsz,nlayer) ! Zonal increment(m/s) due to orographic gravity waves + REAL, intent(out):: d_v(ndomainsz,nlayer) ! Meridional increment(m/s) due to orographic gravity waves + + ! 0.3 Variables locales: + INTEGER i, k + REAL inv_pplev(ndomainsz,nlayer+1) ! Inversed (by inverse the pplev pressure) pressure at 1/2 levels + REAL inv_pplay(ndomainsz,nlayer) ! Inversed (by inverse the pplay pressure) pressure at full levels + REAL zgeom(ndomainsz,nlayer) ! Geopotetial height (Inversed ??) + REAL inv_pt(ndomainsz,nlayer) ! Inversed (by inverse the t) temperature (K) at full levels + REAL inv_pu(ndomainsz,nlayer) ! Inversed (by inverse the u) zonal wind (m/s) at full levels + REAL inv_pv(ndomainsz,nlayer) ! Inversed (by inverse the v) meridional wind (m/s) at full levels + REAL pdtdt(ndomainsz,nlayer) ! Temperature tendency outputs from main routine ORODRAG + REAL pdudt(ndomainsz,nlayer) ! Zonal winds tendency outputs from main routine ORODRAG + REAL pdvdt(ndomainsz,nlayer) ! Meridional winds tendency outputs from main routine ORODRAG + +!----------------------------------------------------------------------------------------------------------------------- +! 1. INITIALISATIONS +!----------------------------------------------------------------------------------------------------------------------- + ! 1.1 Initialize the input/output variables (for security purpose) + DO i = 1,ngrid + pulow(i) = 0.0 + pvlow(i) = 0.0 + pustr(i) = 0.0 + pvstr(i) = 0.0 + ENDDO + + DO k = 1, nlayer + DO i = 1, ngrid + d_t(i,k) = 0.0 + d_u(i,k) = 0.0 + d_v(i,k) = 0.0 + pdudt(i,k)=0.0 + pdvdt(i,k)=0.0 + pdtdt(i,k)=0.0 + ENDDO + ENDDO + + ! 1.2 Prepare the ininv_put varibales (Attention: The order of vertical levels increases from top to bottom ) + ! Here the levels are overturned, that is, the surface becomes to the top and the top becomes to the surface + ! and so forth + DO k = 1, nlayer + DO i = 1, ngrid + inv_pt(i,k) = t(i,nlayer-k+1) + inv_pu(i,k) = u(i,nlayer-k+1) + inv_pv(i,k) = v(i,nlayer-k+1) + inv_pplay(i,k) = pplay(i,nlayer-k+1) + inv_pplev(i,k) = pplev(i,nlayer+1-k+1) + ENDDO + endDO + + DO i = 1, ngrid + inv_pplev(i,nlayer+1) = pplev(i,1) + ENDDO + + !calculate g*dz by R*T*[LN(P(z)/P(z+dz))] + DO i = 1, ngrid + zgeom(i,nlayer) = r * inv_pt(i,nlayer)* LOG(inv_pplev(i,nlayer+1)/inv_pplay(i,nlayer)) + ENDDO + + ! sum g*dz from surface to top to get geopotential height + DO k = nlayer-1, 1, -1 + DO i = 1, ngrid + zgeom(i,k) = zgeom(i,k+1) + r * (inv_pt(i,k)+inv_pt(i,k+1))/2.0 * LOG(inv_pplay(i,k+1)/inv_pplay(i,k)) + ENDDO + endDO + +!----------------------------------------------------------------------------------------------------------------------- +! 2. CALL the Main Rountine OORDRAG +!----------------------------------------------------------------------------------------------------------------------- + ! 2.1 call the main rountine to get the tendencies + CALL ORODRAG(ngrid,nlayer,kgwd,kgwdim,kdx,ktest,ptimestep,inv_pplev,inv_pplay,zgeom,& + inv_pt, inv_pu, inv_pv, pvar, psig, pgam, pthe, & + ! output: + pulow,pvlow,pdudt,pdvdt,pdtdt) + + ! 2.2 Transfer tendency into increment by multiply the physical time steps + ! maybe in the future here cancel multiply the ptimestep thus it will not divide ptimestep again in the main rountine + ! calldrag_noro_mod.F90 + DO k = 1, nlayer + DO i = 1, ngrid + d_u(i,nlayer+1-k) = ptimestep*pdudt(i,k) + d_v(i,nlayer+1-k) = ptimestep*pdvdt(i,k) + d_t(i,nlayer+1-k) = ptimestep*pdtdt(i,k) + pustr(i) = pustr(i)+g*pdudt(i,k)*(inv_pplev(i,k+1)-inv_pplev(i,k)) + pvstr(i) = pvstr(i)+g*pdvdt(i,k)*(inv_pplev(i,k+1)-inv_pplev(i,k)) + ENDDO + ENDDO + + END SUBROUTINE drag_noro + +END MODULE drag_noro_mod Index: /trunk/LMDZ.PLUTO/libf/phypluto/gwprofil_mod.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/gwprofil_mod.F90 (revision 3754) +++ /trunk/LMDZ.PLUTO/libf/phypluto/gwprofil_mod.F90 (revision 3754) @@ -0,0 +1,164 @@ +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 dimphy, 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 Index: /trunk/LMDZ.PLUTO/libf/phypluto/gwstress_mod.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/gwstress_mod.F90 (revision 3754) +++ /trunk/LMDZ.PLUTO/libf/phypluto/gwstress_mod.F90 (revision 3754) @@ -0,0 +1,101 @@ +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 dimphy, 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 Index: /trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 (revision 3753) +++ /trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 (revision 3754) @@ -355,4 +355,8 @@ if (is_master) write(*,*) trim(rname)//": callsoil = ",callsoil + calllott=.true. ! default value + call getin_p("calllott",calllott) + write(*,*)" calllott = ",calllott + if (is_master) write(*,*)trim(rname)//& ": Rad transfer is computed every iradia", & Index: /trunk/LMDZ.PLUTO/libf/phypluto/orodrag_mod.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/orodrag_mod.F90 (revision 3754) +++ /trunk/LMDZ.PLUTO/libf/phypluto/orodrag_mod.F90 (revision 3754) @@ -0,0 +1,251 @@ +MODULE orodrag_mod + +IMPLICIT NONE + +CONTAINS + + SUBROUTINE ORODRAG( ngrid,nlayer, kgwd, kgwdim, kdx, ktest, ptimestep, & + pplev,pplay,zgeom,pt,pu, pv, pvar, psig, pgam, pthe, & + !OUTPUTS + PULOW,PVLOW, pdudt,pdvdt,pdtdt) + + !-------------------------------------------------------------------------------- + ! The main routine ORODRAG that does the orographic gravity wave parameterization. + ! This routine computes the physical tendencies of the prognostic variables U,V, and + ! T due to vertical transport by subgridscale orographically excited gravity waves + ! M.MILLER + B.RITTER E.C.M.W.F. 15/06/86. + ! F.LOTT + M. MILLER E.C.M.W.F. 22/11/94 + !-- + ! The purpose of this routine is: + ! 1) call OROSETUP to get the parameters such as the top altitude(level) of the blocked + ! flow and other critical levels. + ! 2) call GWSTRESS and GWPROFILE to get the gw stress. + ! 3) calculate the zonal/meridional wind tendency based on the GW stress and mountain drag + ! on the flow. + ! 4) calculate the temperature tendency based on the differtial of square of the wind + ! between t and t+dt. + ! Rewrited by J.LIU LMDZ.COMMON/phymars/libf 29/01/2022 + !-- + ! 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. + !----------------------------------------------------------------------- + + use dimphy, only: ndomainsz + USE gwstress_mod, ONLY: gwstress + USE gwprofil_mod, ONLY: gwprofil + USE comcstfi_mod, ONLY: g, cpp + USE yoegwd_h, ONLY: gkwake + USE orosetup_mod, ONLY: orosetup + + 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):: kgwd ! Number of points at which the scheme is called + INTEGER, intent(in):: kgwdim ! kgwdim=max(1,kgwd) + INTEGER, intent(in):: kdx(ndomainsz) ! Points at which to call the scheme. + INTEGER, intent(in):: ktest(ndomainsz) ! Map of calling points + + 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):: pvar(ndomainsz) ! Sub-grid scale standard deviation + REAL, intent(in):: psig(ndomainsz) ! Sub-grid scale slope + REAL, intent(inout):: pgam(ndomainsz) ! Sub-grid scale anisotropy + REAL, intent(in):: pthe(ndomainsz) ! Sub-grid scale principal axes angle + REAL, intent(in):: zgeom(ndomainsz,nlayer) ! Geopotential height of full levels + REAL, intent(in):: pplay(ndomainsz,nlayer) ! Pressure at full levels(Pa) (has been inversed by DRAG_NORO=inv_pplay) + REAL, intent(in):: pplev(ndomainsz,nlayer+1) ! Pressure at 1/2 levels(Pa) (has been inversed by DRAG_NORO=inv_pplev) + REAL, intent(in):: ptimestep ! Time step of the Physics(s) + + ! 0.2 OUTPUTS + REAL, intent(out):: pdtdt(ndomainsz,nlayer) ! Tendency on temperature (K/s) due to orographic gravity waves + REAL, intent(out):: pdvdt(ndomainsz,nlayer) ! Tendency on meridional wind (m/s/s) due to orographic gravity waves + REAL, intent(out):: pdudt(ndomainsz,nlayer) ! Tendency on zonal wind (m/s/s) due to orographic gravity waves + REAL, intent(out):: PULOW(ndomainsz) ! Low level zonal wind + REAL, intent(out):: PVLOW(ndomainsz) ! Low level meridional wind + + !0.3 LOCAL ARRAYS +! INTEGER ISECT(ndomainsz) ! not used ? + INTEGER ICRIT(ndomainsz) ! Critical layer where orographic GW breaks + INTEGER IKCRITH(ndomainsz) ! Dynamical mixing height for the breaking of gravity waves + INTEGER IKenvh(ndomainsz) ! Top of the blocked layer + INTEGER IKNU(ndomainsz) ! 4*pvar layer + INTEGER IKNU2(ndomainsz) ! 3*pvar layer + INTEGER IKCRIT(ndomainsz) ! Top of low level flow height +! INTEGER IKHLIM(ndomainsz) ! not used? + + integer ji,jk,jl + integer ilevp1 !=nlayer+1 just used once. Maybe directly use nlayer+1 to replace in the future + +! real ztauf(ndomainsz,nlayer+1) ! not used + real zdelp ! =pplev(nlayer+1)-pplev(nlayer) dp of two 1/2 levels + real zb ! B(pgam), is a seceond order polynomial fit of SS ANISOTRPY pgam(Reference 2) + real zc ! C(pgam), is a seceond order polynomial fit of SS ANISOTRPY pgam(Reference 2) + real zbet ! Equation (16) blocked flow drag + real zconb ! Middle part of the equation (16) + real zabsv ! Tail of equation (16),U*|U|/2+V*|V|/2 + real zzd1 ! Second tail part of the equation (16) + real ratio ! Aspect ratio of the obstacle (mountain),equation (14) + real zust ! U(t+dt) + real zvst ! V(t+dt) + real zdis ! =1/2*[U(t)*|U(t)|+V(t)*|V(t)|-U(t+dt)*|U(t+dt)|-V(t+dt)*|V(t+dt)|] + real ztemp ! T=-g*dtao/U_fistlevel*dp temperature due to GW stress? + REAL ZTAU(ndomainsz,nlayer+1) ! Output of subrountine GWPROFILE: Gravtiy wave stress + REAL BV(ndomainsz,nlayer+1) ! Brunt–Väisälä frequency at 1/2 level (output from OROSETUP input for GWSTRESS and GWPROFIL) + REAL ZVPH(ndomainsz,nlayer+1) ! Low level wind speed U_H in Ref.2 + REAL ZRHO(ndomainsz,nlayer+1) ! density calculated (output) by the OROSETUP rountines(but not used by ORODRAG) + ! used as input for GWSTRESS and GWPROFIL + REAL PRI(ndomainsz,nlayer+1) ! Mean flow richardson number at 1/2 level + REAL ZpsI(ndomainsz,nlayer+1) ! The angle between the incident flow direction and the normal ridge direction of the mountain + REAL Zzdep(ndomainsz,nlayer) ! dp by full level + REAL ZDUDT(ndomainsz) ! du/dt due to oro-gw + REAL ZDVDT(ndomainsz) ! dv/dt due to oro-gw + REAL ZDTDT(ndomainsz) ! dT/dt due to oro-gw + REAL ZDEDT(ndomainsz) ! zdis/ptimestemp, zdedt/Cp=dT/dt + REAL ZVIDIS(ndomainsz) ! in an invalid line, not used + REAL Znu(ndomainsz) ! A critical value see equation 9(Since it only compute inside OROSETUP,Maybe needs delete in future) + REAL Zd1(ndomainsz) ! Bcos^2(psi)-Csin^2(psi) see equation 17 or 18 in Ref.2 + REAL Zd2(ndomainsz) ! (B-C)sin(psi)cos(psi) see equation 17 or 18 in Ref.2 + REAL Zdmod(ndomainsz) ! sqrt(zd1^2+zd2^2) + +!------------------------------------------------------------------------------------ +! 1.INITIALIZATION (NOT USED ) +!------------------------------------------------------------------------------------ +! 100 CONTINUE ! continue tag without source, maybe need delete in future +!* 1.1 Computional constants +! 110 CONTINUE ! continue tag without source, maybe need delete in future +! kfdia=ndomainsz +! ZTMST=TWODT +! IF(NSTEP.EQ.NSTART) ZTMST=0.5*TWODT +! nlayerM1=nlayer-1 +! ZTMST=ptimestep +! ZRTMST=1./ptimestep +! 120 CONTINUE ! continue tag without source, maybe need delete in future +!* 1.3 Check whether row contains point for printing +! 130 CONTINUE ! continue tag without source, maybe need delete in future + +!------------------------------------------------------------------------------------ +! 2. Precompute basic state variables . +!------------------------------------------------------------------------------------ +! 200 CONTINUE ! continue tag without source, maybe need delete in future + ! 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 + CALL 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) + +!------------------------------------------------------------------------------------ +! 3. Computes low level stresses using subcritical and super critical forms. +! Computes anisotropy coefficient as measure of orographic two-dimensionality +!------------------------------------------------------------------------------------ +! 300 CONTINUE ! continue tag without source, maybe need delete in future + ! Computes low level stresses + CALL 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 ) + +!------------------------------------------------------------------------------------ +! 4. Compute stress profile. +!------------------------------------------------------------------------------------ +! 400 CONTINUE ! continue tag without source, maybe need delete in future + ! Compute stress profile. + CALL 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 ) + +!------------------------------------------------------------------------------------ +! 5. Compute tendencies. +!------------------------------------------------------------------------------------ +! 500 CONTINUE ! continue tag without source, maybe need delete in future +! EXplicit solution at all levels for gravity wave +! Implicit solution for the blocked levels +! DO JL=KIDIA,KFDIA + ! Initialization the tendencies + DO JL=1,ndomainsz + ZVIDIS(JL)=0.0 !An invalid lines contain this variable + ZDUDT(JL)=0.0 + ZDVDT(JL)=0.0 + ZDTDT(JL)=0.0 + ENDDO + + ! all the calculation of equation (16-18) are in the flowing loop + ! 1) if the wind flow over the mountain [Either the mountain is too low or the + ! flow higher than the blocked level], then calculate the wind tendencies by oro-GW stress directly + ! 2) if the flow level lower than the blocked level, then calculate the wind tendencies by obstacle drag + ! + ILEVP1=nlayer+1 + DO JK=1,nlayer + DO JL=1,kgwd + ! Pick up the position of each calling points + JI=kdx(JL) + ! dp + ZDELP=pplev(Ji,JK+1)-pplev(Ji,JK) + ! T=-g*dtao/U_fistlevel*dp temperature due to GW stress? + ZTEMP=-g*(ZTAU(Ji,JK+1)-ZTAU(Ji,JK))/(ZVPH(Ji,ILEVP1)*ZDELP) + ! du/dt=(U_low*tao1-V_low*tao2)*T/(tao1^2+tao2^2)*0.5 + ZDUDT(JI)=(PULOW(JI)*Zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji) + ! dv/dt=(V_low*tao1+U_low*tao2)*T/(tao1^2+tao2^2)*0.5 + ZDVDT(JI)=(PVLOW(JI)*Zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji) + if(jk.ge.ikenvh(ji)) then ! IF the level lower than(notice here the level is inversed thus it is lower) + ! the top of the blocked layer (commonlly is the highest top of a mountain.) + ! Then use equation (16) to calculate the mountain drag on the incident flow + ! Coefficient B and C Ref.2 + zb=1.0-0.18*pgam(ji)-0.04*pgam(ji)**2 + zc=0.48*pgam(ji)+0.3*pgam(ji)**2 + ! Middle part of the equation (16) + zconb=2.*ptimestep*GKWAKE*psig(ji)/(4.*pvar(ji)) + ! Tail part of the equation (16),U*|U|/2+V*|V|/2 + zabsv=sqrt(pu(JI,JK)**2+pv(JI,JK)**2)/2. + ! Second tail part of the equation (16) + zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2 + ! Aspect ratio of the obstacle(topography), equation(14) + ratio=(cos(zpsi(ji,jk))**2+pgam(ji)*sin(zpsi(ji,jk))**2)/(pgam(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2) + ! Equation (16) + zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv + ! du/dt and dv/dt due to the mountain drag + zdudt(ji)=-pu(ji,jk)/ptimestep + zdvdt(ji)=-pv(ji,jk)/ptimestep + zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet)) + zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet)) + end if + ! wind tendencies due to the orographic gravity wave (due to stress(flow over) and drag (blocked)) + pdudt(JI,JK)=ZDUDT(JI) + pdvdt(JI,JK)=ZDVDT(JI) + ! winds at t+dt (oro-gw induced increaments are added) + ZUST=pu(JI,JK)+ptimestep*ZDUDT(JI) + ZVST=pv(JI,JK)+ptimestep*ZDVDT(JI) + ! =1/2*[U(t)*|U(t)|+V(t)*|V(t)|-U(t+dt)*|U(t+dt)|-V(t+dt)*|V(t+dt)|] + ZDIS=0.5*(pu(JI,JK)**2+pv(JI,JK)**2-ZUST**2-ZVST**2) + ZDEDT(JI)=ZDIS/ptimestep + ! This line no use maybe delete in the future + ZVIDIS(JI)=ZVIDIS(JI)+ZDIS*ZDELP + ! Temperature tendencies due to the orographic gravity wave (due to stress(flow over) and drag (blocked)) + ZDTDT(JI)=ZDEDT(JI)/cpp + pdtdt(JI,JK)=ZDTDT(JI) + ENDDO !JL=1,kgwd + end DO !JK=1,nlayer + + END SUBROUTINE ORODRAG + +END MODULE orodrag_mod Index: /trunk/LMDZ.PLUTO/libf/phypluto/orosetup.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/orosetup.F90 (revision 3754) +++ /trunk/LMDZ.PLUTO/libf/phypluto/orosetup.F90 (revision 3754) @@ -0,0 +1,476 @@ +MODULE orosetup_mod + +IMPLICIT NONE + +CONTAINS + +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 dimphy, only: ndomainsz + use comcstfi_mod, 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 +!-------------------------------------------------------------------------------- + + !* 1.1 COMPUTATIONAL CONSTANTS + kidia=1 + kfdia=ngrid + 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 +!----------------------------------------------------------------------------------------------------- + + ! 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 + + ! 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 + + 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 + + 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 + +END SUBROUTINE OROSETUP + +END MODULE orosetup_mod Index: /trunk/LMDZ.PLUTO/libf/phypluto/phys_state_var_mod.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/phys_state_var_mod.F90 (revision 3753) +++ /trunk/LMDZ.PLUTO/libf/phypluto/phys_state_var_mod.F90 (revision 3754) @@ -9,5 +9,5 @@ !====================================================================== ! Declaration des variables - USE dimphy, only : klon,klev + USE dimphy, only : klon,klev,ndomainsz USE comsoil_h, only : nsoilmx use comsaison_h, only: mu0, fract @@ -117,4 +117,5 @@ !rugoro(:) ! longueur de rugosite de l'OESM + ndomainsz=(klon-1)/20 + 1 ALLOCATE(phisfi(klon)) ALLOCATE(tsurf(klon)) Index: /trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3753) +++ /trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3754) @@ -44,7 +44,8 @@ obliquit, z0, adjust, tpal use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp + use calldrag_noro_mod, only: calldrag_noro use time_phylmdz_mod, only: daysec use callkeys_mod, only: albedo_spectral_mode, calladj, calldifv, & - callrad, callsoil, nosurf, & + calllott, callrad, callsoil, nosurf, & callconduct,callmolvis,callmoldiff, & corrk, & @@ -117,4 +118,6 @@ ! II.2.b Option 2 : Atmosphere has no radiative effect. ! +! II.3 Gravity wave and subgrid scale topography drag : +! ! III. Vertical diffusion (turbulent mixing) ! @@ -385,4 +388,6 @@ real zdhdif(ngrid,nlayer) ! Turbdiff/vdifc routines. real zdhadj(ngrid,nlayer) ! Convadj routine. + REAL zdtgw(ngrid,nlayer) ! Gravity waves (K/s) + REAL zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer) ! Gravity waves (m.s-2) REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer) ! condense_n2 routine. @@ -1143,4 +1148,23 @@ endif + +!----------------------------------------------------------------------- +! II.3 Gravity wave and subgrid scale topography drag : +! ------------------------------------------------- + + IF(calllott)THEN + CALL calldrag_noro(ngrid,nlayer,ptimestep, & + zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw) + + DO l=1,nlayer + DO ig=1,ngrid + pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l) + pdu(ig,l)=pdu(ig,l)+zdugw(ig,l) + pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l) + ENDDO + ENDDO + ENDIF + + ! -------------------------------------------- ! III. Vertical diffusion (turbulent mixing) : Index: /trunk/LMDZ.PLUTO/libf/phypluto/sugwd.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/sugwd.F90 (revision 3754) +++ /trunk/LMDZ.PLUTO/libf/phypluto/sugwd.F90 (revision 3754) @@ -0,0 +1,73 @@ + SUBROUTINE SUGWD(nlayer,sigtest) +! ============================================================================== +! Initialize common variables in yoegwd.h to control the orographic +! gravity wave drag parameterization. That means, all the tunable parameters +! for oro-GW scheme are in this subroutine. +! MARTIN MILLER *ECMWF* ORIGINAL : 90-01-01 +! Update: Jiandong Liu 2022/03/15 Rewirite into .F90 and +! comment. +! REFERENCE. +! ---------- +! ECMWF Research Department documentation of the IFS +!=============================================================================== + + USE yoegwd_h, ONLY: GFRCRIT, GRCRIT, GVCRIT + USE yoegwd_h, ONLY: GKDRAG, GKDRAGL, GHMAX + USE yoegwd_h, ONLY: GRAHILO, GSIGCR, GSSEC + USE yoegwd_h, ONLY: GTSEC, GVSEC, GKWAKE + USE yoegwd_h, ONLY: NKTOPG + + implicit none + + ! 0.1 Inputs: + integer,intent(in):: nlayer ! Number of model levels + REAL,intent(in):: sigtest(nlayer+1) ! Vertical coordinate table + + ! 0.2 Outputs: + ! None. + + ! 0.3 Local variables + real zsigt ! Top of the sigma coordinates? + real zpr ! Reference pressure ? + real zpm1r ! Pressure at full layer ? + integer jk + +!------------------------------------------------------------------------------- +! 1. Set the values of the parameters +!------------------------------------------------------------------------------- +! PRINT *,' Dans sugwd nlayer=',nlayer,' SIG=',sigtest + GHMAX=10000. + +! old ZSIGT=0.94 +! old ZPR=80000. + ZSIGT=0.85 ! Sigmal levels + ZPR=100000. ! Surface (Reference) Pressure? + + ! ! Condition to find NKTOPG layer, which NKTOPG is a condition to set + ! 1*pvar and 2*pvar layers (OROSETUP) + DO JK=1,nlayer-1 + ZPM1R=0.5*ZPR*(sigtest(JK)+sigtest(JK+1)) + IF((ZPM1R/ZPR).GE.ZSIGT)THEN + NKTOPG=JK + ENDIF + ENDDO + WRITE(*,*) 'In sugwd NKTOPG=',NKTOPG + + GSIGCR=0.80 ! Sigmal levels to found the top of low level flow height (OROSETUP) + GKDRAG= 0.1 ! used to be 0.1 for mcd Version 1 and 2 (before 10/2000) (OROSETUP) + + GFRCRIT=1.0 + GKWAKE=1.0 ! The G in equation (16) + GRCRIT=0.25 ! Critical value for Mean flow richardson number(OROSETUP) + GKDRAGL=4.*GKDRAG + GRAHILO=1. + GVCRIT =0.0 +!------------------------------------------------------------------------------- +! 2. Set values of security parameters +!------------------------------------------------------------------------------- + GVSEC=0.10 ! Security values for For normal wind (pu^2+pv^2)^0.5(OROSETUP,GWSTRESS) + GSSEC=1.E-12 ! Security values for Brunt–Väisälä frequency N^2 (OROSETUP) + GTSEC=1.E-07 ! Security values for Sub-grid scale anisotropy(OROSETUP,GWSTRESS) + + RETURN + END Index: /trunk/LMDZ.PLUTO/libf/phypluto/yoegwd_h.F90 =================================================================== --- /trunk/LMDZ.PLUTO/libf/phypluto/yoegwd_h.F90 (revision 3754) +++ /trunk/LMDZ.PLUTO/libf/phypluto/yoegwd_h.F90 (revision 3754) @@ -0,0 +1,27 @@ +module yoegwd_h +! ----------------------------------------------------------------- +!* *COMMON* *YOEGWD* - PARAMETERS FOR GRAVITY WAVE DRAG CALCULATIONS +! ----------------------------------------------------------------- +implicit none + + real,save :: GFRCRIT + real,save :: GRCRIT + real,save :: GVCRIT +!$OMP THREADPRIVATE(GFRCRIT,GRCRIT,GVCRIT) + real,save :: GKDRAG + real,save :: GKDRAGL + real,save :: GHMAX +!$OMP THREADPRIVATE(GKDRAG,GKDRAGL,GHMAX) + real,save :: GRAHILO + real,save :: GSIGCR + real,save :: GSSEC +!$OMP THREADPRIVATE(GRAHILO,GSIGCR,GSSEC) + real,save :: GTSEC + real,save :: GVSEC + real,save :: GKWAKE +!$OMP THREADPRIVATE(GTSEC,GVSEC,GKWAKE) + integer :: NKTOPG +!$OMP THREADPRIVATE(NKTOPG) + +end module yoegwd_h +