source: flexpart.git/src/calcpv.f90 @ 02095e3

devrelease-10univie
Last change on this file since 02095e3 was db712a8, checked in by Espen Sollum ATMOS <eso@…>, 4 years ago

Completed handling of nested wind fields with cloud water (for wet deposition).

  • Property mode set to 100644
File size: 11.1 KB
RevLine 
[e200b7a]1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine calcpv(n,uuh,vvh,pvh)
[4fbe7a5]23  !               i  i   i   o
[e200b7a]24  !*****************************************************************************
25  !                                                                            *
26  !  Calculation of potential vorticity on 3-d grid.                           *
27  !                                                                            *
28  !     Author: P. James                                                       *
29  !     3 February 2000                                                        *
30  !                                                                            *
31  !     Adaptation to FLEXPART, A. Stohl, 1 May 2000                           *
32  !                                                                            *
33  !*****************************************************************************
34  !                                                                            *
35  ! Variables:                                                                 *
36  ! n                  temporal index for meteorological fields (1 to 2)       *
37  !                                                                            *
38  ! Constants:                                                                 *
39  !                                                                            *
40  !*****************************************************************************
41
42  use par_mod
43  use com_mod
44
45  implicit none
46
47  integer :: n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch
48  integer :: jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr
49  integer :: nlck
50  real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
[4fbe7a5]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)
[e200b7a]53  real :: thup,thdn
54  real,parameter :: eps=1.e-5, p0=101325
55  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
56  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
57  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
58
59  ! Set number of levels to check for adjacent theta
60  nlck=nuvz/3
61  !
62  ! Loop over entire grid
63  !**********************
[4fbe7a5]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
[8a65cb0]71
72!  ppmk(:,:,1:nuvz)=(100000./ppml(:,:,1:nuvz))**kappa
[db712a8]73  ppmk(0:nxmin1,0:nymin1,1:nuvz)=(100000./ppml(0:nxmin1,0:nymin1,1:nuvz))**kappa
[4fbe7a5]74
[e200b7a]75  do jy=0,nymin1
76    if (sglobal.and.jy.eq.0) goto 10
77    if (nglobal.and.jy.eq.nymin1) goto 10
78    phi = (ylat0 + jy * dy) * pi / 180.
79    f = 0.00014585 * sin(phi)
80    tanphi = tan(phi)
81    cosphi = cos(phi)
82  ! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat)
83      jyvp=jy+1
84      jyvm=jy-1
85      if (jy.eq.0) jyvm=0
86      if (jy.eq.nymin1) jyvp=nymin1
87  ! Define absolute gap length
88      jumpy=2
89      if (jy.eq.0.or.jy.eq.nymin1) jumpy=1
90      if (sglobal.and.jy.eq.1) then
91         jyvm=1
92         jumpy=1
93      end if
94      if (nglobal.and.jy.eq.ny-2) then
95         jyvp=ny-2
96         jumpy=1
97      end if
98      juy=jumpy
99  !
100    do ix=0,nxmin1
101  ! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long)
102      ixvp=ix+1
103      ixvm=ix-1
104      jumpx=2
105      if (xglobal) then
106         ivrp=ixvp
107         ivrm=ixvm
108         if (ixvm.lt.0) ivrm=ixvm+nxmin1
109         if (ixvp.ge.nx) ivrp=ixvp-nx+1
110      else
111        if (ix.eq.0) ixvm=0
112        if (ix.eq.nxmin1) ixvp=nxmin1
113        ivrp=ixvp
114        ivrm=ixvm
115  ! Define absolute gap length
116        if (ix.eq.0.or.ix.eq.nxmin1) jumpx=1
117      end if
118      jux=jumpx
119  !
120  ! Loop over the vertical
121  !***********************
122
123      do kl=1,nuvz
[4fbe7a5]124        theta=tth(ix,jy,kl,n)*ppmk(ix,jy,kl)
[e200b7a]125        klvrp=kl+1
126        klvrm=kl-1
127        klpt=kl
128  ! If top or bottom level, dthetadp is evaluated between the current
129  ! level and the level inside, otherwise between level+1 and level-1
130  !
131        if (klvrp.gt.nuvz) klvrp=nuvz
132        if (klvrm.lt.1) klvrm=1
[4fbe7a5]133        thetap=tth(ix,jy,klvrp,n)*ppmk(ix,jy,klvrp)
134        thetam=tth(ix,jy,klvrm,n)*ppmk(ix,jy,klvrm)
135        dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm))
[e200b7a]136
137  ! Compute vertical position at pot. temperature surface on subgrid
138  ! and the wind at that position
139  !*****************************************************************
140  ! a) in x direction
141        ii=0
142        do i=ixvm,ixvp,jumpx
143          ivr=i
144          if (xglobal) then
145             if (i.lt.0) ivr=ivr+nxmin1
146             if (i.ge.nx) ivr=ivr-nx+1
147          end if
148          ii=ii+1
149  ! Search adjacent levels for current theta value
150  ! Spiral out from current level for efficiency
151          kup=klpt-1
152          kdn=klpt
153          kch=0
15440        continue
155  ! Upward branch
156          kup=kup+1
157          if (kch.ge.nlck) goto 21     ! No more levels to check,
158  !                                       ! and no values found
159          if (kup.ge.nuvz) goto 41
160          kch=kch+1
161          k=kup
[4fbe7a5]162          thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k)
163          thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1)
[e200b7a]164
165
[4fbe7a5]166          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
167          ((thdn.le.theta).and.(thup.ge.theta))) then
168            dt1=abs(theta-thdn)
169            dt2=abs(theta-thup)
170            dt=dt1+dt2
171            if (dt.lt.eps) then   ! Avoid division by zero error
172              dt1=0.5             ! G.W., 10.4.1996
173              dt2=0.5
174              dt=1.0
[e200b7a]175            endif
[4fbe7a5]176            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
177            goto 20
178          endif
[e200b7a]17941        continue
180  ! Downward branch
181          kdn=kdn-1
182          if (kdn.lt.1) goto 40
183          kch=kch+1
184          k=kdn
[4fbe7a5]185          thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k)
186          thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1)
[e200b7a]187
[4fbe7a5]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
[e200b7a]197            endif
[4fbe7a5]198            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
199            goto 20
200          endif
201          goto 40
[e200b7a]202  ! This section used when no values were found
20321      continue
204  ! Must use vv at current level and long. jux becomes smaller by 1
205        vx(ii)=vvh(ix,jy,kl)
206        jux=jux-1
207  ! Otherwise OK
20820        continue
209        end do
210      if (jux.gt.0) then
211      dvdx=(vx(2)-vx(1))/real(jux)/(dx*pi/180.)
212      else
213      dvdx=vvh(ivrp,jy,kl)-vvh(ivrm,jy,kl)
214      dvdx=dvdx/real(jumpx)/(dx*pi/180.)
215  ! Only happens if no equivalent theta value
216  ! can be found on either side, hence must use values
217  ! from either side, same pressure level.
218      end if
219
220  ! b) in y direction
221
222        jj=0
223        do j=jyvm,jyvp,jumpy
224          jj=jj+1
225  ! Search adjacent levels for current theta value
226  ! Spiral out from current level for efficiency
227          kup=klpt-1
228          kdn=klpt
229          kch=0
23070        continue
231  ! Upward branch
232          kup=kup+1
233          if (kch.ge.nlck) goto 51     ! No more levels to check,
234  !                                     ! and no values found
235          if (kup.ge.nuvz) goto 71
236          kch=kch+1
237          k=kup
[4fbe7a5]238          thdn=tth(ix,j,k,n)*ppmk(ix,j,k)
239          thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1)
240          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
241          ((thdn.le.theta).and.(thup.ge.theta))) then
242            dt1=abs(theta-thdn)
243            dt2=abs(theta-thup)
244            dt=dt1+dt2
245            if (dt.lt.eps) then   ! Avoid division by zero error
246              dt1=0.5             ! G.W., 10.4.1996
247              dt2=0.5
248              dt=1.0
[e200b7a]249            endif
[4fbe7a5]250            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
251            goto 50
252          endif
[e200b7a]25371        continue
254  ! Downward branch
255          kdn=kdn-1
256          if (kdn.lt.1) goto 70
257          kch=kch+1
258          k=kdn
[4fbe7a5]259          thdn=tth(ix,j,k,n)*ppmk(ix,j,k)
260          thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1)
261          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
262          ((thdn.le.theta).and.(thup.ge.theta))) then
263            dt1=abs(theta-thdn)
264            dt2=abs(theta-thup)
265            dt=dt1+dt2
266            if (dt.lt.eps) then   ! Avoid division by zero error
267              dt1=0.5             ! G.W., 10.4.1996
268              dt2=0.5
269              dt=1.0
[e200b7a]270            endif
[4fbe7a5]271            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
272            goto 50
273          endif
274          goto 70
[e200b7a]275  ! This section used when no values were found
27651      continue
277  ! Must use uu at current level and lat. juy becomes smaller by 1
[4fbe7a5]278          uy(jj)=uuh(ix,jy,kl)
279          juy=juy-1
[e200b7a]280  ! Otherwise OK
28150        continue
282        end do
283      if (juy.gt.0) then
284      dudy=(uy(2)-uy(1))/real(juy)/(dy*pi/180.)
285      else
286      dudy=uuh(ix,jyvp,kl)-uuh(ix,jyvm,kl)
287      dudy=dudy/real(jumpy)/(dy*pi/180.)
288      end if
289  !
290      pvh(ix,jy,kl)=dthetadp*(f+(dvdx/cosphi-dudy &
291           +uuh(ix,jy,kl)*tanphi)/r_earth)*(-1.e6)*9.81
292
293
294  !
295  ! Resest jux and juy
296      jux=jumpx
297      juy=jumpy
298      end do
299    end do
30010  continue
301  end do
302  !
303  ! Fill in missing PV values on poles, if present
304  ! Use mean PV of surrounding latitude ring
305  !
306      if (sglobal) then
307         do kl=1,nuvz
308            pvavr=0.
309            do ix=0,nxmin1
310               pvavr=pvavr+pvh(ix,1,kl)
311            end do
312            pvavr=pvavr/real(nx)
313            jy=0
314            do ix=0,nxmin1
315               pvh(ix,jy,kl)=pvavr
316            end do
317         end do
318      end if
319      if (nglobal) then
320         do kl=1,nuvz
321            pvavr=0.
322            do ix=0,nxmin1
323               pvavr=pvavr+pvh(ix,ny-2,kl)
324            end do
325            pvavr=pvavr/real(nx)
326            jy=nymin1
327            do ix=0,nxmin1
328               pvh(ix,jy,kl)=pvavr
329            end do
330         end do
331      end if
332
333end subroutine calcpv
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG