source: flexpart.git/src/verttransform_gfs.f90 @ 6ecb30a

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 6ecb30a was 6ecb30a, checked in by Espen Sollum ATMOS <eso@…>, 7 years ago

Merged changes from CTBTO project

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