Ignore:
Timestamp:
May 23, 2014, 11:48:41 AM (8 years ago)
Author:
igpis
Message:

version 9.2 beta. Changes from HH, AST, MC, NIK, IP. Changes in vert transform. New SPECIES input includes scavenging coefficients

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/verttransform.f90

    r20 r24  
    2121
    2222subroutine verttransform(n,uuh,vvh,wwh,pvh)
    23   !                         i  i   i   i   i
     23  !                      i  i   i   i   i
    2424  !*****************************************************************************
    2525  !                                                                            *
     
    4949  ! Sabine Eckhardt, March 2007
    5050  ! added the variable cloud for use with scavenging - descr. in com_mod
    51   ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
    52   ! note that also other subroutines are affected by the fix
    5351  !*****************************************************************************
    5452  !                                                                            *
     
    7270
    7371  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
    7673  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)
    8076  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
    8177  real :: xlon,ylat,xlonr,dzdx,dzdy
    82   real :: dzdx1,dzdx2,dzdy1,dzdy2, precmin
     78  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
    8379  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
    8480  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
     
    8783  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
    8884  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
    89   logical lconvectprec
    9085  real,parameter :: const=r_air/ga
    91   parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
    9286
    9387  logical :: init = .true.
     
    9791  ! If verttransform is called the first time, initialize heights of the   *
    9892  ! 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.                                                *
    10494  !*************************************************************************
    105 
    106 
    107   !      do 897 kz=1,nuvz
    108   !       write (*,*) 'akz: ',akz(kz),'bkz',bkz(kz)
    109   !897     continue
    11095
    11196  if (init) then
     
    11398  ! Search for a point with high surface pressure (i.e. not above significant topography)
    11499  ! Then, use this point to construct a reference z profile, to be used at all times
    115   !*****************************************************************************
     100  !**************************************************************************************
    116101
    117102    do jy=0,nymin1
     
    128113
    129114    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
    130          ps(ixm,jym,1,n))
     115    ps(ixm,jym,1,n))
    131116    pold=ps(ixm,jym,1,n)
    132117    height(1)=0.
     
    136121      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
    137122
    138 
    139   ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
    140   ! upon the transformation to z levels. In order to save computer memory, this is
    141   ! not done anymore in the standard version. However, this option can still be
    142   ! switched on by replacing the following lines with those below, that are
    143   ! 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 
    148123      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)
    152126      else
    153         height(kz)=height(kz-1)+ &
    154              const*log(pold/pint)*tv
     127        height(kz)=height(kz-1)+const*log(pold/pint)*tv
    155128      endif
    156 
    157   ! Switch on following lines to use doubled vertical resolution
    158   !*************************************************************
    159   !    if (abs(tv-tvold).gt.0.2) then
    160   !      height((kz-1)*2)=
    161   !    +      height(max((kz-2)*2,1))+const*log(pold/pint)*
    162   !    +      (tv-tvold)/log(tv/tvold)
    163   !    else
    164   !      height((kz-1)*2)=height(max((kz-2)*2,1))+
    165   !    +      const*log(pold/pint)*tv
    166   !    endif
    167   ! End doubled vertical resolution
    168129
    169130      tvold=tv
    170131      pold=pint
    171132    end do
    172 
    173 
    174   ! Switch on following lines to use doubled vertical resolution
    175   !*************************************************************
    176   !  do 7 kz=3,nz-1,2
    177   !    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 resolution
    180133
    181134
     
    204157  do jy=0,nymin1
    205158    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))
    208160      pold=ps(ix,jy,1,n)
    209       uvzlev(1)=0.
     161      uvwzlev(ix,jy,1)=0.
    210162      wzlev(1)=0.
    211163      rhoh(1)=pold/(r_air*tvold)
     
    221173
    222174        if (abs(tv-tvold).gt.0.2) then
    223           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
    224                (tv-tvold)/log(tv/tvold)
     175          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
     176          (tv-tvold)/log(tv/tvold)
    225177        else
    226           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
     178          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
    227179        endif
    228180
     
    233185
    234186      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)
    252190
    253191  ! pinmconv=(h2-h1)/(p2-p1)
    254192
    255193      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
    256            ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
    257            (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
     194      ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
     195      (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
    258196      do kz=2,nz-1
    259197        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
    260              ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
    261              (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
     198        ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
     199        (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
    262200      end do
    263201      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
    264            ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
    265            (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
     202      ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
     203      (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
    266204
    267205  ! Levels, where u,v,t and q are given
     
    283221      do iz=2,nz-1
    284222        do kz=kmin,nuvz
    285           if(height(iz).gt.uvzlev(nuvz)) then
     223          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
    286224            uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
    287225            vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
     
    292230            goto 30
    293231          endif
    294           if ((height(iz).gt.uvzlev(kz-1)).and. &
    295                (height(iz).le.uvzlev(kz))) then
    296            dz1=height(iz)-uvzlev(kz-1)
    297            dz2=uvzlev(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)
    298236           dz=dz1+dz2
    299237           uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
     
    339277
    340278      drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
    341            (height(2)-height(1))
     279      (height(2)-height(1))
    342280      do kz=2,nz-1
    343281        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
    344              (height(kz+1)-height(kz-1))
     282        (height(kz+1)-height(kz-1))
    345283      end do
    346284      drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
     
    356294
    357295  do jy=1,ny-2
     296    cosf=cos((real(jy)*dy+ylat0)*pi180)
    358297    do ix=1,nx-2
    359298
     
    361300      do iz=2,nz-1
    362301
    363         ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
     302        ui=uu(ix,jy,iz,n)*dxconst/cosf
    364303        vi=vv(ix,jy,iz,n)*dyconst
    365304
    366305        do kz=kmin,nz
    367306          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
    368                (height(iz).le.uvwzlev(ix,jy,kz))) then
     307          (height(iz).le.uvwzlev(ix,jy,kz))) then
    369308            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
    370309            dz2=uvwzlev(ix,jy,kz)-height(iz)
     
    409348        do iz=1,nz
    410349          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))
    413351        end do
    414352      end do
     
    424362      xlon=xlon0+real(nx/2-1)*dx
    425363      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)
    428365      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
    431367      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
    432368        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
    433              vv(nx/2-1,nymin1,iz,n))-xlonr
     369        vv(nx/2-1,nymin1,iz,n))-xlonr
    434370      else
    435371        ddpol=pi/2-xlonr
     
    444380      uuaux=-ffpol*sin(xlonr+ddpol)
    445381      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)
    449383      jy=nymin1
    450384      do ix=0,nxmin1
     
    485419        do iz=1,nz
    486420          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))
    489422        end do
    490423      end do
     
    499432      xlon=xlon0+real(nx/2-1)*dx
    500433      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)
    503435      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
    506437      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
    509439      else
    510440        ddpol=pi/2-xlonr
     
    519449      uuaux=+ffpol*sin(xlonr-ddpol)
    520450      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)
    523452
    524453      jy=0
     
    551480  !   create a cloud and rainout/washout field, clouds occur where rh>80%
    552481  !   total cloudheight is stored at level 0
    553 
    554 
    555 
    556482  do jy=0,nymin1
    557483    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
    620505            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
    643513            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
    657516    end do
    658517  end do
    659518
    660   !do 102 kz=1,nuvz
    661   !write(an,'(i02)') kz+10
    662   !write(*,*) nuvz,nymin1,nxmin1,'--',an,'--'
    663   !open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')
    664   !do 101 jy=0,nymin1
    665   !    write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)
    666   !101   continue
    667   ! close(4)
    668   !102   continue
    669 
    670   ! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')
    671   ! do 103 jy=0,nymin1
    672   !     write (4,*)
    673   !+       (height(kz),kz=1,nuvz)
    674   !103   continue
    675   ! close(4)
    676 
    677   !open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')
    678   ! do 104 jy=0,nymin1
    679   !     write (4,*)
    680   !+       (r_air*tt(ix,jy,1,n)*rho(ix,jy,1,n),ix=1,nxmin1)
    681   !104   continue
    682   ! close(4)
    683 
    684 
    685519end subroutine verttransform
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG