source: branches/flexpart91_hasod/src_parallel/calcpar.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: 10.0 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(n,uuh,vvh,pvh)
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  !                                                                            *
30  !     Author: A. Stohl                                                       *
31  !                                                                            *
32  !     21 May 1995                                                            *
33  !                                                                            *
34  ! ------------------------------------------------------------------         *
35  !     Petra Seibert, Feb 2000:                                               *
36  !     convection scheme:                                                     *
37  !     new variables in call to richardson                                    *
38  !                                                                            *
39  !*****************************************************************************
40  !  Changes, Bernd C. Krueger, Feb. 2001:
41  !   Variables tth and qvh (on eta coordinates) in common block
42  !*****************************************************************************
43  !                                                                            *
44  ! Variables:                                                                 *
45  ! n                  temporal index for meteorological fields (1 to 3)       *
46  !                                                                            *
47  ! Constants:                                                                 *
48  !                                                                            *
49  !                                                                            *
50  ! Functions:                                                                 *
51  ! scalev             computation of ustar                                    *
52  ! obukhov            computatio of Obukhov length                            *
53  !                                                                            *
54  !*****************************************************************************
55
56  use par_mod
57  use com_mod
58
59  implicit none
60
61  integer :: n,ix,jy,i,kz,lz,kzmin
62  real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
63  real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
64  real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
65  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
66  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
67  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
68  real,parameter :: const=r_air/ga
69
70  !write(*,*) 'in calcpar writting snowheight'
71  !***********************************
72  ! for test: write out snow depths
73
74  ! open(4,file='slandusetest',form='formatted')
75  ! do 5 ix=0,nxmin1
76  !5       write (4,*) (sd(ix,jy,1,n),jy=0,nymin1)
77  !  close(4)
78
79
80  ! Loop over entire grid
81  !**********************
82
83!$OMP PARALLEL DEFAULT(none) &
84!$OMP PRIVATE(jy, ix, ylat, altmin, ol, ulev, i, vlev, ttlev, qvlev, &
85!$OMP   hmixplus, kz, tvold, pold, zold, tv, pint, zlev, subsceff, rh, vd, &
86!$OMP   kzmin) &
87!$OMP SHARED(nymin1, nxmin1, dx, dy, ylat0, n, ustar, surfstr, td2, tt2, ps, &
88!$OMP   akm, bkm, sshf, oli, nuvz, tth, uuh, vvh, qvh, hmix, wstar, akz, bkz, &
89!$OMP   lsubgrid, excessoro, DRYDEP, lsprec, convprec, vdep, sd, ssr, tropopause, &
90!$OMP   nspec) &
91!$OMP COPYIN(z0)
92
93#if (defined STATIC_SCHED)
94!$OMP DO SCHEDULE(static)
95#else
96!$OMP DO SCHEDULE(dynamic, max(1,nymin1/10))
97#endif
98
99  latloop : do jy=0,nymin1
100
101  ! Set minimum height for tropopause
102  !**********************************
103
104    ylat=ylat0+real(jy)*dy
105    if ((ylat.ge.-20.).and.(ylat.le.20.)) then
106      altmin = 5000.
107    else
108      if ((ylat.gt.20.).and.(ylat.lt.40.)) then
109        altmin=2500.+(40.-ylat)*125.
110      else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
111        altmin=2500.+(40.+ylat)*125.
112      else
113        altmin=2500.
114      endif
115    endif
116
117    lonloop : do ix=0,nxmin1
118
119  ! 1) Calculation of friction velocity
120  !************************************
121
122      ustar(ix,jy,1,n)=scalev(ps(ix,jy,1,n),tt2(ix,jy,1,n), &
123           td2(ix,jy,1,n),surfstr(ix,jy,1,n))
124      if (ustar(ix,jy,1,n).le.1.e-8) ustar(ix,jy,1,n)=1.e-8
125
126  ! 2) Calculation of inverse Obukhov length scale
127  !***********************************************
128
129      ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
130           tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm)
131      if (ol.ne.0.) then
132        oli(ix,jy,1,n)=1./ol
133      else
134        oli(ix,jy,1,n)=99999.
135      endif
136
137
138  ! 3) Calculation of convective velocity scale and mixing height
139  !**************************************************************
140
141      do i=1,nuvz
142        ulev(i)=uuh(ix,jy,i)
143        vlev(i)=vvh(ix,jy,i)
144        ttlev(i)=tth(ix,jy,i,n)
145        qvlev(i)=qvh(ix,jy,i,n)
146      end do
147
148      call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
149           ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), &
150           td2(ix,jy,1,n),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus)
151
152      if(lsubgrid.eq.1) then
153        subsceff=min(excessoro(ix,jy),hmixplus)
154      else
155        subsceff=0.0
156      endif
157  !
158  ! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
159  !
160      hmix(ix,jy,1,n)=hmix(ix,jy,1,n)+subsceff
161      hmix(ix,jy,1,n)=max(hmixmin,hmix(ix,jy,1,n)) ! set minimum PBL height
162      hmix(ix,jy,1,n)=min(hmixmax,hmix(ix,jy,1,n)) ! set maximum PBL height
163
164  ! 4) Calculation of dry deposition velocities
165  !********************************************
166
167      if (DRYDEP) then
168  ! Sabine Eckhardt, Dec 06: use new index for z0 for water depending on
169  ! windspeed
170        z0(7)=0.016*ustar(ix,jy,1,n)*ustar(ix,jy,1,n)/ga
171
172  ! Calculate relative humidity at surface
173  !***************************************
174        rh=ew(td2(ix,jy,1,n))/ew(tt2(ix,jy,1,n))
175
176        call getvdep(n,ix,jy,ustar(ix,jy,1,n), &
177             tt2(ix,jy,1,n),ps(ix,jy,1,n),1./oli(ix,jy,1,n), &
178             ssr(ix,jy,1,n),rh,lsprec(ix,jy,1,n)+convprec(ix,jy,1,n), &
179             sd(ix,jy,1,n),vd)
180        do i=1,nspec
181          vdep(ix,jy,i,n)=vd(i)
182        end do
183
184      endif
185
186  !******************************************************
187  ! Calculate height of thermal tropopause (Hoinka, 1997)
188  !******************************************************
189
190  ! 1) Calculate altitudes of ECMWF model levels
191  !*********************************************
192
193      tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
194           ps(ix,jy,1,n))
195      pold=ps(ix,jy,1,n)
196      zold=0.
197      do kz=2,nuvz
198        pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)  ! pressure on model layers
199        tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
200
201        if (abs(tv-tvold).gt.0.2) then
202         zlev(kz)=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
203        else
204          zlev(kz)=zold+const*log(pold/pint)*tv
205        endif
206        tvold=tv
207        pold=pint
208        zold=zlev(kz)
209      end do
210
211  ! 2) Define a minimum level kzmin, from which upward the tropopause is
212  !    searched for. This is to avoid inversions in the lower troposphere
213  !    to be identified as the tropopause
214  !************************************************************************
215
216      do kz=1,nuvz
217        if (zlev(kz).ge.altmin) then
218          kzmin=kz
219          goto 45
220        endif
221      end do
22245    continue
223
224  ! 3) Search for first stable layer above minimum height that fulfills the
225  !    thermal tropopause criterion
226  !************************************************************************
227
228      do kz=kzmin,nuvz
229        do lz=kz+1,nuvz
230          if ((zlev(lz)-zlev(kz)).gt.2000.) then
231            if (((tth(ix,jy,kz,n)-tth(ix,jy,lz,n))/ &
232                 (zlev(lz)-zlev(kz))).lt.0.002) then
233              tropopause(ix,jy,1,n)=zlev(kz)
234              goto 51
235            endif
236            goto 50
237          endif
238        end do
23950      continue
240      end do
24151    continue
242
243
244    end do lonloop
245
246  end do latloop
247
248!$OMP END DO
249!$OMP END PARALLEL
250
251!  call dump_field("vdep", vdep(0:nxmin1,0:nymin1,1:nspec,n), nx, ny, nspec)
252!  call dump_field("oli", oli(0:nxmin1,0:nymin1,1,n), nx, ny, 1)
253!  call dump_field("hmix", hmix(0:nxmin1,0:nymin1,1,n), nx, ny, 1)
254!  call dump_field("xlanduse", xlanduse(0:nxmin1,0:nymin1,1:numclass), nx, ny, numclass)
255!   stop
256
257  ! Calculation of potential vorticity on 3-d grid
258  !***********************************************
259
260  call calcpv(n,uuh,vvh,pvh)
261
262
263end subroutine calcpar
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG