Changeset 24 for trunk/src/verttransform.f90
- Timestamp:
- May 23, 2014, 11:48:41 AM (10 years ago)
- File:
-
- 1 edited
-
trunk/src/verttransform.f90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/verttransform.f90
r20 r24 21 21 22 22 subroutine verttransform(n,uuh,vvh,wwh,pvh) 23 ! 23 ! i i i i i 24 24 !***************************************************************************** 25 25 ! * … … 49 49 ! Sabine Eckhardt, March 2007 50 50 ! added the variable cloud for use with scavenging - descr. in com_mod 51 ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification52 ! note that also other subroutines are affected by the fix53 51 !***************************************************************************** 54 52 ! * … … 72 70 73 71 integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym 74 integer :: rain_cloud_above,kz_inv !SE 75 integer icloudtop !PS 72 integer :: rain_cloud_above,kz_inv 76 73 real :: f_qvsat,pressure 77 !real :: rh,lsp,convp 78 real :: rh,lsp,convp,prec,rhmin 79 real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax) 74 real :: rh,lsp,convp 75 real :: rhoh(nuvzmax),pinmconv(nzmax) 80 76 real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi 81 77 real :: xlon,ylat,xlonr,dzdx,dzdy 82 real :: dzdx1,dzdx2,dzdy1,dzdy2, precmin78 real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf 83 79 real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy 84 80 real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax) … … 87 83 real :: wwh(0:nxmax-1,0:nymax-1,nwzmax) 88 84 real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax) 89 logical lconvectprec90 85 real,parameter :: const=r_air/ga 91 parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics92 86 93 87 logical :: init = .true. … … 97 91 ! If verttransform is called the first time, initialize heights of the * 98 92 ! z levels in meter. The heights are the heights of model levels, where * 99 ! u,v,T and qv are given, and of the interfaces, where w is given. So, * 100 ! the vertical resolution in the z system is doubled. As reference point,* 101 ! the lower left corner of the grid is used. * 102 ! Unlike in the eta system, no difference between heights for u,v and * 103 ! heights for w exists. * 93 ! u,v,T and qv are given. * 104 94 !************************************************************************* 105 106 107 ! do 897 kz=1,nuvz108 ! write (*,*) 'akz: ',akz(kz),'bkz',bkz(kz)109 !897 continue110 95 111 96 if (init) then … … 113 98 ! Search for a point with high surface pressure (i.e. not above significant topography) 114 99 ! Then, use this point to construct a reference z profile, to be used at all times 115 !***************************************************************************** 100 !************************************************************************************** 116 101 117 102 do jy=0,nymin1 … … 128 113 129 114 tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ & 130 115 ps(ixm,jym,1,n)) 131 116 pold=ps(ixm,jym,1,n) 132 117 height(1)=0. … … 136 121 tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n)) 137 122 138 139 ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled140 ! upon the transformation to z levels. In order to save computer memory, this is141 ! not done anymore in the standard version. However, this option can still be142 ! switched on by replacing the following lines with those below, that are143 ! currently commented out.144 ! Note that two more changes are necessary in this subroutine below.145 ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.146 !*****************************************************************************147 148 123 if (abs(tv-tvold).gt.0.2) then 149 height(kz)= & 150 height(kz-1)+const*log(pold/pint)* & 151 (tv-tvold)/log(tv/tvold) 124 height(kz)=height(kz-1)+const*log(pold/pint)* & 125 (tv-tvold)/log(tv/tvold) 152 126 else 153 height(kz)=height(kz-1)+ & 154 const*log(pold/pint)*tv 127 height(kz)=height(kz-1)+const*log(pold/pint)*tv 155 128 endif 156 157 ! Switch on following lines to use doubled vertical resolution158 !*************************************************************159 ! if (abs(tv-tvold).gt.0.2) then160 ! height((kz-1)*2)=161 ! + height(max((kz-2)*2,1))+const*log(pold/pint)*162 ! + (tv-tvold)/log(tv/tvold)163 ! else164 ! height((kz-1)*2)=height(max((kz-2)*2,1))+165 ! + const*log(pold/pint)*tv166 ! endif167 ! End doubled vertical resolution168 129 169 130 tvold=tv 170 131 pold=pint 171 132 end do 172 173 174 ! Switch on following lines to use doubled vertical resolution175 !*************************************************************176 ! do 7 kz=3,nz-1,2177 ! height(kz)=0.5*(height(kz-1)+height(kz+1))178 ! height(nz)=height(nz-1)+height(nz-1)-height(nz-2)179 ! End doubled vertical resolution180 133 181 134 … … 204 157 do jy=0,nymin1 205 158 do ix=0,nxmin1 206 tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ & 207 ps(ix,jy,1,n)) 159 tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ps(ix,jy,1,n)) 208 160 pold=ps(ix,jy,1,n) 209 uv zlev(1)=0.161 uvwzlev(ix,jy,1)=0. 210 162 wzlev(1)=0. 211 163 rhoh(1)=pold/(r_air*tvold) … … 221 173 222 174 if (abs(tv-tvold).gt.0.2) then 223 uv zlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &224 175 uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* & 176 (tv-tvold)/log(tv/tvold) 225 177 else 226 uv zlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv178 uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv 227 179 endif 228 180 … … 233 185 234 186 do kz=2,nwz-1 235 wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2. 236 end do 237 wzlev(nwz)=wzlev(nwz-1)+ & 238 uvzlev(nuvz)-uvzlev(nuvz-1) 239 240 uvwzlev(ix,jy,1)=0.0 241 do kz=2,nuvz 242 uvwzlev(ix,jy,kz)=uvzlev(kz) 243 end do 244 245 ! Switch on following lines to use doubled vertical resolution 246 ! Switch off the three lines above. 247 !************************************************************* 248 !22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz) 249 ! do 23 kz=2,nwz 250 !23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz) 251 ! End doubled vertical resolution 187 wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2. 188 end do 189 wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1) 252 190 253 191 ! pinmconv=(h2-h1)/(p2-p1) 254 192 255 193 pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ & 256 257 194 ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- & 195 (aknew(1)+bknew(1)*ps(ix,jy,1,n))) 258 196 do kz=2,nz-1 259 197 pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ & 260 261 198 ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- & 199 (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n))) 262 200 end do 263 201 pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ & 264 265 202 ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- & 203 (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n))) 266 204 267 205 ! Levels, where u,v,t and q are given … … 283 221 do iz=2,nz-1 284 222 do kz=kmin,nuvz 285 if(height(iz).gt.uv zlev(nuvz)) then223 if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then 286 224 uu(ix,jy,iz,n)=uu(ix,jy,nz,n) 287 225 vv(ix,jy,iz,n)=vv(ix,jy,nz,n) … … 292 230 goto 30 293 231 endif 294 if ((height(iz).gt.uv zlev(kz-1)).and. &295 (height(iz).le.uvzlev(kz))) then296 dz1=height(iz)-uv zlev(kz-1)297 dz2=uv zlev(kz)-height(iz)232 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 233 (height(iz).le.uvwzlev(ix,jy,kz))) then 234 dz1=height(iz)-uvwzlev(ix,jy,kz-1) 235 dz2=uvwzlev(ix,jy,kz)-height(iz) 298 236 dz=dz1+dz2 299 237 uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz … … 339 277 340 278 drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ & 341 279 (height(2)-height(1)) 342 280 do kz=2,nz-1 343 281 drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ & 344 282 (height(kz+1)-height(kz-1)) 345 283 end do 346 284 drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n) … … 356 294 357 295 do jy=1,ny-2 296 cosf=cos((real(jy)*dy+ylat0)*pi180) 358 297 do ix=1,nx-2 359 298 … … 361 300 do iz=2,nz-1 362 301 363 ui=uu(ix,jy,iz,n)*dxconst/cos ((real(jy)*dy+ylat0)*pi180)302 ui=uu(ix,jy,iz,n)*dxconst/cosf 364 303 vi=vv(ix,jy,iz,n)*dyconst 365 304 366 305 do kz=kmin,nz 367 306 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 368 307 (height(iz).le.uvwzlev(ix,jy,kz))) then 369 308 dz1=height(iz)-uvwzlev(ix,jy,kz-1) 370 309 dz2=uvwzlev(ix,jy,kz)-height(iz) … … 409 348 do iz=1,nz 410 349 call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), & 411 vv(ix,jy,iz,n),uupol(ix,jy,iz,n), & 412 vvpol(ix,jy,iz,n)) 350 vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n)) 413 351 end do 414 352 end do … … 424 362 xlon=xlon0+real(nx/2-1)*dx 425 363 xlonr=xlon*pi/180. 426 ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ & 427 vv(nx/2-1,nymin1,iz,n)**2) 364 ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2) 428 365 if (vv(nx/2-1,nymin1,iz,n).lt.0.) then 429 ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ & 430 vv(nx/2-1,nymin1,iz,n))-xlonr 366 ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr 431 367 else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then 432 368 ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ & 433 369 vv(nx/2-1,nymin1,iz,n))-xlonr 434 370 else 435 371 ddpol=pi/2-xlonr … … 444 380 uuaux=-ffpol*sin(xlonr+ddpol) 445 381 vvaux=-ffpol*cos(xlonr+ddpol) 446 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, & 447 vvpolaux) 448 382 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux) 449 383 jy=nymin1 450 384 do ix=0,nxmin1 … … 485 419 do iz=1,nz 486 420 call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), & 487 vv(ix,jy,iz,n),uupol(ix,jy,iz,n), & 488 vvpol(ix,jy,iz,n)) 421 vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n)) 489 422 end do 490 423 end do … … 499 432 xlon=xlon0+real(nx/2-1)*dx 500 433 xlonr=xlon*pi/180. 501 ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ & 502 vv(nx/2-1,0,iz,n)**2) 434 ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2) 503 435 if (vv(nx/2-1,0,iz,n).lt.0.) then 504 ddpol=atan(uu(nx/2-1,0,iz,n)/ & 505 vv(nx/2-1,0,iz,n))+xlonr 436 ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr 506 437 else if (vv(nx/2-1,0,iz,n).gt.0.) then 507 ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ & 508 vv(nx/2-1,0,iz,n))+xlonr 438 ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr 509 439 else 510 440 ddpol=pi/2-xlonr … … 519 449 uuaux=+ffpol*sin(xlonr-ddpol) 520 450 vvaux=-ffpol*cos(xlonr-ddpol) 521 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, & 522 vvpolaux) 451 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux) 523 452 524 453 jy=0 … … 551 480 ! create a cloud and rainout/washout field, clouds occur where rh>80% 552 481 ! total cloudheight is stored at level 0 553 554 555 556 482 do jy=0,nymin1 557 483 do ix=0,nxmin1 558 559 560 561 ! rain_cloud_above=0 562 ! lsp=lsprec(ix,jy,1,n) 563 ! convp=convprec(ix,jy,1,n) 564 ! cloudsh(ix,jy,n)=0 565 ! do kz_inv=1,nz-1 566 ! kz=nz-kz_inv+1 567 ! pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n) 568 ! rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n)) 569 ! clouds(ix,jy,kz,n)=0 570 ! if (rh.gt.0.8) then ! in cloud 571 ! if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation 572 ! rain_cloud_above=1 573 ! cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ & 574 ! height(kz)-height(kz-1) 575 ! if (lsp.ge.convp) then 576 ! clouds(ix,jy,kz,n)=3 ! lsp dominated rainout 577 ! else 578 ! clouds(ix,jy,kz,n)=2 ! convp dominated rainout 579 ! endif 580 ! else ! no precipitation 581 ! clouds(ix,jy,kz,n)=1 ! cloud 582 ! endif 583 ! else ! no cloud 584 ! if (rain_cloud_above.eq.1) then ! scavenging 585 ! if (lsp.ge.convp) then 586 ! clouds(ix,jy,kz,n)=5 ! lsp dominated washout 587 ! else 588 ! clouds(ix,jy,kz,n)=4 ! convp dominated washout 589 ! endif 590 ! endif 591 ! endif 592 ! end do 593 594 595 ! PS 3012 596 597 lsp=lsprec(ix,jy,1,n) 598 convp=convprec(ix,jy,1,n) 599 prec=lsp+convp 600 if (lsp.gt.convp) then ! prectype='lsp' 601 lconvectprec = .false. 602 else ! prectype='cp ' 603 lconvectprec = .true. 604 endif 605 rhmin = 0.90 ! standard condition for presence of clouds 606 !PS note that original by Sabine Eckhart was 80% 607 !PS however, for T<-20 C we consider saturation over ice 608 !PS so I think 90% should be enough 609 icloudbot(ix,jy,n)=icmv 610 icloudtop=icmv ! this is just a local variable 611 98 do kz=1,nz 612 pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n) 613 rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n)) 614 !ps if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz) 615 if (rh .gt. rhmin) then 616 if (icloudbot(ix,jy,n) .eq. icmv) then 617 icloudbot(ix,jy,n)=nint(height(kz)) 618 endif 619 icloudtop=nint(height(kz)) ! use int to save memory 484 rain_cloud_above=0 485 lsp=lsprec(ix,jy,1,n) 486 convp=convprec(ix,jy,1,n) 487 cloudsh(ix,jy,n)=0 488 do kz_inv=1,nz-1 489 kz=nz-kz_inv+1 490 pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n) 491 rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n)) 492 clouds(ix,jy,kz,n)=0 493 if (rh.gt.0.8) then ! in cloud 494 if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation 495 rain_cloud_above=1 496 cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ & 497 height(kz)-height(kz-1) 498 if (lsp.ge.convp) then 499 clouds(ix,jy,kz,n)=3 ! lsp dominated rainout 500 else 501 clouds(ix,jy,kz,n)=2 ! convp dominated rainout 502 endif 503 else ! no precipitation 504 clouds(ix,jy,kz,n)=1 ! cloud 620 505 endif 621 enddo 622 623 !PS try to get a cloud thicker than 50 m 624 !PS if there is at least .01 mm/h - changed to 0.002 and put into 625 !PS parameter precpmin 626 if ((icloudbot(ix,jy,n) .eq. icmv .or. & 627 icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. & 628 prec .gt. precmin) then 629 rhmin = rhmin - 0.05 630 if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum. 631 endif 632 !PS implement a rough fix for badly represented convection 633 !PS is based on looking at a limited set of comparison data 634 if (lconvectprec .and. icloudtop .lt. 6000 .and. & 635 prec .gt. precmin) then 636 637 if (convp .lt. 0.1) then 638 icloudbot(ix,jy,n) = 500 639 icloudtop = 8000 640 else 641 icloudbot(ix,jy,n) = 0 642 icloudtop = 10000 506 else ! no cloud 507 if (rain_cloud_above.eq.1) then ! scavenging 508 if (lsp.ge.convp) then 509 clouds(ix,jy,kz,n)=5 ! lsp dominated washout 510 else 511 clouds(ix,jy,kz,n)=4 ! convp dominated washout 512 endif 643 513 endif 644 endif 645 if (icloudtop .ne. icmv) then 646 icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n) 647 else 648 icloudthck(ix,jy,n) = icmv 649 endif 650 !PS get rid of too thin clouds 651 if (icloudthck(ix,jy,n) .lt. 50) then 652 icloudbot(ix,jy,n)=icmv 653 icloudthck(ix,jy,n)=icmv 654 endif 655 656 514 endif 515 end do 657 516 end do 658 517 end do 659 518 660 !do 102 kz=1,nuvz661 !write(an,'(i02)') kz+10662 !write(*,*) nuvz,nymin1,nxmin1,'--',an,'--'663 !open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')664 !do 101 jy=0,nymin1665 ! write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)666 !101 continue667 ! close(4)668 !102 continue669 670 ! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')671 ! do 103 jy=0,nymin1672 ! write (4,*)673 !+ (height(kz),kz=1,nuvz)674 !103 continue675 ! close(4)676 677 !open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')678 ! do 104 jy=0,nymin1679 ! write (4,*)680 !+ (r_air*tt(ix,jy,1,n)*rho(ix,jy,1,n),ix=1,nxmin1)681 !104 continue682 ! close(4)683 684 685 519 end subroutine verttransform
Note: See TracChangeset
for help on using the changeset viewer.
