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