source: trunk/LMDZ.MARS/libf/phymars/gwprofil_mod.F90 @ 2937

Last change on this file since 2937 was 2651, checked in by emillour, 4 years ago

Mars GCM:

  • turn sugwd.F to sugwd.F90 with extra comments
  • turn yoegwd.h into a module

JL+EM
File size: 6.7 KB
RevLine 
[2642]1MODULE gwprofil_mod
[1913]2     
[2642]3IMPLICIT NONE
[1913]4     
[2642]5CONTAINS
[1913]6     
[2642]7      SUBROUTINE GWPROFIL( ngrid, nlayer, kgwd ,kdx, ktest, &
8                 IKCRITH, ICRIT, IKENVH, IKNU,IKNU2,pplev,  &
9                 ZRHO, BV, ZVPH, PRI,zdmod, psig , pvar,    &
10                 !not used variables
11                 !IKCRIT,ZTAUf,ZTFR,znu,pgam,
12                 ! in-output
13                 ZTAU      )
[38]14
[2642]15        !---------------------------------------------------------------------------
16        !     THe stress profile for gravity waves is computed as follows:
17        !     It is constant (no GWD) at the levels between the ground and
18        !     the top of the  blocked layer (IKENVH). It decreases linearly
19        !     with heights from the top of the blocked layer to 3*varor(IKNU)
20        !     to simulates lee waves or nonlinear gravity wave breaking.
21        !     Above it is constant, except when the wave encounters a critical
22        !     level(ICRIT) or when it breaks.!
23        !     PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93)
24        !     REFERENCE.
25        !     SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
26        !     MODIFICATIONS.
27        !     Rewrited by J.Liu 03/03/2022
28        !---------------------------------------------------------------------------
29       
30      use dimradmars_mod, only: ndomainsz 
[2651]31      use yoegwd_h, only: gkdrag, grcrit, gssec, gtsec
32     
[38]33      implicit none
[2642]34 
35      ! 0. DECLARATIONS:
36      ! 0.1   ARGUMENTS
37      integer, intent(in):: ngrid    ! Number of atmospheric columns
38      integer, intent(in):: nlayer   ! Number of atmospheric layers
39      INTEGER, intent(in) :: kgwd    ! Number of points at which the scheme is called     
40!      INTEGER, INTENT(IN) :: IKCRIT(ndomainsz  ) ! Top of low level flow height (not used)
41      INTEGER, INTENT(IN) :: IKCRITH(ndomainsz  ) ! dynamical mixing height for the breaking of gravity waves
42      INTEGER, INTENT(IN) :: ICRIT(ndomainsz  )   ! Critical layer where orographic GW breaks
43      INTEGER, INTENT(IN) :: kdx(ndomainsz  )     ! Points at which to call the scheme   
44      INTEGER, INTENT(IN) :: ktest(ndomainsz  )   ! Map of calling points
45      INTEGER, INTENT(IN) :: IKENVH(ndomainsz  )  ! Top of the  blocked layer
46      INTEGER, INTENT(IN) :: IKNU(ndomainsz  )    ! 4*pvar layer
47      INTEGER, INTENT(IN) :: IKNU2(ndomainsz  )   ! 3*pvar layer
[38]48
[2642]49      REAL, INTENT(IN):: pplev(ndomainsz  ,nlayer+1)   ! Pressure at 1/2 levels(Pa)
50      REAL, INTENT(IN):: BV(ndomainsz  ,nlayer+1)      ! Brunt–Väisälä frequency at 1/2 level (output from OROSETUP input for GWSTRESS and GWPROFIL)
51      REAL, INTENT(IN):: ZRHO  (ndomainsz  ,nlayer+1)  ! Atmospheric density at 1/2 level
52      REAL, INTENT(IN):: ZVPH (ndomainsz  ,nlayer+1)   ! SUB-GRID SCALE STANDARD DEVIATION at 1/2 level
53      REAL, INTENT(IN):: PRI   (ndomainsz  ,nlayer+1)  ! Mean flow richardson number at 1/2 level
54!      REAL, INTENT(IN):: ZTFR (ndomainsz  )           ! not used
55!      REAL, INTENT(IN):: ZTAUf (ndomainsz  ,nlayer+1) ! not used
56!     
57      REAL, INTENT(IN):: zdmod (ndomainsz  )
58!      REAL, INTENT(IN):: znu (ndomainsz  )            ! not used
59      REAL, INTENT(IN):: psig(ndomainsz  )             ! Sub-grid scale slope
60!      REAL, INTENT(IN):: pgam(ndomainsz  )            ! Sub-grid scale anisotropy(not used)
61      REAL, INTENT(IN):: pvar(ndomainsz  )             ! Sub-grid scale standard deviation
62      REAL, INTENT(inout):: ZTAU(ndomainsz  ,nlayer+1) ! Gravity waves induced stress
[38]63
[2642]64      ! 0.2   LOCAL ARRAYS
[38]65      real zsqr,zalfa,zriw,zalpha,zb,zdel,zdz2n,zdelp,zdelpt
66      integer ji,jk,jl,ilevh
[2642]67      REAL ZDZ2 (ndomainsz  ,nlayer) , ZNORM(ndomainsz  ) , zoro(ndomainsz  )
68      REAL Z_TAU (ndomainsz  ,nlayer+1)
[38]69
[2642]70!-----------------------------------------------------------------------
71!*         1.    Initialization
72!-----------------------------------------------------------------------
73!      kidia=1
74!      kfdia=ngrid
75 100  CONTINUE
76      ! 1.1 Computional constants .
77      ilevh=nlayer/3
[38]78
[2642]79      DO ji=1,kgwd
80        jl=kdx(ji)
81        Zoro(JL)=psig(JL)*zdmod(JL)/4./max(pvar(jl),1.0)
82        Z_TAU(JL,IKNU(JL)+1)=ZTAU(JL,IKNU(JL)+1)
83        Z_TAU(JL,nlayer+1)=ZTAU(JL,nlayer+1)
84      ENDDO
85     
86      ! all start here with this long loop
87      DO JK=nlayer,2,-1   
88        ! 4.1   Constant wave stress until top of the blocking layer
[38]89  410 CONTINUE
[2642]90        DO ji=1,kgwd
91           jl=kdx(ji)
92           IF(JK.GE.IKNU2(JL)) THEN
93           ZTAU(JL,JK)=Z_TAU(JL,nlayer+1)
[38]94           ENDIF
[2642]95        ENDDO           
96
97      !*         4.15  Constant shear stress until the top of the low level flow layer
[38]98 415  CONTINUE
[2642]99      ! 4.2    Wave displacement at next level.
[38]100  420 CONTINUE
[2642]101
102        DO ji=1,kgwd
103           jl=kdx(ji)
104           IF(JK.LT.IKNU2(JL)) THEN
105                ZNORM(JL)=gkdrag*ZRHO(JL,JK)*SQRT(BV(JL,JK))*ZVPH(JL,JK)*zoro(jl)
106                ZDZ2(JL,JK)=ZTAU(JL,JK+1)/max(ZNORM(JL),gssec)
107           ENDIF
108        ENDDO
109
110       ! 4.3    Wave Richardson number, new wave displacement and stress:
111       ! Breaking evaluation and critical level C
112        DO ji=1,kgwd
113           jl=kdx(ji)
114           IF(JK.LT.IKNU2(JL)) THEN
115                IF((ZTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.ICRIT(JL))) THEN
116                   ZTAU(JL,JK)=0.0
117                ELSE
118                   ZSQR=SQRT(PRI(JL,JK))
119                   ZALFA=SQRT(BV(JL,JK)*ZDZ2(JL,JK))/ZVPH(JL,JK)
120                   ZRIW=PRI(JL,JK)*(1.-ZALFA)/(1+ZALFA*ZSQR)**2
121                   IF(ZRIW.LT.GRCRIT) THEN
122                     ZDEL=4./ZSQR/GRCRIT+1./GRCRIT**2+4./GRCRIT
123                     ZB=1./GRCRIT+2./ZSQR
124                     ZALPHA=0.5*(-ZB+SQRT(ZDEL))
125                     ZDZ2N=(ZVPH(JL,JK)*ZALPHA)**2/BV(JL,JK)
126                     ZTAU(JL,JK)=ZNORM(JL)*ZDZ2N
127                   ELSE
128                     ZTAU(JL,JK)=ZNORM(JL)*ZDZ2(JL,JK)
129                   ENDIF
130                   ZTAU(JL,JK)=MIN(ZTAU(JL,JK),ZTAU(JL,JK+1))
131                ENDIF
132           ENDIF
133        ENDDO 
134      end DO !JK=nlayer,2,-1
135!-------------------------------------------------------------------------------------------
136     
[38]137  440 CONTINUE
138 
[2642]139!     write(*,*) 'ZTAU'
140!     write(*,99) ((ji,ilevh,ZTAU(ji,ilevh),ji=1,ndomainsz  ),
141!    .                  ilevh=1,nlayer+1)
142 99   FORMAT(i3,i3,f15.5)  ! not used
[38]143
[2642]144      !  Proganisation of the stress profile if breaking occurs at low level
145      DO ji=1,kgwd
146         jl=kdx(ji)
147         Z_TAU(JL,IKENVH(JL))=ZTAU(JL,IKENVH(JL))
148         Z_TAU(JL,IKCRITH(JL))=ZTAU(JL,IKCRITH(JL))
149      ENDDO     
[38]150
[2642]151      DO JK=1,nlayer     
152        DO ji=1,kgwd
153           jl=kdx(ji)               
154           IF(JK.GT.IKCRITH(JL).AND.JK.LT.IKENVH(JL))THEN
155             ZDELP=pplev(JL,JK)-pplev(JL,IKENVH(JL))
156             ZDELPT=pplev(JL,IKCRITH(JL))-pplev(JL,IKENVH(JL))
157             ZTAU(JL,JK)=Z_TAU(JL,IKENVH(JL))+(Z_TAU(JL,IKCRITH(JL))-Z_TAU(JL,IKENVH(JL)))*ZDELP/ZDELPT   
158          ENDIF           
159          ENDDO
160        end DO       
[38]161
[1913]162      END SUBROUTINE GWPROFIL
163
[2642]164END MODULE gwprofil_mod
Note: See TracBrowser for help on using the repository browser.