calcpv_nests.f90

Go to the documentation of this file.

!**********************************************************************
! 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)
  !                        i i  i    i    o
  !*****************************************************************************
  !                                                                            *
  !  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
  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
  real :: ppml(nuvzmax)
  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
  !**********************
  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
  ! Precalculate pressure values for efficiency
      do kl=1,nuvz
        ppml(kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
      end do
  !
  ! Loop over the vertical
  !***********************

      do kl=1,nuvz
        ppmk=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
        theta=tthn(ix,jy,kl,n,l)*(100000./ppmk)**kappa
        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
        ppmk=akz(klvrp)+bkz(klvrp)*psn(ix,jy,1,n,l)
        thetap=tthn(ix,jy,klvrp,n,l)*(100000./ppmk)**kappa
        ppmk=akz(klvrm)+bkz(klvrm)*psn(ix,jy,1,n,l)
        thetam=tthn(ix,jy,klvrm,n,l)*(100000./ppmk)**kappa
        dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))

  ! 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
          ppmk=akz(k)+bkz(k)*psn(ivr,jy,1,n,l)
          thdn=tthn(ivr,jy,k,n,l)*(100000./ppmk)**kappa
          ppmk=akz(k+1)+bkz(k+1)*psn(ivr,jy,1,n,l)
          thup=tthn(ivr,jy,k+1,n,l)*(100000./ppmk)**kappa

      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
          ppmk=akz(k)+bkz(k)*psn(ivr,jy,1,n,l)
          thdn=tthn(ivr,jy,k,n,l)*(100000./ppmk)**kappa
          ppmk=akz(k+1)+bkz(k+1)*psn(ivr,jy,1,n,l)
          thup=tthn(ivr,jy,k+1,n,l)*(100000./ppmk)**kappa
      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
          ppmk=akz(k)+bkz(k)*psn(ix,j,1,n,l)
          thdn=tthn(ix,j,k,n,l)*(100000./ppmk)**kappa
          ppmk=akz(k+1)+bkz(k+1)*psn(ix,j,1,n,l)
          thup=tthn(ix,j,k+1,n,l)*(100000./ppmk)**kappa
      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
          ppmk=akz(k)+bkz(k)*psn(ix,j,1,n,l)
          thdn=tthn(ix,j,k,n,l)*(100000./ppmk)**kappa
          ppmk=akz(k+1)+bkz(k+1)*psn(ix,j,1,n,l)
          thup=tthn(ix,j,k+1,n,l)*(100000./ppmk)**kappa
      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