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_nests.f90

    r4 r24  
    2121
    2222subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn)
    23   !                               i   i    i    i   i
     23  !                            i   i    i    i   i
    2424  !*****************************************************************************
    2525  !                                                                            *
     
    7272  integer :: rain_cloud_above,kz_inv
    7373  real :: f_qvsat,pressure,rh,lsp,convp
    74   real :: uvzlev(nuvzmax),wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)
     74  real :: wzlev(nwzmax),rhoh(nuvzmax),pinmconv(nzmax)
    7575  real :: uvwzlev(0:nxmaxn-1,0:nymaxn-1,nzmax)
    76   real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
     76  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi,cosf
    7777  real :: dzdx,dzdy
    7878  real :: dzdx1,dzdx2,dzdy1,dzdy2
     
    9696
    9797      tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l))/ &
    98            psn(ix,jy,1,n,l))
     98      psn(ix,jy,1,n,l))
    9999      pold=psn(ix,jy,1,n,l)
    100       uvzlev(1)=0.
     100      uvwzlev(ix,jy,1)=0.
    101101      wzlev(1)=0.
    102102      rhoh(1)=pold/(r_air*tvold)
     
    112112
    113113        if (abs(tv-tvold).gt.0.2) then
    114           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
    115                (tv-tvold)/log(tv/tvold)
     114          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
     115          (tv-tvold)/log(tv/tvold)
    116116        else
    117           uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
     117          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
    118118        endif
    119119
     
    124124
    125125      do kz=2,nwz-1
    126         wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
    127       end do
    128       wzlev(nwz)=wzlev(nwz-1)+ &
    129            uvzlev(nuvz)-uvzlev(nuvz-1)
    130 
    131   ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
    132   ! upon the transformation to z levels. In order to save computer memory, this is
    133   ! not done anymore in the standard version. However, this option can still be
    134   ! switched on by replacing the following lines with those below, that are
    135   ! currently commented out.
    136   ! Note that one change is also necessary in gridcheck.f,
    137   ! and three changes in verttransform.f
    138   !*****************************************************************************
    139       uvwzlev(ix,jy,1)=0.0
    140       do kz=2,nuvz
    141         uvwzlev(ix,jy,kz)=uvzlev(kz)
    142       end do
    143 
    144   ! Switch on following lines to use doubled vertical resolution
    145   ! Switch off the three lines above.
    146   !*************************************************************
    147   !22          uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
    148   !     do 23 kz=2,nwz
    149   !23          uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
    150   ! End doubled vertical resolution
     126        wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2.
     127      end do
     128      wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
    151129
    152130  ! pinmconv=(h2-h1)/(p2-p1)
    153131
    154132      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
    155            ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
    156            (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
     133      ((aknew(2)+bknew(2)*psn(ix,jy,1,n,l))- &
     134      (aknew(1)+bknew(1)*psn(ix,jy,1,n,l)))
    157135      do kz=2,nz-1
    158136        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
    159              ((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- &
    160              (aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l)))
     137        ((aknew(kz+1)+bknew(kz+1)*psn(ix,jy,1,n,l))- &
     138        (aknew(kz-1)+bknew(kz-1)*psn(ix,jy,1,n,l)))
    161139      end do
    162140      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
    163            ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
    164            (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
     141      ((aknew(nz)+bknew(nz)*psn(ix,jy,1,n,l))- &
     142      (aknew(nz-1)+bknew(nz-1)*psn(ix,jy,1,n,l)))
    165143
    166144
     
    183161      do iz=2,nz-1
    184162        do kz=kmin,nuvz
    185           if(height(iz).gt.uvzlev(nuvz)) then
     163          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
    186164            uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
    187165            vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
     
    192170            goto 30
    193171          endif
    194           if ((height(iz).gt.uvzlev(kz-1)).and. &
    195                (height(iz).le.uvzlev(kz))) then
    196            dz1=height(iz)-uvzlev(kz-1)
    197            dz2=uvzlev(kz)-height(iz)
     172          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
     173          (height(iz).le.uvwzlev(ix,jy,kz))) then
     174           dz1=height(iz)-uvwzlev(ix,jy,kz-1)
     175           dz2=uvwzlev(ix,jy,kz)-height(iz)
    198176           dz=dz1+dz2
    199177           uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
    200                 uuhn(ix,jy,kz,l)*dz1)/dz
     178           uuhn(ix,jy,kz,l)*dz1)/dz
    201179           vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
    202                 vvhn(ix,jy,kz,l)*dz1)/dz
     180           vvhn(ix,jy,kz,l)*dz1)/dz
    203181           ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
    204                 tthn(ix,jy,kz,n,l)*dz1)/dz
     182           tthn(ix,jy,kz,n,l)*dz1)/dz
    205183           qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
    206                 qvhn(ix,jy,kz,n,l)*dz1)/dz
     184           qvhn(ix,jy,kz,n,l)*dz1)/dz
    207185           pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
    208                 pvhn(ix,jy,kz,l)*dz1)/dz
     186           pvhn(ix,jy,kz,l)*dz1)/dz
    209187           rhon(ix,jy,iz,n,l)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
    210188           kmin=kz
     
    225203        do kz=kmin,nwz
    226204          if ((height(iz).gt.wzlev(kz-1)).and. &
    227                (height(iz).le.wzlev(kz))) then
    228            dz1=height(iz)-wzlev(kz-1)
    229            dz2=wzlev(kz)-height(iz)
    230            dz=dz1+dz2
    231            wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
    232                 +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
    233            kmin=kz
    234            goto 40
     205          (height(iz).le.wzlev(kz))) then
     206            dz1=height(iz)-wzlev(kz-1)
     207            dz2=wzlev(kz)-height(iz)
     208            dz=dz1+dz2
     209            wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)*pinmconv(kz-1)*dz2 &
     210            +wwhn(ix,jy,kz,l)*pinmconv(kz)*dz1)/dz
     211            kmin=kz
     212            goto 40
    235213          endif
    236214        end do
     
    242220
    243221      drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
    244            (height(2)-height(1))
     222      (height(2)-height(1))
    245223      do kz=2,nz-1
    246224        drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
    247              rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
     225        rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
    248226      end do
    249227      drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
     
    259237
    260238  do jy=1,nyn(l)-2
     239    cosf=cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
    261240    do ix=1,nxn(l)-2
    262241
     
    264243      do iz=2,nz-1
    265244
    266         ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/ &
    267              cos((real(jy)*dyn(l)+ylat0n(l))*pi180)
     245        ui=uun(ix,jy,iz,n,l)*dxconst*xresoln(l)/cosf
    268246        vi=vvn(ix,jy,iz,n,l)*dyconst*yresoln(l)
    269247
     
    319297               rain_cloud_above=1
    320298               cloudsnh(ix,jy,n,l)=cloudsnh(ix,jy,n,l)+ &
    321                     height(kz)-height(kz-1)
     299               height(kz)-height(kz-1)
    322300               if (lsp.ge.convp) then
    323301                  cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG