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

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

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/Sources/phylmd/Orography/gwprofil.f revision 158 by guez, Tue Jul 21 14:44:45 2015 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 ilevh, 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      ilevh = klev/3
51        IMPLICIT NONE  
52        DO jl = 1, klon
53           IF (ktest(jl)==1) THEN
54              zoro(jl) = psig(jl)*pdmod(jl)/4./max(pvar(jl), 1.0)
55              ztau(jl, klev+1) = ptau(jl, klev+1)
56           END IF
57  !-----------------------------------------------------------------------      end DO
58    
59  !*       0.1   ARGUMENTS      DO jk = klev, 2, -1
60  !              ---------         ! 4.1 CONSTANT WAVE STRESS UNTIL TOP OF THE
61           ! BLOCKING LAYER.
62        INTEGER nlon, nlev         DO jl = 1, klon
63        INTEGER kkcrith(nlon), kcrit(nlon), kdx(nlon), ktest(nlon)            IF (ktest(jl)==1) THEN
64                 IF (jk>kkcrith(jl)) THEN
65                    ptau(jl, jk) = ztau(jl, klev+1)
66        REAL paphm1(nlon,nlev+1), pstab(nlon,nlev+1), prho(nlon,nlev+1), &               ELSE
67          pvph(nlon,nlev+1), pri(nlon,nlev+1), ptau(nlon,nlev+1)                  ptau(jl, jk) = grahilo*ztau(jl, klev+1)
68                 END IF
69        REAL pdmod(nlon)            END IF
70        REAL, INTENT (IN) :: psig(nlon)         end DO
71        REAL, INTENT (IN) :: pvar(nlon)  
72           ! 4.2 WAVE DISPLACEMENT AT NEXT LEVEL.
73  !-----------------------------------------------------------------------         DO jl = 1, klon
74              IF (ktest(jl)==1) THEN
75  !*       0.2   LOCAL ARRAYS               IF (jk<kkcrith(jl)) THEN
76  !              ------------                  znorm(jl) = gkdrag * prho(jl, jk) * sqrt(pstab(jl, jk)) &
77                         * pvph(jl, jk)* zoro(jl)
78        INTEGER ilevh, ji, kgwd, jl, jk                  zdz2(jl, jk) = ptau(jl, jk+1)/max(znorm(jl), gssec)
79        REAL zsqr, zalfa, zriw, zdel, zb, zalpha, zdz2n               END IF
80        REAL zdelp, zdelpt            END IF
81        REAL zdz2(klon,klev), znorm(klon), zoro(klon)         end DO
82        REAL ztau(klon,klev+1)  
83           ! 4.3 WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
84  !-----------------------------------------------------------------------         ! AND STRESS: BREAKING EVALUATION AND CRITICAL
85           ! LEVEL
86  !*         1.    INITIALIZATION         DO jl = 1, klon
87  !                --------------            IF (ktest(jl)==1) THEN
88                 IF (jk<kkcrith(jl)) THEN
89  !      print *,' entree gwprofil'                  IF ((ptau(jl, jk+1)<gtsec) .OR. (jk<=kcrit(jl))) THEN
90  100   CONTINUE                     ptau(jl, jk) = 0.0
   
   
 !*    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  
91                  ELSE                  ELSE
92                    ptau(jl,jk) = znorm(jl)*zdz2(jl,jk)                     zsqr = sqrt(pri(jl, jk))
93                       zalfa = sqrt(pstab(jl, jk)*zdz2(jl, jk))/pvph(jl, jk)
94                       zriw = pri(jl, jk)*(1.-zalfa)/(1+zalfa*zsqr)**2
95                       IF (zriw<grcrit) THEN
96                          zdel = 4./zsqr/grcrit + 1./grcrit**2 + 4./grcrit
97                          zb = 1./grcrit + 2./zsqr
98                          zalpha = 0.5*(-zb+sqrt(zdel))
99                          zdz2n = (pvph(jl, jk)*zalpha)**2/pstab(jl, jk)
100                          ptau(jl, jk) = znorm(jl)*zdz2n
101                       ELSE
102                          ptau(jl, jk) = znorm(jl)*zdz2(jl, jk)
103                       END IF
104                       ptau(jl, jk) = min(ptau(jl, jk), ptau(jl, jk+1))
105                  END IF                  END IF
106                  ptau(jl,jk) = min(ptau(jl,jk),ptau(jl,jk+1))               END IF
               END IF  
             END IF  
107            END IF            END IF
108  431     CONTINUE         end DO
109        end DO
110    
111  430   CONTINUE      ! REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL
 440   CONTINUE  
112    
113  !  REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL      DO jl = 1, klon
114           IF (ktest(jl)==1) THEN
115              ztau(jl, kkcrith(jl)) = ptau(jl, kkcrith(jl))
116              ztau(jl, nstra) = ptau(jl, nstra)
117           END IF
118        end DO
119    
120  !     DO 530 ji=1,kgwd      DO jk = 1, klev
121  !     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  
122            IF (ktest(jl)==1) THEN            IF (ktest(jl)==1) THEN
123                 IF (jk>kkcrith(jl)) THEN
124                    zdelp = paphm1(jl, jk) - paphm1(jl, klev+1)
125              IF (jk>kkcrith(jl)) THEN                  zdelpt = paphm1(jl, kkcrith(jl)) - paphm1(jl, klev+1)
126                    ptau(jl, jk) = ztau(jl, klev+1) &
127                zdelp = paphm1(jl,jk) - paphm1(jl,klev+1)                       + (ztau(jl, kkcrith(jl)) - ztau(jl, klev+1))*zdelp/zdelpt
128                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  
   
129            END IF            END IF
130  532     CONTINUE         end DO
   
 !  REORGANISATION IN THE STRATOSPHERE  
131    
132  !     DO 533 ji=1,kgwd         ! REORGANISATION IN THE STRATOSPHERE
133  !     jl=kdx(ji)         DO jl = 1, klon
 !  Modif vectorisation 02/04/2004  
         DO 533 jl = 1, klon  
134            IF (ktest(jl)==1) THEN            IF (ktest(jl)==1) THEN
135                 IF (jk < nstra) THEN
136                    zdelp = paphm1(jl, nstra)
137              IF (jk<nstra) THEN                  zdelpt = paphm1(jl, jk)
138                    ptau(jl, jk) = ztau(jl, nstra) * zdelpt / zdelp
139                zdelp = paphm1(jl,nstra)               END IF
               zdelpt = paphm1(jl,jk)  
               ptau(jl,jk) = ztau(jl,nstra)*zdelpt/zdelp  
   
             END IF  
   
140            END IF            END IF
141  533     CONTINUE         end DO
   
 ! REORGANISATION IN THE TROPOSPHERE  
142    
143  !      DO 534 ji=1,kgwd         ! REORGANISATION IN THE TROPOSPHERE
144  !      jl=kdx(ji)         DO jl = 1, klon
 !  Modif vectorisation 02/04/2004  
         DO 534 jl = 1, klon  
145            IF (ktest(jl)==1) THEN            IF (ktest(jl)==1) THEN
146                 IF (jk<kkcrith(jl) .AND. jk > nstra) THEN
147                    zdelp = paphm1(jl, jk) - paphm1(jl, kkcrith(jl))
148              IF (jk<kkcrith(jl) .AND. jk>nstra) THEN                  zdelpt = paphm1(jl, nstra) - paphm1(jl, kkcrith(jl))
149                    ptau(jl, jk) = ztau(jl, kkcrith(jl)) &
150                zdelp = paphm1(jl,jk) - paphm1(jl,kkcrith(jl))                       + (ztau(jl, nstra) - ztau(jl, kkcrith(jl)))*zdelp/zdelpt
151                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  
152            END IF            END IF
153  534     CONTINUE         end DO
154        end DO
   
 531   CONTINUE  
155    
156      END SUBROUTINE gwprofil
157    
158        RETURN  end module gwprofil_m
     END  
Legend:
Removed from v.23  
changed lines
  Added in v.158
  ViewVC Help
Powered by ViewVC 1.1.21