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

revision 82 by guez, Wed Mar 5 14:57:53 2014 UTC revision 254 by guez, Mon Feb 5 10:39:38 2018 UTC
# Line 1  Line 1 
1      SUBROUTINE orodrag(nlon,nlev,kgwd,kdx,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 kgwd, 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 kdx(nlon), ktest(nlon)      INTEGER ktest(nlon)
 !-----------------------------------------------------------------------  
69    
70  !*       0.2   local arrays      !* 0.2 local arrays
 !              ------------  
       INTEGER isect(klon), icrit(klon), ikcrith(klon), ikenvh(klon), &  
         iknu(klon), iknu2(klon), ikcrit(klon), ikhlim(klon)  
71    
72        REAL ztau(klon,klev+1), ztauf(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  100   CONTINUE      !* 1. initialization
86    
87  !     ------------------------------------------------------------------      !* 1.1 computational constants
88    
89  !*         1.1   computational constants      ztmst = ptsphy
 !                -----------------------  
90    
91  110   CONTINUE      !* 1.3 check whether row contains point for printing
92    
93  !     ztmst=twodt      !* 2. precompute basic state variables.
 !     if(nstep.eq.nstart) ztmst=0.5*twodt  
       klevm1 = klev - 1  
       ztmst = ptsphy  
       zrtmst = 1./ztmst  
 !     ------------------------------------------------------------------  
94    
95  120   CONTINUE      !* 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  !*         1.3   check whether row contains point for printing      !* 3. compute low level stresses using subcritical and
105  !                ---------------------------------------------      !* supercritical forms.computes anisotropy coefficient
106        !* as measure of orographic twodimensionality.
107    
108  130   CONTINUE      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  !*         2.     precompute basic state variables.      CALL gwprofil(nlon, nlev, ktest, ikcrith, icrit, paphm1, zrho, zstab, &
114  !*                ---------- ----- ----- ----------           zvph, zri, ztau, zdmod, psig, pstd)
 !*                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.  
115    
116  200   CONTINUE      !* 5. compute tendencies.
117    
118        ! explicit solution at all levels for the gravity wave
119        ! implicit solution for the blocked levels
120    
121        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 orosetup(nlon,ktest,ikcrit,ikcrith,icrit,ikenvh,iknu,iknu2,paphm1, &      ilevp1 = klev + 1
         papm1,pum1,pvm1,ptm1,pgeom1,pstd,zrho,zri,zstab,ztau,zvph,zpsi,zzdep, &  
         pulow,pvlow,ptheta,pgamma,pmea,ppic,pval,znu,zd1,zd2,zdmod)  
128    
129        DO jk = 1, klev
130    
131           ! Modif vectorisation 02/04/2004
132  !***********************************************************         DO ji = 1, klon
   
   
 !*         3.      compute low level stresses using subcritical and  
 !*                 supercritical forms.computes anisotropy coefficient  
 !*                 as measure of orographic twodimensionality.  
   
 300   CONTINUE  
   
       CALL gwstress(nlon,nlev,ktest,icrit,ikenvh,iknu,zrho,zstab,zvph,pstd, &  
         psig,pmea,ppic,ztau,pgeom1,zdmod)  
   
   
 !*         4.      compute stress profile.  
 !*                 ------- ------ --------  
   
 400   CONTINUE  
   
   
       CALL gwprofil(nlon,nlev,kgwd,kdx,ktest,ikcrith,icrit,paphm1,zrho,zstab, &  
         zvph,zri,ztau,zdmod,psig,pstd)  
   
   
 !*         5.      compute tendencies.  
 !*                 -------------------  
   
 500   CONTINUE  
   
 !  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  
   
   
 !     do 523 jl=1,kgwd  
 !     ji=kdx(jl)  
 !  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.82  
changed lines
  Added in v.254
  ViewVC Help
Powered by ViewVC 1.1.21