source: flexpart.git/src/calcpv.f90

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 4 years ago

add SPDX-License-Identifier to all .f90 files

  • Property mode set to 100644
File size: 9.8 KB
RevLine 
[92fab65]1! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2! SPDX-License-Identifier: GPL-3.0-or-later
[332fbbd]3
[e200b7a]4subroutine calcpv(n,uuh,vvh,pvh)
[4fbe7a5]5  !               i  i   i   o
[e200b7a]6  !*****************************************************************************
7  !                                                                            *
8  !  Calculation of potential vorticity on 3-d grid.                           *
9  !                                                                            *
10  !     Author: P. James                                                       *
11  !     3 February 2000                                                        *
12  !                                                                            *
13  !     Adaptation to FLEXPART, A. Stohl, 1 May 2000                           *
14  !                                                                            *
15  !*****************************************************************************
16  !                                                                            *
17  ! Variables:                                                                 *
18  ! n                  temporal index for meteorological fields (1 to 2)       *
19  !                                                                            *
20  ! Constants:                                                                 *
21  !                                                                            *
22  !*****************************************************************************
23
24  use par_mod
25  use com_mod
26
27  implicit none
28
29  integer :: n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch
30  integer :: jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr
31  integer :: nlck
32  real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
[4fbe7a5]33  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt
34  real :: pvavr,ppml(0:nxmax-1,0:nymax-1,nuvzmax),ppmk(0:nxmax-1,0:nymax-1,nuvzmax)
[e200b7a]35  real :: thup,thdn
36  real,parameter :: eps=1.e-5, p0=101325
37  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
38  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
39  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
40
41  ! Set number of levels to check for adjacent theta
42  nlck=nuvz/3
43  !
44  ! Loop over entire grid
45  !**********************
[4fbe7a5]46  do kl=1,nuvz
47    do jy=0,nymin1
48      do ix=0,nxmin1
49         ppml(ix,jy,kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
50      enddo
51    enddo
52  enddo
[8a65cb0]53
54!  ppmk(:,:,1:nuvz)=(100000./ppml(:,:,1:nuvz))**kappa
[db712a8]55  ppmk(0:nxmin1,0:nymin1,1:nuvz)=(100000./ppml(0:nxmin1,0:nymin1,1:nuvz))**kappa
[4fbe7a5]56
[e200b7a]57  do jy=0,nymin1
58    if (sglobal.and.jy.eq.0) goto 10
59    if (nglobal.and.jy.eq.nymin1) goto 10
60    phi = (ylat0 + jy * dy) * pi / 180.
61    f = 0.00014585 * sin(phi)
62    tanphi = tan(phi)
63    cosphi = cos(phi)
64  ! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat)
65      jyvp=jy+1
66      jyvm=jy-1
67      if (jy.eq.0) jyvm=0
68      if (jy.eq.nymin1) jyvp=nymin1
69  ! Define absolute gap length
70      jumpy=2
71      if (jy.eq.0.or.jy.eq.nymin1) jumpy=1
72      if (sglobal.and.jy.eq.1) then
73         jyvm=1
74         jumpy=1
75      end if
76      if (nglobal.and.jy.eq.ny-2) then
77         jyvp=ny-2
78         jumpy=1
79      end if
80      juy=jumpy
81  !
82    do ix=0,nxmin1
83  ! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long)
84      ixvp=ix+1
85      ixvm=ix-1
86      jumpx=2
87      if (xglobal) then
88         ivrp=ixvp
89         ivrm=ixvm
90         if (ixvm.lt.0) ivrm=ixvm+nxmin1
91         if (ixvp.ge.nx) ivrp=ixvp-nx+1
92      else
93        if (ix.eq.0) ixvm=0
94        if (ix.eq.nxmin1) ixvp=nxmin1
95        ivrp=ixvp
96        ivrm=ixvm
97  ! Define absolute gap length
98        if (ix.eq.0.or.ix.eq.nxmin1) jumpx=1
99      end if
100      jux=jumpx
101  !
102  ! Loop over the vertical
103  !***********************
104
105      do kl=1,nuvz
[4fbe7a5]106        theta=tth(ix,jy,kl,n)*ppmk(ix,jy,kl)
[e200b7a]107        klvrp=kl+1
108        klvrm=kl-1
109        klpt=kl
110  ! If top or bottom level, dthetadp is evaluated between the current
111  ! level and the level inside, otherwise between level+1 and level-1
112  !
113        if (klvrp.gt.nuvz) klvrp=nuvz
114        if (klvrm.lt.1) klvrm=1
[4fbe7a5]115        thetap=tth(ix,jy,klvrp,n)*ppmk(ix,jy,klvrp)
116        thetam=tth(ix,jy,klvrm,n)*ppmk(ix,jy,klvrm)
117        dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm))
[e200b7a]118
119  ! Compute vertical position at pot. temperature surface on subgrid
120  ! and the wind at that position
121  !*****************************************************************
122  ! a) in x direction
123        ii=0
124        do i=ixvm,ixvp,jumpx
125          ivr=i
126          if (xglobal) then
127             if (i.lt.0) ivr=ivr+nxmin1
128             if (i.ge.nx) ivr=ivr-nx+1
129          end if
130          ii=ii+1
131  ! Search adjacent levels for current theta value
132  ! Spiral out from current level for efficiency
133          kup=klpt-1
134          kdn=klpt
135          kch=0
13640        continue
137  ! Upward branch
138          kup=kup+1
139          if (kch.ge.nlck) goto 21     ! No more levels to check,
140  !                                       ! and no values found
141          if (kup.ge.nuvz) goto 41
142          kch=kch+1
143          k=kup
[4fbe7a5]144          thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k)
145          thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1)
[e200b7a]146
147
[4fbe7a5]148          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
149          ((thdn.le.theta).and.(thup.ge.theta))) then
150            dt1=abs(theta-thdn)
151            dt2=abs(theta-thup)
152            dt=dt1+dt2
153            if (dt.lt.eps) then   ! Avoid division by zero error
154              dt1=0.5             ! G.W., 10.4.1996
155              dt2=0.5
156              dt=1.0
[e200b7a]157            endif
[4fbe7a5]158            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
159            goto 20
160          endif
[e200b7a]16141        continue
162  ! Downward branch
163          kdn=kdn-1
164          if (kdn.lt.1) goto 40
165          kch=kch+1
166          k=kdn
[4fbe7a5]167          thdn=tth(ivr,jy,k,n)*ppmk(ivr,jy,k)
168          thup=tth(ivr,jy,k+1,n)*ppmk(ivr,jy,k+1)
[e200b7a]169
[4fbe7a5]170          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
171          ((thdn.le.theta).and.(thup.ge.theta))) then
172            dt1=abs(theta-thdn)
173            dt2=abs(theta-thup)
174            dt=dt1+dt2
175            if (dt.lt.eps) then   ! Avoid division by zero error
176              dt1=0.5             ! G.W., 10.4.1996
177              dt2=0.5
178              dt=1.0
[e200b7a]179            endif
[4fbe7a5]180            vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
181            goto 20
182          endif
183          goto 40
[e200b7a]184  ! This section used when no values were found
18521      continue
186  ! Must use vv at current level and long. jux becomes smaller by 1
187        vx(ii)=vvh(ix,jy,kl)
188        jux=jux-1
189  ! Otherwise OK
19020        continue
191        end do
192      if (jux.gt.0) then
193      dvdx=(vx(2)-vx(1))/real(jux)/(dx*pi/180.)
194      else
195      dvdx=vvh(ivrp,jy,kl)-vvh(ivrm,jy,kl)
196      dvdx=dvdx/real(jumpx)/(dx*pi/180.)
197  ! Only happens if no equivalent theta value
198  ! can be found on either side, hence must use values
199  ! from either side, same pressure level.
200      end if
201
202  ! b) in y direction
203
204        jj=0
205        do j=jyvm,jyvp,jumpy
206          jj=jj+1
207  ! Search adjacent levels for current theta value
208  ! Spiral out from current level for efficiency
209          kup=klpt-1
210          kdn=klpt
211          kch=0
21270        continue
213  ! Upward branch
214          kup=kup+1
215          if (kch.ge.nlck) goto 51     ! No more levels to check,
216  !                                     ! and no values found
217          if (kup.ge.nuvz) goto 71
218          kch=kch+1
219          k=kup
[4fbe7a5]220          thdn=tth(ix,j,k,n)*ppmk(ix,j,k)
221          thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1)
222          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
223          ((thdn.le.theta).and.(thup.ge.theta))) then
224            dt1=abs(theta-thdn)
225            dt2=abs(theta-thup)
226            dt=dt1+dt2
227            if (dt.lt.eps) then   ! Avoid division by zero error
228              dt1=0.5             ! G.W., 10.4.1996
229              dt2=0.5
230              dt=1.0
[e200b7a]231            endif
[4fbe7a5]232            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
233            goto 50
234          endif
[e200b7a]23571        continue
236  ! Downward branch
237          kdn=kdn-1
238          if (kdn.lt.1) goto 70
239          kch=kch+1
240          k=kdn
[4fbe7a5]241          thdn=tth(ix,j,k,n)*ppmk(ix,j,k)
242          thup=tth(ix,j,k+1,n)*ppmk(ix,j,k+1)
243          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
244          ((thdn.le.theta).and.(thup.ge.theta))) then
245            dt1=abs(theta-thdn)
246            dt2=abs(theta-thup)
247            dt=dt1+dt2
248            if (dt.lt.eps) then   ! Avoid division by zero error
249              dt1=0.5             ! G.W., 10.4.1996
250              dt2=0.5
251              dt=1.0
[e200b7a]252            endif
[4fbe7a5]253            uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
254            goto 50
255          endif
256          goto 70
[e200b7a]257  ! This section used when no values were found
25851      continue
259  ! Must use uu at current level and lat. juy becomes smaller by 1
[4fbe7a5]260          uy(jj)=uuh(ix,jy,kl)
261          juy=juy-1
[e200b7a]262  ! Otherwise OK
26350        continue
264        end do
265      if (juy.gt.0) then
266      dudy=(uy(2)-uy(1))/real(juy)/(dy*pi/180.)
267      else
268      dudy=uuh(ix,jyvp,kl)-uuh(ix,jyvm,kl)
269      dudy=dudy/real(jumpy)/(dy*pi/180.)
270      end if
271  !
272      pvh(ix,jy,kl)=dthetadp*(f+(dvdx/cosphi-dudy &
273           +uuh(ix,jy,kl)*tanphi)/r_earth)*(-1.e6)*9.81
274
275
276  !
277  ! Resest jux and juy
278      jux=jumpx
279      juy=jumpy
280      end do
281    end do
28210  continue
283  end do
284  !
285  ! Fill in missing PV values on poles, if present
286  ! Use mean PV of surrounding latitude ring
287  !
288      if (sglobal) then
289         do kl=1,nuvz
290            pvavr=0.
291            do ix=0,nxmin1
292               pvavr=pvavr+pvh(ix,1,kl)
293            end do
294            pvavr=pvavr/real(nx)
295            jy=0
296            do ix=0,nxmin1
297               pvh(ix,jy,kl)=pvavr
298            end do
299         end do
300      end if
301      if (nglobal) then
302         do kl=1,nuvz
303            pvavr=0.
304            do ix=0,nxmin1
305               pvavr=pvavr+pvh(ix,ny-2,kl)
306            end do
307            pvavr=pvavr/real(nx)
308            jy=nymin1
309            do ix=0,nxmin1
310               pvh(ix,jy,kl)=pvavr
311            end do
312         end do
313      end if
314
315end subroutine calcpv
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG