!*********************************************************************** !* Copyright 2012,2013 * !* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, * !* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,* !* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, * !* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, * !* * !* This file is part of FLEXPART WRF * !* * !* 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 . * !*********************************************************************** subroutine calcpv_nests(l,n,uuhn,vvhn,pvhn) ! i i i i o !******************************************************************************* ! * ! Note: This is the FLEXPART_WRF version of subroutine calcpv_nests. * ! The computational grid is the WRF x-y grid rather than lat-lon. * ! * ! Calculation of potential vorticity on 3-d nested grid * ! * ! Author: P. James * ! 22 February 2000 * ! * ! 11 Nov 2005, R. Easter - changes associated with WRF horizontal grid. * ! For pressure use pph instead of (akz + bkz*ps) * ! *** Note -- see ??? comments below regarding the pvh calculation. * ! * !******************************************************************************* ! * ! Variables: * ! n temporal index for meteorological fields (1 to 2) * ! l index of current nest * ! * ! Constants: * ! * !******************************************************************************* ! include 'includepar' ! include 'includecom' ! ! 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 altit(nuvzmax),ppml(nuvzmax) ! real thup,thdn,eps,p0 ! parameter(eps=1.e-5,p0=101325) ! real thhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) ! 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) 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(kind=4) :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) real(kind=4) :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) real :: dumlon,dumlat real :: thhn(0:nxmax-1,0:nymax-1,nuvzmax,maxnests) real :: altit(nuvzmax) ! Set number of levels to check for adjacent theta nlck=nuvz/3 ! FLEXPART_WRF -- altit is never used, so don't calculate it ! do 5 k=1,nuvz ! altit(k)=akz(k)/p0+bkz(k) !5 continue ! *** Precalculate all theta levels for efficiency do jy=0,nyn(l)-1 do kl=1,nuvz do ix=0,nxn(l)-1 ! FLEXPART_WRF -- use pph here ! ppmk=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l) ppmk=pphn(ix,jy,kl,n,l) thhn(ix,jy,kl,l)=tthn(ix,jy,kl,n,l)*(100000./ppmk)**kappa enddo enddo enddo ! ! Loop over entire grid !********************** do jy=0,nyn(l)-1 ! for FLEXPART_WRF, x & y coords are in meters ! and true latitude varies with both i and j ! 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 ! for FLEXPART_WRF, x & y coords are in meters ! and true latitude varies with both i and j call xyindex_to_ll_wrf( & l, real(ix), real(jy), dumlon, dumlat ) phi = dumlat * pi / 180. f = 0.00014585 * sin(phi) tanphi = tan(phi) cosphi = cos(phi) ! 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 ! FLEXPART_WRF -- use pph here ! ppml(kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l) ppml(kl)=pphn(ix,jy,kl,n,l) enddo ! ! Loop over the vertical !*********************** do kl=1,nuvz theta=thhn(ix,jy,kl,l) 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 thetap=thhn(ix,jy,klvrp,l) thetam=thhn(ix,jy,klvrm,l) 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 thdn=thhn(ivr,jy,k,l) thup=thhn(ivr,jy,k+1,l) 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 thdn=thhn(ivr,jy,k,l) thup=thhn(ivr,jy,k+1,l) 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 ! for FLEXPART_WRF, dx & dy are in meters. ! dvdx=(vx(2)-vx(1))/real(jux)/(dxn(l)*pi/180.) dvdx=(vx(2)-vx(1))/real(jux)/dxn(l) else dvdx=vvhn(ivrp,jy,kl,l)-vvhn(ivrm,jy,kl,l) ! dvdx=dvdx/real(jumpx)/(dxn(l)*pi/180.) dvdx=dvdx/real(jumpx)/dxn(l) ! 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 thdn=thhn(ix,j,k,l) thup=thhn(ix,j,k+1,l) 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 thdn=thhn(ix,j,k,l) thup=thhn(ix,j,k+1,l) 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 ! for FLEXPART_WRF, dx & dy are in meters. ! dudy=(uy(2)-uy(1))/real(juy)/(dyn(l)*pi/180.) dudy=(uy(2)-uy(1))/real(juy)/dyn(l) else dudy=uuhn(ix,jyvp,kl,l)-uuhn(ix,jyvm,kl,l) ! dudy=dudy/real(jumpy)/(dyn(l)*pi/180.) dudy=dudy/real(jumpy)/dyn(l) end if ! ! for FLEXPART_WRF, dx & dy are in meters. ! don't need to divide by r_earth when doing d/dy ! don't need to divide by r_earth*cosphi when doing d/dx ! ??? I don't understand the uuhn*tanphi term, but leave it in for now ??? ! ??? What is the "-1.e6" factor ??? ! ! pvhn(ix,jy,kl,l)=dthetadp*(f+(dvdx/cosphi-dudy ! + +uuhn(ix,jy,kl,l)*tanphi)/r_earth)*(-1.e6)*9.81 pvhn(ix,jy,kl,l)=dthetadp*( f + dvdx - 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