/[lmdze]/trunk/phylmd/Orography/gwprofil.f90
ViewVC logotype

Diff of /trunk/phylmd/Orography/gwprofil.f90

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/Orography/gwprofil.f90 revision 23 by guez, Mon Dec 14 15:25:16 2009 UTC trunk/phylmd/Orography/gwprofil.f90 revision 328 by guez, Thu Jun 13 14:40:06 2019 UTC
# Line 1  Line 1 
1      SUBROUTINE gwprofil(nlon,nlev,kgwd,kdx,ktest,kkcrith,kcrit,paphm1,prho, &  module gwprofil_m
         pstab,pvph,pri,ptau,pdmod,psig,pvar)  
2    
3  !**** *GWPROFIL*    IMPLICIT NONE
4    
5  !     PURPOSE.  contains
 !     --------  
6    
7  !**   INTERFACE.    SUBROUTINE gwprofil(nlon, nlev, ktest, kkcrith, kcrit, paphm1, prho, pstab, &
8  !     ----------         pvph, pri, ptau, pdmod, psig, pvar)
 !          FROM *GWDRAG*  
9    
10  !        EXPLICIT ARGUMENTS :      ! Method. The stress profile for gravity waves is computed as
11  !        --------------------      ! follows: it is constant (no gwd) at the levels between the
12  !     ==== INPUTS ===      ! ground and the top of the blocked layer (kkenvh).  It decreases
13  !     ==== OUTPUTS ===      ! linearly with height from the top of the blocked layer to 3 *
14        ! varor (kknu), to simulate lee waves or nonlinear gravity wave
15  !        IMPLICIT ARGUMENTS :   NONE      ! breaking. Above it is constant, except when the wave encounters
16  !        --------------------      ! a critical level (kcrit) or when it breaks.
17    
18  !     METHOD:      ! Reference. See ECMWF research department documentation of the
19  !     -------      ! "I.F.S."
20  !     THE STRESS PROFILE FOR GRAVITY WAVES IS COMPUTED AS FOLLOWS:  
21  !     IT IS CONSTANT (NO GWD) AT THE LEVELS BETWEEN THE GROUND      ! Modifications. Passage of the new gwdrag TO I.F.S. (F. LOTT,
22  !     AND THE TOP OF THE BLOCKED LAYER (KKENVH).      ! 22/11/93)
23  !     IT DECREASES LINEARLY WITH HEIGHTS FROM THE TOP OF THE  
24  !     BLOCKED LAYER TO 3*VAROR (kKNU), TO SIMULATES LEE WAVES OR      USE dimphy, ONLY : klev, klon
25  !     NONLINEAR GRAVITY WAVE BREAKING.      USE yoegwd, ONLY : gkdrag, grahilo, grcrit, gssec, gtsec, nstra
26  !     ABOVE IT IS CONSTANT, EXCEPT WHEN THE WAVE ENCOUNTERS A CRITICAL  
27  !     LEVEL (KCRIT) OR WHEN IT BREAKS.      INTEGER, intent(in):: nlon, nlev
28        INTEGER, intent(in):: ktest(nlon), kkcrith(nlon), kcrit(nlon)
29        REAL, intent(in):: paphm1(nlon, nlev+1), prho(nlon, nlev+1)
30        REAL, intent(in):: pstab(nlon, nlev+1)
31  !     EXTERNALS.      real, intent(in):: pvph(nlon, nlev+1), pri(nlon, nlev+1)
32  !     ----------      real ptau(nlon, nlev+1)
33        REAL, intent(in):: pdmod(nlon)
34        REAL, INTENT (IN) :: psig(nlon)
35  !     REFERENCE.      REAL, INTENT (IN) :: pvar(nlon)
36  !     ----------  
37        ! Local:
38  !        SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."      INTEGER jl, jk
39        REAL zsqr, zalfa, zriw, zdel, zb, zalpha, zdz2n
40  !     AUTHOR.      REAL zdelp, zdelpt
41  !     -------      REAL zdz2(klon, klev), znorm(klon), zoro(klon)
42        REAL ztau(klon, klev+1)
43  !     MODIFICATIONS.  
44  !     --------------      !-----------------------------------------------------------------------
45  !     PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93)  
46  !-----------------------------------------------------------------------      ! 1. INITIALIZATION
47        USE dimens_m  
48        USE dimphy      ! COMPUTATIONAL CONSTANTS.
49        USE yomcst  
50        USE yoegwd      DO jl = 1, klon
51        IMPLICIT NONE         IF (ktest(jl)==1) THEN
52              zoro(jl) = psig(jl)*pdmod(jl)/4./max(pvar(jl), 1.0)
53              ztau(jl, klev+1) = ptau(jl, klev+1)
54           END IF
55        end DO
56    
57  !-----------------------------------------------------------------------      DO jk = klev, 2, -1
58           ! 4.1 CONSTANT WAVE STRESS UNTIL TOP OF THE
59  !*       0.1   ARGUMENTS         ! BLOCKING LAYER.
60  !              ---------         DO jl = 1, klon
61              IF (ktest(jl)==1) THEN
62        INTEGER nlon, nlev               IF (jk>kkcrith(jl)) THEN
63        INTEGER kkcrith(nlon), kcrit(nlon), kdx(nlon), ktest(nlon)                  ptau(jl, jk) = ztau(jl, klev+1)
64                 ELSE
65                    ptau(jl, jk) = grahilo*ztau(jl, klev+1)
66        REAL paphm1(nlon,nlev+1), pstab(nlon,nlev+1), prho(nlon,nlev+1), &               END IF
67          pvph(nlon,nlev+1), pri(nlon,nlev+1), ptau(nlon,nlev+1)            END IF
68           end DO
69        REAL pdmod(nlon)  
70        REAL, INTENT (IN) :: psig(nlon)         ! 4.2 WAVE DISPLACEMENT AT NEXT LEVEL.
71        REAL, INTENT (IN) :: pvar(nlon)         DO jl = 1, klon
72              IF (ktest(jl)==1) THEN
73  !-----------------------------------------------------------------------               IF (jk<kkcrith(jl)) THEN
74                    znorm(jl) = gkdrag * prho(jl, jk) * sqrt(pstab(jl, jk)) &
75  !*       0.2   LOCAL ARRAYS                       * pvph(jl, jk)* zoro(jl)
76  !              ------------                  zdz2(jl, jk) = ptau(jl, jk+1)/max(znorm(jl), gssec)
77                 END IF
78        INTEGER ilevh, ji, kgwd, jl, jk            END IF
79        REAL zsqr, zalfa, zriw, zdel, zb, zalpha, zdz2n         end DO
80        REAL zdelp, zdelpt  
81        REAL zdz2(klon,klev), znorm(klon), zoro(klon)         ! 4.3 WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
82        REAL ztau(klon,klev+1)         ! AND STRESS: BREAKING EVALUATION AND CRITICAL
83           ! LEVEL
84  !-----------------------------------------------------------------------         DO jl = 1, klon
85              IF (ktest(jl)==1) THEN
86  !*         1.    INITIALIZATION               IF (jk<kkcrith(jl)) THEN
87  !                --------------                  IF ((ptau(jl, jk+1)<gtsec) .OR. (jk<=kcrit(jl))) THEN
88                       ptau(jl, jk) = 0.0
 !      print *,' entree gwprofil'  
 100   CONTINUE  
   
   
 !*    COMPUTATIONAL CONSTANTS.  
 !     ------------- ----------  
   
       ilevh = klev/3  
   
 !     DO 400 ji=1,kgwd  
 !     jl=kdx(ji)  
 !  Modif vectorisation 02/04/2004  
       DO 400 jl = 1, klon  
         IF (ktest(jl)==1) THEN  
           zoro(jl) = psig(jl)*pdmod(jl)/4./max(pvar(jl),1.0)  
           ztau(jl,klev+1) = ptau(jl,klev+1)  
         END IF  
 400   CONTINUE  
   
   
       DO 430 jk = klev, 2, -1  
   
   
 !*         4.1    CONSTANT WAVE STRESS UNTIL TOP OF THE  
 !                 BLOCKING LAYER.  
 410     CONTINUE  
   
 !     DO 411 ji=1,kgwd  
 !     jl=kdx(ji)  
 !  Modif vectorisation 02/04/2004  
         DO 411 jl = 1, klon  
           IF (ktest(jl)==1) THEN  
             IF (jk>kkcrith(jl)) THEN  
               ptau(jl,jk) = ztau(jl,klev+1)  
 !          ENDIF  
 !          IF(JK.EQ.KKCRITH(JL)) THEN  
             ELSE  
               ptau(jl,jk) = grahilo*ztau(jl,klev+1)  
             END IF  
           END IF  
 411     CONTINUE  
   
 !*         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 421 ji=1,kgwd  
 !     jl=kdx(ji)  
 !  Modif vectorisation 02/04/2004  
         DO 421 jl = 1, klon  
           IF (ktest(jl)==1) THEN  
             IF (jk<kkcrith(jl)) THEN  
               znorm(jl) = gkdrag*prho(jl,jk)*sqrt(pstab(jl,jk))*pvph(jl,jk)* &  
                 zoro(jl)  
               zdz2(jl,jk) = ptau(jl,jk+1)/max(znorm(jl),gssec)  
             END IF  
           END IF  
 421     CONTINUE  
   
 !*         4.3    WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT  
 !*                AND STRESS:  BREAKING EVALUATION AND CRITICAL  
 !                 LEVEL  
   
   
 !     DO 431 ji=1,kgwd  
 !     jl=Kdx(ji)  
 !  Modif vectorisation 02/04/2004  
         DO 431 jl = 1, klon  
           IF (ktest(jl)==1) THEN  
   
             IF (jk<kkcrith(jl)) THEN  
               IF ((ptau(jl,jk+1)<gtsec) .OR. (jk<=kcrit(jl))) THEN  
                 ptau(jl,jk) = 0.0  
               ELSE  
                 zsqr = sqrt(pri(jl,jk))  
                 zalfa = sqrt(pstab(jl,jk)*zdz2(jl,jk))/pvph(jl,jk)  
                 zriw = pri(jl,jk)*(1.-zalfa)/(1+zalfa*zsqr)**2  
                 IF (zriw<grcrit) THEN  
                   zdel = 4./zsqr/grcrit + 1./grcrit**2 + 4./grcrit  
                   zb = 1./grcrit + 2./zsqr  
                   zalpha = 0.5*(-zb+sqrt(zdel))  
                   zdz2n = (pvph(jl,jk)*zalpha)**2/pstab(jl,jk)  
                   ptau(jl,jk) = znorm(jl)*zdz2n  
89                  ELSE                  ELSE
90                    ptau(jl,jk) = znorm(jl)*zdz2(jl,jk)                     zsqr = sqrt(pri(jl, jk))
91                       zalfa = sqrt(pstab(jl, jk)*zdz2(jl, jk))/pvph(jl, jk)
92                       zriw = pri(jl, jk)*(1.-zalfa)/(1+zalfa*zsqr)**2
93                       IF (zriw<grcrit) THEN
94                          zdel = 4./zsqr/grcrit + 1./grcrit**2 + 4./grcrit
95                          zb = 1./grcrit + 2./zsqr
96                          zalpha = 0.5*(-zb+sqrt(zdel))
97                          zdz2n = (pvph(jl, jk)*zalpha)**2/pstab(jl, jk)
98                          ptau(jl, jk) = znorm(jl)*zdz2n
99                       ELSE
100                          ptau(jl, jk) = znorm(jl)*zdz2(jl, jk)
101                       END IF
102                       ptau(jl, jk) = min(ptau(jl, jk), ptau(jl, jk+1))
103                  END IF                  END IF
104                  ptau(jl,jk) = min(ptau(jl,jk),ptau(jl,jk+1))               END IF
               END IF  
             END IF  
105            END IF            END IF
106  431     CONTINUE         end DO
107        end DO
108    
109  430   CONTINUE      ! REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL
 440   CONTINUE  
110    
111  !  REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL      DO jl = 1, klon
112           IF (ktest(jl)==1) THEN
113              ztau(jl, kkcrith(jl)) = ptau(jl, kkcrith(jl))
114              ztau(jl, nstra) = ptau(jl, nstra)
115           END IF
116        end DO
117    
118  !     DO 530 ji=1,kgwd      DO jk = 1, klev
119  !     jl=kdx(ji)         DO jl = 1, klon
 !  Modif vectorisation 02/04/2004  
       DO 530 jl = 1, klon  
         IF (ktest(jl)==1) THEN  
           ztau(jl,kkcrith(jl)) = ptau(jl,kkcrith(jl))  
           ztau(jl,nstra) = ptau(jl,nstra)  
         END IF  
 530   CONTINUE  
   
       DO 531 jk = 1, klev  
   
 !     DO 532 ji=1,kgwd  
 !     jl=kdx(ji)  
 !  Modif vectorisation 02/04/2004  
         DO 532 jl = 1, klon  
120            IF (ktest(jl)==1) THEN            IF (ktest(jl)==1) THEN
121                 IF (jk>kkcrith(jl)) THEN
122                    zdelp = paphm1(jl, jk) - paphm1(jl, klev+1)
123              IF (jk>kkcrith(jl)) THEN                  zdelpt = paphm1(jl, kkcrith(jl)) - paphm1(jl, klev+1)
124                    ptau(jl, jk) = ztau(jl, klev+1) &
125                zdelp = paphm1(jl,jk) - paphm1(jl,klev+1)                       + (ztau(jl, kkcrith(jl)) - ztau(jl, klev+1))*zdelp/zdelpt
126                zdelpt = paphm1(jl,kkcrith(jl)) - paphm1(jl,klev+1)               END IF
               ptau(jl,jk) = ztau(jl,klev+1) + (ztau(jl,kkcrith(jl))-ztau(jl, &  
                 klev+1))*zdelp/zdelpt  
   
             END IF  
   
127            END IF            END IF
128  532     CONTINUE         end DO
   
 !  REORGANISATION IN THE STRATOSPHERE  
129    
130  !     DO 533 ji=1,kgwd         ! REORGANISATION IN THE STRATOSPHERE
131  !     jl=kdx(ji)         DO jl = 1, klon
 !  Modif vectorisation 02/04/2004  
         DO 533 jl = 1, klon  
132            IF (ktest(jl)==1) THEN            IF (ktest(jl)==1) THEN
133                 IF (jk < nstra) THEN
134                    zdelp = paphm1(jl, nstra)
135              IF (jk<nstra) THEN                  zdelpt = paphm1(jl, jk)
136                    ptau(jl, jk) = ztau(jl, nstra) * zdelpt / zdelp
137                zdelp = paphm1(jl,nstra)               END IF
               zdelpt = paphm1(jl,jk)  
               ptau(jl,jk) = ztau(jl,nstra)*zdelpt/zdelp  
   
             END IF  
   
138            END IF            END IF
139  533     CONTINUE         end DO
   
 ! REORGANISATION IN THE TROPOSPHERE  
140    
141  !      DO 534 ji=1,kgwd         ! REORGANISATION IN THE TROPOSPHERE
142  !      jl=kdx(ji)         DO jl = 1, klon
 !  Modif vectorisation 02/04/2004  
         DO 534 jl = 1, klon  
143            IF (ktest(jl)==1) THEN            IF (ktest(jl)==1) THEN
144                 IF (jk<kkcrith(jl) .AND. jk > nstra) THEN
145                    zdelp = paphm1(jl, jk) - paphm1(jl, kkcrith(jl))
146              IF (jk<kkcrith(jl) .AND. jk>nstra) THEN                  zdelpt = paphm1(jl, nstra) - paphm1(jl, kkcrith(jl))
147                    ptau(jl, jk) = ztau(jl, kkcrith(jl)) &
148                zdelp = paphm1(jl,jk) - paphm1(jl,kkcrith(jl))                       + (ztau(jl, nstra) - ztau(jl, kkcrith(jl)))*zdelp/zdelpt
149                zdelpt = paphm1(jl,nstra) - paphm1(jl,kkcrith(jl))               END IF
               ptau(jl,jk) = ztau(jl,kkcrith(jl)) + (ztau(jl,nstra)-ztau(jl, &  
                 kkcrith(jl)))*zdelp/zdelpt  
   
             END IF  
150            END IF            END IF
151  534     CONTINUE         end DO
152        end DO
   
 531   CONTINUE  
153    
154      END SUBROUTINE gwprofil
155    
156        RETURN  end module gwprofil_m
     END  
Legend:
Removed from v.23  
changed lines
  Added in v.328
  ViewVC Help
Powered by ViewVC 1.1.21