source: trunk/src/verttransform.f90 @ 28

Last change on this file since 28 was 24, checked in by igpis, 10 years ago

version 9.2 beta. Changes from HH, AST, MC, NIK, IP. Changes in vert transform. New SPECIES input includes scavenging coefficients

File size: 17.3 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 verttransform(n,uuh,vvh,wwh,pvh)
23  !                      i  i   i   i   i
24  !*****************************************************************************
25  !                                                                            *
26  !     This subroutine transforms temperature, dew point temperature and      *
27  !     wind components from eta to meter coordinates.                         *
28  !     The vertical wind component is transformed from Pa/s to m/s using      *
29  !     the conversion factor pinmconv.                                        *
30  !     In addition, this routine calculates vertical density gradients        *
31  !     needed for the parameterization of the turbulent velocities.           *
32  !                                                                            *
33  !     Author: A. Stohl, G. Wotawa                                            *
34  !                                                                            *
35  !     12 August 1996                                                         *
36  !     Update: 16 January 1998                                                *
37  !                                                                            *
38  !     Major update: 17 February 1999                                         *
39  !     by G. Wotawa                                                           *
40  !                                                                            *
41  !     - Vertical levels for u, v and w are put together                      *
42  !     - Slope correction for vertical velocity: Modification of calculation  *
43  !       procedure                                                            *
44  !                                                                            *
45  !*****************************************************************************
46  !  Changes, Bernd C. Krueger, Feb. 2001:
47  !   Variables tth and qvh (on eta coordinates) from common block
48  !*****************************************************************************
49  ! Sabine Eckhardt, March 2007
50  ! added the variable cloud for use with scavenging - descr. in com_mod
51  !*****************************************************************************
52  !                                                                            *
53  ! Variables:                                                                 *
54  ! nx,ny,nz                        field dimensions in x,y and z direction    *
55  ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition           *
56  ! uu(0:nxmax,0:nymax,nzmax,2)     wind components in x-direction [m/s]       *
57  ! vv(0:nxmax,0:nymax,nzmax,2)     wind components in y-direction [m/s]       *
58  ! ww(0:nxmax,0:nymax,nzmax,2)     wind components in z-direction [deltaeta/s]*
59  ! tt(0:nxmax,0:nymax,nzmax,2)     temperature [K]                            *
60  ! pv(0:nxmax,0:nymax,nzmax,2)     potential voriticity (pvu)                 *
61  ! ps(0:nxmax,0:nymax,2)           surface pressure [Pa]                      *
62  !                                                                            *
63  !*****************************************************************************
64
65  use par_mod
66  use com_mod
67  use cmapf_mod, only: cc2gll
68
69  implicit none
70
71  integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
72  integer :: rain_cloud_above,kz_inv
73  real :: f_qvsat,pressure
74  real :: rh,lsp,convp
75  real :: rhoh(nuvzmax),pinmconv(nzmax)
76  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
77  real :: xlon,ylat,xlonr,dzdx,dzdy
78  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
79  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
80  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
81  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
82  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
83  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
84  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
85  real,parameter :: const=r_air/ga
86
87  logical :: init = .true.
88
89
90  !*************************************************************************
91  ! If verttransform is called the first time, initialize heights of the   *
92  ! z levels in meter. The heights are the heights of model levels, where  *
93  ! u,v,T and qv are given.                                                *
94  !*************************************************************************
95
96  if (init) then
97
98  ! Search for a point with high surface pressure (i.e. not above significant topography)
99  ! Then, use this point to construct a reference z profile, to be used at all times
100  !**************************************************************************************
101
102    do jy=0,nymin1
103      do ix=0,nxmin1
104        if (ps(ix,jy,1,n).gt.100000.) then
105          ixm=ix
106          jym=jy
107          goto 3
108        endif
109      end do
110    end do
1113   continue
112
113
114    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
115    ps(ixm,jym,1,n))
116    pold=ps(ixm,jym,1,n)
117    height(1)=0.
118
119    do kz=2,nuvz
120      pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
121      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
122
123      if (abs(tv-tvold).gt.0.2) then
124        height(kz)=height(kz-1)+const*log(pold/pint)* &
125        (tv-tvold)/log(tv/tvold)
126      else
127        height(kz)=height(kz-1)+const*log(pold/pint)*tv
128      endif
129
130      tvold=tv
131      pold=pint
132    end do
133
134
135  ! Determine highest levels that can be within PBL
136  !************************************************
137
138    do kz=1,nz
139      if (height(kz).gt.hmixmax) then
140        nmixz=kz
141        goto 9
142      endif
143    end do
1449   continue
145
146  ! Do not repeat initialization of the Cartesian z grid
147  !*****************************************************
148
149    init=.false.
150
151  endif
152
153
154  ! Loop over the whole grid
155  !*************************
156
157  do jy=0,nymin1
158    do ix=0,nxmin1
159      tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ps(ix,jy,1,n))
160      pold=ps(ix,jy,1,n)
161      uvwzlev(ix,jy,1)=0.
162      wzlev(1)=0.
163      rhoh(1)=pold/(r_air*tvold)
164
165
166  ! Compute heights of eta levels
167  !******************************
168
169      do kz=2,nuvz
170        pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
171        tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
172        rhoh(kz)=pint/(r_air*tv)
173
174        if (abs(tv-tvold).gt.0.2) then
175          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
176          (tv-tvold)/log(tv/tvold)
177        else
178          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
179        endif
180
181        tvold=tv
182        pold=pint
183      end do
184
185
186      do kz=2,nwz-1
187        wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2.
188      end do
189      wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
190
191  ! pinmconv=(h2-h1)/(p2-p1)
192
193      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
194      ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
195      (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
196      do kz=2,nz-1
197        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
198        ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
199        (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
200      end do
201      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
202      ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
203      (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
204
205  ! Levels, where u,v,t and q are given
206  !************************************
207
208      uu(ix,jy,1,n)=uuh(ix,jy,1)
209      vv(ix,jy,1,n)=vvh(ix,jy,1)
210      tt(ix,jy,1,n)=tth(ix,jy,1,n)
211      qv(ix,jy,1,n)=qvh(ix,jy,1,n)
212      pv(ix,jy,1,n)=pvh(ix,jy,1)
213      rho(ix,jy,1,n)=rhoh(1)
214      uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
215      vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
216      tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
217      qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
218      pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
219      rho(ix,jy,nz,n)=rhoh(nuvz)
220      kmin=2
221      do iz=2,nz-1
222        do kz=kmin,nuvz
223          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
224            uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
225            vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
226            tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
227            qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
228            pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
229            rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
230            goto 30
231          endif
232          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
233          (height(iz).le.uvwzlev(ix,jy,kz))) then
234           dz1=height(iz)-uvwzlev(ix,jy,kz-1)
235           dz2=uvwzlev(ix,jy,kz)-height(iz)
236           dz=dz1+dz2
237           uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
238           vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
239           tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
240                +tth(ix,jy,kz,n)*dz1)/dz
241           qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
242                +qvh(ix,jy,kz,n)*dz1)/dz
243           pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
244           rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
245           kmin=kz
246           goto 30
247          endif
248        end do
24930      continue
250      end do
251
252
253  ! Levels, where w is given
254  !*************************
255
256      ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
257      ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
258      kmin=2
259      do iz=2,nz
260        do kz=kmin,nwz
261          if ((height(iz).gt.wzlev(kz-1)).and. &
262               (height(iz).le.wzlev(kz))) then
263           dz1=height(iz)-wzlev(kz-1)
264           dz2=wzlev(kz)-height(iz)
265           dz=dz1+dz2
266           ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
267                +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
268           kmin=kz
269           goto 40
270          endif
271        end do
27240      continue
273      end do
274
275  ! Compute density gradients at intermediate levels
276  !*************************************************
277
278      drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
279      (height(2)-height(1))
280      do kz=2,nz-1
281        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
282        (height(kz+1)-height(kz-1))
283      end do
284      drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
285
286    end do
287  end do
288
289
290  !****************************************************************
291  ! Compute slope of eta levels in windward direction and resulting
292  ! vertical wind correction
293  !****************************************************************
294
295  do jy=1,ny-2
296    cosf=cos((real(jy)*dy+ylat0)*pi180)
297    do ix=1,nx-2
298
299      kmin=2
300      do iz=2,nz-1
301
302        ui=uu(ix,jy,iz,n)*dxconst/cosf
303        vi=vv(ix,jy,iz,n)*dyconst
304
305        do kz=kmin,nz
306          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
307          (height(iz).le.uvwzlev(ix,jy,kz))) then
308            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
309            dz2=uvwzlev(ix,jy,kz)-height(iz)
310            dz=dz1+dz2
311            kl=kz-1
312            klp=kz
313            kmin=kz
314            goto 47
315          endif
316        end do
317
31847      ix1=ix-1
319        jy1=jy-1
320        ixp=ix+1
321        jyp=jy+1
322
323        dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
324        dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
325        dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
326
327        dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
328        dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
329        dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
330
331        ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*ui+dzdy*vi)
332
333      end do
334
335    end do
336  end do
337
338
339  ! If north pole is in the domain, calculate wind velocities in polar
340  ! stereographic coordinates
341  !*******************************************************************
342
343  if (nglobal) then
344    do jy=int(switchnorthg)-2,nymin1
345      ylat=ylat0+real(jy)*dy
346      do ix=0,nxmin1
347        xlon=xlon0+real(ix)*dx
348        do iz=1,nz
349          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
350          vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
351        end do
352      end do
353    end do
354
355
356    do iz=1,nz
357
358  ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
359  !
360  !   AMSnauffer Nov 18 2004 Added check for case vv=0
361  !
362      xlon=xlon0+real(nx/2-1)*dx
363      xlonr=xlon*pi/180.
364      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2)
365      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
366        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
367      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
368        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
369        vv(nx/2-1,nymin1,iz,n))-xlonr
370      else
371        ddpol=pi/2-xlonr
372      endif
373      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
374      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
375
376  ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
377      xlon=180.0
378      xlonr=xlon*pi/180.
379      ylat=90.0
380      uuaux=-ffpol*sin(xlonr+ddpol)
381      vvaux=-ffpol*cos(xlonr+ddpol)
382      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
383      jy=nymin1
384      do ix=0,nxmin1
385        uupol(ix,jy,iz,n)=uupolaux
386        vvpol(ix,jy,iz,n)=vvpolaux
387      end do
388    end do
389
390
391  ! Fix: Set W at pole to the zonally averaged W of the next equator-
392  ! ward parallel of latitude
393
394  do iz=1,nz
395      wdummy=0.
396      jy=ny-2
397      do ix=0,nxmin1
398        wdummy=wdummy+ww(ix,jy,iz,n)
399      end do
400      wdummy=wdummy/real(nx)
401      jy=nymin1
402      do ix=0,nxmin1
403        ww(ix,jy,iz,n)=wdummy
404      end do
405  end do
406
407  endif
408
409
410  ! If south pole is in the domain, calculate wind velocities in polar
411  ! stereographic coordinates
412  !*******************************************************************
413
414  if (sglobal) then
415    do jy=0,int(switchsouthg)+3
416      ylat=ylat0+real(jy)*dy
417      do ix=0,nxmin1
418        xlon=xlon0+real(ix)*dx
419        do iz=1,nz
420          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
421          vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
422        end do
423      end do
424    end do
425
426    do iz=1,nz
427
428  ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
429  !
430  !   AMSnauffer Nov 18 2004 Added check for case vv=0
431  !
432      xlon=xlon0+real(nx/2-1)*dx
433      xlonr=xlon*pi/180.
434      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2)
435      if (vv(nx/2-1,0,iz,n).lt.0.) then
436        ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
437      else if (vv(nx/2-1,0,iz,n).gt.0.) then
438        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
439      else
440        ddpol=pi/2-xlonr
441      endif
442      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
443      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
444
445  ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
446      xlon=180.0
447      xlonr=xlon*pi/180.
448      ylat=-90.0
449      uuaux=+ffpol*sin(xlonr-ddpol)
450      vvaux=-ffpol*cos(xlonr-ddpol)
451      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
452
453      jy=0
454      do ix=0,nxmin1
455        uupol(ix,jy,iz,n)=uupolaux
456        vvpol(ix,jy,iz,n)=vvpolaux
457      end do
458    end do
459
460
461  ! Fix: Set W at pole to the zonally averaged W of the next equator-
462  ! ward parallel of latitude
463
464    do iz=1,nz
465      wdummy=0.
466      jy=1
467      do ix=0,nxmin1
468        wdummy=wdummy+ww(ix,jy,iz,n)
469      end do
470      wdummy=wdummy/real(nx)
471      jy=0
472      do ix=0,nxmin1
473        ww(ix,jy,iz,n)=wdummy
474      end do
475    end do
476  endif
477
478
479  !write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
480  !   create a cloud and rainout/washout field, clouds occur where rh>80%
481  !   total cloudheight is stored at level 0
482  do jy=0,nymin1
483    do ix=0,nxmin1
484      rain_cloud_above=0
485      lsp=lsprec(ix,jy,1,n)
486      convp=convprec(ix,jy,1,n)
487      cloudsh(ix,jy,n)=0
488      do kz_inv=1,nz-1
489         kz=nz-kz_inv+1
490         pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
491         rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
492         clouds(ix,jy,kz,n)=0
493         if (rh.gt.0.8) then ! in cloud
494            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
495               rain_cloud_above=1
496               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
497               height(kz)-height(kz-1)
498               if (lsp.ge.convp) then
499                  clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
500               else
501                  clouds(ix,jy,kz,n)=2 ! convp dominated rainout
502               endif
503            else ! no precipitation
504                  clouds(ix,jy,kz,n)=1 ! cloud
505            endif
506         else ! no cloud
507            if (rain_cloud_above.eq.1) then ! scavenging
508               if (lsp.ge.convp) then
509                  clouds(ix,jy,kz,n)=5 ! lsp dominated washout
510               else
511                  clouds(ix,jy,kz,n)=4 ! convp dominated washout
512               endif
513            endif
514         endif
515      end do
516    end do
517  end do
518
519end subroutine verttransform
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG