source: trunk/src/calcpv.f90 @ 28

Last change on this file since 28 was 24, checked in by igpis, 8 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