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

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

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

trunk/Sources/phylmd/Orography/orodrag.f revision 150 by guez, Thu Jun 18 13:49:26 2015 UTC trunk/phylmd/Orography/orodrag.f revision 254 by guez, Mon Feb 5 10:39:38 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 suphec_m  
       USE yoegwd  
       use gwprofil_m, only: gwprofil  
       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  !     externals.      ! implicit arguments : none
 !     ----------  
       INTEGER ismin, ismax  
       EXTERNAL ismin, ismax  
40    
41  !     reference.      ! implicit logical (l)
 !     ----------  
42    
43  !     author.      ! method.
 !     -------  
 !     m.miller + b.ritter   e.c.m.w.f.     15/06/86.  
44    
45  !     f.lott + m. miller    e.c.m.w.f.     22/11/94      ! reference.
 !-----------------------------------------------------------------------  
46    
47        ! author.
48    
49  !-----------------------------------------------------------------------      ! m.miller + b.ritter e.c.m.w.f. 15/06/86.
50    
51  !*       0.1   arguments      ! f.lott + m. miller e.c.m.w.f. 22/11/94
 !              ---------  
52    
53        !* 0.1 arguments
54    
55        INTEGER nlon, nlev, klevm1      INTEGER nlon, nlev
56        INTEGER jl, ilevp1, jk, ji      INTEGER jl, ilevp1, jk, ji
57        REAL zdelp, ztemp, zforc, ztend      REAL zdelp, ztemp, zforc, ztend
58        REAL rover, zb, zc, zconb, zabsv      REAL rover, zb, zc, zconb, zabsv
59        REAL zzd1, ratio, zbet, zust, zvst, zdis      REAL zzd1, ratio, zbet, zust, zvst, zdis
60        REAL pte(nlon,nlev), pvol(nlon,nlev), pvom(nlon,nlev), pulow(klon), &      REAL pte(nlon, nlev), pvol(nlon, nlev), pvom(nlon, nlev), pulow(klon), &
61          pvlow(klon)           pvlow(klon)
62        REAL pum1(nlon,nlev), pvm1(nlon,nlev), ptm1(nlon,nlev), pmea(nlon)      REAL pum1(nlon, nlev), pvm1(nlon, nlev), ptm1(nlon, nlev), pmea(nlon)
63        REAL, INTENT (IN) :: pstd(nlon)      REAL, INTENT (IN) :: pstd(nlon)
64        REAL, INTENT (IN) :: psig(nlon)      REAL, INTENT (IN) :: psig(nlon)
65        REAL pgamma(nlon), ptheta(nlon), ppic(nlon), pval(nlon), &      REAL pgamma(nlon), ptheta(nlon), ppic(nlon), pval(nlon), &
66          pgeom1(nlon,nlev), papm1(nlon,nlev), paphm1(nlon,nlev+1)           pgeom1(nlon, nlev), papm1(nlon, nlev), paphm1(nlon, nlev+1)
67    
68        INTEGER ktest(nlon)      INTEGER ktest(nlon)
 !-----------------------------------------------------------------------  
69    
70  !*       0.2   local arrays      !* 0.2 local arrays
 !              ------------  
       INTEGER icrit(klon), ikcrith(klon), ikenvh(klon), &  
         iknu(klon), iknu2(klon), ikcrit(klon)  
71    
72        REAL ztau(klon,klev+1), zstab(klon,klev+1), &      INTEGER icrit(klon), ikcrith(klon), ikenvh(klon), &
73          zvph(klon,klev+1), zrho(klon,klev+1), zri(klon,klev+1), &           iknu(klon), iknu2(klon), ikcrit(klon)
         zpsi(klon,klev+1), zzdep(klon,klev)  
       REAL zdudt(klon), zdvdt(klon), zdtdt(klon), zdedt(klon), zvidis(klon), &  
         znu(klon), zd1(klon), zd2(klon), zdmod(klon)  
       REAL ztmst, zrtmst  
       REAL, INTENT (IN) :: ptsphy  
74    
75  !------------------------------------------------------------------      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.    initialization      !------------------------------------------------------------------
 !                --------------  
84    
85  !*         1.1   computational constants      !* 1. initialization
 !                -----------------------  
86    
87        klevm1 = klev - 1      !* 1.1 computational constants
       ztmst = ptsphy  
       zrtmst = 1./ztmst  
 !     ------------------------------------------------------------------  
88    
89  !*         1.3   check whether row contains point for printing      ztmst = ptsphy
 !                ---------------------------------------------  
90    
91  !*         2.     precompute basic state variables.      !* 1.3 check whether row contains point for printing
 !*                ---------- ----- ----- ----------  
 !*                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.  
92    
93        !* 2. precompute basic state variables.
94    
95        CALL orosetup(nlon,ktest,ikcrit,ikcrith,icrit,ikenvh,iknu,iknu2,paphm1, &      !* define low level wind, project winds in plane of
96          papm1,pum1,pvm1,ptm1,pgeom1,pstd,zrho,zri,zstab,ztau,zvph,zpsi,zzdep, &      !* low level wind, determine sector in which to take
97          pulow,pvlow,ptheta,pgamma,pmea,ppic,pval,znu,zd1,zd2,zdmod)      !* 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        !* 4. compute stress profile.
112    
113  !*         3.      compute low level stresses using subcritical and      CALL gwprofil(nlon, nlev, ktest, ikcrith, icrit, paphm1, zrho, zstab, &
114  !*                 supercritical forms.computes anisotropy coefficient           zvph, zri, ztau, zdmod, psig, pstd)
 !*                 as measure of orographic twodimensionality.  
115    
116        CALL gwstress(nlon,nlev,ktest,icrit,ikenvh,iknu,zrho,zstab,zvph,pstd, &      !* 5. compute tendencies.
         psig,pmea,ppic,ztau,pgeom1,zdmod)  
117    
118        ! explicit solution at all levels for the gravity wave
119        ! implicit solution for the blocked levels
120    
121  !*         4.      compute stress profile.      DO jl = 1, klon
122  !*                 ------- ------ --------         zvidis(jl) = 0.0
123           zdudt(jl) = 0.0
124           zdvdt(jl) = 0.0
125        end DO
126    
127        CALL gwprofil(nlon,nlev,ktest,ikcrith,icrit,paphm1,zrho,zstab, &      ilevp1 = klev + 1
         zvph,zri,ztau,zdmod,psig,pstd)  
128    
129        DO jk = 1, klev
130    
131  !*         5.      compute tendencies.         ! Modif vectorisation 02/04/2004
132  !*                 -------------------         DO ji = 1, klon
   
 !  explicit solution at all levels for the gravity wave  
 !  implicit solution for the blocked levels  
   
       DO 510 jl = 1, klon  
         zvidis(jl) = 0.0  
         zdudt(jl) = 0.0  
         zdvdt(jl) = 0.0  
         zdtdt(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              zdedt(ji) = zdis/ztmst               zvidis(ji) = zvidis(ji) + zdis*zdelp
             zvidis(ji) = zvidis(ji) + zdis*zdelp  
             zdtdt(ji) = zdedt(ji)/rcpd  
 !     pte(ji,jk)=zdtdt(ji)  
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.150  
changed lines
  Added in v.254
  ViewVC Help
Powered by ViewVC 1.1.21