source: trunk/src/verttransform.f90 @ 28

Last change on this file since 28 was 24, checked in by igpis, 8 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