Changeset 4fbe7a5 in flexpart.git for src/calcpv.f90


Ignore:
Timestamp:
May 23, 2014, 11:48:41 AM (7 years ago)
Author:
Ignacio Pisso <ip@…>
Branches:
master, 10.4.1_pesei, FPv9.3.1, FPv9.3.1b_testing, FPv9.3.2, NetCDF, deposition, dev, fp9.3.1-20161214-nc4, grib2nc4_repair, inputlist, laptop, release-10, release-10.4.1, svn-petra, svn-trunk, univie
Children:
0aded10
Parents:
f13406c
Message:

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

git-svn-id: http://flexpart.flexpart.eu:8088/svn/FlexPart90/trunk@24 ef8cc7e1-21b7-489e-abab-c1baa636049d

File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/calcpv.f90

    re200b7a r4fbe7a5  
    2121
    2222subroutine calcpv(n,uuh,vvh,pvh)
    23   !                  i  i   i   o
     23  !               i  i   i   o
    2424  !*****************************************************************************
    2525  !                                                                            *
     
    4949  integer :: nlck
    5050  real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
    51   real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
    52   real :: pvavr,ppml(nuvzmax)
     51  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt
     52  real :: pvavr,ppml(0:nxmax-1,0:nymax-1,nuvzmax),ppmk(0:nxmax-1,0:nymax-1,nuvzmax)
    5353  real :: thup,thdn
    5454  real,parameter :: eps=1.e-5, p0=101325
     
    6262  ! Loop over entire grid
    6363  !**********************
     64  do kl=1,nuvz
     65    do jy=0,nymin1
     66      do ix=0,nxmin1
     67         ppml(ix,jy,kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
     68      enddo
     69    enddo
     70  enddo
     71  ppmk=(100000./ppml)**kappa
     72
    6473  do jy=0,nymin1
    6574    if (sglobal.and.jy.eq.0) goto 10
     
    106115      end if
    107116      jux=jumpx
    108   ! Precalculate pressure values for efficiency
    109       do kl=1,nuvz
    110         ppml(kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
    111       end do
    112117  !
    113118  ! Loop over the vertical
     
    115120
    116121      do kl=1,nuvz
    117         ppmk=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
    118         theta=tth(ix,jy,kl,n)*(100000./ppmk)**kappa
     122        theta=tth(ix,jy,kl,n)*ppmk(ix,jy,kl)
    119123        klvrp=kl+1
    120124        klvrm=kl-1
     
    125129        if (klvrp.gt.nuvz) klvrp=nuvz
    126130        if (klvrm.lt.1) klvrm=1
    127         ppmk=akz(klvrp)+bkz(klvrp)*ps(ix,jy,1,n)
    128         thetap=tth(ix,jy,klvrp,n)*(100000./ppmk)**kappa
    129         ppmk=akz(klvrm)+bkz(klvrm)*ps(ix,jy,1,n)
    130         thetam=tth(ix,jy,klvrm,n)*(100000./ppmk)**kappa
    131         dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))
     131        thetap=tth(ix,jy,klvrp,n)*ppmk(ix,jy,klvrp)
     132        thetam=tth(ix,jy,klvrm,n)*ppmk(ix,jy,klvrm)
     133        dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm))
    132134
    133135  ! Compute vertical position at pot. temperature surface on subgrid
     
    156158          kch=kch+1
    157159          k=kup
    158           ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n)
    159           thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa
    160           ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n)
    161           thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa
    162 
    163 
    164       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    165            ((thdn.le.theta).and.(thup.ge.theta))) then
    166               dt1=abs(theta-thdn)
    167               dt2=abs(theta-thup)
    168               dt=dt1+dt2
    169               if (dt.lt.eps) then   ! Avoid division by zero error
    170                 dt1=0.5             ! G.W., 10.4.1996
    171                 dt2=0.5
    172                 dt=1.0
    173               endif
    174           vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
    175               goto 20
     160          thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k)
     161          thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1)
     162
     163
     164          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
     165          ((thdn.le.theta).and.(thup.ge.theta))) then
     166            dt1=abs(theta-thdn)
     167            dt2=abs(theta-thup)
     168            dt=dt1+dt2
     169            if (dt.lt.eps) then   ! Avoid division by zero error
     170              dt1=0.5             ! G.W., 10.4.1996
     171              dt2=0.5
     172              dt=1.0
    176173            endif
     174            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
     175            goto 20
     176          endif
    17717741        continue
    178178  ! Downward branch
     
    181181          kch=kch+1
    182182          k=kdn
    183           ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n)
    184           thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa
    185           ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n)
    186           thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa
    187 
    188       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    189            ((thdn.le.theta).and.(thup.ge.theta))) then
    190               dt1=abs(theta-thdn)
    191               dt2=abs(theta-thup)
    192               dt=dt1+dt2
    193               if (dt.lt.eps) then   ! Avoid division by zero error
    194                 dt1=0.5             ! G.W., 10.4.1996
    195                 dt2=0.5
    196                 dt=1.0
    197               endif
    198           vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
    199               goto 20
     183          thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k)
     184          thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1)
     185
     186          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
     187          ((thdn.le.theta).and.(thup.ge.theta))) then
     188            dt1=abs(theta-thdn)
     189            dt2=abs(theta-thup)
     190            dt=dt1+dt2
     191            if (dt.lt.eps) then   ! Avoid division by zero error
     192              dt1=0.5             ! G.W., 10.4.1996
     193              dt2=0.5
     194              dt=1.0
    200195            endif
    201             goto 40
     196            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
     197            goto 20
     198          endif
     199          goto 40
    202200  ! This section used when no values were found
    20320121      continue
     
    236234          kch=kch+1
    237235          k=kup
    238           ppmk=akz(k)+bkz(k)*ps(ix,j,1,n)
    239           thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa
    240           ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n)
    241           thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa
    242       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    243            ((thdn.le.theta).and.(thup.ge.theta))) then
    244               dt1=abs(theta-thdn)
    245               dt2=abs(theta-thup)
    246               dt=dt1+dt2
    247               if (dt.lt.eps) then   ! Avoid division by zero error
    248                 dt1=0.5             ! G.W., 10.4.1996
    249                 dt2=0.5
    250                 dt=1.0
    251               endif
    252               uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
    253               goto 50
     236          thdn=tth(ix,j,k,n)*ppmk(ix,j,k)
     237          thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1)
     238          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
     239          ((thdn.le.theta).and.(thup.ge.theta))) then
     240            dt1=abs(theta-thdn)
     241            dt2=abs(theta-thup)
     242            dt=dt1+dt2
     243            if (dt.lt.eps) then   ! Avoid division by zero error
     244              dt1=0.5             ! G.W., 10.4.1996
     245              dt2=0.5
     246              dt=1.0
    254247            endif
     248            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
     249            goto 50
     250          endif
    25525171        continue
    256252  ! Downward branch
     
    259255          kch=kch+1
    260256          k=kdn
    261           ppmk=akz(k)+bkz(k)*ps(ix,j,1,n)
    262           thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa
    263           ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n)
    264           thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa
    265       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
    266            ((thdn.le.theta).and.(thup.ge.theta))) then
    267               dt1=abs(theta-thdn)
    268               dt2=abs(theta-thup)
    269               dt=dt1+dt2
    270               if (dt.lt.eps) then   ! Avoid division by zero error
    271                 dt1=0.5             ! G.W., 10.4.1996
    272                 dt2=0.5
    273                 dt=1.0
    274               endif
    275               uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
    276               goto 50
     257          thdn=tth(ix,j,k,n)*ppmk(ix,j,k)
     258          thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1)
     259          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
     260          ((thdn.le.theta).and.(thup.ge.theta))) then
     261            dt1=abs(theta-thdn)
     262            dt2=abs(theta-thup)
     263            dt=dt1+dt2
     264            if (dt.lt.eps) then   ! Avoid division by zero error
     265              dt1=0.5             ! G.W., 10.4.1996
     266              dt2=0.5
     267              dt=1.0
    277268            endif
    278             goto 70
     269            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
     270            goto 50
     271          endif
     272          goto 70
    279273  ! This section used when no values were found
    28027451      continue
    281275  ! Must use uu at current level and lat. juy becomes smaller by 1
    282         uy(jj)=uuh(ix,jy,kl)
    283         juy=juy-1
     276          uy(jj)=uuh(ix,jy,kl)
     277          juy=juy-1
    284278  ! Otherwise OK
    28527950        continue
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG