source: trunk/src/verttransform_gfs.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.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