/[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

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