source: branches/ignacio/FLEXPART_9.1.8/src/calcpv.f90 @ 18

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

add verbose mode to version 9.1.7.1

File size: 11.5 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,ppmk
52  real :: pvavr,ppml(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 jy=0,nymin1
65    if (sglobal.and.jy.eq.0) goto 10
66    if (nglobal.and.jy.eq.nymin1) goto 10
67    phi = (ylat0 + jy * dy) * pi / 180.
68    f = 0.00014585 * sin(phi)
69    tanphi = tan(phi)
70    cosphi = cos(phi)
71  ! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat)
72      jyvp=jy+1
73      jyvm=jy-1
74      if (jy.eq.0) jyvm=0
75      if (jy.eq.nymin1) jyvp=nymin1
76  ! Define absolute gap length
77      jumpy=2
78      if (jy.eq.0.or.jy.eq.nymin1) jumpy=1
79      if (sglobal.and.jy.eq.1) then
80         jyvm=1
81         jumpy=1
82      end if
83      if (nglobal.and.jy.eq.ny-2) then
84         jyvp=ny-2
85         jumpy=1
86      end if
87      juy=jumpy
88  !
89    do ix=0,nxmin1
90  ! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long)
91      ixvp=ix+1
92      ixvm=ix-1
93      jumpx=2
94      if (xglobal) then
95         ivrp=ixvp
96         ivrm=ixvm
97         if (ixvm.lt.0) ivrm=ixvm+nxmin1
98         if (ixvp.ge.nx) ivrp=ixvp-nx+1
99      else
100        if (ix.eq.0) ixvm=0
101        if (ix.eq.nxmin1) ixvp=nxmin1
102        ivrp=ixvp
103        ivrm=ixvm
104  ! Define absolute gap length
105        if (ix.eq.0.or.ix.eq.nxmin1) jumpx=1
106      end if
107      jux=jumpx
108  ! Precalculate pressure values for efficiency
109      do kl=1,nuvz
110        ppml(kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
111      end do
112  !
113  ! Loop over the vertical
114  !***********************
115
116      do kl=1,nuvz
117        ppmk=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
118        theta=tth(ix,jy,kl,n)*(100000./ppmk)**kappa
119        klvrp=kl+1
120        klvrm=kl-1
121        klpt=kl
122  ! If top or bottom level, dthetadp is evaluated between the current
123  ! level and the level inside, otherwise between level+1 and level-1
124  !
125        if (klvrp.gt.nuvz) klvrp=nuvz
126        if (klvrm.lt.1) klvrm=1
127        ppmk=akz(klvrp)+bkz(klvrp)*ps(ix,jy,1,n)
128        thetap=tth(ix,jy,klvrp,n)*(100000./ppmk)**kappa
129        ppmk=akz(klvrm)+bkz(klvrm)*ps(ix,jy,1,n)
130        thetam=tth(ix,jy,klvrm,n)*(100000./ppmk)**kappa
131        dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))
132
133  ! Compute vertical position at pot. temperature surface on subgrid
134  ! and the wind at that position
135  !*****************************************************************
136  ! a) in x direction
137        ii=0
138        do i=ixvm,ixvp,jumpx
139          ivr=i
140          if (xglobal) then
141             if (i.lt.0) ivr=ivr+nxmin1
142             if (i.ge.nx) ivr=ivr-nx+1
143          end if
144          ii=ii+1
145  ! Search adjacent levels for current theta value
146  ! Spiral out from current level for efficiency
147          kup=klpt-1
148          kdn=klpt
149          kch=0
15040        continue
151  ! Upward branch
152          kup=kup+1
153          if (kch.ge.nlck) goto 21     ! No more levels to check,
154  !                                       ! and no values found
155          if (kup.ge.nuvz) goto 41
156          kch=kch+1
157          k=kup
158          ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n)
159          thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa
160          ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n)
161          thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa
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          ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n)
184          thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa
185          ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n)
186          thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa
187
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
197              endif
198          vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
199              goto 20
200            endif
201            goto 40
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
238          ppmk=akz(k)+bkz(k)*ps(ix,j,1,n)
239          thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa
240          ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n)
241          thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa
242      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
243           ((thdn.le.theta).and.(thup.ge.theta))) then
244              dt1=abs(theta-thdn)
245              dt2=abs(theta-thup)
246              dt=dt1+dt2
247              if (dt.lt.eps) then   ! Avoid division by zero error
248                dt1=0.5             ! G.W., 10.4.1996
249                dt2=0.5
250                dt=1.0
251              endif
252              uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
253              goto 50
254            endif
25571        continue
256  ! Downward branch
257          kdn=kdn-1
258          if (kdn.lt.1) goto 70
259          kch=kch+1
260          k=kdn
261          ppmk=akz(k)+bkz(k)*ps(ix,j,1,n)
262          thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa
263          ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n)
264          thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa
265      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
266           ((thdn.le.theta).and.(thup.ge.theta))) then
267              dt1=abs(theta-thdn)
268              dt2=abs(theta-thup)
269              dt=dt1+dt2
270              if (dt.lt.eps) then   ! Avoid division by zero error
271                dt1=0.5             ! G.W., 10.4.1996
272                dt2=0.5
273                dt=1.0
274              endif
275              uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
276              goto 50
277            endif
278            goto 70
279  ! This section used when no values were found
28051      continue
281  ! Must use uu at current level and lat. juy becomes smaller by 1
282        uy(jj)=uuh(ix,jy,kl)
283        juy=juy-1
284  ! Otherwise OK
28550        continue
286        end do
287      if (juy.gt.0) then
288      dudy=(uy(2)-uy(1))/real(juy)/(dy*pi/180.)
289      else
290      dudy=uuh(ix,jyvp,kl)-uuh(ix,jyvm,kl)
291      dudy=dudy/real(jumpy)/(dy*pi/180.)
292      end if
293  !
294      pvh(ix,jy,kl)=dthetadp*(f+(dvdx/cosphi-dudy &
295           +uuh(ix,jy,kl)*tanphi)/r_earth)*(-1.e6)*9.81
296
297
298  !
299  ! Resest jux and juy
300      jux=jumpx
301      juy=jumpy
302      end do
303    end do
30410  continue
305  end do
306  !
307  ! Fill in missing PV values on poles, if present
308  ! Use mean PV of surrounding latitude ring
309  !
310      if (sglobal) then
311         do kl=1,nuvz
312            pvavr=0.
313            do ix=0,nxmin1
314               pvavr=pvavr+pvh(ix,1,kl)
315            end do
316            pvavr=pvavr/real(nx)
317            jy=0
318            do ix=0,nxmin1
319               pvh(ix,jy,kl)=pvavr
320            end do
321         end do
322      end if
323      if (nglobal) then
324         do kl=1,nuvz
325            pvavr=0.
326            do ix=0,nxmin1
327               pvavr=pvavr+pvh(ix,ny-2,kl)
328            end do
329            pvavr=pvavr/real(nx)
330            jy=nymin1
331            do ix=0,nxmin1
332               pvh(ix,jy,kl)=pvavr
333            end do
334         end do
335      end if
336
337end subroutine calcpv
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG