source: branches/jerome/src_flexwrf_v3.1/calcpar_nests.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 10 years ago

sources for flexwrf v3.1

File size: 12.6 KB
Line 
1!***********************************************************************
2!* Copyright 2012,2013                                                *
3!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine,       *
4!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso,          *
6!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann,   *
7!*                                                                     *
8!* This file is part of FLEXPART WRF                                   *
9!*                                                                     *
10!* FLEXPART is free software: you can redistribute it and/or modify    *
11!* it under the terms of the GNU General Public License as published by*
12!* the Free Software Foundation, either version 3 of the License, or   *
13!* (at your option) any later version.                                 *
14!*                                                                     *
15!* FLEXPART is distributed in the hope that it will be useful,         *
16!* but WITHOUT ANY WARRANTY; without even the implied warranty of      *
17!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
18!* GNU General Public License for more details.                        *
19!*                                                                     *
20!* You should have received a copy of the GNU General Public License   *
21!* along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
22!***********************************************************************
23      subroutine calcpar_nests(n,uuhn,vvhn,pvhn)
24!                              i  i    i    o
25!*******************************************************************************
26!                                                                              *
27!     Computation of several boundary layer parameters needed for the          *
28!     dispersion calculation and calculation of dry deposition velocities.     *
29!     All parameters are calculated over the entire grid.                      *
30!     This routine is similar to calcpar, but is used for the nested grids.    *
31!                                                                              *
32!     Note:  This is the FLEXPART_WRF version of subroutine calcpar.           *
33!                                                                              *
34!     Author: A. Stohl                                                         *
35!                                                                              *
36!     8 February 1999                                                          *
37!                                                                              *
38! ------------------------------------------------------------------           *
39!     Petra Seibert, Feb 2000:                                                 *
40!     convection scheme:                                                       *
41!     new variables in call to richardson                                      *
42!                                                                              *
43!     Changes, Bernd C. Krueger, Feb. 2001:                                    *
44!        Variables tth and qvh (on eta coordinates) in common block            *
45!                                                                              *
46!     14 Nov 2005 - R. Easter -                                                *
47!          use xyindex_to_ll_wrf to get latitude                               *
48!          limit ustar to < 5.0 m/s                                            *
49!          added ierr in call to richardson                                    *
50!          use pph for calculating zlev                                        *
51!          pass level-2 pph directly to obukhov                                *
52!     15 Nov 2005 - R. Easter - pass pplev to richardson instead of akz,bkz    *
53!    Jul 2012: J. Brioude: coded in fortran90 and parallelized                 *
54!*******************************************************************************
55!                                                                              *
56! Variables:                                                                   *
57! n                  temporal index for meteorological fields (1 to 3)         *
58!                                                                              *
59! Constants:                                                                   *
60!                                                                              *
61!                                                                              *
62! Functions:                                                                   *
63! scalev             computation of ustar                                      *
64! obukhov            computatio of Obukhov length                              *
65!                                                                              *
66!*******************************************************************************
67
68  use par_mod
69  use com_mod
70
71  implicit none
72
73  integer :: n,ix,jy,i,l,kz,lz,kzmin,ierr
74  real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
75  real :: pplev(nuvzmax),xlon
76  real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
77  real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
78  real(kind=4) :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
79  real(kind=4) :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
80  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
81  real,parameter :: const=r_air/ga
82
83
84! Loop over all nests
85!********************
86!     ientry = ientry + 1
87
88      do l=1,numbnests
89
90! Loop over entire grid
91!**********************
92!$OMP PARALLEL DEFAULT(SHARED) &
93!$OMP PRIVATE(i,ix,jy,kz,lz,kzmin,tvold,pold,zold,zlev,tv,pint, &
94!$OMP rh,ierr,subsceff,ulev,vlev,pplev,ttlev,qvlev,ol,altmin,xlon,ylat )
95!$OMP DO
96      do jy=0,nyn(l)-1
97
98        do ix=0,nxn(l)-1
99
100! Set minimum height for tropopause
101!**********************************
102
103! FLEXPART_WRF - use this routine to get lat,lon
104!       ylat=ylat0n(l)+real(jy)*dyn(l)
105        call xyindex_to_ll_wrf( l, real(ix), real(jy), xlon, ylat )
106
107!       if ( ((ix.eq.0) .or. (ix.eq.(nxn(l)-1)) .or.
108!    &                       (ix.eq.(nxn(l)-1)/2)) .and.
109!    &       ((jy.eq.0) .or. (jy.eq.(nyn(l)-1)) .or.
110!    &                       (jy.eq.(nyn(l)-1)/2)) ) then
111!           if (ientry .eq. 1) then
112!               write(*,'(a,3i4,2f12.5)')
113!    &              'calcpar_nests l,i,j, xlon,ylat',
114!    &              l, ix, jy, xlon, ylat
115!               write(*,'(a,12x,2f12.5)')
116!    &              '                     dlon,dlat',
117!    &              (xlon-xlon2dn(ix,jy,l)), (ylat-ylat2dn(ix,jy,l))
118!           end if
119!       end if
120
121        if ((ylat.ge.-20.).and.(ylat.le.20.)) then
122          altmin = 5000.
123        else
124          if ((ylat.gt.20.).and.(ylat.lt.40.)) then
125            altmin=2500.+(40.-ylat)*125.
126          else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
127            altmin=2500.+(40.+ylat)*125.
128          else
129            altmin=2500.
130          endif
131        endif
132
133! 1) Calculation of friction velocity
134!************************************
135          if ( (.not.strswitch)) then
136
137          ustarn(ix,jy,1,n,l)=scalev(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
138          td2n(ix,jy,1,n,l),surfstrn(ix,jy,1,n,l))
139          endif
140! FLEXPART_WRF - limit ustar
141          if (ustarn(ix,jy,1,n,l).le.1.e-8) ustarn(ix,jy,1,n,l)=1.e-8
142          if (ustarn(ix,jy,1,n,l).ge.5.0)   ustarn(ix,jy,1,n,l)=5.0
143
144! 2) Calculation of inverse Obukhov length scale
145!***********************************************
146
147! FLEXPART_WRF - pass k=2 pressure directly
148!         ol=obukhov(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l),
149!    +    td2n(ix,jy,1,n,l),tthn(ix,jy,2,n,l),ustarn(ix,jy,1,n,l),
150!    +    sshfn(ix,jy,1,n,l),akm,bkm)
151          ol=obukhov(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
152          td2n(ix,jy,1,n,l),tthn(ix,jy,2,n,l),ustarn(ix,jy,1,n,l), &
153          sshfn(ix,jy,1,n,l),pphn(ix,jy,2,n,l) )
154          if (ol.ne.0.) then
155            olin(ix,jy,1,n,l)=1./ol
156          else
157            olin(ix,jy,1,n,l)=99999.
158          endif
159
160
161! 3) Calculation of convective velocity scale and mixing height
162!**************************************************************
163
164          do i=1,nuvz
165            ulev(i) =uuhn(ix,jy,i,l)
166            vlev(i) =vvhn(ix,jy,i,l)
167            pplev(i)=pphn(ix,jy,i,n,l)
168            ttlev(i)=tthn(ix,jy,i,n,l)
169            qvlev(i)=qvhn(ix,jy,i,n,l)
170      end do
171
172! FLEXPART_WRF - use  & check ierr argument
173! FLEXPART_WRF - pass pplev instead of akz,bkz
174!         call richardson(psn(ix,jy,1,n,l),ustarn(ix,jy,1,n,l),ttlev,
175!    +    qvlev,ulev,vlev,nuvz,akz,bkz,sshfn(ix,jy,1,n,l),
176!    +    tt2n(ix,jy,1,n,l),td2n(ix,jy,1,n,l),hmixn(ix,jy,1,n,l),
177!    +    wstarn(ix,jy,1,n,l),hmixplus)
178          call richardson(psn(ix,jy,1,n,l),ustarn(ix,jy,1,n,l),ttlev, &
179          qvlev,ulev,vlev,nuvz,  pplev,sshfn(ix,jy,1,n,l), &
180          tt2n(ix,jy,1,n,l),td2n(ix,jy,1,n,l),hmixn(ix,jy,1,n,l), &
181          wstarn(ix,jy,1,n,l),hmixplus,ierr,sfc_option)
182
183          if (ierr .gt. 0) then
184              write(*,9500) 'warning', l, ix, jy
185          else if (ierr .lt. 0) then
186              write(*,9500) 'failure', l, ix, jy
187              stop
188          end if
1899500      format( 'calcpar_nests - richardson ', a, ' - l,ix,jy=', 3i5 )
190
191          if(lsubgrid.eq.1) then
192            subsceff=min(excessoron(ix,jy,l),hmixplus)
193          else
194            subsceff=0
195          endif
196!
197! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
198!
199          hmixn(ix,jy,1,n,l)=hmixn(ix,jy,1,n,l)+subsceff
200          hmixn(ix,jy,1,n,l)=max(hmixmin,hmixn(ix,jy,1,n,l)) ! minim PBL height
201          hmixn(ix,jy,1,n,l)=min(hmixmax,hmixn(ix,jy,1,n,l)) ! maxim PBL height
202
203
204! 4) Calculation of dry deposition velocities
205!********************************************
206
207          if (DRYDEP) then
208            z0(4)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
209            z0(9)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
210
211! Calculate relative humidity at surface
212!***************************************
213            rh=ew(td2n(ix,jy,1,n,l))/ew(tt2n(ix,jy,1,n,l))
214
215            call getvdep_nests(n,ix,jy,ustarn(ix,jy,1,n,l), &
216            tt2n(ix,jy,1,n,l),psn(ix,jy,1,n,l),1./olin(ix,jy,1,n,l), &
217            ssrn(ix,jy,1,n,l),rh,lsprecn(ix,jy,1,n,l)+ &
218            convprecn(ix,jy,1,n,l),vd,l)
219
220            do i=1,nspec
221              vdepn(ix,jy,i,n,l)=vd(i)
222            enddo
223          endif
224
225!******************************************************
226! Calculate height of thermal tropopause (Hoinka, 1997)
227!******************************************************
228
229! 1) Calculate altitudes of ECMWF model levels
230!*********************************************
231
232          tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l))/ &
233                                         psn(ix,jy,1,n,l))
234          pold=psn(ix,jy,1,n,l)
235          zold=0.
236! FLEXPART_WRF - set zlev(1)
237          zlev(1)=zold
238      do kz=2,nuvz
239        pint=akz(kz)+bkz(kz)*psn(ix,jy,1,n,l)  ! pressure on model layers
240        tv=tthn(ix,jy,kz,n,l)*(1.+0.608*qvhn(ix,jy,kz,n,l))
241
242        if (abs(tv-tvold).gt.0.2) then
243         zlev(kz)=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
244        else
245          zlev(kz)=zold+const*log(pold/pint)*tv
246        endif
247        tvold=tv
248        pold=pint
249        zold=zlev(kz)
250      end do
251
252  ! 2) Define a minimum level kzmin, from which upward the tropopause is
253  !    searched for. This is to avoid inversions in the lower troposphere
254  !    to be identified as the tropopause
255  !************************************************************************
256
257      do kz=1,nuvz
258        if (zlev(kz).ge.altmin) then
259          kzmin=kz
260          goto 45
261        endif
262      end do
26345    continue
264
265! 3) Search for first stable layer above minimum height that fulfills the
266!    thermal tropopause criterion
267!************************************************************************
268
269      do kz=kzmin,nuvz
270        do lz=kz+1,nuvz
271          if ((zlev(lz)-zlev(kz)).gt.2000.) then
272            if (((tthn(ix,jy,kz,n,l)-tthn(ix,jy,lz,n,l))/ &
273                 (zlev(lz)-zlev(kz))).lt.0.002) then
274              tropopausen(ix,jy,1,n,l)=zlev(kz)
275              goto 51
276            endif
277            goto 50
278          endif
279        end do
28050      continue
281      end do
28251    continue
283
284
285    end do
286  end do
287!$OMP END DO
288!$OMP END PARALLEL
289
290! Calculation of potential vorticity on 3-d grid, if plume trajectory mode is used
291!*********************************************************************************
292
293        if ((iout.eq.4).or.(iout.eq.5)) then
294          call calcpv_nests(l,n,uuhn,vvhn,pvhn)
295        endif
296
297        enddo
298
299
300end subroutine calcpar_nests
301
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG