source: flexpart.git/src/calcpar_nests.f90 @ 02095e3

devrelease-10univie
Last change on this file since 02095e3 was 6ecb30a, checked in by Espen Sollum ATMOS <eso@…>, 2 years ago

Merged changes from CTBTO project

  • Property mode set to 100644
File size: 9.9 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 calcpar_nests(n,uuhn,vvhn,pvhn,metdata_format)
23  !                         i  i    i    o
24  !*****************************************************************************
25  !                                                                            *
26  !     Computation of several boundary layer parameters needed for the        *
27  !     dispersion calculation and calculation of dry deposition velocities.   *
28  !     All parameters are calculated over the entire grid.                    *
29  !     This routine is similar to calcpar, but is used for the nested grids.  *
30  !                                                                            *
31  !     Author: A. Stohl                                                       *
32  !                                                                            *
33  !     8 February 1999                                                        *
34  !                                                                            *
35  ! ------------------------------------------------------------------         *
36  !     Petra Seibert, Feb 2000:                                               *
37  !     convection scheme:                                                     *
38  !     new variables in call to richardson                                    *
39  !                                                                            *
40  !*****************************************************************************
41  !  Changes, Bernd C. Krueger, Feb. 2001:                                     *
42  !   Variables tth and qvh (on eta coordinates) in common block               *
43  !                                                                            *
44  !   Unified ECMWF and GFS builds                                             *
45  !   Marian Harustak, 12.5.2017                                               *
46  !     - Added passing of metdata_format as it was needed by called routines  *
47  !*****************************************************************************
48  !                                                                            *
49  ! Variables:                                                                 *
50  ! n                  temporal index for meteorological fields (1 to 3)       *
51  ! metdata_format     format of metdata (ecmwf/gfs)                           *
52  !                                                                            *
53  ! Constants:                                                                 *
54  !                                                                            *
55  !                                                                            *
56  ! Functions:                                                                 *
57  ! scalev             computation of ustar                                    *
58  ! obukhov            computatio of Obukhov length                            *
59  !                                                                            *
60  !*****************************************************************************
61
62  use par_mod
63  use com_mod
64
65  implicit none
66
67  integer :: metdata_format
68  integer :: n,ix,jy,i,l,kz,lz,kzmin
69  real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus,dummyakzllev
70  real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
71  real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
72  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
73  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
74  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
75  real,parameter :: const=r_air/ga
76
77
78  ! Loop over all nests
79  !********************
80
81  do l=1,numbnests
82
83  ! Loop over entire grid
84  !**********************
85
86  do jy=0,nyn(l)-1
87
88  ! Set minimum height for tropopause
89  !**********************************
90
91    ylat=ylat0n(l)+real(jy)*dyn(l)
92    if ((ylat.ge.-20.).and.(ylat.le.20.)) then
93      altmin = 5000.
94    else
95      if ((ylat.gt.20.).and.(ylat.lt.40.)) then
96        altmin=2500.+(40.-ylat)*125.
97      else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
98        altmin=2500.+(40.+ylat)*125.
99      else
100        altmin=2500.
101      endif
102    endif
103
104    do ix=0,nxn(l)-1
105
106  ! 1) Calculation of friction velocity
107  !************************************
108
109      ustarn(ix,jy,1,n,l)=scalev(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
110           td2n(ix,jy,1,n,l),surfstrn(ix,jy,1,n,l))
111      if (ustarn(ix,jy,1,n,l).le.1.e-8) ustarn(ix,jy,1,n,l)=1.e-8
112
113  ! 2) Calculation of inverse Obukhov length scale
114  !***********************************************
115
116      ol=obukhov(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
117           td2n(ix,jy,1,n,l),tthn(ix,jy,2,n,l),ustarn(ix,jy,1,n,l), &
118           sshfn(ix,jy,1,n,l),akm,bkm,dummyakzllev,metdata_format)
119      if (ol.ne.0.) then
120        olin(ix,jy,1,n,l)=1./ol
121      else
122        olin(ix,jy,1,n,l)=99999.
123      endif
124
125
126  ! 3) Calculation of convective velocity scale and mixing height
127  !**************************************************************
128
129      do i=1,nuvz
130        ulev(i)=uuhn(ix,jy,i,l)
131        vlev(i)=vvhn(ix,jy,i,l)
132        ttlev(i)=tthn(ix,jy,i,n,l)
133        qvlev(i)=qvhn(ix,jy,i,n,l)
134      end do
135
136      call richardson(psn(ix,jy,1,n,l),ustarn(ix,jy,1,n,l),ttlev, &
137           qvlev,ulev,vlev,nuvz,akz,bkz,sshfn(ix,jy,1,n,l), &
138           tt2n(ix,jy,1,n,l),td2n(ix,jy,1,n,l),hmixn(ix,jy,1,n,l), &
139           wstarn(ix,jy,1,n,l),hmixplus,metdata_format)
140
141      if(lsubgrid.eq.1) then
142        subsceff=min(excessoron(ix,jy,l),hmixplus)
143      else
144        subsceff=0.0
145      endif
146  !
147  ! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
148  !
149      hmixn(ix,jy,1,n,l)=hmixn(ix,jy,1,n,l)+subsceff
150      hmixn(ix,jy,1,n,l)=max(hmixmin,hmixn(ix,jy,1,n,l)) ! minim PBL height
151      hmixn(ix,jy,1,n,l)=min(hmixmax,hmixn(ix,jy,1,n,l)) ! maxim PBL height
152
153
154  ! 4) Calculation of dry deposition velocities
155  !********************************************
156
157      if (DRYDEP) then
158        ! z0(4)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
159        ! z0(9)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
160        z0(7)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
161
162  ! Calculate relative humidity at surface
163  !***************************************
164        rh=ew(td2n(ix,jy,1,n,l))/ew(tt2n(ix,jy,1,n,l))
165
166        call getvdep_nests(n,ix,jy,ustarn(ix,jy,1,n,l), &
167             tt2n(ix,jy,1,n,l),psn(ix,jy,1,n,l),1./olin(ix,jy,1,n,l), &
168             ssrn(ix,jy,1,n,l),rh,lsprecn(ix,jy,1,n,l)+ &
169             convprecn(ix,jy,1,n,l),sdn(ix,jy,1,n,l),vd,l)
170
171        do i=1,nspec
172          vdepn(ix,jy,i,n,l)=vd(i)
173        end do
174
175      endif
176
177  !******************************************************
178  ! Calculate height of thermal tropopause (Hoinka, 1997)
179  !******************************************************
180
181  ! 1) Calculate altitudes of ECMWF model levels
182  !*********************************************
183
184      tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l))/ &
185           psn(ix,jy,1,n,l))
186      pold=psn(ix,jy,1,n,l)
187      zold=0.
188      do kz=2,nuvz
189        pint=akz(kz)+bkz(kz)*psn(ix,jy,1,n,l)  ! pressure on model layers
190        tv=tthn(ix,jy,kz,n,l)*(1.+0.608*qvhn(ix,jy,kz,n,l))
191
192        if (abs(tv-tvold).gt.0.2) then
193         zlev(kz)=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
194        else
195          zlev(kz)=zold+const*log(pold/pint)*tv
196        endif
197        tvold=tv
198        pold=pint
199        zold=zlev(kz)
200      end do
201
202  ! 2) Define a minimum level kzmin, from which upward the tropopause is
203  !    searched for. This is to avoid inversions in the lower troposphere
204  !    to be identified as the tropopause
205  !************************************************************************
206
207      do kz=1,nuvz
208        if (zlev(kz).ge.altmin) then
209          kzmin=kz
210          goto 45
211        endif
212      end do
21345    continue
214
215  ! 3) Search for first stable layer above minimum height that fulfills the
216  !    thermal tropopause criterion
217  !************************************************************************
218
219      do kz=kzmin,nuvz
220        do lz=kz+1,nuvz
221          if ((zlev(lz)-zlev(kz)).gt.2000.) then
222            if (((tthn(ix,jy,kz,n,l)-tthn(ix,jy,lz,n,l))/ &
223                 (zlev(lz)-zlev(kz))).lt.0.002) then
224              tropopausen(ix,jy,1,n,l)=zlev(kz)
225              goto 51
226            endif
227            goto 50
228          endif
229        end do
23050      continue
231      end do
23251    continue
233
234
235    end do
236  end do
237
238  ! Calculation of potential vorticity on 3-d grid
239  !***********************************************
240
241  call calcpv_nests(l,n,uuhn,vvhn,pvhn)
242
243  end do
244
245
246end subroutine calcpar_nests
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG