source: trunk/src/calcpv.f90 @ 28

Last change on this file since 28 was 24, checked in by igpis, 10 years ago

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

File size: 11.0 KB
Line 
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)
23  !               i  i   i   o
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
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)
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  !**********************
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
73  do jy=0,nymin1
74    if (sglobal.and.jy.eq.0) goto 10
75    if (nglobal.and.jy.eq.nymin1) goto 10
76    phi = (ylat0 + jy * dy) * pi / 180.
77    f = 0.00014585 * sin(phi)
78    tanphi = tan(phi)
79    cosphi = cos(phi)
80  ! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat)
81      jyvp=jy+1
82      jyvm=jy-1
83      if (jy.eq.0) jyvm=0
84      if (jy.eq.nymin1) jyvp=nymin1
85  ! Define absolute gap length
86      jumpy=2
87      if (jy.eq.0.or.jy.eq.nymin1) jumpy=1
88      if (sglobal.and.jy.eq.1) then
89         jyvm=1
90         jumpy=1
91      end if
92      if (nglobal.and.jy.eq.ny-2) then
93         jyvp=ny-2
94         jumpy=1
95      end if
96      juy=jumpy
97  !
98    do ix=0,nxmin1
99  ! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long)
100      ixvp=ix+1
101      ixvm=ix-1
102      jumpx=2
103      if (xglobal) then
104         ivrp=ixvp
105         ivrm=ixvm
106         if (ixvm.lt.0) ivrm=ixvm+nxmin1
107         if (ixvp.ge.nx) ivrp=ixvp-nx+1
108      else
109        if (ix.eq.0) ixvm=0
110        if (ix.eq.nxmin1) ixvp=nxmin1
111        ivrp=ixvp
112        ivrm=ixvm
113  ! Define absolute gap length
114        if (ix.eq.0.or.ix.eq.nxmin1) jumpx=1
115      end if
116      jux=jumpx
117  !
118  ! Loop over the vertical
119  !***********************
120
121      do kl=1,nuvz
122        theta=tth(ix,jy,kl,n)*ppmk(ix,jy,kl)
123        klvrp=kl+1
124        klvrm=kl-1
125        klpt=kl
126  ! If top or bottom level, dthetadp is evaluated between the current
127  ! level and the level inside, otherwise between level+1 and level-1
128  !
129        if (klvrp.gt.nuvz) klvrp=nuvz
130        if (klvrm.lt.1) klvrm=1
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))
134
135  ! Compute vertical position at pot. temperature surface on subgrid
136  ! and the wind at that position
137  !*****************************************************************
138  ! a) in x direction
139        ii=0
140        do i=ixvm,ixvp,jumpx
141          ivr=i
142          if (xglobal) then
143             if (i.lt.0) ivr=ivr+nxmin1
144             if (i.ge.nx) ivr=ivr-nx+1
145          end if
146          ii=ii+1
147  ! Search adjacent levels for current theta value
148  ! Spiral out from current level for efficiency
149          kup=klpt-1
150          kdn=klpt
151          kch=0
15240        continue
153  ! Upward branch
154          kup=kup+1
155          if (kch.ge.nlck) goto 21     ! No more levels to check,
156  !                                       ! and no values found
157          if (kup.ge.nuvz) goto 41
158          kch=kch+1
159          k=kup
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
173            endif
174            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
175            goto 20
176          endif
17741        continue
178  ! Downward branch
179          kdn=kdn-1
180          if (kdn.lt.1) goto 40
181          kch=kch+1
182          k=kdn
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
195            endif
196            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
197            goto 20
198          endif
199          goto 40
200  ! This section used when no values were found
20121      continue
202  ! Must use vv at current level and long. jux becomes smaller by 1
203        vx(ii)=vvh(ix,jy,kl)
204        jux=jux-1
205  ! Otherwise OK
20620        continue
207        end do
208      if (jux.gt.0) then
209      dvdx=(vx(2)-vx(1))/real(jux)/(dx*pi/180.)
210      else
211      dvdx=vvh(ivrp,jy,kl)-vvh(ivrm,jy,kl)
212      dvdx=dvdx/real(jumpx)/(dx*pi/180.)
213  ! Only happens if no equivalent theta value
214  ! can be found on either side, hence must use values
215  ! from either side, same pressure level.
216      end if
217
218  ! b) in y direction
219
220        jj=0
221        do j=jyvm,jyvp,jumpy
222          jj=jj+1
223  ! Search adjacent levels for current theta value
224  ! Spiral out from current level for efficiency
225          kup=klpt-1
226          kdn=klpt
227          kch=0
22870        continue
229  ! Upward branch
230          kup=kup+1
231          if (kch.ge.nlck) goto 51     ! No more levels to check,
232  !                                     ! and no values found
233          if (kup.ge.nuvz) goto 71
234          kch=kch+1
235          k=kup
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
247            endif
248            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
249            goto 50
250          endif
25171        continue
252  ! Downward branch
253          kdn=kdn-1
254          if (kdn.lt.1) goto 70
255          kch=kch+1
256          k=kdn
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
268            endif
269            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
270            goto 50
271          endif
272          goto 70
273  ! This section used when no values were found
27451      continue
275  ! Must use uu at current level and lat. juy becomes smaller by 1
276          uy(jj)=uuh(ix,jy,kl)
277          juy=juy-1
278  ! Otherwise OK
27950        continue
280        end do
281      if (juy.gt.0) then
282      dudy=(uy(2)-uy(1))/real(juy)/(dy*pi/180.)
283      else
284      dudy=uuh(ix,jyvp,kl)-uuh(ix,jyvm,kl)
285      dudy=dudy/real(jumpy)/(dy*pi/180.)
286      end if
287  !
288      pvh(ix,jy,kl)=dthetadp*(f+(dvdx/cosphi-dudy &
289           +uuh(ix,jy,kl)*tanphi)/r_earth)*(-1.e6)*9.81
290
291
292  !
293  ! Resest jux and juy
294      jux=jumpx
295      juy=jumpy
296      end do
297    end do
29810  continue
299  end do
300  !
301  ! Fill in missing PV values on poles, if present
302  ! Use mean PV of surrounding latitude ring
303  !
304      if (sglobal) then
305         do kl=1,nuvz
306            pvavr=0.
307            do ix=0,nxmin1
308               pvavr=pvavr+pvh(ix,1,kl)
309            end do
310            pvavr=pvavr/real(nx)
311            jy=0
312            do ix=0,nxmin1
313               pvh(ix,jy,kl)=pvavr
314            end do
315         end do
316      end if
317      if (nglobal) then
318         do kl=1,nuvz
319            pvavr=0.
320            do ix=0,nxmin1
321               pvavr=pvavr+pvh(ix,ny-2,kl)
322            end do
323            pvavr=pvavr/real(nx)
324            jy=nymin1
325            do ix=0,nxmin1
326               pvh(ix,jy,kl)=pvavr
327            end do
328         end do
329      end if
330
331end subroutine calcpv
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG