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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 178 - (show annotations)
Fri Mar 11 18:47:26 2016 UTC (8 years, 2 months ago) by guez
File size: 6887 byte(s)
Moved variables date0, deltat, datasz_max, ncvar_ids, point, buff_pos,
buffer, regular from module histcom_var to modules where they are
defined.

Removed procedure ioipslmpp, useless for a sequential program.

Added argument datasz_max to histwrite_real (to avoid circular
dependency with histwrite).

Removed useless variables and computations everywhere.

Changed real litteral constants from default kind to double precision
in lwb, lwu, lwvn, sw1s, swtt, swtt1, swu.

Removed unused arguments: paer of sw, sw1s, sw2s, swclr; pcldsw of
sw1s, sw2s; pdsig, prayl of swr; co2_ppm of clmain, clqh; tsol of
transp_lay; nsrf of screenp; kcrit and kknu of gwstress; pstd of
orosetup.

Added output of relative humidity.

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

  ViewVC Help
Powered by ViewVC 1.1.21