calcpv_nests.f90 9.98 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
!                                                                     *
! This file is part of FLEXPART.                                      *
!                                                                     *
! FLEXPART is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or   *
! (at your option) any later version.                                 *
!                                                                     *
! FLEXPART is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
! GNU General Public License for more details.                        *
!                                                                     *
! You should have received a copy of the GNU General Public License   *
! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
!**********************************************************************

subroutine calcpv_nests(l,n,uuhn,vvhn,pvhn)
23
  !                     i i  i    i    o
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
  !*****************************************************************************
  !                                                                            *
  !  Calculation of potential vorticity on 3-d nested grid                     *
  !                                                                            *
  !     Author: P. James                                                       *
  !     22 February 2000                                                       *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! n                  temporal index for meteorological fields (1 to 2)       *
  ! l                  index of current nest                                   *
  !                                                                            *
  ! Constants:                                                                 *
  !                                                                            *
  !*****************************************************************************

  use par_mod
  use com_mod

  implicit none

  integer :: n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch
  integer :: jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr
  integer :: nlck,l
  real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
50 51
  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt
  real :: ppml(0:nxmaxn-1,0:nymaxn-1,nuvzmax),ppmk(0:nxmaxn-1,0:nymaxn-1,nuvzmax)
52 53 54 55 56 57 58 59 60 61 62
  real :: thup,thdn
  real,parameter :: eps=1.e-5,p0=101325
  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)

  ! Set number of levels to check for adjacent theta
  nlck=nuvz/3
  !
  ! Loop over entire grid
  !**********************
63 64 65 66 67 68 69
  do kl=1,nuvz
    do jy=0,nyn(l)-1
      do ix=0,nxn(l)-1
         ppml(ix,jy,kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
      enddo
    enddo
  enddo
70 71
!  ppmk=(100000./ppml)**kappa
  ppmk(0:nxn(l)-1,0:nyn(l)-1,1:nuvz)=(100000./ppml(0:nxn(l)-1,0:nyn(l)-1,1:nuvz))**kappa
72

73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104
  do jy=0,nyn(l)-1
    phi = (ylat0n(l) + jy * dyn(l)) * pi / 180.
    f = 0.00014585 * sin(phi)
    tanphi = tan(phi)
    cosphi = cos(phi)
  ! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat)
      jyvp=jy+1
      jyvm=jy-1
      if (jy.eq.0) jyvm=0
      if (jy.eq.nyn(l)-1) jyvp=nyn(l)-1
  ! Define absolute gap length
      jumpy=2
      if (jy.eq.0.or.jy.eq.nyn(l)-1) jumpy=1
      juy=jumpy
  !
    do ix=0,nxn(l)-1
  ! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long)
      ixvp=ix+1
      ixvm=ix-1
      jumpx=2
      if (ix.eq.0) ixvm=0
      if (ix.eq.nxn(l)-1) ixvp=nxn(l)-1
      ivrp=ixvp
      ivrm=ixvm
  ! Define absolute gap length
      if (ix.eq.0.or.ix.eq.nxn(l)-1) jumpx=1
      jux=jumpx
  !
  ! Loop over the vertical
  !***********************

      do kl=1,nuvz
105
        theta=tthn(ix,jy,kl,n,l)*ppmk(ix,jy,kl)
106 107 108 109 110 111 112 113
        klvrp=kl+1
        klvrm=kl-1
        klpt=kl
  ! If top or bottom level, dthetadp is evaluated between the current
  ! level and the level inside, otherwise between level+1 and level-1
  !
        if (klvrp.gt.nuvz) klvrp=nuvz
        if (klvrm.lt.1) klvrm=1
114 115 116
        thetap=tthn(ix,jy,klvrp,n,l)*ppmk(ix,jy,klvrp)
        thetam=tthn(ix,jy,klvrm,n,l)*ppmk(ix,jy,klvrm)
        dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm))
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138

  ! Compute vertical position at pot. temperature surface on subgrid
  ! and the wind at that position
  !*****************************************************************
  ! a) in x direction
        ii=0
        do i=ixvm,ixvp,jumpx
          ivr=i
          ii=ii+1
  ! Search adjacent levels for current theta value
  ! Spiral out from current level for efficiency
          kup=klpt-1
          kdn=klpt
          kch=0
40        continue
  ! Upward branch
          kup=kup+1
          if (kch.ge.nlck) goto 21     ! No more levels to check,
  !                                       ! and no values found
          if (kup.ge.nuvz) goto 41
          kch=kch+1
          k=kup
139 140
          thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k)
          thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1)
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160

      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
           ((thdn.le.theta).and.(thup.ge.theta))) then
              dt1=abs(theta-thdn)
              dt2=abs(theta-thup)
              dt=dt1+dt2
              if (dt.lt.eps) then   ! Avoid division by zero error
                dt1=0.5             ! G.W., 10.4.1996
                dt2=0.5
                dt=1.0
              endif
    vx(ii)=(vvhn(ivr,jy,k,l)*dt2+vvhn(ivr,jy,k+1,l)*dt1)/dt
              goto 20
            endif
41        continue
  ! Downward branch
          kdn=kdn-1
          if (kdn.lt.1) goto 40
          kch=kch+1
          k=kdn
161 162
          thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k)
          thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1)
163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212
      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
           ((thdn.le.theta).and.(thup.ge.theta))) then
              dt1=abs(theta-thdn)
              dt2=abs(theta-thup)
              dt=dt1+dt2
              if (dt.lt.eps) then   ! Avoid division by zero error
                dt1=0.5             ! G.W., 10.4.1996
                dt2=0.5
                dt=1.0
              endif
    vx(ii)=(vvhn(ivr,jy,k,l)*dt2+vvhn(ivr,jy,k+1,l)*dt1)/dt
              goto 20
            endif
            goto 40
  ! This section used when no values were found
21      continue
  ! Must use vv at current level and long. jux becomes smaller by 1
        vx(ii)=vvhn(ix,jy,kl,l)
        jux=jux-1
  ! Otherwise OK
20        continue
        end do
      if (jux.gt.0) then
      dvdx=(vx(2)-vx(1))/real(jux)/(dxn(l)*pi/180.)
      else
      dvdx=vvhn(ivrp,jy,kl,l)-vvhn(ivrm,jy,kl,l)
      dvdx=dvdx/real(jumpx)/(dxn(l)*pi/180.)
  ! Only happens if no equivalent theta value
  ! can be found on either side, hence must use values
  ! from either side, same pressure level.
      end if

  ! b) in y direction

        jj=0
        do j=jyvm,jyvp,jumpy
          jj=jj+1
  ! Search adjacent levels for current theta value
  ! Spiral out from current level for efficiency
          kup=klpt-1
          kdn=klpt
          kch=0
70        continue
  ! Upward branch
          kup=kup+1
          if (kch.ge.nlck) goto 51     ! No more levels to check,
  !                                     ! and no values found
          if (kup.ge.nuvz) goto 71
          kch=kch+1
          k=kup
213 214
          thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k)
          thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1)
215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233
      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
           ((thdn.le.theta).and.(thup.ge.theta))) then
              dt1=abs(theta-thdn)
              dt2=abs(theta-thup)
              dt=dt1+dt2
              if (dt.lt.eps) then   ! Avoid division by zero error
                dt1=0.5             ! G.W., 10.4.1996
                dt2=0.5
                dt=1.0
              endif
        uy(jj)=(uuhn(ix,j,k,l)*dt2+uuhn(ix,j,k+1,l)*dt1)/dt
              goto 50
            endif
71        continue
  ! Downward branch
          kdn=kdn-1
          if (kdn.lt.1) goto 70
          kch=kch+1
          k=kdn
234 235
          thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k)
          thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1)
236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276
      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
           ((thdn.le.theta).and.(thup.ge.theta))) then
              dt1=abs(theta-thdn)
              dt2=abs(theta-thup)
              dt=dt1+dt2
              if (dt.lt.eps) then   ! Avoid division by zero error
                dt1=0.5             ! G.W., 10.4.1996
                dt2=0.5
                dt=1.0
              endif
        uy(jj)=(uuhn(ix,j,k,l)*dt2+uuhn(ix,j,k+1,l)*dt1)/dt
              goto 50
            endif
            goto 70
  ! This section used when no values were found
51      continue
  ! Must use uu at current level and lat. juy becomes smaller by 1
        uy(jj)=uuhn(ix,jy,kl,l)
        juy=juy-1
  ! Otherwise OK
50        continue
        end do
      if (juy.gt.0) then
      dudy=(uy(2)-uy(1))/real(juy)/(dyn(l)*pi/180.)
      else
      dudy=uuhn(ix,jyvp,kl,l)-uuhn(ix,jyvm,kl,l)
      dudy=dudy/real(jumpy)/(dyn(l)*pi/180.)
      end if
  !
      pvhn(ix,jy,kl,l)=dthetadp*(f+(dvdx/cosphi-dudy &
           +uuhn(ix,jy,kl,l)*tanphi)/r_earth)*(-1.e6)*9.81

  !
  ! Resest jux and juy
      jux=jumpx
      juy=jumpy
      end do
    end do
  end do
  !
end subroutine calcpv_nests