Changeset 24 for trunk/src/verttransform_gfs.f90
- Timestamp:
- May 23, 2014, 11:48:41 AM (10 years ago)
- File:
-
- 1 edited
-
trunk/src/verttransform_gfs.f90 (modified) (22 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/verttransform_gfs.f90
r4 r24 21 21 22 22 subroutine verttransform(n,uuh,vvh,wwh,pvh) 23 ! 23 ! i i i i i 24 24 !***************************************************************************** 25 25 ! * … … 71 71 real :: f_qvsat,pressure 72 72 real :: rh,lsp,convp 73 real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)73 real :: rhoh(nuvzmax),pinmconv(nzmax) 74 74 real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi 75 75 real :: xlon,ylat,xlonr,dzdx,dzdy 76 real :: dzdx1,dzdx2,dzdy1,dzdy2 76 real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf 77 77 real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy 78 78 real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax) … … 92 92 ! If verttransform is called the first time, initialize heights of the * 93 93 ! z levels in meter. The heights are the heights of model levels, where * 94 ! u,v,T and qv are given, and of the interfaces, where w is given. So, * 95 ! the vertical resolution in the z system is doubled. As reference point,* 96 ! the lower left corner of the grid is used. * 97 ! Unlike in the eta system, no difference between heights for u,v and * 98 ! heights for w exists. * 94 ! u,v,T and qv are given. * 99 95 !************************************************************************* 100 96 … … 118 114 119 115 tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ & 120 116 ps(ixm,jym,1,n)) 121 117 pold=ps(ixm,jym,1,n) 122 118 height(1)=0. … … 126 122 tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n)) 127 123 128 129 ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled130 ! upon the transformation to z levels. In order to save computer memory, this is131 ! not done anymore in the standard version. However, this option can still be132 ! switched on by replacing the following lines with those below, that are133 ! currently commented out.134 ! Note that two more changes are necessary in this subroutine below.135 ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.136 !*****************************************************************************137 138 124 if (abs(tv-tvold).gt.0.2) then 139 height(kz)= & 140 height(kz-1)+const*log(pold/pint)* & 141 (tv-tvold)/log(tv/tvold) 125 height(kz)=height(kz-1)+const*log(pold/pint)* & 126 (tv-tvold)/log(tv/tvold) 142 127 else 143 height(kz)=height(kz-1)+ & 144 const*log(pold/pint)*tv 128 height(kz)=height(kz-1)+const*log(pold/pint)*tv 145 129 endif 146 147 ! Switch on following lines to use doubled vertical resolution148 !*************************************************************149 ! if (abs(tv-tvold).gt.0.2) then150 ! height((kz-1)*2)=151 ! + height(max((kz-2)*2,1))+const*log(pold/pint)*152 ! + (tv-tvold)/log(tv/tvold)153 ! else154 ! height((kz-1)*2)=height(max((kz-2)*2,1))+155 ! + const*log(pold/pint)*tv156 ! endif157 ! End doubled vertical resolution158 130 159 131 tvold=tv 160 132 pold=pint 161 133 end do 162 163 164 ! Switch on following lines to use doubled vertical resolution165 !*************************************************************166 ! do 7 kz=3,nz-1,2167 ! height(kz)=0.5*(height(kz-1)+height(kz+1))168 ! height(nz)=height(nz-1)+height(nz-1)-height(nz-2)169 ! End doubled vertical resolution170 134 171 135 … … 198 162 llev = 0 199 163 do i=1,nuvz 200 if (ps(ix,jy,1,n).lt.akz(i)) llev=i164 if (ps(ix,jy,1,n).lt.akz(i)) llev=i 201 165 end do 202 166 llev = llev+1 … … 212 176 tvold=tth(ix,jy,llev,n)*(1.+0.608*qvh(ix,jy,llev,n)) 213 177 pold=akz(llev) 214 uvzlev(llev)=0.215 178 wzlev(llev)=0. 216 179 uvwzlev(ix,jy,llev)=0. … … 223 186 224 187 if (abs(tv-tvold).gt.0.2) then 225 uv zlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &226 188 uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* & 189 (tv-tvold)/log(tv/tvold) 227 190 else 228 uv zlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv191 uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv 229 192 endif 230 wzlev(kz)=uvzlev(kz) 231 uvwzlev(ix,jy,kz)=uvzlev(kz) 193 wzlev(kz)=uvwzlev(ix,jy,kz) 232 194 233 195 tvold=tv 234 196 pold=pint 235 197 end do 236 237 238 ! Switch on following lines to use doubled vertical resolution239 ! Switch off the three lines above.240 !*************************************************************241 !22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)242 ! do 23 kz=2,nwz243 !23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)244 ! End doubled vertical resolution245 198 246 199 ! pinmconv=(h2-h1)/(p2-p1) … … 279 232 do iz=2,nz-1 280 233 do kz=kmin,nuvz 281 if(height(iz).gt.uv zlev(nuvz)) then234 if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then 282 235 uu(ix,jy,iz,n)=uu(ix,jy,nz,n) 283 236 vv(ix,jy,iz,n)=vv(ix,jy,nz,n) … … 289 242 goto 30 290 243 endif 291 if ((height(iz).gt.uv zlev(kz-1)).and. &292 (height(iz).le.uvzlev(kz))) then293 dz1=height(iz)-uvzlev(kz-1)294 dz2=uvzlev(kz)-height(iz)295 dz=dz1+dz2296 uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz297 vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz298 tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &299 300 qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &301 302 pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz303 rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz304 pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz244 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 245 (height(iz).le.uvwzlev(ix,jy,kz))) then 246 dz1=height(iz)-uvwzlev(ix,jy,kz-1) 247 dz2=uvwzlev(ix,jy,kz)-height(iz) 248 dz=dz1+dz2 249 uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz 250 vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz 251 tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 & 252 +tth(ix,jy,kz,n)*dz1)/dz 253 qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 & 254 +qvh(ix,jy,kz,n)*dz1)/dz 255 pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz 256 rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz 257 pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz 305 258 endif 306 259 end do … … 318 271 do kz=kmin,nwz 319 272 if ((height(iz).gt.wzlev(kz-1)).and. & 320 (height(iz).le.wzlev(kz))) then 321 dz1=height(iz)-wzlev(kz-1) 322 dz2=wzlev(kz)-height(iz) 323 dz=dz1+dz2 324 ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 & 325 +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz 326 273 (height(iz).le.wzlev(kz))) then 274 dz1=height(iz)-wzlev(kz-1) 275 dz2=wzlev(kz)-height(iz) 276 dz=dz1+dz2 277 ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 & 278 +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz 327 279 endif 328 280 end do … … 337 289 do kz=2,nz-1 338 290 drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ & 339 291 (height(kz+1)-height(kz-1)) 340 292 end do 341 293 drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n) … … 351 303 352 304 do jy=1,ny-2 305 cosf=cos((real(jy)*dy+ylat0)*pi180) 353 306 do ix=1,nx-2 354 307 … … 367 320 do iz=2,nz-1 368 321 369 ui=uu(ix,jy,iz,n)*dxconst/cos ((real(jy)*dy+ylat0)*pi180)322 ui=uu(ix,jy,iz,n)*dxconst/cosf 370 323 vi=vv(ix,jy,iz,n)*dyconst 371 324 372 325 do kz=kmin,nz 373 326 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 374 327 (height(iz).le.uvwzlev(ix,jy,kz))) then 375 328 dz1=height(iz)-uvwzlev(ix,jy,kz-1) 376 329 dz2=uvwzlev(ix,jy,kz)-height(iz) … … 426 379 xlon=xlon0+real(nx/2-1)*dx 427 380 xlonr=xlon*pi/180. 428 ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ & 429 vv(nx/2-1,nymin1,iz,n)**2) 430 if(vv(nx/2-1,nymin1,iz,n).lt.0.) then 431 ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ & 432 vv(nx/2-1,nymin1,iz,n))-xlonr 381 ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2) 382 if (vv(nx/2-1,nymin1,iz,n).lt.0.) then 383 ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr 433 384 elseif (vv(nx/2-1,nymin1,iz,n).gt.0.) then 434 385 ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ & 435 386 vv(nx/2-1,nymin1,iz,n))-xlonr 436 387 else 437 388 ddpol=pi/2-xlonr … … 446 397 uuaux=-ffpol*sin(xlonr+ddpol) 447 398 vvaux=-ffpol*cos(xlonr+ddpol) 448 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, & 449 vvpolaux) 450 399 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux) 451 400 jy=nymin1 452 401 do ix=0,nxmin1 … … 460 409 ! ward parallel of latitude 461 410 462 do iz=1,nz411 do iz=1,nz 463 412 wdummy=0. 464 413 jy=ny-2 … … 471 420 ww(ix,jy,iz,n)=wdummy 472 421 end do 473 end do422 end do 474 423 475 424 endif … … 487 436 do iz=1,nz 488 437 call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), & 489 vv(ix,jy,iz,n),uupol(ix,jy,iz,n), & 490 vvpol(ix,jy,iz,n)) 438 vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n)) 491 439 end do 492 440 end do … … 498 446 xlon=xlon0+real(nx/2-1)*dx 499 447 xlonr=xlon*pi/180. 500 ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ & 501 vv(nx/2-1,0,iz,n)**2) 448 ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2) 502 449 if(vv(nx/2-1,0,iz,n).lt.0.) then 503 ddpol=atan(uu(nx/2-1,0,iz,n)/ & 504 vv(nx/2-1,0,iz,n))+xlonr 450 ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr 505 451 elseif (vv(nx/2-1,0,iz,n).gt.0.) then 506 ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ & 507 vv(nx/2-1,0,iz,n))-xlonr 452 ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))-xlonr 508 453 else 509 454 ddpol=pi/2-xlonr … … 518 463 uuaux=+ffpol*sin(xlonr-ddpol) 519 464 vvaux=-ffpol*cos(xlonr-ddpol) 520 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, & 521 vvpolaux) 465 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux) 522 466 523 467 jy=0 … … 562 506 clouds(ix,jy,kz,n)=0 563 507 if (rh.gt.0.8) then ! in cloud 564 if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation 565 rain_cloud_above=1 566 cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ & 567 height(kz)-height(kz-1) 568 if (lsp.ge.convp) then 569 clouds(ix,jy,kz,n)=3 ! lsp dominated rainout 570 else 571 clouds(ix,jy,kz,n)=2 ! convp dominated rainout 572 endif 573 else ! no precipitation 574 clouds(ix,jy,kz,n)=1 ! cloud 575 endif 508 if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation 509 rain_cloud_above=1 510 cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1) 511 if (lsp.ge.convp) then 512 clouds(ix,jy,kz,n)=3 ! lsp dominated rainout 513 else 514 clouds(ix,jy,kz,n)=2 ! convp dominated rainout 515 endif 516 else ! no precipitation 517 clouds(ix,jy,kz,n)=1 ! cloud 518 endif 576 519 else ! no cloud 577 578 579 580 581 582 583 520 if (rain_cloud_above.eq.1) then ! scavenging 521 if (lsp.ge.convp) then 522 clouds(ix,jy,kz,n)=5 ! lsp dominated washout 523 else 524 clouds(ix,jy,kz,n)=4 ! convp dominated washout 525 endif 526 endif 584 527 endif 585 528 end do
Note: See TracChangeset
for help on using the changeset viewer.
