/[lmdze]/trunk/libf/phylmd/Orography/orodrag.f90
ViewVC logotype

Contents of /trunk/libf/phylmd/Orography/orodrag.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 54 - (show annotations)
Tue Dec 6 15:07:04 2011 UTC (12 years, 5 months ago) by guez
File size: 7668 byte(s) 
Removed Numerical Recipes procedure "ran1". Replaced calls to "ran1"
in "inidissip" by calls to intrinsic procedures.

Split file "interface_surf.f90" into a file with a module containing
only variables, "interface_surf", and single-procedure files. Gathered
files into directory "Interface_surf".

Added argument "cdivu" to "gradiv" and "gradiv2", "cdivh" to
"divgrad2" and "divgrad", and "crot" to "nxgraro2" and
"nxgrarot". "dissip" now uses variables "cdivu", "cdivh" and "crot"
from module "inidissip_m", so it can pass them to "gradiv2",
etc. Thanks to this modification, we avoid a circular dependency
betwwen "inidissip.f90" and "gradiv2.f90", etc. The value -1. used by
"gradiv2", for instance, during computation of eigenvalues is not the
value "cdivu" computed by "inidissip".

Extracted procedure "start_inter_3d" from module "startdyn", to its
own module.

In "inidissip", unrolled loop on "ii". I find it clearer now.

Moved variables "matriceun", "matriceus", "matricevn", "matricevs",
"matrinvn" and "matrinvs" from module "parafilt" to module
"inifilr_m". Moved variables "jfiltnu", "jfiltnv", "jfiltsu",
"jfiltsv" from module "coefils" to module "inifilr_m".
1 SUBROUTINE orodrag(nlon,nlev,kgwd,kdx,ktest,ptsphy,paphm1,papm1,pgeom1, &
2 ptm1,pum1,pvm1,pmea,pstd,psig,pgamma,ptheta,ppic,pval &
3 ,pulow,pvlow,pvom,pvol,pte)
4
5 USE dimens_m
6 USE dimphy
7 USE suphec_m
8 USE yoegwd
9 use gwprofil_m, only: gwprofil
10 IMPLICIT NONE
11
12
13
14 !**** *gwdrag* - does the gravity wave parametrization.
15
16 ! purpose.
17 ! --------
18
19 ! this routine computes the physical tendencies of the
20 ! prognostic variables u,v and t due to vertical transports by
21 ! subgridscale orographically excited gravity waves
22
23 !** interface.
24 ! ----------
25 ! called from *callpar*.
26
27 ! the routine takes its input from the long-term storage:
28 ! u,v,t and p at t-1.
29
30 ! explicit arguments :
31 ! --------------------
32 ! ==== inputs ===
33 ! ==== outputs ===
34
35 ! implicit arguments : none
36 ! --------------------
37
38 ! implicit logical (l)
39
40 ! method.
41 ! -------
42
43 ! externals.
44 ! ----------
45 INTEGER ismin, ismax
46 EXTERNAL ismin, ismax
47
48 ! reference.
49 ! ----------
50
51 ! author.
52 ! -------
53 ! m.miller + b.ritter e.c.m.w.f. 15/06/86.
54
55 ! f.lott + m. miller e.c.m.w.f. 22/11/94
56 !-----------------------------------------------------------------------
57
58
59 !-----------------------------------------------------------------------
60
61 !* 0.1 arguments
62 ! ---------
63
64
65 INTEGER nlon, nlev, klevm1
66 INTEGER kgwd, jl, ilevp1, jk, ji
67 REAL zdelp, ztemp, zforc, ztend
68 REAL rover, zb, zc, zconb, zabsv
69 REAL zzd1, ratio, zbet, zust, zvst, zdis
70 REAL pte(nlon,nlev), pvol(nlon,nlev), pvom(nlon,nlev), pulow(klon), &
71 pvlow(klon)
72 REAL pum1(nlon,nlev), pvm1(nlon,nlev), ptm1(nlon,nlev), pmea(nlon)
73 REAL, INTENT (IN) :: pstd(nlon)
74 REAL, INTENT (IN) :: psig(nlon)
75 REAL pgamma(nlon), ptheta(nlon), ppic(nlon), pval(nlon), &
76 pgeom1(nlon,nlev), papm1(nlon,nlev), paphm1(nlon,nlev+1)
77
78 INTEGER kdx(nlon), ktest(nlon)
79 !-----------------------------------------------------------------------
80
81 !* 0.2 local arrays
82 ! ------------
83 INTEGER isect(klon), icrit(klon), ikcrith(klon), ikenvh(klon), &
84 iknu(klon), iknu2(klon), ikcrit(klon), ikhlim(klon)
85
86 REAL ztau(klon,klev+1), ztauf(klon,klev+1), zstab(klon,klev+1), &
87 zvph(klon,klev+1), zrho(klon,klev+1), zri(klon,klev+1), &
88 zpsi(klon,klev+1), zzdep(klon,klev)
89 REAL zdudt(klon), zdvdt(klon), zdtdt(klon), zdedt(klon), zvidis(klon), &
90 znu(klon), zd1(klon), zd2(klon), zdmod(klon)
91 REAL ztmst, zrtmst
92 REAL, INTENT (IN) :: ptsphy
93
94 !------------------------------------------------------------------
95
96 !* 1. initialization
97 ! --------------
98
99 100 CONTINUE
100
101 ! ------------------------------------------------------------------
102
103 !* 1.1 computational constants
104 ! -----------------------
105
106 110 CONTINUE
107
108 ! ztmst=twodt
109 ! if(nstep.eq.nstart) ztmst=0.5*twodt
110 klevm1 = klev - 1
111 ztmst = ptsphy
112 zrtmst = 1./ztmst
113 ! ------------------------------------------------------------------
114
115 120 CONTINUE
116
117 ! ------------------------------------------------------------------
118
119 !* 1.3 check whether row contains point for printing
120 ! ---------------------------------------------
121
122 130 CONTINUE
123
124 ! ------------------------------------------------------------------
125
126 !* 2. precompute basic state variables.
127 !* ---------- ----- ----- ----------
128 !* define low level wind, project winds in plane of
129 !* low level wind, determine sector in which to take
130 !* the variance and set indicator for critical levels.
131
132 200 CONTINUE
133
134
135
136 CALL orosetup(nlon,ktest,ikcrit,ikcrith,icrit,ikenvh,iknu,iknu2,paphm1, &
137 papm1,pum1,pvm1,ptm1,pgeom1,pstd,zrho,zri,zstab,ztau,zvph,zpsi,zzdep, &
138 pulow,pvlow,ptheta,pgamma,pmea,ppic,pval,znu,zd1,zd2,zdmod)
139
140
141
142 !***********************************************************
143
144
145 !* 3. compute low level stresses using subcritical and
146 !* supercritical forms.computes anisotropy coefficient
147 !* as measure of orographic twodimensionality.
148
149 300 CONTINUE
150
151 CALL gwstress(nlon,nlev,ktest,icrit,ikenvh,iknu,zrho,zstab,zvph,pstd, &
152 psig,pmea,ppic,ztau,pgeom1,zdmod)
153
154
155 !* 4. compute stress profile.
156 !* ------- ------ --------
157
158 400 CONTINUE
159
160
161 CALL gwprofil(nlon,nlev,kgwd,kdx,ktest,ikcrith,icrit,paphm1,zrho,zstab, &
162 zvph,zri,ztau,zdmod,psig,pstd)
163
164
165 !* 5. compute tendencies.
166 !* -------------------
167
168 500 CONTINUE
169
170 ! explicit solution at all levels for the gravity wave
171 ! implicit solution for the blocked levels
172
173 DO 510 jl = 1, klon
174 zvidis(jl) = 0.0
175 zdudt(jl) = 0.0
176 zdvdt(jl) = 0.0
177 zdtdt(jl) = 0.0
178 510 CONTINUE
179
180 ilevp1 = klev + 1
181
182
183 DO 524 jk = 1, klev
184
185
186 ! do 523 jl=1,kgwd
187 ! ji=kdx(jl)
188 ! Modif vectorisation 02/04/2004
189 DO 523 ji = 1, klon
190 IF (ktest(ji)==1) THEN
191
192 zdelp = paphm1(ji,jk+1) - paphm1(ji,jk)
193 ztemp = -rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)*zdelp)
194 zdudt(ji) = (pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
195 zdvdt(ji) = (pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
196
197 ! controle des overshoots:
198
199 zforc = sqrt(zdudt(ji)**2+zdvdt(ji)**2) + 1.E-12
200 ztend = sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst + 1.E-12
201 rover = 0.25
202 IF (zforc>=rover*ztend) THEN
203 zdudt(ji) = rover*ztend/zforc*zdudt(ji)
204 zdvdt(ji) = rover*ztend/zforc*zdvdt(ji)
205 END IF
206
207 ! fin du controle des overshoots
208
209 IF (jk>=ikenvh(ji)) THEN
210 zb = 1.0 - 0.18*pgamma(ji) - 0.04*pgamma(ji)**2
211 zc = 0.48*pgamma(ji) + 0.3*pgamma(ji)**2
212 zconb = 2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji))
213 zabsv = sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2.
214 zzd1 = zb*cos(zpsi(ji,jk))**2 + zc*sin(zpsi(ji,jk))**2
215 ratio = (cos(zpsi(ji,jk))**2+pgamma(ji)*sin(zpsi(ji, &
216 jk))**2)/(pgamma(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)
217 zbet = max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv
218
219 ! simplement oppose au vent
220
221 zdudt(ji) = -pum1(ji,jk)/ztmst
222 zdvdt(ji) = -pvm1(ji,jk)/ztmst
223
224 ! projection dans la direction de l'axe principal de l'orographie
225 !mod zdudt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.)
226 !mod * +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.))
227 !mod * *cos(ptheta(ji)*rpi/180.)/ztmst
228 !mod zdvdt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.)
229 !mod * +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.))
230 !mod * *sin(ptheta(ji)*rpi/180.)/ztmst
231 zdudt(ji) = zdudt(ji)*(zbet/(1.+zbet))
232 zdvdt(ji) = zdvdt(ji)*(zbet/(1.+zbet))
233 END IF
234 pvom(ji,jk) = zdudt(ji)
235 pvol(ji,jk) = zdvdt(ji)
236 zust = pum1(ji,jk) + ztmst*zdudt(ji)
237 zvst = pvm1(ji,jk) + ztmst*zdvdt(ji)
238 zdis = 0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2)
239 zdedt(ji) = zdis/ztmst
240 zvidis(ji) = zvidis(ji) + zdis*zdelp
241 zdtdt(ji) = zdedt(ji)/rcpd
242 ! pte(ji,jk)=zdtdt(ji)
243
244 ! ENCORE UN TRUC POUR EVITER LES EXPLOSIONS
245
246 pte(ji,jk) = 0.0
247
248 END IF
249 523 CONTINUE
250
251 524 CONTINUE
252
253
254 RETURN
255 END
  ViewVC Help
Powered by ViewVC 1.1.21