/[lmdze]/trunk/Sources/phylmd/Orography/orodrag.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/Orography/orodrag.f

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

revision 246 by guez, Fri Mar 11 18:47:26 2016 UTC revision 247 by guez, Fri Jan 5 14:45:45 2018 UTC
# Line 1  Line 1 
1      SUBROUTINE orodrag(nlon,nlev,ktest,ptsphy,paphm1,papm1,pgeom1, &  module orodrag_m
         ptm1,pum1,pvm1,pmea,pstd,psig,pgamma,ptheta,ppic,pval &  
         ,pulow,pvlow,pvom,pvol,pte)  
2    
3        USE dimens_m    IMPLICIT NONE
       USE dimphy  
       use gwstress_m, only: gwstress  
       USE suphec_m  
       USE yoegwd  
       use gwprofil_m, only: gwprofil  
       use orosetup_m, only: orosetup  
       IMPLICIT NONE  
4    
5    contains
6    
7      SUBROUTINE orodrag(nlon, nlev, ktest, ptsphy, paphm1, papm1, pgeom1, ptm1, &
8           pum1, pvm1, pmea, pstd, psig, pgamma, ptheta, ppic, pval, pulow, &
9           pvlow, pvom, pvol, pte)
10    
11  !**** *gwdrag* - does the gravity wave parametrization.      USE dimens_m
12        USE dimphy
13        use gwstress_m, only: gwstress
14        USE suphec_m
15        USE yoegwd
16        use gwprofil_m, only: gwprofil
17        use orosetup_m, only: orosetup
18    
19  !     purpose.      !**** *gwdrag* - does the gravity wave parametrization.
 !     --------  
20    
21  !          this routine computes the physical tendencies of the      ! purpose.
 !     prognostic variables u,v  and t due to  vertical transports by  
 !     subgridscale orographically excited gravity waves  
22    
23  !**   interface.      ! this routine computes the physical tendencies of the
24  !     ----------      ! prognostic variables u, v and t due to vertical transports by
25  !          called from *callpar*.      ! subgridscale orographically excited gravity waves
26    
27  !          the routine takes its input from the long-term storage:      !** interface.
 !          u,v,t and p at t-1.  
28    
29  !        explicit arguments :      ! called from *callpar*.
 !        --------------------  
 !     ==== inputs ===  
 !     ==== outputs ===  
30    
31  !        implicit arguments :   none      ! the routine takes its input from the long-term storage:
32  !        --------------------      ! u, v, t and p at t-1.
33    
34  !      implicit logical (l)      ! explicit arguments :
35    
36  !     method.      ! ==== inputs ===
37  !     -------      ! ==== outputs ===
38    
39  !     reference.      ! implicit arguments : none
 !     ----------  
40    
41  !     author.      ! implicit logical (l)
 !     -------  
 !     m.miller + b.ritter   e.c.m.w.f.     15/06/86.  
42    
43  !     f.lott + m. miller    e.c.m.w.f.     22/11/94      ! method.
 !-----------------------------------------------------------------------  
44    
45        ! reference.
46    
47  !-----------------------------------------------------------------------      ! author.
48    
49  !*       0.1   arguments      ! m.miller + b.ritter e.c.m.w.f. 15/06/86.
 !              ---------  
50    
51        ! f.lott + m. miller e.c.m.w.f. 22/11/94
52    
53        INTEGER nlon, nlev      !* 0.1 arguments
       INTEGER jl, ilevp1, jk, ji  
       REAL zdelp, ztemp, zforc, ztend  
       REAL rover, zb, zc, zconb, zabsv  
       REAL zzd1, ratio, zbet, zust, zvst, zdis  
       REAL pte(nlon,nlev), pvol(nlon,nlev), pvom(nlon,nlev), pulow(klon), &  
         pvlow(klon)  
       REAL pum1(nlon,nlev), pvm1(nlon,nlev), ptm1(nlon,nlev), pmea(nlon)  
       REAL, INTENT (IN) :: pstd(nlon)  
       REAL, INTENT (IN) :: psig(nlon)  
       REAL pgamma(nlon), ptheta(nlon), ppic(nlon), pval(nlon), &  
         pgeom1(nlon,nlev), papm1(nlon,nlev), paphm1(nlon,nlev+1)  
54    
55        INTEGER ktest(nlon)      INTEGER nlon, nlev
56  !-----------------------------------------------------------------------      INTEGER jl, ilevp1, jk, ji
57        REAL zdelp, ztemp, zforc, ztend
58        REAL rover, zb, zc, zconb, zabsv
59        REAL zzd1, ratio, zbet, zust, zvst, zdis
60        REAL pte(nlon, nlev), pvol(nlon, nlev), pvom(nlon, nlev), pulow(klon), &
61             pvlow(klon)
62        REAL pum1(nlon, nlev), pvm1(nlon, nlev), ptm1(nlon, nlev), pmea(nlon)
63        REAL, INTENT (IN) :: pstd(nlon)
64        REAL, INTENT (IN) :: psig(nlon)
65        REAL pgamma(nlon), ptheta(nlon), ppic(nlon), pval(nlon), &
66             pgeom1(nlon, nlev), papm1(nlon, nlev), paphm1(nlon, nlev+1)
67    
68  !*       0.2   local arrays      INTEGER ktest(nlon)
 !              ------------  
       INTEGER icrit(klon), ikcrith(klon), ikenvh(klon), &  
         iknu(klon), iknu2(klon), ikcrit(klon)  
69    
70        REAL ztau(klon,klev+1), zstab(klon,klev+1), &      !* 0.2 local arrays
         zvph(klon,klev+1), zrho(klon,klev+1), zri(klon,klev+1), &  
         zpsi(klon,klev+1), zzdep(klon,klev)  
       REAL zdudt(klon), zdvdt(klon), zvidis(klon), &  
         znu(klon), zd1(klon), zd2(klon), zdmod(klon)  
       REAL ztmst  
       REAL, INTENT (IN) :: ptsphy  
71    
72  !------------------------------------------------------------------      INTEGER icrit(klon), ikcrith(klon), ikenvh(klon), &
73             iknu(klon), iknu2(klon), ikcrit(klon)
74    
75  !*         1.    initialization      REAL ztau(klon, klev+1), zstab(klon, klev+1), &
76  !                --------------           zvph(klon, klev+1), zrho(klon, klev+1), zri(klon, klev+1), &
77             zpsi(klon, klev+1), zzdep(klon, klev)
78        REAL zdudt(klon), zdvdt(klon), zvidis(klon), &
79             znu(klon), zd1(klon), zd2(klon), zdmod(klon)
80        REAL ztmst
81        REAL, INTENT (IN) :: ptsphy
82    
83  !*         1.1   computational constants      !------------------------------------------------------------------
 !                -----------------------  
84    
85        ztmst = ptsphy      !* 1. initialization
 !     ------------------------------------------------------------------  
86    
87  !*         1.3   check whether row contains point for printing      !* 1.1 computational constants
 !                ---------------------------------------------  
88    
89  !*         2.     precompute basic state variables.      ztmst = ptsphy
 !*                ---------- ----- ----- ----------  
 !*                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.  
90    
91        !* 1.3 check whether row contains point for printing
92    
93        CALL orosetup(nlon,ktest,ikcrit,ikcrith,icrit,ikenvh,iknu,iknu2,paphm1, &      !* 2. precompute basic state variables.
         papm1,pum1,pvm1,ptm1,pgeom1,zrho,zri,zstab,ztau,zvph,zpsi,zzdep, &  
         pulow,pvlow,ptheta,pgamma,pmea,ppic,pval,znu,zd1,zd2,zdmod)  
94    
95        !* define low level wind, project winds in plane of
96        !* low level wind, determine sector in which to take
97        !* the variance and set indicator for critical levels.
98    
99        CALL orosetup(nlon, ktest, ikcrit, ikcrith, icrit, ikenvh, iknu, iknu2, &
100             paphm1, papm1, pum1, pvm1, ptm1, pgeom1, zrho, zri, zstab, ztau, &
101             zvph, zpsi, zzdep, pulow, pvlow, ptheta, pgamma, pmea, ppic, pval, &
102             znu, zd1, zd2, zdmod)
103    
104  !***********************************************************      !* 3. compute low level stresses using subcritical and
105        !* supercritical forms.computes anisotropy coefficient
106        !* as measure of orographic twodimensionality.
107    
108        CALL gwstress(nlon, nlev, ktest, ikenvh, zrho, zstab, zvph, pstd, &
109             psig, pmea, ppic, ztau, pgeom1, zdmod)
110    
111  !*         3.      compute low level stresses using subcritical and      !* 4. compute stress profile.
 !*                 supercritical forms.computes anisotropy coefficient  
 !*                 as measure of orographic twodimensionality.  
112    
113        CALL gwstress(nlon,nlev,ktest,ikenvh,zrho,zstab,zvph,pstd, &      CALL gwprofil(nlon, nlev, ktest, ikcrith, icrit, paphm1, zrho, zstab, &
114          psig,pmea,ppic,ztau,pgeom1,zdmod)           zvph, zri, ztau, zdmod, psig, pstd)
115    
116        !* 5. compute tendencies.
117    
118  !*         4.      compute stress profile.      ! explicit solution at all levels for the gravity wave
119  !*                 ------- ------ --------      ! implicit solution for the blocked levels
120    
121        CALL gwprofil(nlon,nlev,ktest,ikcrith,icrit,paphm1,zrho,zstab, &      DO jl = 1, klon
122          zvph,zri,ztau,zdmod,psig,pstd)         zvidis(jl) = 0.0
123           zdudt(jl) = 0.0
124           zdvdt(jl) = 0.0
125        end DO
126    
127        ilevp1 = klev + 1
128    
129  !*         5.      compute tendencies.      DO jk = 1, klev
 !*                 -------------------  
130    
131  !  explicit solution at all levels for the gravity wave         ! Modif vectorisation 02/04/2004
132  !  implicit solution for the blocked levels         DO ji = 1, klon
   
       DO 510 jl = 1, klon  
         zvidis(jl) = 0.0  
         zdudt(jl) = 0.0  
         zdvdt(jl) = 0.0  
 510   CONTINUE  
   
       ilevp1 = klev + 1  
   
   
       DO 524 jk = 1, klev  
   
   
 !  Modif vectorisation 02/04/2004  
         DO 523 ji = 1, klon  
133            IF (ktest(ji)==1) THEN            IF (ktest(ji)==1) THEN
134    
135              zdelp = paphm1(ji,jk+1) - paphm1(ji,jk)               zdelp = paphm1(ji, jk+1) - paphm1(ji, jk)
136              ztemp = -rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)*zdelp)               ztemp = -rg*(ztau(ji, jk+1)-ztau(ji, jk))/(zvph(ji, ilevp1)*zdelp)
137              zdudt(ji) = (pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)               zdudt(ji) = (pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
138              zdvdt(ji) = (pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)               zdvdt(ji) = (pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
139    
140  ! controle des overshoots:               ! controle des overshoots:
141    
142              zforc = sqrt(zdudt(ji)**2+zdvdt(ji)**2) + 1.E-12               zforc = sqrt(zdudt(ji)**2+zdvdt(ji)**2) + 1.E-12
143              ztend = sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst + 1.E-12               ztend = sqrt(pum1(ji, jk)**2+pvm1(ji, jk)**2)/ztmst + 1.E-12
144              rover = 0.25               rover = 0.25
145              IF (zforc>=rover*ztend) THEN               IF (zforc>=rover*ztend) THEN
146                zdudt(ji) = rover*ztend/zforc*zdudt(ji)                  zdudt(ji) = rover*ztend/zforc*zdudt(ji)
147                zdvdt(ji) = rover*ztend/zforc*zdvdt(ji)                  zdvdt(ji) = rover*ztend/zforc*zdvdt(ji)
148              END IF               END IF
149    
150  ! fin du controle des overshoots               ! fin du controle des overshoots
151    
152              IF (jk>=ikenvh(ji)) THEN               IF (jk>=ikenvh(ji)) THEN
153                zb = 1.0 - 0.18*pgamma(ji) - 0.04*pgamma(ji)**2                  zb = 1.0 - 0.18*pgamma(ji) - 0.04*pgamma(ji)**2
154                zc = 0.48*pgamma(ji) + 0.3*pgamma(ji)**2                  zc = 0.48*pgamma(ji) + 0.3*pgamma(ji)**2
155                zconb = 2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji))                  zconb = 2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji))
156                zabsv = sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2.                  zabsv = sqrt(pum1(ji, jk)**2+pvm1(ji, jk)**2)/2.
157                zzd1 = zb*cos(zpsi(ji,jk))**2 + zc*sin(zpsi(ji,jk))**2                  zzd1 = zb*cos(zpsi(ji, jk))**2 + zc*sin(zpsi(ji, jk))**2
158                ratio = (cos(zpsi(ji,jk))**2+pgamma(ji)*sin(zpsi(ji, &                  ratio = (cos(zpsi(ji, jk))**2+pgamma(ji)*sin(zpsi(ji, &
159                  jk))**2)/(pgamma(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)                       jk))**2)/(pgamma(ji)*cos(zpsi(ji, jk))**2+sin(zpsi(ji, jk))**2)
160                zbet = max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv                  zbet = max(0., 2.-1./ratio)*zconb*zzdep(ji, jk)*zzd1*zabsv
161    
162  ! simplement oppose au vent                  ! simplement oppose au vent
163    
164                zdudt(ji) = -pum1(ji,jk)/ztmst                  zdudt(ji) = -pum1(ji, jk)/ztmst
165                zdvdt(ji) = -pvm1(ji,jk)/ztmst                  zdvdt(ji) = -pvm1(ji, jk)/ztmst
166    
167  !  projection dans la direction de l'axe principal de l'orographie                  ! projection dans la direction de l'axe principal de l'orographie
168  !mod     zdudt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.)                  !mod zdudt(ji)=-(pum1(ji, jk)*cos(ptheta(ji)*rpi/180.)
169  !mod *              +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.))                  !mod * +pvm1(ji, jk)*sin(ptheta(ji)*rpi/180.))
170  !mod *              *cos(ptheta(ji)*rpi/180.)/ztmst                  !mod * *cos(ptheta(ji)*rpi/180.)/ztmst
171  !mod     zdvdt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.)                  !mod zdvdt(ji)=-(pum1(ji, jk)*cos(ptheta(ji)*rpi/180.)
172  !mod *              +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.))                  !mod * +pvm1(ji, jk)*sin(ptheta(ji)*rpi/180.))
173  !mod *              *sin(ptheta(ji)*rpi/180.)/ztmst                  !mod * *sin(ptheta(ji)*rpi/180.)/ztmst
174                zdudt(ji) = zdudt(ji)*(zbet/(1.+zbet))                  zdudt(ji) = zdudt(ji)*(zbet/(1.+zbet))
175                zdvdt(ji) = zdvdt(ji)*(zbet/(1.+zbet))                  zdvdt(ji) = zdvdt(ji)*(zbet/(1.+zbet))
176              END IF               END IF
177              pvom(ji,jk) = zdudt(ji)               pvom(ji, jk) = zdudt(ji)
178              pvol(ji,jk) = zdvdt(ji)               pvol(ji, jk) = zdvdt(ji)
179              zust = pum1(ji,jk) + ztmst*zdudt(ji)               zust = pum1(ji, jk) + ztmst*zdudt(ji)
180              zvst = pvm1(ji,jk) + ztmst*zdvdt(ji)               zvst = pvm1(ji, jk) + ztmst*zdvdt(ji)
181              zdis = 0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2)               zdis = 0.5*(pum1(ji, jk)**2+pvm1(ji, jk)**2-zust**2-zvst**2)
182              zvidis(ji) = zvidis(ji) + zdis*zdelp               zvidis(ji) = zvidis(ji) + zdis*zdelp
183    
184  !  ENCORE UN TRUC POUR EVITER LES EXPLOSIONS               ! ENCORE UN TRUC POUR EVITER LES EXPLOSIONS
185    
186              pte(ji,jk) = 0.0               pte(ji, jk) = 0.0
187    
188            END IF            END IF
189  523     CONTINUE         end DO
190    
191  524   CONTINUE      end DO
192    
193        RETURN
194      END SUBROUTINE orodrag
195    
196        RETURN  end module orodrag_m
     END  
Legend:
Removed from v.246  
changed lines
  Added in v.247
  ViewVC Help
Powered by ViewVC 1.1.21