source: branches/flexpart91_hasod/src_parallel/calcpv.f90 @ 8

Last change on this file since 8 was 8, checked in by hasod, 11 years ago

Added parallel version of Flexpart91

File size: 12.3 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, only: pi, nxmax, nymax, kappa, r_earth, nuvzmax
43  use com_mod, only: tth, dx, dy, nglobal, ps, akz, bkz, ylat0, xlon0, xglobal, &
44                     sglobal, nglobal, nymin1, nxmin1, nx, ny, nuvz
45
46  implicit none
47
48  integer :: n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch
49  integer :: jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr
50  integer :: nlck
51  real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
52  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
53  real :: pvavr,ppml(nuvzmax)
54  real :: thup,thdn
55  real,parameter :: eps=1.e-5, p0=101325
56  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
57  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
58  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
59
60  ! Set number of levels to check for adjacent theta
61  nlck=nuvz/3
62
63!$OMP PARALLEL DEFAULT(none) &
64!$OMP SHARED(nlck, nuvz, nxmin1, nymin1, sglobal, nglobal, dx, dy, xlon0, ylat0,  &
65!$OMP   nx, ny, xglobal, n, akz, bkz, ps, tth, uuh, vvh, pvh) &
66!$OMP PRIVATE(ix, jy, ii, i, phi, f, cosphi, tanphi, jyvp, jyvm, jumpy, juy, &
67!$OMP   ixvp, ixvm, jumpx, ivrp, ivrm, jux, ppml, ppmk, theta, klvrp, klvrm, klpt, &
68!$OMP   thetap, thetam, dthetadp, ivr, kup, kdn, kch, k, thdn, thup, dt1, dt2, dt, &
69!$OMP   vx, dvdx, jj, uy, dudy)
70
71#if (defined STATIC_SCHED)
72!$OMP DO SCHEDULE(static)
73#else
74!$OMP DO SCHEDULE(dynamic, max(1,nymin1/10))
75#endif
76  !
77  ! Loop over entire grid
78  !**********************
79  latloop : do jy=0,nymin1
80    if (sglobal.and.jy.eq.0) goto 10
81    if (nglobal.and.jy.eq.nymin1) goto 10
82    phi = (ylat0 + jy * dy) * pi / 180.
83    f = 0.00014585 * sin(phi)
84    tanphi = tan(phi)
85    cosphi = cos(phi)
86  ! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat)
87      jyvp=jy+1
88      jyvm=jy-1
89      if (jy.eq.0) jyvm=0
90      if (jy.eq.nymin1) jyvp=nymin1
91  ! Define absolute gap length
92      jumpy=2
93      if (jy.eq.0.or.jy.eq.nymin1) jumpy=1
94      if (sglobal.and.jy.eq.1) then
95         jyvm=1
96         jumpy=1
97      end if
98      if (nglobal.and.jy.eq.ny-2) then
99         jyvp=ny-2
100         jumpy=1
101      end if
102      juy=jumpy
103  !
104    lonloop : do ix=0,nxmin1
105  ! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long)
106      ixvp=ix+1
107      ixvm=ix-1
108      jumpx=2
109      if (xglobal) then
110         ivrp=ixvp
111         ivrm=ixvm
112         if (ixvm.lt.0) ivrm=ixvm+nxmin1
113         if (ixvp.ge.nx) ivrp=ixvp-nx+1
114      else
115        if (ix.eq.0) ixvm=0
116        if (ix.eq.nxmin1) ixvp=nxmin1
117        ivrp=ixvp
118        ivrm=ixvm
119  ! Define absolute gap length
120        if (ix.eq.0.or.ix.eq.nxmin1) jumpx=1
121      end if
122      jux=jumpx
123  ! Precalculate pressure values for efficiency
124      do kl=1,nuvz
125        ppml(kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
126      end do
127  !
128  ! Loop over the vertical
129  !***********************
130
131      do kl=1,nuvz
132        ppmk=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
133        theta=tth(ix,jy,kl,n)*(100000./ppmk)**kappa
134        klvrp=kl+1
135        klvrm=kl-1
136        klpt=kl
137  ! If top or bottom level, dthetadp is evaluated between the current
138  ! level and the level inside, otherwise between level+1 and level-1
139  !
140        if (klvrp.gt.nuvz) klvrp=nuvz
141        if (klvrm.lt.1) klvrm=1
142        ppmk=akz(klvrp)+bkz(klvrp)*ps(ix,jy,1,n)
143        thetap=tth(ix,jy,klvrp,n)*(100000./ppmk)**kappa
144        ppmk=akz(klvrm)+bkz(klvrm)*ps(ix,jy,1,n)
145        thetam=tth(ix,jy,klvrm,n)*(100000./ppmk)**kappa
146        dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))
147
148  ! Compute vertical position at pot. temperature surface on subgrid
149  ! and the wind at that position
150  !*****************************************************************
151  ! a) in x direction
152        ii=0
153        do i=ixvm,ixvp,jumpx
154          ivr=i
155          if (xglobal) then
156             if (i.lt.0) ivr=ivr+nxmin1
157             if (i.ge.nx) ivr=ivr-nx+1
158          end if
159          ii=ii+1
160  ! Search adjacent levels for current theta value
161  ! Spiral out from current level for efficiency
162          kup=klpt-1
163          kdn=klpt
164          kch=0
16540        continue
166  ! Upward branch
167          kup=kup+1
168          if (kch.ge.nlck) goto 21     ! No more levels to check,
169  !                                       ! and no values found
170          if (kup.ge.nuvz) goto 41
171          kch=kch+1
172          k=kup
173          ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n)
174          thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa
175          ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n)
176          thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa
177
178
179      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
180           ((thdn.le.theta).and.(thup.ge.theta))) then
181              dt1=abs(theta-thdn)
182              dt2=abs(theta-thup)
183              dt=dt1+dt2
184              if (dt.lt.eps) then   ! Avoid division by zero error
185                dt1=0.5             ! G.W., 10.4.1996
186                dt2=0.5
187                dt=1.0
188              endif
189          vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
190              goto 20
191            endif
19241        continue
193  ! Downward branch
194          kdn=kdn-1
195          if (kdn.lt.1) goto 40
196          kch=kch+1
197          k=kdn
198          ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n)
199          thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa
200          ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n)
201          thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa
202
203      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
204           ((thdn.le.theta).and.(thup.ge.theta))) then
205              dt1=abs(theta-thdn)
206              dt2=abs(theta-thup)
207              dt=dt1+dt2
208              if (dt.lt.eps) then   ! Avoid division by zero error
209                dt1=0.5             ! G.W., 10.4.1996
210                dt2=0.5
211                dt=1.0
212              endif
213          vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
214              goto 20
215            endif
216            goto 40
217  ! This section used when no values were found
21821      continue
219  ! Must use vv at current level and long. jux becomes smaller by 1
220        vx(ii)=vvh(ix,jy,kl)
221        jux=jux-1
222  ! Otherwise OK
22320        continue
224        end do
225      if (jux.gt.0) then
226      dvdx=(vx(2)-vx(1))/real(jux)/(dx*pi/180.)
227      else
228      dvdx=vvh(ivrp,jy,kl)-vvh(ivrm,jy,kl)
229      dvdx=dvdx/real(jumpx)/(dx*pi/180.)
230  ! Only happens if no equivalent theta value
231  ! can be found on either side, hence must use values
232  ! from either side, same pressure level.
233      end if
234
235  ! b) in y direction
236
237        jj=0
238        do j=jyvm,jyvp,jumpy
239          jj=jj+1
240  ! Search adjacent levels for current theta value
241  ! Spiral out from current level for efficiency
242          kup=klpt-1
243          kdn=klpt
244          kch=0
24570        continue
246  ! Upward branch
247          kup=kup+1
248          if (kch.ge.nlck) goto 51     ! No more levels to check,
249  !                                     ! and no values found
250          if (kup.ge.nuvz) goto 71
251          kch=kch+1
252          k=kup
253          ppmk=akz(k)+bkz(k)*ps(ix,j,1,n)
254          thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa
255          ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n)
256          thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa
257      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
258           ((thdn.le.theta).and.(thup.ge.theta))) then
259              dt1=abs(theta-thdn)
260              dt2=abs(theta-thup)
261              dt=dt1+dt2
262              if (dt.lt.eps) then   ! Avoid division by zero error
263                dt1=0.5             ! G.W., 10.4.1996
264                dt2=0.5
265                dt=1.0
266              endif
267              uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
268              goto 50
269            endif
27071        continue
271  ! Downward branch
272          kdn=kdn-1
273          if (kdn.lt.1) goto 70
274          kch=kch+1
275          k=kdn
276          ppmk=akz(k)+bkz(k)*ps(ix,j,1,n)
277          thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa
278          ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n)
279          thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa
280      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
281           ((thdn.le.theta).and.(thup.ge.theta))) then
282              dt1=abs(theta-thdn)
283              dt2=abs(theta-thup)
284              dt=dt1+dt2
285              if (dt.lt.eps) then   ! Avoid division by zero error
286                dt1=0.5             ! G.W., 10.4.1996
287                dt2=0.5
288                dt=1.0
289              endif
290              uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
291              goto 50
292            endif
293            goto 70
294  ! This section used when no values were found
29551      continue
296  ! Must use uu at current level and lat. juy becomes smaller by 1
297        uy(jj)=uuh(ix,jy,kl)
298        juy=juy-1
299  ! Otherwise OK
30050        continue
301        end do
302      if (juy.gt.0) then
303      dudy=(uy(2)-uy(1))/real(juy)/(dy*pi/180.)
304      else
305      dudy=uuh(ix,jyvp,kl)-uuh(ix,jyvm,kl)
306      dudy=dudy/real(jumpy)/(dy*pi/180.)
307      end if
308  !
309      pvh(ix,jy,kl)=dthetadp*(f+(dvdx/cosphi-dudy &
310           +uuh(ix,jy,kl)*tanphi)/r_earth)*(-1.e6)*9.81
311
312
313  !
314  ! Resest jux and juy
315      jux=jumpx
316      juy=jumpy
317      end do
318    end do lonloop
31910  continue
320  end do latloop
321
322!$OMP END DO
323!$OMP END PARALLEL
324
325  !
326  ! Fill in missing PV values on poles, if present
327  ! Use mean PV of surrounding latitude ring
328  !
329      if (sglobal) then
330         do kl=1,nuvz
331            pvavr=0.
332            do ix=0,nxmin1
333               pvavr=pvavr+pvh(ix,1,kl)
334            end do
335            pvavr=pvavr/real(nx)
336            jy=0
337            do ix=0,nxmin1
338               pvh(ix,jy,kl)=pvavr
339            end do
340         end do
341      end if
342      if (nglobal) then
343         do kl=1,nuvz
344            pvavr=0.
345            do ix=0,nxmin1
346               pvavr=pvavr+pvh(ix,ny-2,kl)
347            end do
348            pvavr=pvavr/real(nx)
349            jy=nymin1
350            do ix=0,nxmin1
351               pvh(ix,jy,kl)=pvavr
352            end do
353         end do
354      end if
355
356end subroutine calcpv
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG