source: flexpart.git/src/calcpv_nests.f90 @ 5f2e8f6

flexpart-noresm
Last change on this file since 5f2e8f6 was 5f2e8f6, checked in by Ignacio Pisso <Ignacio.Pisso@…>, 8 years ago

new flexpart noresm code in src

  • Property mode set to 100755
File size: 10.4 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_nests(l,n,uuhn,vvhn,pvhn)
23  !                        i i  i    i    o
24  !*****************************************************************************
25  !                                                                            *
26  !  Calculation of potential vorticity on 3-d nested grid                     *
27  !                                                                            *
28  !     Author: P. James                                                       *
29  !     22 February 2000                                                       *
30  !                                                                            *
31  !*****************************************************************************
32  !                                                                            *
33  ! Variables:                                                                 *
34  ! n                  temporal index for meteorological fields (1 to 2)       *
35  ! l                  index of current nest                                   *
36  !                                                                            *
37  ! Constants:                                                                 *
38  !                                                                            *
39  !*****************************************************************************
40
41  use par_mod
42  use com_mod
43
44  implicit none
45
46  integer :: n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch
47  integer :: jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr
48  integer :: nlck,l
49  real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
50  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
51  real :: ppml(nuvzmax)
52  real :: thup,thdn
53  real,parameter :: eps=1.e-5,p0=101325
54  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
55  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
56  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
57
58  ! Set number of levels to check for adjacent theta
59  nlck=nuvz/3
60  !
61  ! Loop over entire grid
62  !**********************
63  do jy=0,nyn(l)-1
64    phi = (ylat0n(l) + jy * dyn(l)) * pi / 180.
65    f = 0.00014585 * sin(phi)
66    tanphi = tan(phi)
67    cosphi = cos(phi)
68  ! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat)
69      jyvp=jy+1
70      jyvm=jy-1
71      if (jy.eq.0) jyvm=0
72      if (jy.eq.nyn(l)-1) jyvp=nyn(l)-1
73  ! Define absolute gap length
74      jumpy=2
75      if (jy.eq.0.or.jy.eq.nyn(l)-1) jumpy=1
76      juy=jumpy
77  !
78    do ix=0,nxn(l)-1
79  ! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long)
80      ixvp=ix+1
81      ixvm=ix-1
82      jumpx=2
83      if (ix.eq.0) ixvm=0
84      if (ix.eq.nxn(l)-1) ixvp=nxn(l)-1
85      ivrp=ixvp
86      ivrm=ixvm
87  ! Define absolute gap length
88      if (ix.eq.0.or.ix.eq.nxn(l)-1) jumpx=1
89      jux=jumpx
90  ! Precalculate pressure values for efficiency
91      do kl=1,nuvz
92        ppml(kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
93      end do
94  !
95  ! Loop over the vertical
96  !***********************
97
98      do kl=1,nuvz
99        ppmk=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
100        theta=tthn(ix,jy,kl,n,l)*(100000./ppmk)**kappa
101        klvrp=kl+1
102        klvrm=kl-1
103        klpt=kl
104  ! If top or bottom level, dthetadp is evaluated between the current
105  ! level and the level inside, otherwise between level+1 and level-1
106  !
107        if (klvrp.gt.nuvz) klvrp=nuvz
108        if (klvrm.lt.1) klvrm=1
109        ppmk=akz(klvrp)+bkz(klvrp)*psn(ix,jy,1,n,l)
110        thetap=tthn(ix,jy,klvrp,n,l)*(100000./ppmk)**kappa
111        ppmk=akz(klvrm)+bkz(klvrm)*psn(ix,jy,1,n,l)
112        thetam=tthn(ix,jy,klvrm,n,l)*(100000./ppmk)**kappa
113        dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))
114
115  ! Compute vertical position at pot. temperature surface on subgrid
116  ! and the wind at that position
117  !*****************************************************************
118  ! a) in x direction
119        ii=0
120        do i=ixvm,ixvp,jumpx
121          ivr=i
122          ii=ii+1
123  ! Search adjacent levels for current theta value
124  ! Spiral out from current level for efficiency
125          kup=klpt-1
126          kdn=klpt
127          kch=0
12840        continue
129  ! Upward branch
130          kup=kup+1
131          if (kch.ge.nlck) goto 21     ! No more levels to check,
132  !                                       ! and no values found
133          if (kup.ge.nuvz) goto 41
134          kch=kch+1
135          k=kup
136          ppmk=akz(k)+bkz(k)*psn(ivr,jy,1,n,l)
137          thdn=tthn(ivr,jy,k,n,l)*(100000./ppmk)**kappa
138          ppmk=akz(k+1)+bkz(k+1)*psn(ivr,jy,1,n,l)
139          thup=tthn(ivr,jy,k+1,n,l)*(100000./ppmk)**kappa
140
141      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
142           ((thdn.le.theta).and.(thup.ge.theta))) then
143              dt1=abs(theta-thdn)
144              dt2=abs(theta-thup)
145              dt=dt1+dt2
146              if (dt.lt.eps) then   ! Avoid division by zero error
147                dt1=0.5             ! G.W., 10.4.1996
148                dt2=0.5
149                dt=1.0
150              endif
151    vx(ii)=(vvhn(ivr,jy,k,l)*dt2+vvhn(ivr,jy,k+1,l)*dt1)/dt
152              goto 20
153            endif
15441        continue
155  ! Downward branch
156          kdn=kdn-1
157          if (kdn.lt.1) goto 40
158          kch=kch+1
159          k=kdn
160          ppmk=akz(k)+bkz(k)*psn(ivr,jy,1,n,l)
161          thdn=tthn(ivr,jy,k,n,l)*(100000./ppmk)**kappa
162          ppmk=akz(k+1)+bkz(k+1)*psn(ivr,jy,1,n,l)
163          thup=tthn(ivr,jy,k+1,n,l)*(100000./ppmk)**kappa
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)=(vvhn(ivr,jy,k,l)*dt2+vvhn(ivr,jy,k+1,l)*dt1)/dt
175              goto 20
176            endif
177            goto 40
178  ! This section used when no values were found
17921      continue
180  ! Must use vv at current level and long. jux becomes smaller by 1
181        vx(ii)=vvhn(ix,jy,kl,l)
182        jux=jux-1
183  ! Otherwise OK
18420        continue
185        end do
186      if (jux.gt.0) then
187      dvdx=(vx(2)-vx(1))/real(jux)/(dxn(l)*pi/180.)
188      else
189      dvdx=vvhn(ivrp,jy,kl,l)-vvhn(ivrm,jy,kl,l)
190      dvdx=dvdx/real(jumpx)/(dxn(l)*pi/180.)
191  ! Only happens if no equivalent theta value
192  ! can be found on either side, hence must use values
193  ! from either side, same pressure level.
194      end if
195
196  ! b) in y direction
197
198        jj=0
199        do j=jyvm,jyvp,jumpy
200          jj=jj+1
201  ! Search adjacent levels for current theta value
202  ! Spiral out from current level for efficiency
203          kup=klpt-1
204          kdn=klpt
205          kch=0
20670        continue
207  ! Upward branch
208          kup=kup+1
209          if (kch.ge.nlck) goto 51     ! No more levels to check,
210  !                                     ! and no values found
211          if (kup.ge.nuvz) goto 71
212          kch=kch+1
213          k=kup
214          ppmk=akz(k)+bkz(k)*psn(ix,j,1,n,l)
215          thdn=tthn(ix,j,k,n,l)*(100000./ppmk)**kappa
216          ppmk=akz(k+1)+bkz(k+1)*psn(ix,j,1,n,l)
217          thup=tthn(ix,j,k+1,n,l)*(100000./ppmk)**kappa
218      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
219           ((thdn.le.theta).and.(thup.ge.theta))) then
220              dt1=abs(theta-thdn)
221              dt2=abs(theta-thup)
222              dt=dt1+dt2
223              if (dt.lt.eps) then   ! Avoid division by zero error
224                dt1=0.5             ! G.W., 10.4.1996
225                dt2=0.5
226                dt=1.0
227              endif
228        uy(jj)=(uuhn(ix,j,k,l)*dt2+uuhn(ix,j,k+1,l)*dt1)/dt
229              goto 50
230            endif
23171        continue
232  ! Downward branch
233          kdn=kdn-1
234          if (kdn.lt.1) goto 70
235          kch=kch+1
236          k=kdn
237          ppmk=akz(k)+bkz(k)*psn(ix,j,1,n,l)
238          thdn=tthn(ix,j,k,n,l)*(100000./ppmk)**kappa
239          ppmk=akz(k+1)+bkz(k+1)*psn(ix,j,1,n,l)
240          thup=tthn(ix,j,k+1,n,l)*(100000./ppmk)**kappa
241      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
242           ((thdn.le.theta).and.(thup.ge.theta))) then
243              dt1=abs(theta-thdn)
244              dt2=abs(theta-thup)
245              dt=dt1+dt2
246              if (dt.lt.eps) then   ! Avoid division by zero error
247                dt1=0.5             ! G.W., 10.4.1996
248                dt2=0.5
249                dt=1.0
250              endif
251        uy(jj)=(uuhn(ix,j,k,l)*dt2+uuhn(ix,j,k+1,l)*dt1)/dt
252              goto 50
253            endif
254            goto 70
255  ! This section used when no values were found
25651      continue
257  ! Must use uu at current level and lat. juy becomes smaller by 1
258        uy(jj)=uuhn(ix,jy,kl,l)
259        juy=juy-1
260  ! Otherwise OK
26150        continue
262        end do
263      if (juy.gt.0) then
264      dudy=(uy(2)-uy(1))/real(juy)/(dyn(l)*pi/180.)
265      else
266      dudy=uuhn(ix,jyvp,kl,l)-uuhn(ix,jyvm,kl,l)
267      dudy=dudy/real(jumpy)/(dyn(l)*pi/180.)
268      end if
269  !
270      pvhn(ix,jy,kl,l)=dthetadp*(f+(dvdx/cosphi-dudy &
271           +uuhn(ix,jy,kl,l)*tanphi)/r_earth)*(-1.e6)*9.81
272
273  !
274  ! Resest jux and juy
275      jux=jumpx
276      juy=jumpy
277      end do
278    end do
279  end do
280  !
281end subroutine calcpv_nests
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG