source: flexpart.git/src/calcpv.f90 @ b0434e1

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since b0434e1 was 8a65cb0, checked in by Espen Sollum ATMOS <espen@…>, 9 years ago

Added code, makefile for dev branch

  • Property mode set to 100644
File size: 11.2 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! :DBG: halts with SIGFPE if compiling with -ffpe-trap
73!  where(ppml.le.0) ppml=10000.0
74!  ppmk(:,:,1:nuvz)=(100000./ppml(:,:,1:nuvz))**kappa
[4fbe7a5]75  ppmk=(100000./ppml)**kappa
76
[e200b7a]77  do jy=0,nymin1
78    if (sglobal.and.jy.eq.0) goto 10
79    if (nglobal.and.jy.eq.nymin1) goto 10
80    phi = (ylat0 + jy * dy) * pi / 180.
81    f = 0.00014585 * sin(phi)
82    tanphi = tan(phi)
83    cosphi = cos(phi)
84  ! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat)
85      jyvp=jy+1
86      jyvm=jy-1
87      if (jy.eq.0) jyvm=0
88      if (jy.eq.nymin1) jyvp=nymin1
89  ! Define absolute gap length
90      jumpy=2
91      if (jy.eq.0.or.jy.eq.nymin1) jumpy=1
92      if (sglobal.and.jy.eq.1) then
93         jyvm=1
94         jumpy=1
95      end if
96      if (nglobal.and.jy.eq.ny-2) then
97         jyvp=ny-2
98         jumpy=1
99      end if
100      juy=jumpy
101  !
102    do ix=0,nxmin1
103  ! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long)
104      ixvp=ix+1
105      ixvm=ix-1
106      jumpx=2
107      if (xglobal) then
108         ivrp=ixvp
109         ivrm=ixvm
110         if (ixvm.lt.0) ivrm=ixvm+nxmin1
111         if (ixvp.ge.nx) ivrp=ixvp-nx+1
112      else
113        if (ix.eq.0) ixvm=0
114        if (ix.eq.nxmin1) ixvp=nxmin1
115        ivrp=ixvp
116        ivrm=ixvm
117  ! Define absolute gap length
118        if (ix.eq.0.or.ix.eq.nxmin1) jumpx=1
119      end if
120      jux=jumpx
121  !
122  ! Loop over the vertical
123  !***********************
124
125      do kl=1,nuvz
[4fbe7a5]126        theta=tth(ix,jy,kl,n)*ppmk(ix,jy,kl)
[e200b7a]127        klvrp=kl+1
128        klvrm=kl-1
129        klpt=kl
130  ! If top or bottom level, dthetadp is evaluated between the current
131  ! level and the level inside, otherwise between level+1 and level-1
132  !
133        if (klvrp.gt.nuvz) klvrp=nuvz
134        if (klvrm.lt.1) klvrm=1
[4fbe7a5]135        thetap=tth(ix,jy,klvrp,n)*ppmk(ix,jy,klvrp)
136        thetam=tth(ix,jy,klvrm,n)*ppmk(ix,jy,klvrm)
137        dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm))
[e200b7a]138
139  ! Compute vertical position at pot. temperature surface on subgrid
140  ! and the wind at that position
141  !*****************************************************************
142  ! a) in x direction
143        ii=0
144        do i=ixvm,ixvp,jumpx
145          ivr=i
146          if (xglobal) then
147             if (i.lt.0) ivr=ivr+nxmin1
148             if (i.ge.nx) ivr=ivr-nx+1
149          end if
150          ii=ii+1
151  ! Search adjacent levels for current theta value
152  ! Spiral out from current level for efficiency
153          kup=klpt-1
154          kdn=klpt
155          kch=0
15640        continue
157  ! Upward branch
158          kup=kup+1
159          if (kch.ge.nlck) goto 21     ! No more levels to check,
160  !                                       ! and no values found
161          if (kup.ge.nuvz) goto 41
162          kch=kch+1
163          k=kup
[4fbe7a5]164          thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k)
165          thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1)
[e200b7a]166
167
[4fbe7a5]168          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
169          ((thdn.le.theta).and.(thup.ge.theta))) then
170            dt1=abs(theta-thdn)
171            dt2=abs(theta-thup)
172            dt=dt1+dt2
173            if (dt.lt.eps) then   ! Avoid division by zero error
174              dt1=0.5             ! G.W., 10.4.1996
175              dt2=0.5
176              dt=1.0
[e200b7a]177            endif
[4fbe7a5]178            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
179            goto 20
180          endif
[e200b7a]18141        continue
182  ! Downward branch
183          kdn=kdn-1
184          if (kdn.lt.1) goto 40
185          kch=kch+1
186          k=kdn
[4fbe7a5]187          thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k)
188          thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1)
[e200b7a]189
[4fbe7a5]190          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
191          ((thdn.le.theta).and.(thup.ge.theta))) then
192            dt1=abs(theta-thdn)
193            dt2=abs(theta-thup)
194            dt=dt1+dt2
195            if (dt.lt.eps) then   ! Avoid division by zero error
196              dt1=0.5             ! G.W., 10.4.1996
197              dt2=0.5
198              dt=1.0
[e200b7a]199            endif
[4fbe7a5]200            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
201            goto 20
202          endif
203          goto 40
[e200b7a]204  ! This section used when no values were found
20521      continue
206  ! Must use vv at current level and long. jux becomes smaller by 1
207        vx(ii)=vvh(ix,jy,kl)
208        jux=jux-1
209  ! Otherwise OK
21020        continue
211        end do
212      if (jux.gt.0) then
213      dvdx=(vx(2)-vx(1))/real(jux)/(dx*pi/180.)
214      else
215      dvdx=vvh(ivrp,jy,kl)-vvh(ivrm,jy,kl)
216      dvdx=dvdx/real(jumpx)/(dx*pi/180.)
217  ! Only happens if no equivalent theta value
218  ! can be found on either side, hence must use values
219  ! from either side, same pressure level.
220      end if
221
222  ! b) in y direction
223
224        jj=0
225        do j=jyvm,jyvp,jumpy
226          jj=jj+1
227  ! Search adjacent levels for current theta value
228  ! Spiral out from current level for efficiency
229          kup=klpt-1
230          kdn=klpt
231          kch=0
23270        continue
233  ! Upward branch
234          kup=kup+1
235          if (kch.ge.nlck) goto 51     ! No more levels to check,
236  !                                     ! and no values found
237          if (kup.ge.nuvz) goto 71
238          kch=kch+1
239          k=kup
[4fbe7a5]240          thdn=tth(ix,j,k,n)*ppmk(ix,j,k)
241          thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1)
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
[e200b7a]251            endif
[4fbe7a5]252            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
253            goto 50
254          endif
[e200b7a]25571        continue
256  ! Downward branch
257          kdn=kdn-1
258          if (kdn.lt.1) goto 70
259          kch=kch+1
260          k=kdn
[4fbe7a5]261          thdn=tth(ix,j,k,n)*ppmk(ix,j,k)
262          thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1)
263          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
264          ((thdn.le.theta).and.(thup.ge.theta))) then
265            dt1=abs(theta-thdn)
266            dt2=abs(theta-thup)
267            dt=dt1+dt2
268            if (dt.lt.eps) then   ! Avoid division by zero error
269              dt1=0.5             ! G.W., 10.4.1996
270              dt2=0.5
271              dt=1.0
[e200b7a]272            endif
[4fbe7a5]273            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
274            goto 50
275          endif
276          goto 70
[e200b7a]277  ! This section used when no values were found
27851      continue
279  ! Must use uu at current level and lat. juy becomes smaller by 1
[4fbe7a5]280          uy(jj)=uuh(ix,jy,kl)
281          juy=juy-1
[e200b7a]282  ! Otherwise OK
28350        continue
284        end do
285      if (juy.gt.0) then
286      dudy=(uy(2)-uy(1))/real(juy)/(dy*pi/180.)
287      else
288      dudy=uuh(ix,jyvp,kl)-uuh(ix,jyvm,kl)
289      dudy=dudy/real(jumpy)/(dy*pi/180.)
290      end if
291  !
292      pvh(ix,jy,kl)=dthetadp*(f+(dvdx/cosphi-dudy &
293           +uuh(ix,jy,kl)*tanphi)/r_earth)*(-1.e6)*9.81
294
295
296  !
297  ! Resest jux and juy
298      jux=jumpx
299      juy=jumpy
300      end do
301    end do
30210  continue
303  end do
304  !
305  ! Fill in missing PV values on poles, if present
306  ! Use mean PV of surrounding latitude ring
307  !
308      if (sglobal) then
309         do kl=1,nuvz
310            pvavr=0.
311            do ix=0,nxmin1
312               pvavr=pvavr+pvh(ix,1,kl)
313            end do
314            pvavr=pvavr/real(nx)
315            jy=0
316            do ix=0,nxmin1
317               pvh(ix,jy,kl)=pvavr
318            end do
319         end do
320      end if
321      if (nglobal) then
322         do kl=1,nuvz
323            pvavr=0.
324            do ix=0,nxmin1
325               pvavr=pvavr+pvh(ix,ny-2,kl)
326            end do
327            pvavr=pvavr/real(nx)
328            jy=nymin1
329            do ix=0,nxmin1
330               pvh(ix,jy,kl)=pvavr
331            end do
332         end do
333      end if
334
335end subroutine calcpv
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG