| 1 |
! |
| 2 |
! $Header: /home/cvsroot/LMDZ4/libf/phylmd/orografi.F,v 1.4 2005/12/01 11:27:29 fairhead Exp $ |
| 3 |
! |
| 4 |
SUBROUTINE drag_noro (nlon,nlev,dtime,paprs,pplay, |
| 5 |
e pmea,pstd, psig, pgam, pthe,ppic,pval, |
| 6 |
e kgwd,kdx,ktest, |
| 7 |
e t, u, v, |
| 8 |
s pulow, pvlow, pustr, pvstr, |
| 9 |
s d_t, d_u, d_v) |
| 10 |
c |
| 11 |
use dimens_m |
| 12 |
use dimphy |
| 13 |
use YOMCST |
| 14 |
IMPLICIT none |
| 15 |
c====================================================================== |
| 16 |
c Auteur(s): F.Lott (LMD/CNRS) date: 19950201 |
| 17 |
c Objet: Frottement de la montagne Interface |
| 18 |
c====================================================================== |
| 19 |
c Arguments: |
| 20 |
c dtime---input-R- pas d'integration (s) |
| 21 |
c paprs---input-R-pression pour chaque inter-couche (en Pa) |
| 22 |
c pplay---input-R-pression pour le mileu de chaque couche (en Pa) |
| 23 |
c t-------input-R-temperature (K) |
| 24 |
c u-------input-R-vitesse horizontale (m/s) |
| 25 |
c v-------input-R-vitesse horizontale (m/s) |
| 26 |
c |
| 27 |
c d_t-----output-R-increment de la temperature |
| 28 |
c d_u-----output-R-increment de la vitesse u |
| 29 |
c d_v-----output-R-increment de la vitesse v |
| 30 |
c====================================================================== |
| 31 |
c |
| 32 |
c ARGUMENTS |
| 33 |
c |
| 34 |
INTEGER nlon,nlev |
| 35 |
REAL, intent(in):: dtime |
| 36 |
REAL, intent(in):: paprs(klon,klev+1) |
| 37 |
REAL, intent(in):: pplay(klon,klev) |
| 38 |
REAL pmea(nlon),pstd(nlon),psig(nlon),pgam(nlon),pthe(nlon) |
| 39 |
REAL ppic(nlon),pval(nlon) |
| 40 |
REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon) |
| 41 |
REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev) |
| 42 |
REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev) |
| 43 |
c |
| 44 |
INTEGER i, k, kgwd, kdx(nlon), ktest(nlon) |
| 45 |
c |
| 46 |
c Variables locales: |
| 47 |
c |
| 48 |
REAL zgeom(klon,klev) |
| 49 |
REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev) |
| 50 |
REAL pt(klon,klev), pu(klon,klev), pv(klon,klev) |
| 51 |
REAL papmf(klon,klev),papmh(klon,klev+1) |
| 52 |
c |
| 53 |
c initialiser les variables de sortie (pour securite) |
| 54 |
c |
| 55 |
DO i = 1,klon |
| 56 |
pulow(i) = 0.0 |
| 57 |
pvlow(i) = 0.0 |
| 58 |
pustr(i) = 0.0 |
| 59 |
pvstr(i) = 0.0 |
| 60 |
ENDDO |
| 61 |
DO k = 1, klev |
| 62 |
DO i = 1, klon |
| 63 |
d_t(i,k) = 0.0 |
| 64 |
d_u(i,k) = 0.0 |
| 65 |
d_v(i,k) = 0.0 |
| 66 |
pdudt(i,k)=0.0 |
| 67 |
pdvdt(i,k)=0.0 |
| 68 |
pdtdt(i,k)=0.0 |
| 69 |
ENDDO |
| 70 |
ENDDO |
| 71 |
c |
| 72 |
c preparer les variables d'entree (attention: l'ordre des niveaux |
| 73 |
c verticaux augmente du haut vers le bas) |
| 74 |
c |
| 75 |
DO k = 1, klev |
| 76 |
DO i = 1, klon |
| 77 |
pt(i,k) = t(i,klev-k+1) |
| 78 |
pu(i,k) = u(i,klev-k+1) |
| 79 |
pv(i,k) = v(i,klev-k+1) |
| 80 |
papmf(i,k) = pplay(i,klev-k+1) |
| 81 |
ENDDO |
| 82 |
ENDDO |
| 83 |
DO k = 1, klev+1 |
| 84 |
DO i = 1, klon |
| 85 |
papmh(i,k) = paprs(i,klev-k+2) |
| 86 |
ENDDO |
| 87 |
ENDDO |
| 88 |
DO i = 1, klon |
| 89 |
zgeom(i,klev) = RD * pt(i,klev) |
| 90 |
. * LOG(papmh(i,klev+1)/papmf(i,klev)) |
| 91 |
ENDDO |
| 92 |
DO k = klev-1, 1, -1 |
| 93 |
DO i = 1, klon |
| 94 |
zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0 |
| 95 |
. * LOG(papmf(i,k+1)/papmf(i,k)) |
| 96 |
ENDDO |
| 97 |
ENDDO |
| 98 |
c |
| 99 |
c appeler la routine principale |
| 100 |
c |
| 101 |
CALL orodrag(klon,klev,kgwd,kdx,ktest, |
| 102 |
. dtime, |
| 103 |
. papmh, papmf, zgeom, |
| 104 |
. pt, pu, pv, |
| 105 |
. pmea, pstd, psig, pgam, pthe, ppic,pval, |
| 106 |
. pulow,pvlow, |
| 107 |
. pdudt,pdvdt,pdtdt) |
| 108 |
C |
| 109 |
DO k = 1, klev |
| 110 |
DO i = 1, klon |
| 111 |
d_u(i,klev+1-k) = dtime*pdudt(i,k) |
| 112 |
d_v(i,klev+1-k) = dtime*pdvdt(i,k) |
| 113 |
d_t(i,klev+1-k) = dtime*pdtdt(i,k) |
| 114 |
pustr(i) = pustr(i) |
| 115 |
cIM BUG . +rg*pdudt(i,k)*(papmh(i,k+1)-papmh(i,k)) |
| 116 |
. +pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG |
| 117 |
pvstr(i) = pvstr(i) |
| 118 |
cIM BUG . +rg*pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k)) |
| 119 |
. +pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG |
| 120 |
ENDDO |
| 121 |
ENDDO |
| 122 |
c |
| 123 |
RETURN |
| 124 |
END |
| 125 |
SUBROUTINE orodrag( nlon,nlev |
| 126 |
i , kgwd, kdx, ktest |
| 127 |
r , ptsphy |
| 128 |
r , paphm1,papm1,pgeom1,ptm1,pum1,pvm1 |
| 129 |
r , pmea, pstd, psig, pgamma, ptheta, ppic, pval |
| 130 |
c outputs |
| 131 |
r , pulow,pvlow |
| 132 |
r , pvom,pvol,pte ) |
| 133 |
|
| 134 |
use dimens_m |
| 135 |
use dimphy |
| 136 |
use YOMCST |
| 137 |
use yoegwd |
| 138 |
implicit none |
| 139 |
|
| 140 |
c |
| 141 |
c |
| 142 |
c**** *gwdrag* - does the gravity wave parametrization. |
| 143 |
c |
| 144 |
c purpose. |
| 145 |
c -------- |
| 146 |
c |
| 147 |
c this routine computes the physical tendencies of the |
| 148 |
c prognostic variables u,v and t due to vertical transports by |
| 149 |
c subgridscale orographically excited gravity waves |
| 150 |
c |
| 151 |
c** interface. |
| 152 |
c ---------- |
| 153 |
c called from *callpar*. |
| 154 |
c |
| 155 |
c the routine takes its input from the long-term storage: |
| 156 |
c u,v,t and p at t-1. |
| 157 |
c |
| 158 |
c explicit arguments : |
| 159 |
c -------------------- |
| 160 |
c ==== inputs === |
| 161 |
c ==== outputs === |
| 162 |
c |
| 163 |
c implicit arguments : none |
| 164 |
c -------------------- |
| 165 |
c |
| 166 |
c implicit logical (l) |
| 167 |
c |
| 168 |
c method. |
| 169 |
c ------- |
| 170 |
c |
| 171 |
c externals. |
| 172 |
c ---------- |
| 173 |
integer ismin, ismax |
| 174 |
external ismin, ismax |
| 175 |
c |
| 176 |
c reference. |
| 177 |
c ---------- |
| 178 |
c |
| 179 |
c author. |
| 180 |
c ------- |
| 181 |
c m.miller + b.ritter e.c.m.w.f. 15/06/86. |
| 182 |
c |
| 183 |
c f.lott + m. miller e.c.m.w.f. 22/11/94 |
| 184 |
c----------------------------------------------------------------------- |
| 185 |
c |
| 186 |
c |
| 187 |
c----------------------------------------------------------------------- |
| 188 |
c |
| 189 |
c* 0.1 arguments |
| 190 |
c --------- |
| 191 |
c |
| 192 |
c |
| 193 |
integer nlon, nlev, klevm1 |
| 194 |
integer kgwd, jl, ilevp1, jk, ji |
| 195 |
real zdelp, ztemp, zforc, ztend |
| 196 |
real rover, zb, zc, zconb, zabsv |
| 197 |
real zzd1, ratio, zbet, zust,zvst, zdis |
| 198 |
real pte(nlon,nlev), |
| 199 |
* pvol(nlon,nlev), |
| 200 |
* pvom(nlon,nlev), |
| 201 |
* pulow(klon), |
| 202 |
* pvlow(klon) |
| 203 |
real pum1(nlon,nlev), |
| 204 |
* pvm1(nlon,nlev), |
| 205 |
* ptm1(nlon,nlev), |
| 206 |
* pmea(nlon),pstd(nlon),psig(nlon), |
| 207 |
* pgamma(nlon),ptheta(nlon),ppic(nlon),pval(nlon), |
| 208 |
* pgeom1(nlon,nlev), |
| 209 |
* papm1(nlon,nlev), |
| 210 |
* paphm1(nlon,nlev+1) |
| 211 |
c |
| 212 |
integer kdx(nlon),ktest(nlon) |
| 213 |
c----------------------------------------------------------------------- |
| 214 |
c |
| 215 |
c* 0.2 local arrays |
| 216 |
c ------------ |
| 217 |
integer isect(klon), |
| 218 |
* icrit(klon), |
| 219 |
* ikcrith(klon), |
| 220 |
* ikenvh(klon), |
| 221 |
* iknu(klon), |
| 222 |
* iknu2(klon), |
| 223 |
* ikcrit(klon), |
| 224 |
* ikhlim(klon) |
| 225 |
c |
| 226 |
real ztau(klon,klev+1), |
| 227 |
$ ztauf(klon,klev+1), |
| 228 |
* zstab(klon,klev+1), |
| 229 |
* zvph(klon,klev+1), |
| 230 |
* zrho(klon,klev+1), |
| 231 |
* zri(klon,klev+1), |
| 232 |
* zpsi(klon,klev+1), |
| 233 |
* zzdep(klon,klev) |
| 234 |
real zdudt(klon), |
| 235 |
* zdvdt(klon), |
| 236 |
* zdtdt(klon), |
| 237 |
* zdedt(klon), |
| 238 |
* zvidis(klon), |
| 239 |
* znu(klon), |
| 240 |
* zd1(klon), |
| 241 |
* zd2(klon), |
| 242 |
* zdmod(klon) |
| 243 |
real ztmst, zrtmst |
| 244 |
real, intent(in):: ptsphy |
| 245 |
c |
| 246 |
c------------------------------------------------------------------ |
| 247 |
c |
| 248 |
c* 1. initialization |
| 249 |
c -------------- |
| 250 |
c |
| 251 |
100 continue |
| 252 |
c |
| 253 |
c ------------------------------------------------------------------ |
| 254 |
c |
| 255 |
c* 1.1 computational constants |
| 256 |
c ----------------------- |
| 257 |
c |
| 258 |
110 continue |
| 259 |
c |
| 260 |
c ztmst=twodt |
| 261 |
c if(nstep.eq.nstart) ztmst=0.5*twodt |
| 262 |
klevm1=klev-1 |
| 263 |
ztmst=ptsphy |
| 264 |
zrtmst=1./ztmst |
| 265 |
c ------------------------------------------------------------------ |
| 266 |
c |
| 267 |
120 continue |
| 268 |
c |
| 269 |
c ------------------------------------------------------------------ |
| 270 |
c |
| 271 |
c* 1.3 check whether row contains point for printing |
| 272 |
c --------------------------------------------- |
| 273 |
c |
| 274 |
130 continue |
| 275 |
c |
| 276 |
c ------------------------------------------------------------------ |
| 277 |
c |
| 278 |
c* 2. precompute basic state variables. |
| 279 |
c* ---------- ----- ----- ---------- |
| 280 |
c* define low level wind, project winds in plane of |
| 281 |
c* low level wind, determine sector in which to take |
| 282 |
c* the variance and set indicator for critical levels. |
| 283 |
c |
| 284 |
200 continue |
| 285 |
c |
| 286 |
c |
| 287 |
c |
| 288 |
call orosetup |
| 289 |
* ( nlon, ktest |
| 290 |
* , ikcrit, ikcrith, icrit, ikenvh,iknu,iknu2 |
| 291 |
* , paphm1, papm1 , pum1 , pvm1 , ptm1 , pgeom1, pstd |
| 292 |
* , zrho , zri , zstab , ztau , zvph , zpsi, zzdep |
| 293 |
* , pulow, pvlow |
| 294 |
* , ptheta,pgamma,pmea,ppic,pval,znu ,zd1, zd2, zdmod ) |
| 295 |
c |
| 296 |
c |
| 297 |
c |
| 298 |
c*********************************************************** |
| 299 |
c |
| 300 |
c |
| 301 |
c* 3. compute low level stresses using subcritical and |
| 302 |
c* supercritical forms.computes anisotropy coefficient |
| 303 |
c* as measure of orographic twodimensionality. |
| 304 |
c |
| 305 |
300 continue |
| 306 |
c |
| 307 |
call gwstress |
| 308 |
* ( nlon , nlev |
| 309 |
* , ktest , icrit, ikenvh, iknu |
| 310 |
* , zrho , zstab, zvph , pstd, psig, pmea, ppic |
| 311 |
* , ztau |
| 312 |
* , pgeom1,zdmod) |
| 313 |
c |
| 314 |
c |
| 315 |
c* 4. compute stress profile. |
| 316 |
c* ------- ------ -------- |
| 317 |
c |
| 318 |
400 continue |
| 319 |
c |
| 320 |
c |
| 321 |
call gwprofil |
| 322 |
* ( nlon , nlev |
| 323 |
* , kgwd , kdx , ktest |
| 324 |
* , ikcrith, icrit |
| 325 |
* , paphm1, zrho , zstab , zvph |
| 326 |
* , zri , ztau |
| 327 |
* , zdmod , psig , pstd) |
| 328 |
c |
| 329 |
c |
| 330 |
c* 5. compute tendencies. |
| 331 |
c* ------------------- |
| 332 |
c |
| 333 |
500 continue |
| 334 |
c |
| 335 |
c explicit solution at all levels for the gravity wave |
| 336 |
c implicit solution for the blocked levels |
| 337 |
|
| 338 |
do 510 jl=1,klon |
| 339 |
zvidis(jl)=0.0 |
| 340 |
zdudt(jl)=0.0 |
| 341 |
zdvdt(jl)=0.0 |
| 342 |
zdtdt(jl)=0.0 |
| 343 |
510 continue |
| 344 |
c |
| 345 |
ilevp1=klev+1 |
| 346 |
c |
| 347 |
c |
| 348 |
do 524 jk=1,klev |
| 349 |
c |
| 350 |
c |
| 351 |
c do 523 jl=1,kgwd |
| 352 |
c ji=kdx(jl) |
| 353 |
c Modif vectorisation 02/04/2004 |
| 354 |
do 523 ji=1,klon |
| 355 |
if(ktest(ji).eq.1) then |
| 356 |
|
| 357 |
zdelp=paphm1(ji,jk+1)-paphm1(ji,jk) |
| 358 |
ztemp=-rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)*zdelp) |
| 359 |
zdudt(ji)=(pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji) |
| 360 |
zdvdt(ji)=(pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji) |
| 361 |
c |
| 362 |
c controle des overshoots: |
| 363 |
c |
| 364 |
zforc=sqrt(zdudt(ji)**2+zdvdt(ji)**2)+1.E-12 |
| 365 |
ztend=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst+1.E-12 |
| 366 |
rover=0.25 |
| 367 |
if(zforc.ge.rover*ztend)then |
| 368 |
zdudt(ji)=rover*ztend/zforc*zdudt(ji) |
| 369 |
zdvdt(ji)=rover*ztend/zforc*zdvdt(ji) |
| 370 |
endif |
| 371 |
c |
| 372 |
c fin du controle des overshoots |
| 373 |
c |
| 374 |
if(jk.ge.ikenvh(ji)) then |
| 375 |
zb=1.0-0.18*pgamma(ji)-0.04*pgamma(ji)**2 |
| 376 |
zc=0.48*pgamma(ji)+0.3*pgamma(ji)**2 |
| 377 |
zconb=2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji)) |
| 378 |
zabsv=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2. |
| 379 |
zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2 |
| 380 |
ratio=(cos(zpsi(ji,jk))**2+pgamma(ji)*sin(zpsi(ji,jk))**2)/ |
| 381 |
* (pgamma(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2) |
| 382 |
zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv |
| 383 |
c |
| 384 |
c simplement oppose au vent |
| 385 |
c |
| 386 |
zdudt(ji)=-pum1(ji,jk)/ztmst |
| 387 |
zdvdt(ji)=-pvm1(ji,jk)/ztmst |
| 388 |
c |
| 389 |
c projection dans la direction de l'axe principal de l'orographie |
| 390 |
cmod zdudt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.) |
| 391 |
cmod * +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.)) |
| 392 |
cmod * *cos(ptheta(ji)*rpi/180.)/ztmst |
| 393 |
cmod zdvdt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.) |
| 394 |
cmod * +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.)) |
| 395 |
cmod * *sin(ptheta(ji)*rpi/180.)/ztmst |
| 396 |
zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet)) |
| 397 |
zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet)) |
| 398 |
end if |
| 399 |
pvom(ji,jk)=zdudt(ji) |
| 400 |
pvol(ji,jk)=zdvdt(ji) |
| 401 |
zust=pum1(ji,jk)+ztmst*zdudt(ji) |
| 402 |
zvst=pvm1(ji,jk)+ztmst*zdvdt(ji) |
| 403 |
zdis=0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2) |
| 404 |
zdedt(ji)=zdis/ztmst |
| 405 |
zvidis(ji)=zvidis(ji)+zdis*zdelp |
| 406 |
zdtdt(ji)=zdedt(ji)/rcpd |
| 407 |
c pte(ji,jk)=zdtdt(ji) |
| 408 |
c |
| 409 |
c ENCORE UN TRUC POUR EVITER LES EXPLOSIONS |
| 410 |
c |
| 411 |
pte(ji,jk)=0.0 |
| 412 |
|
| 413 |
endif |
| 414 |
523 continue |
| 415 |
|
| 416 |
524 continue |
| 417 |
c |
| 418 |
c |
| 419 |
return |
| 420 |
end |
| 421 |
SUBROUTINE orosetup |
| 422 |
* ( nlon , ktest |
| 423 |
* , kkcrit, kkcrith, kcrit |
| 424 |
* , kkenvh, kknu , kknu2 |
| 425 |
* , paphm1, papm1 , pum1 , pvm1 , ptm1 , pgeom1, pstd |
| 426 |
* , prho , pri , pstab , ptau , pvph ,ppsi, pzdep |
| 427 |
* , pulow , pvlow |
| 428 |
* , ptheta, pgamma, pmea, ppic, pval |
| 429 |
* , pnu , pd1 , pd2 ,pdmod ) |
| 430 |
c |
| 431 |
c**** *gwsetup* |
| 432 |
c |
| 433 |
c purpose. |
| 434 |
c -------- |
| 435 |
c |
| 436 |
c** interface. |
| 437 |
c ---------- |
| 438 |
c from *orodrag* |
| 439 |
c |
| 440 |
c explicit arguments : |
| 441 |
c -------------------- |
| 442 |
c ==== inputs === |
| 443 |
c ==== outputs === |
| 444 |
c |
| 445 |
c implicit arguments : none |
| 446 |
c -------------------- |
| 447 |
c |
| 448 |
c method. |
| 449 |
c ------- |
| 450 |
c |
| 451 |
c |
| 452 |
c externals. |
| 453 |
c ---------- |
| 454 |
c |
| 455 |
c |
| 456 |
c reference. |
| 457 |
c ---------- |
| 458 |
c |
| 459 |
c see ecmwf research department documentation of the "i.f.s." |
| 460 |
c |
| 461 |
c author. |
| 462 |
c ------- |
| 463 |
c |
| 464 |
c modifications. |
| 465 |
c -------------- |
| 466 |
c f.lott for the new-gwdrag scheme november 1993 |
| 467 |
c |
| 468 |
c----------------------------------------------------------------------- |
| 469 |
use dimens_m |
| 470 |
use dimphy |
| 471 |
use YOMCST |
| 472 |
use yoegwd |
| 473 |
implicit none |
| 474 |
c |
| 475 |
|
| 476 |
|
| 477 |
c----------------------------------------------------------------------- |
| 478 |
c |
| 479 |
c* 0.1 arguments |
| 480 |
c --------- |
| 481 |
c |
| 482 |
integer nlon |
| 483 |
integer jl, jk |
| 484 |
real zdelp |
| 485 |
|
| 486 |
integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon), |
| 487 |
* ktest(nlon),kkenvh(nlon) |
| 488 |
|
| 489 |
c |
| 490 |
real paphm1(nlon,klev+1),papm1(nlon,klev),pum1(nlon,klev), |
| 491 |
* pvm1(nlon,klev),ptm1(nlon,klev),pgeom1(nlon,klev), |
| 492 |
* prho(nlon,klev+1),pri(nlon,klev+1),pstab(nlon,klev+1), |
| 493 |
* ptau(nlon,klev+1),pvph(nlon,klev+1),ppsi(nlon,klev+1), |
| 494 |
* pzdep(nlon,klev) |
| 495 |
real pulow(nlon),pvlow(nlon),ptheta(nlon),pgamma(nlon),pnu(nlon), |
| 496 |
* pd1(nlon),pd2(nlon),pdmod(nlon) |
| 497 |
real pstd(nlon),pmea(nlon),ppic(nlon),pval(nlon) |
| 498 |
c |
| 499 |
c----------------------------------------------------------------------- |
| 500 |
c |
| 501 |
c* 0.2 local arrays |
| 502 |
c ------------ |
| 503 |
c |
| 504 |
c |
| 505 |
integer ilevm1, ilevm2, ilevh |
| 506 |
real zcons1, zcons2,zcons3, zhgeo |
| 507 |
real zu, zphi, zvt1,zvt2, zst, zvar, zdwind, zwind |
| 508 |
real zstabm, zstabp, zrhom, zrhop, alpha |
| 509 |
real zggeenv, zggeom1,zgvar |
| 510 |
logical lo |
| 511 |
logical ll1(klon,klev+1) |
| 512 |
integer kknu(klon),kknu2(klon),kknub(klon),kknul(klon), |
| 513 |
* kentp(klon),ncount(klon) |
| 514 |
c |
| 515 |
real zhcrit(klon,klev),zvpf(klon,klev), |
| 516 |
* zdp(klon,klev) |
| 517 |
real znorm(klon),zb(klon),zc(klon), |
| 518 |
* zulow(klon),zvlow(klon),znup(klon),znum(klon) |
| 519 |
c |
| 520 |
c ------------------------------------------------------------------ |
| 521 |
c |
| 522 |
c* 1. initialization |
| 523 |
c -------------- |
| 524 |
c |
| 525 |
c print *,' entree gwsetup' |
| 526 |
100 continue |
| 527 |
c |
| 528 |
c ------------------------------------------------------------------ |
| 529 |
c |
| 530 |
c* 1.1 computational constants |
| 531 |
c ----------------------- |
| 532 |
c |
| 533 |
110 continue |
| 534 |
c |
| 535 |
ilevm1=klev-1 |
| 536 |
ilevm2=klev-2 |
| 537 |
ilevh =klev/3 |
| 538 |
c |
| 539 |
zcons1=1./rd |
| 540 |
cold zcons2=g**2/cpd |
| 541 |
zcons2=rg**2/rcpd |
| 542 |
cold zcons3=1.5*api |
| 543 |
zcons3=1.5*rpi |
| 544 |
c |
| 545 |
c |
| 546 |
c ------------------------------------------------------------------ |
| 547 |
c |
| 548 |
c* 2. |
| 549 |
c -------------- |
| 550 |
c |
| 551 |
200 continue |
| 552 |
c |
| 553 |
c ------------------------------------------------------------------ |
| 554 |
c |
| 555 |
c* 2.1 define low level wind, project winds in plane of |
| 556 |
c* low level wind, determine sector in which to take |
| 557 |
c* the variance and set indicator for critical levels. |
| 558 |
c |
| 559 |
c |
| 560 |
c |
| 561 |
do 2001 jl=1,klon |
| 562 |
kknu(jl) =klev |
| 563 |
kknu2(jl) =klev |
| 564 |
kknub(jl) =klev |
| 565 |
kknul(jl) =klev |
| 566 |
pgamma(jl) =max(pgamma(jl),gtsec) |
| 567 |
ll1(jl,klev+1)=.false. |
| 568 |
2001 continue |
| 569 |
c |
| 570 |
c Ajouter une initialisation (L. Li, le 23fev99): |
| 571 |
c |
| 572 |
do jk=klev,ilevh,-1 |
| 573 |
do jl=1,klon |
| 574 |
ll1(jl,jk)= .FALSE. |
| 575 |
ENDDO |
| 576 |
ENDDO |
| 577 |
c |
| 578 |
c* define top of low level flow |
| 579 |
c ---------------------------- |
| 580 |
do 2002 jk=klev,ilevh,-1 |
| 581 |
do 2003 jl=1,klon |
| 582 |
lo=(paphm1(jl,jk)/paphm1(jl,klev+1)).ge.gsigcr |
| 583 |
if(lo) then |
| 584 |
kkcrit(jl)=jk |
| 585 |
endif |
| 586 |
zhcrit(jl,jk)=ppic(jl) |
| 587 |
zhgeo=pgeom1(jl,jk)/rg |
| 588 |
ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk)) |
| 589 |
if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then |
| 590 |
kknu(jl)=jk |
| 591 |
endif |
| 592 |
if(.not.ll1(jl,ilevh))kknu(jl)=ilevh |
| 593 |
2003 continue |
| 594 |
2002 continue |
| 595 |
do 2004 jk=klev,ilevh,-1 |
| 596 |
do 2005 jl=1,klon |
| 597 |
zhcrit(jl,jk)=ppic(jl)-pval(jl) |
| 598 |
zhgeo=pgeom1(jl,jk)/rg |
| 599 |
ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk)) |
| 600 |
if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then |
| 601 |
kknu2(jl)=jk |
| 602 |
endif |
| 603 |
if(.not.ll1(jl,ilevh))kknu2(jl)=ilevh |
| 604 |
2005 continue |
| 605 |
2004 continue |
| 606 |
do 2006 jk=klev,ilevh,-1 |
| 607 |
do 2007 jl=1,klon |
| 608 |
zhcrit(jl,jk)=amax1(ppic(jl)-pmea(jl),pmea(jl)-pval(jl)) |
| 609 |
zhgeo=pgeom1(jl,jk)/rg |
| 610 |
ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk)) |
| 611 |
if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then |
| 612 |
kknub(jl)=jk |
| 613 |
endif |
| 614 |
if(.not.ll1(jl,ilevh))kknub(jl)=ilevh |
| 615 |
2007 continue |
| 616 |
2006 continue |
| 617 |
c |
| 618 |
do 2010 jl=1,klon |
| 619 |
kknu(jl)=min(kknu(jl),nktopg) |
| 620 |
kknu2(jl)=min(kknu2(jl),nktopg) |
| 621 |
kknub(jl)=min(kknub(jl),nktopg) |
| 622 |
kknul(jl)=klev |
| 623 |
2010 continue |
| 624 |
c |
| 625 |
|
| 626 |
210 continue |
| 627 |
c |
| 628 |
c |
| 629 |
cc* initialize various arrays |
| 630 |
c |
| 631 |
do 2107 jl=1,klon |
| 632 |
prho(jl,klev+1) =0.0 |
| 633 |
pstab(jl,klev+1) =0.0 |
| 634 |
pstab(jl,1) =0.0 |
| 635 |
pri(jl,klev+1) =9999.0 |
| 636 |
ppsi(jl,klev+1) =0.0 |
| 637 |
pri(jl,1) =0.0 |
| 638 |
pvph(jl,1) =0.0 |
| 639 |
pulow(jl) =0.0 |
| 640 |
pvlow(jl) =0.0 |
| 641 |
zulow(jl) =0.0 |
| 642 |
zvlow(jl) =0.0 |
| 643 |
kkcrith(jl) =klev |
| 644 |
kkenvh(jl) =klev |
| 645 |
kentp(jl) =klev |
| 646 |
kcrit(jl) =1 |
| 647 |
ncount(jl) =0 |
| 648 |
ll1(jl,klev+1) =.false. |
| 649 |
2107 continue |
| 650 |
c |
| 651 |
c* define low-level flow |
| 652 |
c --------------------- |
| 653 |
c |
| 654 |
do 223 jk=klev,2,-1 |
| 655 |
do 222 jl=1,klon |
| 656 |
if(ktest(jl).eq.1) then |
| 657 |
zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1) |
| 658 |
prho(jl,jk)=2.*paphm1(jl,jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1)) |
| 659 |
pstab(jl,jk)=2.*zcons2/(ptm1(jl,jk)+ptm1(jl,jk-1))* |
| 660 |
* (1.-rcpd*prho(jl,jk)*(ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk)) |
| 661 |
pstab(jl,jk)=max(pstab(jl,jk),gssec) |
| 662 |
endif |
| 663 |
222 continue |
| 664 |
223 continue |
| 665 |
c |
| 666 |
c******************************************************************** |
| 667 |
c |
| 668 |
c* define blocked flow |
| 669 |
c ------------------- |
| 670 |
do 2115 jk=klev,ilevh,-1 |
| 671 |
do 2116 jl=1,klon |
| 672 |
if(jk.ge.kknub(jl).and.jk.le.kknul(jl)) then |
| 673 |
pulow(jl)=pulow(jl)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk)) |
| 674 |
pvlow(jl)=pvlow(jl)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk)) |
| 675 |
end if |
| 676 |
2116 continue |
| 677 |
2115 continue |
| 678 |
do 2110 jl=1,klon |
| 679 |
pulow(jl)=pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknub(jl))) |
| 680 |
pvlow(jl)=pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknub(jl))) |
| 681 |
znorm(jl)=max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec) |
| 682 |
pvph(jl,klev+1)=znorm(jl) |
| 683 |
2110 continue |
| 684 |
c |
| 685 |
c******* setup orography axes and define plane of profiles ******* |
| 686 |
c |
| 687 |
do 2112 jl=1,klon |
| 688 |
lo=(pulow(jl).lt.gvsec).and.(pulow(jl).ge.-gvsec) |
| 689 |
if(lo) then |
| 690 |
zu=pulow(jl)+2.*gvsec |
| 691 |
else |
| 692 |
zu=pulow(jl) |
| 693 |
endif |
| 694 |
zphi=atan(pvlow(jl)/zu) |
| 695 |
ppsi(jl,klev+1)=ptheta(jl)*rpi/180.-zphi |
| 696 |
zb(jl)=1.-0.18*pgamma(jl)-0.04*pgamma(jl)**2 |
| 697 |
zc(jl)=0.48*pgamma(jl)+0.3*pgamma(jl)**2 |
| 698 |
pd1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ppsi(jl,klev+1))**2) |
| 699 |
pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))*cos(ppsi(jl,klev+1)) |
| 700 |
pdmod(jl)=sqrt(pd1(jl)**2+pd2(jl)**2) |
| 701 |
2112 continue |
| 702 |
c |
| 703 |
c ************ define flow in plane of lowlevel stress ************* |
| 704 |
c |
| 705 |
do 213 jk=1,klev |
| 706 |
do 212 jl=1,klon |
| 707 |
if(ktest(jl).eq.1) then |
| 708 |
zvt1 =pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk) |
| 709 |
zvt2 =-pvlow(jl)*pum1(jl,jk)+pulow(jl)*pvm1(jl,jk) |
| 710 |
zvpf(jl,jk)=(zvt1*pd1(jl)+zvt2*pd2(jl))/(znorm(jl)*pdmod(jl)) |
| 711 |
endif |
| 712 |
ptau(jl,jk) =0.0 |
| 713 |
pzdep(jl,jk) =0.0 |
| 714 |
ppsi(jl,jk) =0.0 |
| 715 |
ll1(jl,jk) =.false. |
| 716 |
212 continue |
| 717 |
213 continue |
| 718 |
do 215 jk=2,klev |
| 719 |
do 214 jl=1,klon |
| 720 |
if(ktest(jl).eq.1) then |
| 721 |
zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1) |
| 722 |
pvph(jl,jk)=((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+ |
| 723 |
* (papm1(jl,jk)-paphm1(jl,jk))*zvpf(jl,jk-1)) |
| 724 |
* /zdp(jl,jk) |
| 725 |
if(pvph(jl,jk).lt.gvsec) then |
| 726 |
pvph(jl,jk)=gvsec |
| 727 |
kcrit(jl)=jk |
| 728 |
endif |
| 729 |
endif |
| 730 |
214 continue |
| 731 |
215 continue |
| 732 |
c |
| 733 |
c |
| 734 |
c* 2.2 brunt-vaisala frequency and density at half levels. |
| 735 |
c |
| 736 |
220 continue |
| 737 |
c |
| 738 |
do 2211 jk=ilevh,klev |
| 739 |
do 221 jl=1,klon |
| 740 |
if(ktest(jl).eq.1) then |
| 741 |
if(jk.ge.(kknub(jl)+1).and.jk.le.kknul(jl)) then |
| 742 |
zst=zcons2/ptm1(jl,jk)*(1.-rcpd*prho(jl,jk)* |
| 743 |
* (ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk)) |
| 744 |
pstab(jl,klev+1)=pstab(jl,klev+1)+zst*zdp(jl,jk) |
| 745 |
pstab(jl,klev+1)=max(pstab(jl,klev+1),gssec) |
| 746 |
prho(jl,klev+1)=prho(jl,klev+1)+paphm1(jl,jk)*2.*zdp(jl,jk) |
| 747 |
* *zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1)) |
| 748 |
endif |
| 749 |
endif |
| 750 |
221 continue |
| 751 |
2211 continue |
| 752 |
c |
| 753 |
do 2212 jl=1,klon |
| 754 |
pstab(jl,klev+1)=pstab(jl,klev+1)/(papm1(jl,kknul(jl)) |
| 755 |
* -papm1(jl,kknub(jl))) |
| 756 |
prho(jl,klev+1)=prho(jl,klev+1)/(papm1(jl,kknul(jl)) |
| 757 |
* -papm1(jl,kknub(jl))) |
| 758 |
zvar=pstd(jl) |
| 759 |
2212 continue |
| 760 |
c |
| 761 |
c* 2.3 mean flow richardson number. |
| 762 |
c* and critical height for froude layer |
| 763 |
c |
| 764 |
230 continue |
| 765 |
c |
| 766 |
do 232 jk=2,klev |
| 767 |
do 231 jl=1,klon |
| 768 |
if(ktest(jl).eq.1) then |
| 769 |
zdwind=max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)),gvsec) |
| 770 |
pri(jl,jk)=pstab(jl,jk)*(zdp(jl,jk) |
| 771 |
* /(rg*prho(jl,jk)*zdwind))**2 |
| 772 |
pri(jl,jk)=max(pri(jl,jk),grcrit) |
| 773 |
endif |
| 774 |
231 continue |
| 775 |
232 continue |
| 776 |
|
| 777 |
c |
| 778 |
c |
| 779 |
c* define top of 'envelope' layer |
| 780 |
c ---------------------------- |
| 781 |
|
| 782 |
do 233 jl=1,klon |
| 783 |
pnu (jl)=0.0 |
| 784 |
znum(jl)=0.0 |
| 785 |
233 continue |
| 786 |
|
| 787 |
do 234 jk=2,klev-1 |
| 788 |
do 234 jl=1,klon |
| 789 |
|
| 790 |
if(ktest(jl).eq.1) then |
| 791 |
|
| 792 |
if (jk.ge.kknub(jl)) then |
| 793 |
|
| 794 |
znum(jl)=pnu(jl) |
| 795 |
zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/ |
| 796 |
* max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec) |
| 797 |
zwind=max(sqrt(zwind**2),gvsec) |
| 798 |
zdelp=paphm1(jl,jk+1)-paphm1(jl,jk) |
| 799 |
zstabm=sqrt(max(pstab(jl,jk ),gssec)) |
| 800 |
zstabp=sqrt(max(pstab(jl,jk+1),gssec)) |
| 801 |
zrhom=prho(jl,jk ) |
| 802 |
zrhop=prho(jl,jk+1) |
| 803 |
pnu(jl) = pnu(jl) + (zdelp/rg)* |
| 804 |
* ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind |
| 805 |
if((znum(jl).le.gfrcrit).and.(pnu(jl).gt.gfrcrit) |
| 806 |
* .and.(kkenvh(jl).eq.klev)) |
| 807 |
* kkenvh(jl)=jk |
| 808 |
|
| 809 |
endif |
| 810 |
|
| 811 |
endif |
| 812 |
|
| 813 |
234 continue |
| 814 |
|
| 815 |
c calculation of a dynamical mixing height for the breaking |
| 816 |
c of gravity waves: |
| 817 |
|
| 818 |
|
| 819 |
do 235 jl=1,klon |
| 820 |
znup(jl)=0.0 |
| 821 |
znum(jl)=0.0 |
| 822 |
235 continue |
| 823 |
|
| 824 |
do 236 jk=klev-1,2,-1 |
| 825 |
do 236 jl=1,klon |
| 826 |
|
| 827 |
if(ktest(jl).eq.1) then |
| 828 |
|
| 829 |
znum(jl)=znup(jl) |
| 830 |
zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/ |
| 831 |
* max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec) |
| 832 |
zwind=max(sqrt(zwind**2),gvsec) |
| 833 |
zdelp=paphm1(jl,jk+1)-paphm1(jl,jk) |
| 834 |
zstabm=sqrt(max(pstab(jl,jk ),gssec)) |
| 835 |
zstabp=sqrt(max(pstab(jl,jk+1),gssec)) |
| 836 |
zrhom=prho(jl,jk ) |
| 837 |
zrhop=prho(jl,jk+1) |
| 838 |
znup(jl) = znup(jl) + (zdelp/rg)* |
| 839 |
* ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind |
| 840 |
if((znum(jl).le.rpi/2.).and.(znup(jl).gt.rpi/2.) |
| 841 |
* .and.(kkcrith(jl).eq.klev)) |
| 842 |
* kkcrith(jl)=jk |
| 843 |
|
| 844 |
endif |
| 845 |
|
| 846 |
236 continue |
| 847 |
|
| 848 |
do 237 jl=1,klon |
| 849 |
kkcrith(jl)=min0(kkcrith(jl),kknu2(jl)) |
| 850 |
kkcrith(jl)=max0(kkcrith(jl),ilevh*2) |
| 851 |
237 continue |
| 852 |
c |
| 853 |
c directional info for flow blocking ************************* |
| 854 |
c |
| 855 |
do 251 jk=ilevh,klev |
| 856 |
do 252 jl=1,klon |
| 857 |
if(jk.ge.kkenvh(jl)) then |
| 858 |
lo=(pum1(jl,jk).lt.gvsec).and.(pum1(jl,jk).ge.-gvsec) |
| 859 |
if(lo) then |
| 860 |
zu=pum1(jl,jk)+2.*gvsec |
| 861 |
else |
| 862 |
zu=pum1(jl,jk) |
| 863 |
endif |
| 864 |
zphi=atan(pvm1(jl,jk)/zu) |
| 865 |
ppsi(jl,jk)=ptheta(jl)*rpi/180.-zphi |
| 866 |
end if |
| 867 |
252 continue |
| 868 |
251 continue |
| 869 |
c forms the vertical 'leakiness' ************************** |
| 870 |
|
| 871 |
alpha=3. |
| 872 |
|
| 873 |
do 254 jk=ilevh,klev |
| 874 |
do 253 jl=1,klon |
| 875 |
if(jk.ge.kkenvh(jl)) then |
| 876 |
zggeenv=amax1(1., |
| 877 |
* (pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)-1))/2.) |
| 878 |
zggeom1=amax1(pgeom1(jl,jk),1.) |
| 879 |
zgvar=amax1(pstd(jl)*rg,1.) |
| 880 |
cmod pzdep(jl,jk)=sqrt((zggeenv-zggeom1)/(zggeom1+zgvar)) |
| 881 |
pzdep(jl,jk)=(pgeom1(jl,kkenvh(jl)-1)-pgeom1(jl, jk))/ |
| 882 |
* (pgeom1(jl,kkenvh(jl)-1)-pgeom1(jl,klev)) |
| 883 |
end if |
| 884 |
253 continue |
| 885 |
254 continue |
| 886 |
|
| 887 |
260 continue |
| 888 |
|
| 889 |
return |
| 890 |
end |
| 891 |
SUBROUTINE gwstress |
| 892 |
* ( nlon , nlev |
| 893 |
* , ktest, kcrit, kkenvh |
| 894 |
* , kknu |
| 895 |
* , prho , pstab , pvph , pstd, psig |
| 896 |
* , pmea , ppic , ptau |
| 897 |
* , pgeom1 , pdmod ) |
| 898 |
c |
| 899 |
c**** *gwstress* |
| 900 |
c |
| 901 |
c purpose. |
| 902 |
c -------- |
| 903 |
c |
| 904 |
c** interface. |
| 905 |
c ---------- |
| 906 |
c call *gwstress* from *gwdrag* |
| 907 |
c |
| 908 |
c explicit arguments : |
| 909 |
c -------------------- |
| 910 |
c ==== inputs === |
| 911 |
c ==== outputs === |
| 912 |
c |
| 913 |
c implicit arguments : none |
| 914 |
c -------------------- |
| 915 |
c |
| 916 |
c method. |
| 917 |
c ------- |
| 918 |
c |
| 919 |
c |
| 920 |
c externals. |
| 921 |
c ---------- |
| 922 |
c |
| 923 |
c |
| 924 |
c reference. |
| 925 |
c ---------- |
| 926 |
c |
| 927 |
c see ecmwf research department documentation of the "i.f.s." |
| 928 |
c |
| 929 |
c author. |
| 930 |
c ------- |
| 931 |
c |
| 932 |
c modifications. |
| 933 |
c -------------- |
| 934 |
c f. lott put the new gwd on ifs 22/11/93 |
| 935 |
c |
| 936 |
c----------------------------------------------------------------------- |
| 937 |
use dimens_m |
| 938 |
use dimphy |
| 939 |
use YOMCST |
| 940 |
use yoegwd |
| 941 |
implicit none |
| 942 |
|
| 943 |
c----------------------------------------------------------------------- |
| 944 |
c |
| 945 |
c* 0.1 arguments |
| 946 |
c --------- |
| 947 |
c |
| 948 |
integer nlon, nlev |
| 949 |
integer kcrit(nlon), |
| 950 |
* ktest(nlon),kkenvh(nlon),kknu(nlon) |
| 951 |
c |
| 952 |
real prho(nlon,nlev+1),pstab(nlon,nlev+1),ptau(nlon,nlev+1), |
| 953 |
* pvph(nlon,nlev+1), |
| 954 |
* pgeom1(nlon,nlev),pstd(nlon) |
| 955 |
c |
| 956 |
real psig(nlon) |
| 957 |
real pmea(nlon),ppic(nlon) |
| 958 |
real pdmod(nlon) |
| 959 |
c |
| 960 |
c----------------------------------------------------------------------- |
| 961 |
c |
| 962 |
c* 0.2 local arrays |
| 963 |
c ------------ |
| 964 |
integer jl |
| 965 |
real zblock, zvar, zeff |
| 966 |
logical lo |
| 967 |
c |
| 968 |
c----------------------------------------------------------------------- |
| 969 |
c |
| 970 |
c* 0.3 functions |
| 971 |
c --------- |
| 972 |
c ------------------------------------------------------------------ |
| 973 |
c |
| 974 |
c* 1. initialization |
| 975 |
c -------------- |
| 976 |
c |
| 977 |
100 continue |
| 978 |
c |
| 979 |
c* 3.1 gravity wave stress. |
| 980 |
c |
| 981 |
300 continue |
| 982 |
c |
| 983 |
c |
| 984 |
do 301 jl=1,klon |
| 985 |
if(ktest(jl).eq.1) then |
| 986 |
|
| 987 |
c effective mountain height above the blocked flow |
| 988 |
|
| 989 |
if(kkenvh(jl).eq.klev)then |
| 990 |
zblock=0.0 |
| 991 |
else |
| 992 |
zblock=(pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)+1))/2./rg |
| 993 |
endif |
| 994 |
|
| 995 |
zvar=ppic(jl)-pmea(jl) |
| 996 |
zeff=amax1(0.,zvar-zblock) |
| 997 |
|
| 998 |
ptau(jl,klev+1)=prho(jl,klev+1)*gkdrag*psig(jl)*zeff**2 |
| 999 |
* /4./pstd(jl)*pvph(jl,klev+1)*pdmod(jl)*sqrt(pstab(jl,klev+1)) |
| 1000 |
|
| 1001 |
c too small value of stress or low level flow include critical level |
| 1002 |
c or low level flow: gravity wave stress nul. |
| 1003 |
|
| 1004 |
lo=(ptau(jl,klev+1).lt.gtsec).or.(kcrit(jl).ge.kknu(jl)) |
| 1005 |
* .or.(pvph(jl,klev+1).lt.gvcrit) |
| 1006 |
c if(lo) ptau(jl,klev+1)=0.0 |
| 1007 |
|
| 1008 |
else |
| 1009 |
|
| 1010 |
ptau(jl,klev+1)=0.0 |
| 1011 |
|
| 1012 |
endif |
| 1013 |
|
| 1014 |
301 continue |
| 1015 |
c |
| 1016 |
return |
| 1017 |
end |
| 1018 |
SUBROUTINE GWPROFIL |
| 1019 |
* ( NLON, NLEV |
| 1020 |
* , kgwd, kdx , ktest |
| 1021 |
* , KKCRITH, KCRIT |
| 1022 |
* , PAPHM1, PRHO , PSTAB , PVPH , PRI , PTAU |
| 1023 |
* , pdmod , psig , pvar) |
| 1024 |
|
| 1025 |
C**** *GWPROFIL* |
| 1026 |
C |
| 1027 |
C PURPOSE. |
| 1028 |
C -------- |
| 1029 |
C |
| 1030 |
C** INTERFACE. |
| 1031 |
C ---------- |
| 1032 |
C FROM *GWDRAG* |
| 1033 |
C |
| 1034 |
C EXPLICIT ARGUMENTS : |
| 1035 |
C -------------------- |
| 1036 |
C ==== INPUTS === |
| 1037 |
C ==== OUTPUTS === |
| 1038 |
C |
| 1039 |
C IMPLICIT ARGUMENTS : NONE |
| 1040 |
C -------------------- |
| 1041 |
C |
| 1042 |
C METHOD: |
| 1043 |
C ------- |
| 1044 |
C THE STRESS PROFILE FOR GRAVITY WAVES IS COMPUTED AS FOLLOWS: |
| 1045 |
C IT IS CONSTANT (NO GWD) AT THE LEVELS BETWEEN THE GROUND |
| 1046 |
C AND THE TOP OF THE BLOCKED LAYER (KKENVH). |
| 1047 |
C IT DECREASES LINEARLY WITH HEIGHTS FROM THE TOP OF THE |
| 1048 |
C BLOCKED LAYER TO 3*VAROR (kKNU), TO SIMULATES LEE WAVES OR |
| 1049 |
C NONLINEAR GRAVITY WAVE BREAKING. |
| 1050 |
C ABOVE IT IS CONSTANT, EXCEPT WHEN THE WAVE ENCOUNTERS A CRITICAL |
| 1051 |
C LEVEL (KCRIT) OR WHEN IT BREAKS. |
| 1052 |
C |
| 1053 |
C |
| 1054 |
C |
| 1055 |
C EXTERNALS. |
| 1056 |
C ---------- |
| 1057 |
C |
| 1058 |
C |
| 1059 |
C REFERENCE. |
| 1060 |
C ---------- |
| 1061 |
C |
| 1062 |
C SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S." |
| 1063 |
C |
| 1064 |
C AUTHOR. |
| 1065 |
C ------- |
| 1066 |
C |
| 1067 |
C MODIFICATIONS. |
| 1068 |
C -------------- |
| 1069 |
C PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93) |
| 1070 |
C----------------------------------------------------------------------- |
| 1071 |
use dimens_m |
| 1072 |
use dimphy |
| 1073 |
use YOMCST |
| 1074 |
use yoegwd |
| 1075 |
implicit none |
| 1076 |
C |
| 1077 |
|
| 1078 |
C |
| 1079 |
|
| 1080 |
|
| 1081 |
C----------------------------------------------------------------------- |
| 1082 |
C |
| 1083 |
C* 0.1 ARGUMENTS |
| 1084 |
C --------- |
| 1085 |
C |
| 1086 |
integer nlon,nlev |
| 1087 |
INTEGER KKCRITH(NLON),KCRIT(NLON) |
| 1088 |
* ,kdx(nlon) , ktest(nlon) |
| 1089 |
|
| 1090 |
C |
| 1091 |
REAL PAPHM1(NLON,NLEV+1), PSTAB(NLON,NLEV+1), |
| 1092 |
* PRHO (NLON,NLEV+1), PVPH (NLON,NLEV+1), |
| 1093 |
* PRI (NLON,NLEV+1), PTAU(NLON,NLEV+1) |
| 1094 |
|
| 1095 |
REAL pdmod (NLON) , psig(NLON), |
| 1096 |
* pvar(NLON) |
| 1097 |
|
| 1098 |
C----------------------------------------------------------------------- |
| 1099 |
C |
| 1100 |
C* 0.2 LOCAL ARRAYS |
| 1101 |
C ------------ |
| 1102 |
C |
| 1103 |
integer ilevh, ji, kgwd, jl, jk |
| 1104 |
real zsqr, zalfa, zriw, zdel, zb, zalpha,zdz2n |
| 1105 |
real zdelp, zdelpt |
| 1106 |
REAL ZDZ2 (KLON,KLEV) , ZNORM(KLON) , zoro(KLON) |
| 1107 |
REAL ZTAU (KLON,KLEV+1) |
| 1108 |
C |
| 1109 |
C----------------------------------------------------------------------- |
| 1110 |
C |
| 1111 |
C* 1. INITIALIZATION |
| 1112 |
C -------------- |
| 1113 |
C |
| 1114 |
c print *,' entree gwprofil' |
| 1115 |
100 CONTINUE |
| 1116 |
C |
| 1117 |
C |
| 1118 |
C* COMPUTATIONAL CONSTANTS. |
| 1119 |
C ------------- ---------- |
| 1120 |
C |
| 1121 |
ilevh=KLEV/3 |
| 1122 |
C |
| 1123 |
c DO 400 ji=1,kgwd |
| 1124 |
c jl=kdx(ji) |
| 1125 |
c Modif vectorisation 02/04/2004 |
| 1126 |
DO 400 jl=1,klon |
| 1127 |
if (ktest(jl).eq.1) then |
| 1128 |
Zoro(JL)=Psig(JL)*Pdmod(JL)/4./max(pvar(jl),1.0) |
| 1129 |
ZTAU(JL,KLEV+1)=PTAU(JL,KLEV+1) |
| 1130 |
endif |
| 1131 |
400 CONTINUE |
| 1132 |
|
| 1133 |
C |
| 1134 |
DO 430 JK=KLEV,2,-1 |
| 1135 |
C |
| 1136 |
C |
| 1137 |
C* 4.1 CONSTANT WAVE STRESS UNTIL TOP OF THE |
| 1138 |
C BLOCKING LAYER. |
| 1139 |
410 CONTINUE |
| 1140 |
C |
| 1141 |
c DO 411 ji=1,kgwd |
| 1142 |
c jl=kdx(ji) |
| 1143 |
c Modif vectorisation 02/04/2004 |
| 1144 |
do 411 jl=1,klon |
| 1145 |
if (ktest(jl).eq.1) then |
| 1146 |
IF(JK.GT.KKCRITH(JL)) THEN |
| 1147 |
PTAU(JL,JK)=ZTAU(JL,KLEV+1) |
| 1148 |
C ENDIF |
| 1149 |
C IF(JK.EQ.KKCRITH(JL)) THEN |
| 1150 |
ELSE |
| 1151 |
PTAU(JL,JK)=GRAHILO*ZTAU(JL,KLEV+1) |
| 1152 |
ENDIF |
| 1153 |
endif |
| 1154 |
411 CONTINUE |
| 1155 |
C |
| 1156 |
C* 4.15 CONSTANT SHEAR STRESS UNTIL THE TOP OF THE |
| 1157 |
C LOW LEVEL FLOW LAYER. |
| 1158 |
415 CONTINUE |
| 1159 |
C |
| 1160 |
C |
| 1161 |
C* 4.2 WAVE DISPLACEMENT AT NEXT LEVEL. |
| 1162 |
C |
| 1163 |
420 CONTINUE |
| 1164 |
C |
| 1165 |
c DO 421 ji=1,kgwd |
| 1166 |
c jl=kdx(ji) |
| 1167 |
c Modif vectorisation 02/04/2004 |
| 1168 |
do 421 jl=1,klon |
| 1169 |
if(ktest(jl).eq.1) then |
| 1170 |
IF(JK.LT.KKCRITH(JL)) THEN |
| 1171 |
ZNORM(JL)=gkdrag*PRHO(JL,JK)*SQRT(PSTAB(JL,JK))*PVPH(JL,JK) |
| 1172 |
* *zoro(jl) |
| 1173 |
ZDZ2(JL,JK)=PTAU(JL,JK+1)/max(ZNORM(JL),gssec) |
| 1174 |
ENDIF |
| 1175 |
endif |
| 1176 |
421 CONTINUE |
| 1177 |
C |
| 1178 |
C* 4.3 WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT |
| 1179 |
C* AND STRESS: BREAKING EVALUATION AND CRITICAL |
| 1180 |
C LEVEL |
| 1181 |
C |
| 1182 |
|
| 1183 |
c DO 431 ji=1,kgwd |
| 1184 |
c jl=Kdx(ji) |
| 1185 |
c Modif vectorisation 02/04/2004 |
| 1186 |
do 431 jl=1,klon |
| 1187 |
if(ktest(jl).eq.1) then |
| 1188 |
|
| 1189 |
IF(JK.LT.KKCRITH(JL)) THEN |
| 1190 |
IF((PTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.KCRIT(JL))) THEN |
| 1191 |
PTAU(JL,JK)=0.0 |
| 1192 |
ELSE |
| 1193 |
ZSQR=SQRT(PRI(JL,JK)) |
| 1194 |
ZALFA=SQRT(PSTAB(JL,JK)*ZDZ2(JL,JK))/PVPH(JL,JK) |
| 1195 |
ZRIW=PRI(JL,JK)*(1.-ZALFA)/(1+ZALFA*ZSQR)**2 |
| 1196 |
IF(ZRIW.LT.GRCRIT) THEN |
| 1197 |
ZDEL=4./ZSQR/GRCRIT+1./GRCRIT**2+4./GRCRIT |
| 1198 |
ZB=1./GRCRIT+2./ZSQR |
| 1199 |
ZALPHA=0.5*(-ZB+SQRT(ZDEL)) |
| 1200 |
ZDZ2N=(PVPH(JL,JK)*ZALPHA)**2/PSTAB(JL,JK) |
| 1201 |
PTAU(JL,JK)=ZNORM(JL)*ZDZ2N |
| 1202 |
ELSE |
| 1203 |
PTAU(JL,JK)=ZNORM(JL)*ZDZ2(JL,JK) |
| 1204 |
ENDIF |
| 1205 |
PTAU(JL,JK)=MIN(PTAU(JL,JK),PTAU(JL,JK+1)) |
| 1206 |
ENDIF |
| 1207 |
ENDIF |
| 1208 |
endif |
| 1209 |
431 CONTINUE |
| 1210 |
|
| 1211 |
430 CONTINUE |
| 1212 |
440 CONTINUE |
| 1213 |
|
| 1214 |
C REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL |
| 1215 |
|
| 1216 |
c DO 530 ji=1,kgwd |
| 1217 |
c jl=kdx(ji) |
| 1218 |
c Modif vectorisation 02/04/2004 |
| 1219 |
do 530 jl=1,klon |
| 1220 |
if(ktest(jl).eq.1) then |
| 1221 |
ZTAU(JL,KKCRITH(JL))=PTAU(JL,KKCRITH(JL)) |
| 1222 |
ZTAU(JL,NSTRA)=PTAU(JL,NSTRA) |
| 1223 |
endif |
| 1224 |
530 CONTINUE |
| 1225 |
|
| 1226 |
DO 531 JK=1,KLEV |
| 1227 |
|
| 1228 |
c DO 532 ji=1,kgwd |
| 1229 |
c jl=kdx(ji) |
| 1230 |
c Modif vectorisation 02/04/2004 |
| 1231 |
do 532 jl=1,klon |
| 1232 |
if(ktest(jl).eq.1) then |
| 1233 |
|
| 1234 |
|
| 1235 |
IF(JK.GT.KKCRITH(JL))THEN |
| 1236 |
|
| 1237 |
ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KLEV+1 ) |
| 1238 |
ZDELPT=PAPHM1(JL,KKCRITH(JL))-PAPHM1(JL,KLEV+1 ) |
| 1239 |
PTAU(JL,JK)=ZTAU(JL,KLEV+1 ) + |
| 1240 |
. (ZTAU(JL,KKCRITH(JL))-ZTAU(JL,KLEV+1 ) )* |
| 1241 |
. ZDELP/ZDELPT |
| 1242 |
|
| 1243 |
ENDIF |
| 1244 |
|
| 1245 |
endif |
| 1246 |
532 CONTINUE |
| 1247 |
|
| 1248 |
C REORGANISATION IN THE STRATOSPHERE |
| 1249 |
|
| 1250 |
c DO 533 ji=1,kgwd |
| 1251 |
c jl=kdx(ji) |
| 1252 |
c Modif vectorisation 02/04/2004 |
| 1253 |
do 533 jl=1,klon |
| 1254 |
if(ktest(jl).eq.1) then |
| 1255 |
|
| 1256 |
|
| 1257 |
IF(JK.LT.NSTRA)THEN |
| 1258 |
|
| 1259 |
ZDELP =PAPHM1(JL,NSTRA) |
| 1260 |
ZDELPT=PAPHM1(JL,JK) |
| 1261 |
PTAU(JL,JK)=ZTAU(JL,NSTRA)*ZDELPT/ZDELP |
| 1262 |
|
| 1263 |
ENDIF |
| 1264 |
|
| 1265 |
endif |
| 1266 |
533 CONTINUE |
| 1267 |
|
| 1268 |
C REORGANISATION IN THE TROPOSPHERE |
| 1269 |
|
| 1270 |
c DO 534 ji=1,kgwd |
| 1271 |
c jl=kdx(ji) |
| 1272 |
c Modif vectorisation 02/04/2004 |
| 1273 |
do 534 jl=1,klon |
| 1274 |
if(ktest(jl).eq.1) then |
| 1275 |
|
| 1276 |
|
| 1277 |
IF(JK.LT.KKCRITH(JL).AND.JK.GT.NSTRA)THEN |
| 1278 |
|
| 1279 |
ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KKCRITH(JL)) |
| 1280 |
ZDELPT=PAPHM1(JL,NSTRA)-PAPHM1(JL,KKCRITH(JL)) |
| 1281 |
PTAU(JL,JK)=ZTAU(JL,KKCRITH(JL)) + |
| 1282 |
* (ZTAU(JL,NSTRA)-ZTAU(JL,KKCRITH(JL)))*ZDELP |
| 1283 |
. /ZDELPT |
| 1284 |
|
| 1285 |
ENDIF |
| 1286 |
endif |
| 1287 |
534 CONTINUE |
| 1288 |
|
| 1289 |
|
| 1290 |
531 CONTINUE |
| 1291 |
|
| 1292 |
|
| 1293 |
RETURN |
| 1294 |
END |
| 1295 |
SUBROUTINE lift_noro (nlon,nlev,dtime,paprs,pplay, |
| 1296 |
e plat,pmea,pstd, ppic, |
| 1297 |
e ktest, |
| 1298 |
e t, u, v, |
| 1299 |
s pulow, pvlow, pustr, pvstr, |
| 1300 |
s d_t, d_u, d_v) |
| 1301 |
c |
| 1302 |
use dimens_m |
| 1303 |
use dimphy |
| 1304 |
use YOMCST |
| 1305 |
IMPLICIT none |
| 1306 |
c====================================================================== |
| 1307 |
c Auteur(s): F.Lott (LMD/CNRS) date: 19950201 |
| 1308 |
c Objet: Frottement de la montagne Interface |
| 1309 |
c====================================================================== |
| 1310 |
c Arguments: |
| 1311 |
c dtime---input-R- pas d'integration (s) |
| 1312 |
c paprs---input-R-pression pour chaque inter-couche (en Pa) |
| 1313 |
c pplay---input-R-pression pour le mileu de chaque couche (en Pa) |
| 1314 |
c t-------input-R-temperature (K) |
| 1315 |
c u-------input-R-vitesse horizontale (m/s) |
| 1316 |
c v-------input-R-vitesse horizontale (m/s) |
| 1317 |
c |
| 1318 |
c d_t-----output-R-increment de la temperature |
| 1319 |
c d_u-----output-R-increment de la vitesse u |
| 1320 |
c d_v-----output-R-increment de la vitesse v |
| 1321 |
c====================================================================== |
| 1322 |
c |
| 1323 |
c ARGUMENTS |
| 1324 |
c |
| 1325 |
INTEGER nlon,nlev |
| 1326 |
REAL, intent(in):: dtime |
| 1327 |
REAL, intent(in):: paprs(klon,klev+1) |
| 1328 |
REAL, intent(in):: pplay(klon,klev) |
| 1329 |
REAL, intent(in):: plat(nlon) |
| 1330 |
real pmea(nlon) |
| 1331 |
REAL pstd(nlon) |
| 1332 |
REAL ppic(nlon) |
| 1333 |
REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon) |
| 1334 |
REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev) |
| 1335 |
REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev) |
| 1336 |
c |
| 1337 |
INTEGER i, k, ktest(nlon) |
| 1338 |
c |
| 1339 |
c Variables locales: |
| 1340 |
c |
| 1341 |
REAL zgeom(klon,klev) |
| 1342 |
REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev) |
| 1343 |
REAL pt(klon,klev), pu(klon,klev), pv(klon,klev) |
| 1344 |
REAL papmf(klon,klev),papmh(klon,klev+1) |
| 1345 |
c |
| 1346 |
c initialiser les variables de sortie (pour securite) |
| 1347 |
c |
| 1348 |
DO i = 1,klon |
| 1349 |
pulow(i) = 0.0 |
| 1350 |
pvlow(i) = 0.0 |
| 1351 |
pustr(i) = 0.0 |
| 1352 |
pvstr(i) = 0.0 |
| 1353 |
ENDDO |
| 1354 |
DO k = 1, klev |
| 1355 |
DO i = 1, klon |
| 1356 |
d_t(i,k) = 0.0 |
| 1357 |
d_u(i,k) = 0.0 |
| 1358 |
d_v(i,k) = 0.0 |
| 1359 |
pdudt(i,k)=0.0 |
| 1360 |
pdvdt(i,k)=0.0 |
| 1361 |
pdtdt(i,k)=0.0 |
| 1362 |
ENDDO |
| 1363 |
ENDDO |
| 1364 |
c |
| 1365 |
c preparer les variables d'entree (attention: l'ordre des niveaux |
| 1366 |
c verticaux augmente du haut vers le bas) |
| 1367 |
c |
| 1368 |
DO k = 1, klev |
| 1369 |
DO i = 1, klon |
| 1370 |
pt(i,k) = t(i,klev-k+1) |
| 1371 |
pu(i,k) = u(i,klev-k+1) |
| 1372 |
pv(i,k) = v(i,klev-k+1) |
| 1373 |
papmf(i,k) = pplay(i,klev-k+1) |
| 1374 |
ENDDO |
| 1375 |
ENDDO |
| 1376 |
DO k = 1, klev+1 |
| 1377 |
DO i = 1, klon |
| 1378 |
papmh(i,k) = paprs(i,klev-k+2) |
| 1379 |
ENDDO |
| 1380 |
ENDDO |
| 1381 |
DO i = 1, klon |
| 1382 |
zgeom(i,klev) = RD * pt(i,klev) |
| 1383 |
. * LOG(papmh(i,klev+1)/papmf(i,klev)) |
| 1384 |
ENDDO |
| 1385 |
DO k = klev-1, 1, -1 |
| 1386 |
DO i = 1, klon |
| 1387 |
zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0 |
| 1388 |
. * LOG(papmf(i,k+1)/papmf(i,k)) |
| 1389 |
ENDDO |
| 1390 |
ENDDO |
| 1391 |
c |
| 1392 |
c appeler la routine principale |
| 1393 |
c |
| 1394 |
CALL OROLIFT(klon,klev,ktest, |
| 1395 |
. dtime, |
| 1396 |
. papmh, zgeom, |
| 1397 |
. pt, pu, pv, |
| 1398 |
. plat,pmea, pstd, ppic, |
| 1399 |
. pulow,pvlow, |
| 1400 |
. pdudt,pdvdt,pdtdt) |
| 1401 |
C |
| 1402 |
DO k = 1, klev |
| 1403 |
DO i = 1, klon |
| 1404 |
d_u(i,klev+1-k) = dtime*pdudt(i,k) |
| 1405 |
d_v(i,klev+1-k) = dtime*pdvdt(i,k) |
| 1406 |
d_t(i,klev+1-k) = dtime*pdtdt(i,k) |
| 1407 |
pustr(i) = pustr(i) |
| 1408 |
cIM BUG . +RG*pdudt(i,k)*(papmh(i,k+1)-papmh(i,k)) |
| 1409 |
. +pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG |
| 1410 |
pvstr(i) = pvstr(i) |
| 1411 |
cIM BUG . +RG*pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k)) |
| 1412 |
. +pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG |
| 1413 |
ENDDO |
| 1414 |
ENDDO |
| 1415 |
c |
| 1416 |
RETURN |
| 1417 |
END |
| 1418 |
SUBROUTINE OROLIFT( NLON,NLEV |
| 1419 |
I , KTEST |
| 1420 |
R , PTSPHY |
| 1421 |
R , PAPHM1,PGEOM1,PTM1,PUM1,PVM1 |
| 1422 |
R , PLAT |
| 1423 |
R , PMEA, PVAROR, ppic |
| 1424 |
C OUTPUTS |
| 1425 |
R , PULOW,PVLOW |
| 1426 |
R , PVOM,PVOL,PTE ) |
| 1427 |
|
| 1428 |
C |
| 1429 |
C**** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT. |
| 1430 |
C |
| 1431 |
C PURPOSE. |
| 1432 |
C -------- |
| 1433 |
C |
| 1434 |
C** INTERFACE. |
| 1435 |
C ---------- |
| 1436 |
C CALLED FROM *lift_noro |
| 1437 |
C ---------- |
| 1438 |
C |
| 1439 |
C AUTHOR. |
| 1440 |
C ------- |
| 1441 |
C F.LOTT LMD 22/11/95 |
| 1442 |
C |
| 1443 |
use dimens_m |
| 1444 |
use dimphy |
| 1445 |
use YOMCST |
| 1446 |
use yoegwd |
| 1447 |
implicit none |
| 1448 |
C |
| 1449 |
C |
| 1450 |
C----------------------------------------------------------------------- |
| 1451 |
C |
| 1452 |
C* 0.1 ARGUMENTS |
| 1453 |
C --------- |
| 1454 |
C |
| 1455 |
C |
| 1456 |
integer nlon, nlev |
| 1457 |
REAL PTE(NLON,NLEV), |
| 1458 |
* PVOL(NLON,NLEV), |
| 1459 |
* PVOM(NLON,NLEV), |
| 1460 |
* PULOW(NLON), |
| 1461 |
* PVLOW(NLON) |
| 1462 |
REAL PUM1(NLON,NLEV), |
| 1463 |
* PVM1(NLON,NLEV), |
| 1464 |
* PTM1(NLON,NLEV) |
| 1465 |
real, intent(in):: PLAT(NLON) |
| 1466 |
real PMEA(NLON), |
| 1467 |
* PVAROR(NLON), |
| 1468 |
* ppic(NLON), |
| 1469 |
* PGEOM1(NLON,NLEV), |
| 1470 |
* PAPHM1(NLON,NLEV+1) |
| 1471 |
C |
| 1472 |
INTEGER KTEST(NLON) |
| 1473 |
real, intent(in):: ptsphy |
| 1474 |
C----------------------------------------------------------------------- |
| 1475 |
C |
| 1476 |
C* 0.2 LOCAL ARRAYS |
| 1477 |
C ------------ |
| 1478 |
logical lifthigh |
| 1479 |
integer klevm1, jl, ilevh, jk |
| 1480 |
real zcons1, ztmst, zrtmst,zpi, zhgeo |
| 1481 |
real zdelp, zslow, zsqua, zscav, zbet |
| 1482 |
INTEGER |
| 1483 |
* IKNUB(klon), |
| 1484 |
* IKNUL(klon) |
| 1485 |
LOGICAL LL1(KLON,KLEV+1) |
| 1486 |
C |
| 1487 |
REAL ZTAU(KLON,KLEV+1), |
| 1488 |
* ZTAV(KLON,KLEV+1), |
| 1489 |
* ZRHO(KLON,KLEV+1) |
| 1490 |
REAL ZDUDT(KLON), |
| 1491 |
* ZDVDT(KLON) |
| 1492 |
REAL ZHCRIT(KLON,KLEV) |
| 1493 |
C----------------------------------------------------------------------- |
| 1494 |
C |
| 1495 |
C* 1.1 INITIALIZATIONS |
| 1496 |
C --------------- |
| 1497 |
|
| 1498 |
LIFTHIGH=.FALSE. |
| 1499 |
|
| 1500 |
IF(NLON.NE.KLON.OR.NLEV.NE.KLEV)STOP |
| 1501 |
ZCONS1=1./RD |
| 1502 |
KLEVM1=KLEV-1 |
| 1503 |
ZTMST=PTSPHY |
| 1504 |
ZRTMST=1./ZTMST |
| 1505 |
ZPI=ACOS(-1.) |
| 1506 |
C |
| 1507 |
DO 1001 JL=1,klon |
| 1508 |
ZRHO(JL,KLEV+1) =0.0 |
| 1509 |
PULOW(JL) =0.0 |
| 1510 |
PVLOW(JL) =0.0 |
| 1511 |
iknub(JL) =klev |
| 1512 |
iknul(JL) =klev |
| 1513 |
ilevh=klev/3 |
| 1514 |
ll1(jl,klev+1)=.false. |
| 1515 |
DO 1000 JK=1,KLEV |
| 1516 |
PVOM(JL,JK)=0.0 |
| 1517 |
PVOL(JL,JK)=0.0 |
| 1518 |
PTE (JL,JK)=0.0 |
| 1519 |
1000 CONTINUE |
| 1520 |
1001 CONTINUE |
| 1521 |
|
| 1522 |
C |
| 1523 |
C* 2.1 DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF |
| 1524 |
C* LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE |
| 1525 |
C* THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS. |
| 1526 |
C |
| 1527 |
C |
| 1528 |
C |
| 1529 |
DO 2006 JK=KLEV,1,-1 |
| 1530 |
DO 2007 JL=1,klon |
| 1531 |
IF(KTEST(JL).EQ.1) THEN |
| 1532 |
ZHCRIT(JL,JK)=amax1(Ppic(JL)-pmea(JL),100.) |
| 1533 |
ZHGEO=PGEOM1(JL,JK)/RG |
| 1534 |
ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK)) |
| 1535 |
IF(ll1(JL,JK).neqv.ll1(JL,JK+1)) THEN |
| 1536 |
iknub(JL)=JK |
| 1537 |
ENDIF |
| 1538 |
ENDIF |
| 1539 |
2007 CONTINUE |
| 1540 |
2006 CONTINUE |
| 1541 |
C |
| 1542 |
do 2010 jl=1,klon |
| 1543 |
IF(KTEST(JL).EQ.1) THEN |
| 1544 |
iknub(jl)=max(iknub(jl),klev/2) |
| 1545 |
iknul(jl)=max(iknul(jl),2*klev/3) |
| 1546 |
if(iknub(jl).gt.nktopg) iknub(jl)=nktopg |
| 1547 |
if(iknub(jl).eq.nktopg) iknul(jl)=klev |
| 1548 |
if(iknub(jl).eq.iknul(jl)) iknub(jl)=iknul(jl)-1 |
| 1549 |
ENDIF |
| 1550 |
2010 continue |
| 1551 |
|
| 1552 |
C do 2011 jl=1,klon |
| 1553 |
C IF(KTEST(JL).EQ.1) THEN |
| 1554 |
C print *,' iknul= ',iknul(jl),' iknub=',iknub(jl) |
| 1555 |
C ENDIF |
| 1556 |
C2011 continue |
| 1557 |
|
| 1558 |
C PRINT *,' DANS OROLIFT: 2010' |
| 1559 |
|
| 1560 |
DO 223 JK=KLEV,2,-1 |
| 1561 |
DO 222 JL=1,klon |
| 1562 |
ZRHO(JL,JK)=2.*PAPHM1(JL,JK)*ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1)) |
| 1563 |
222 CONTINUE |
| 1564 |
223 CONTINUE |
| 1565 |
C PRINT *,' DANS OROLIFT: 223' |
| 1566 |
|
| 1567 |
C******************************************************************** |
| 1568 |
C |
| 1569 |
C* DEFINE LOW LEVEL FLOW |
| 1570 |
C ------------------- |
| 1571 |
DO 2115 JK=klev,1,-1 |
| 1572 |
DO 2116 JL=1,klon |
| 1573 |
IF(KTEST(JL).EQ.1) THEN |
| 1574 |
if(jk.ge.iknub(jl).and.jk.le.iknul(jl)) then |
| 1575 |
pulow(JL)=pulow(JL)+PUM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK)) |
| 1576 |
pvlow(JL)=pvlow(JL)+PVM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK)) |
| 1577 |
zrho(JL,klev+1)=zrho(JL,klev+1) |
| 1578 |
* +zrho(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK)) |
| 1579 |
end if |
| 1580 |
ENDIF |
| 1581 |
2116 CONTINUE |
| 1582 |
2115 CONTINUE |
| 1583 |
DO 2110 JL=1,klon |
| 1584 |
IF(KTEST(JL).EQ.1) THEN |
| 1585 |
pulow(JL)=pulow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl))) |
| 1586 |
pvlow(JL)=pvlow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl))) |
| 1587 |
zrho(JL,klev+1)=zrho(JL,klev+1) |
| 1588 |
* /(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl))) |
| 1589 |
ENDIF |
| 1590 |
2110 CONTINUE |
| 1591 |
|
| 1592 |
|
| 1593 |
200 CONTINUE |
| 1594 |
|
| 1595 |
C*********************************************************** |
| 1596 |
C |
| 1597 |
C* 3. COMPUTE MOUNTAIN LIFT |
| 1598 |
C |
| 1599 |
300 CONTINUE |
| 1600 |
C |
| 1601 |
DO 301 JL=1,klon |
| 1602 |
IF(KTEST(JL).EQ.1) THEN |
| 1603 |
ZTAU(JL,KLEV+1)= - GKLIFT*ZRHO(JL,KLEV+1)*2.*ROMEGA* |
| 1604 |
C * (2*PVAROR(JL)+PMEA(JL))* |
| 1605 |
* 2*PVAROR(JL)* |
| 1606 |
* SIN(ZPI/180.*PLAT(JL))*PVLOW(JL) |
| 1607 |
ZTAV(JL,KLEV+1)= GKLIFT*ZRHO(JL,KLEV+1)*2.*ROMEGA* |
| 1608 |
C * (2*PVAROR(JL)+PMEA(JL))* |
| 1609 |
* 2*PVAROR(JL)* |
| 1610 |
* SIN(ZPI/180.*PLAT(JL))*PULOW(JL) |
| 1611 |
ELSE |
| 1612 |
ZTAU(JL,KLEV+1)=0.0 |
| 1613 |
ZTAV(JL,KLEV+1)=0.0 |
| 1614 |
ENDIF |
| 1615 |
301 CONTINUE |
| 1616 |
|
| 1617 |
C |
| 1618 |
C* 4. COMPUTE LIFT PROFILE |
| 1619 |
C* -------------------- |
| 1620 |
C |
| 1621 |
|
| 1622 |
400 CONTINUE |
| 1623 |
|
| 1624 |
DO 401 JK=1,KLEV |
| 1625 |
DO 401 JL=1,klon |
| 1626 |
IF(KTEST(JL).EQ.1) THEN |
| 1627 |
ZTAU(JL,JK)=ZTAU(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1) |
| 1628 |
ZTAV(JL,JK)=ZTAV(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1) |
| 1629 |
ELSE |
| 1630 |
ZTAU(JL,JK)=0.0 |
| 1631 |
ZTAV(JL,JK)=0.0 |
| 1632 |
ENDIF |
| 1633 |
401 CONTINUE |
| 1634 |
C |
| 1635 |
C |
| 1636 |
C* 5. COMPUTE TENDENCIES. |
| 1637 |
C* ------------------- |
| 1638 |
IF(LIFTHIGH)THEN |
| 1639 |
C |
| 1640 |
500 CONTINUE |
| 1641 |
C PRINT *,' DANS OROLIFT: 500' |
| 1642 |
C |
| 1643 |
C EXPLICIT SOLUTION AT ALL LEVELS |
| 1644 |
C |
| 1645 |
DO 524 JK=1,klev |
| 1646 |
DO 523 JL=1,KLON |
| 1647 |
IF(KTEST(JL).EQ.1) THEN |
| 1648 |
ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK) |
| 1649 |
ZDUDT(JL)=-RG*(ZTAU(JL,JK+1)-ZTAU(JL,JK))/ZDELP |
| 1650 |
ZDVDT(JL)=-RG*(ZTAV(JL,JK+1)-ZTAV(JL,JK))/ZDELP |
| 1651 |
ENDIF |
| 1652 |
523 CONTINUE |
| 1653 |
524 CONTINUE |
| 1654 |
C |
| 1655 |
C PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY |
| 1656 |
C |
| 1657 |
DO 530 JK=1,klev |
| 1658 |
DO 530 JL=1,KLON |
| 1659 |
IF(KTEST(JL).EQ.1) THEN |
| 1660 |
|
| 1661 |
ZSLOW=SQRT(PULOW(JL)**2+PVLOW(JL)**2) |
| 1662 |
ZSQUA=AMAX1(SQRT(PUM1(JL,JK)**2+PVM1(JL,JK)**2),GVSEC) |
| 1663 |
ZSCAV=-ZDUDT(JL)*PVM1(JL,JK)+ZDVDT(JL)*PUM1(JL,JK) |
| 1664 |
IF(ZSQUA.GT.GVSEC)THEN |
| 1665 |
PVOM(JL,JK)=-ZSCAV*PVM1(JL,JK)/ZSQUA**2 |
| 1666 |
PVOL(JL,JK)= ZSCAV*PUM1(JL,JK)/ZSQUA**2 |
| 1667 |
ELSE |
| 1668 |
PVOM(JL,JK)=0.0 |
| 1669 |
PVOL(JL,JK)=0.0 |
| 1670 |
ENDIF |
| 1671 |
ZSQUA=SQRT(PUM1(JL,JK)**2+PUM1(JL,JK)**2) |
| 1672 |
IF(ZSQUA.LT.ZSLOW)THEN |
| 1673 |
PVOM(JL,JK)=ZSQUA/ZSLOW*PVOM(JL,JK) |
| 1674 |
PVOL(JL,JK)=ZSQUA/ZSLOW*PVOL(JL,JK) |
| 1675 |
ENDIF |
| 1676 |
|
| 1677 |
ENDIF |
| 1678 |
530 CONTINUE |
| 1679 |
C |
| 1680 |
C 6. LOW LEVEL LIFT, SEMI IMPLICIT: |
| 1681 |
C ---------------------------------- |
| 1682 |
|
| 1683 |
ELSE |
| 1684 |
|
| 1685 |
DO 601 JL=1,KLON |
| 1686 |
IF(KTEST(JL).EQ.1) THEN |
| 1687 |
DO JK=KLEV,IKNUB(JL),-1 |
| 1688 |
ZBET=GKLIFT*2.*ROMEGA*SIN(ZPI/180.*PLAT(JL))*ztmst* |
| 1689 |
* (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL, JK))/ |
| 1690 |
* (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL,KLEV)) |
| 1691 |
ZDUDT(JL)=-PUM1(JL,JK)/ztmst/(1+ZBET**2) |
| 1692 |
ZDVDT(JL)=-PVM1(JL,JK)/ztmst/(1+ZBET**2) |
| 1693 |
PVOM(JL,JK)= ZBET**2*ZDUDT(JL) - ZBET *ZDVDT(JL) |
| 1694 |
PVOL(JL,JK)= ZBET *ZDUDT(JL) + ZBET**2*ZDVDT(JL) |
| 1695 |
ENDDO |
| 1696 |
ENDIF |
| 1697 |
601 CONTINUE |
| 1698 |
|
| 1699 |
ENDIF |
| 1700 |
|
| 1701 |
RETURN |
| 1702 |
END |
| 1703 |
SUBROUTINE SUGWD(NLON,NLEV,paprs,pplay) |
| 1704 |
C |
| 1705 |
C**** *SUGWD* INITIALIZE COMMON YOEGWD CONTROLLING GRAVITY WAVE DRAG |
| 1706 |
C |
| 1707 |
C PURPOSE. |
| 1708 |
C -------- |
| 1709 |
C INITIALIZE YOEGWD, THE COMMON THAT CONTROLS THE |
| 1710 |
C GRAVITY WAVE DRAG PARAMETRIZATION. |
| 1711 |
C |
| 1712 |
C** INTERFACE. |
| 1713 |
C ---------- |
| 1714 |
C CALL *SUGWD* FROM *SUPHEC* |
| 1715 |
C ----- ------ |
| 1716 |
C |
| 1717 |
C EXPLICIT ARGUMENTS : |
| 1718 |
C -------------------- |
| 1719 |
C PSIG : VERTICAL COORDINATE TABLE |
| 1720 |
C NLEV : NUMBER OF MODEL LEVELS |
| 1721 |
C |
| 1722 |
C IMPLICIT ARGUMENTS : |
| 1723 |
C -------------------- |
| 1724 |
C COMMON YOEGWD |
| 1725 |
C |
| 1726 |
C METHOD. |
| 1727 |
C ------- |
| 1728 |
C SEE DOCUMENTATION |
| 1729 |
C |
| 1730 |
C EXTERNALS. |
| 1731 |
C ---------- |
| 1732 |
C NONE |
| 1733 |
C |
| 1734 |
C REFERENCE. |
| 1735 |
C ---------- |
| 1736 |
C ECMWF Research Department documentation of the IFS |
| 1737 |
C |
| 1738 |
C AUTHOR. |
| 1739 |
C ------- |
| 1740 |
C MARTIN MILLER *ECMWF* |
| 1741 |
C |
| 1742 |
C MODIFICATIONS. |
| 1743 |
C -------------- |
| 1744 |
C ORIGINAL : 90-01-01 |
| 1745 |
C ------------------------------------------------------------------ |
| 1746 |
use yoegwd |
| 1747 |
implicit none |
| 1748 |
C |
| 1749 |
C ----------------------------------------------------------------- |
| 1750 |
C ---------------------------------------------------------------- |
| 1751 |
C |
| 1752 |
integer nlon,nlev, jk |
| 1753 |
REAL, intent(in):: paprs(nlon,nlev+1) |
| 1754 |
REAL, intent(in):: pplay(nlon,nlev) |
| 1755 |
real zpr,zstra,zsigt,zpm1r |
| 1756 |
C |
| 1757 |
C* 1. SET THE VALUES OF THE PARAMETERS |
| 1758 |
C -------------------------------- |
| 1759 |
C |
| 1760 |
100 CONTINUE |
| 1761 |
C |
| 1762 |
PRINT *,' DANS SUGWD NLEV=',NLEV |
| 1763 |
GHMAX=10000. |
| 1764 |
C |
| 1765 |
ZPR=100000. |
| 1766 |
ZSTRA=0.1 |
| 1767 |
ZSIGT=0.94 |
| 1768 |
cold ZPR=80000. |
| 1769 |
cold ZSIGT=0.85 |
| 1770 |
C |
| 1771 |
DO 110 JK=1,NLEV |
| 1772 |
ZPM1R=pplay(nlon/2,jk)/paprs(nlon/2,1) |
| 1773 |
IF(ZPM1R.GE.ZSIGT)THEN |
| 1774 |
nktopg=JK |
| 1775 |
ENDIF |
| 1776 |
ZPM1R=pplay(nlon/2,jk)/paprs(nlon/2,1) |
| 1777 |
IF(ZPM1R.GE.ZSTRA)THEN |
| 1778 |
NSTRA=JK |
| 1779 |
ENDIF |
| 1780 |
110 CONTINUE |
| 1781 |
c |
| 1782 |
c inversion car dans orodrag on compte les niveaux a l'envers |
| 1783 |
nktopg=nlev-nktopg+1 |
| 1784 |
nstra=nlev-nstra |
| 1785 |
print *,' DANS SUGWD nktopg=', nktopg |
| 1786 |
print *,' DANS SUGWD nstra=', nstra |
| 1787 |
C |
| 1788 |
GSIGCR=0.80 |
| 1789 |
C |
| 1790 |
GKDRAG=0.2 |
| 1791 |
GRAHILO=1. |
| 1792 |
GRCRIT=0.01 |
| 1793 |
GFRCRIT=1.0 |
| 1794 |
GKWAKE=0.50 |
| 1795 |
C |
| 1796 |
GKLIFT=0.50 |
| 1797 |
GVCRIT =0.0 |
| 1798 |
C |
| 1799 |
C |
| 1800 |
C ---------------------------------------------------------------- |
| 1801 |
C |
| 1802 |
C* 2. SET VALUES OF SECURITY PARAMETERS |
| 1803 |
C --------------------------------- |
| 1804 |
C |
| 1805 |
200 CONTINUE |
| 1806 |
C |
| 1807 |
GVSEC=0.10 |
| 1808 |
GSSEC=1.E-12 |
| 1809 |
C |
| 1810 |
GTSEC=1.E-07 |
| 1811 |
C |
| 1812 |
C ---------------------------------------------------------------- |
| 1813 |
C |
| 1814 |
RETURN |
| 1815 |
END |
Parent Directory
| 
Revision Log