source: flexpart.git/src/verttransform.f90 @ 7999df47

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

Completed handling of nested wind fields with cloud water (for wet deposition).

  • Property mode set to 100644
File size: 29.1 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,numwfmem) cloud field for wet deposition    *
56! uu(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in x-direction [m/s]*
57! vv(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in y-direction [m/s]*
58! ww(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in z-direction      *
59!                                          [deltaeta/s]                      *
60! tt(0:nxmax,0:nymax,nzmax,numwfmem)     temperature [K]                     *
61! pv(0:nxmax,0:nymax,nzmax,numwfmem)     potential voriticity (pvu)          *
62! ps(0:nxmax,0:nymax,numwfmem)           surface pressure [Pa]               *
63!                                                                            *
64!*****************************************************************************
65
66  use par_mod
67  use com_mod
68  use cmapf_mod, only: cc2gll
69!  use mpi_mod
70
71  implicit none
72
73  real,intent(in),dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: uuh,vvh,pvh
74  real,intent(in),dimension(0:nxmax-1,0:nymax-1,nwzmax) :: wwh
75
76  real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: rhoh,uvzlev,wzlev
77  real,dimension(0:nxmax-1,0:nymax-1,nzmax) :: pinmconv
78  real,dimension(0:nxmax-1,0:nymax-1) ::  tvold,pold,pint,tv
79  real,dimension(0:nymax-1) :: cosf
80
81  integer,dimension(0:nxmax-1,0:nymax-1) :: rain_cloud_above,idx
82
83  integer :: ix,jy,kz,iz,n,kmin,ix1,jy1,ixp,jyp,ixm,jym,kz_inv
84  real :: f_qvsat,pressure,rh,lsp,convp,prec
85  real :: ew,dz1,dz2,dz
86  real :: xlon,ylat,xlonr,dzdx,dzdy
87  real :: dzdx1,dzdx2,dzdy1,dzdy2
88  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
89  real,parameter :: const=r_air/ga
90  real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
91
92  logical :: init = .true.
93
94  !ZHG SEP 2014 tests 
95  ! integer :: cloud_ver,cloud_min, cloud_max
96  ! integer ::teller(5), convpteller=0, lspteller=0
97  ! real :: cloud_col_wat, cloud_water
98  !ZHG 2015 temporary variables for testing
99  ! real :: rcw(0:nxmax-1,0:nymax-1)
100  ! real :: rpc(0:nxmax-1,0:nymax-1)
101  ! character(len=60) :: zhgpath='/xnilu_wrk/flex_wrk/zhg/'
102  ! character(len=60) :: fnameA,fnameB,fnameC,fnameD,fnameE,fnameF,fnameG,fnameH
103  ! CHARACTER(LEN=3)  :: aspec
104  ! integer :: virr=0
105  real :: tot_cloud_h
106  real :: dbg_height(nzmax)
107!ZHG
108
109!*************************************************************************
110! If verttransform is called the first time, initialize heights of the   *
111! z levels in meter. The heights are the heights of model levels, where  *
112! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
113! the vertical resolution in the z system is doubled. As reference point,*
114! the lower left corner of the grid is used.                             *
115! Unlike in the eta system, no difference between heights for u,v and    *
116! heights for w exists.                                                  *
117!*************************************************************************
118
119
120!eso measure CPU time
121!  call mpif_mtime('verttransform',0)
122
123  if (init) then
124
125! Search for a point with high surface pressure (i.e. not above significant topography)
126! Then, use this point to construct a reference z profile, to be used at all times
127!*****************************************************************************
128
129    do jy=0,nymin1
130      do ix=0,nxmin1
131        if (ps(ix,jy,1,n).gt.100000.) then
132          ixm=ix
133          jym=jy
134          goto 3
135        endif
136      end do
137    end do
1383   continue
139
140
141    tvold(ixm,jym)=tt2(ixm,jym,1,n)*(1.+0.378*ew*(td2(ixm,jym,1,n))/ &
142         ps(ixm,jym,1,n))
143    pold(ixm,jym)=ps(ixm,jym,1,n)
144    height(1)=0.
145
146    do kz=2,nuvz
147      pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
148      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
149
150      if (abs(tv(ixm,jym)-tvold(ixm,jym)).gt.0.2) then
151        height(kz)= &
152             height(kz-1)+const*log(pold(ixm,jym)/pint(ixm,jym))* &
153             (tv(ixm,jym)-tvold(ixm,jym))/log(tv(ixm,jym)/tvold(ixm,jym))
154      else
155        height(kz)=height(kz-1)+ &
156             const*log(pold(ixm,jym)/pint(ixm,jym))*tv(ixm,jym)
157      endif
158
159      tvold(ixm,jym)=tv(ixm,jym)
160      pold(ixm,jym)=pint(ixm,jym)
161    end do
162
163
164! Determine highest levels that can be within PBL
165!************************************************
166
167    do kz=1,nz
168      if (height(kz).gt.hmixmax) then
169        nmixz=kz
170        goto 9
171      endif
172    end do
1739   continue
174
175! Do not repeat initialization of the Cartesian z grid
176!*****************************************************
177
178    init=.false.
179
180    dbg_height = height
181
182  endif
183
184
185! Loop over the whole grid
186!*************************
187
188
189  do jy=0,nymin1
190    do ix=0,nxmin1
191      tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew*(td2(ix,jy,1,n))/ &
192           ps(ix,jy,1,n))
193    enddo
194  enddo
195  pold=ps(:,:,1,n)
196  uvzlev(:,:,1)=0.
197  wzlev(:,:,1)=0.
198  rhoh(:,:,1)=pold/(r_air*tvold)
199
200
201! Compute heights of eta levels
202!******************************
203
204  do kz=2,nuvz
205    pint=akz(kz)+bkz(kz)*ps(:,:,1,n)
206    tv=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n))
207    rhoh(:,:,kz)=pint(:,:)/(r_air*tv)
208
209    where (abs(tv-tvold).gt.0.2)
210      uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)* &
211           (tv-tvold)/log(tv/tvold)
212    elsewhere
213      uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)*tv
214    endwhere
215
216    tvold=tv
217    pold=pint
218  end do
219
220
221  do kz=2,nwz-1
222    wzlev(:,:,kz)=(uvzlev(:,:,kz+1)+uvzlev(:,:,kz))/2.
223  end do
224  wzlev(:,:,nwz)=wzlev(:,:,nwz-1)+ &
225       uvzlev(:,:,nuvz)-uvzlev(:,:,nuvz-1)
226
227! pinmconv=(h2-h1)/(p2-p1)
228
229  pinmconv(:,:,1)=(uvzlev(:,:,2))/ &
230       ((aknew(2)+bknew(2)*ps(:,:,1,n))- &
231       (aknew(1)+bknew(1)*ps(:,:,1,n)))
232  do kz=2,nz-1
233    pinmconv(:,:,kz)=(uvzlev(:,:,kz+1)-uvzlev(:,:,kz-1))/ &
234         ((aknew(kz+1)+bknew(kz+1)*ps(:,:,1,n))- &
235         (aknew(kz-1)+bknew(kz-1)*ps(:,:,1,n)))
236  end do
237  pinmconv(:,:,nz)=(uvzlev(:,:,nz)-uvzlev(:,:,nz-1))/ &
238       ((aknew(nz)+bknew(nz)*ps(:,:,1,n))- &
239       (aknew(nz-1)+bknew(nz-1)*ps(:,:,1,n)))
240
241! Levels, where u,v,t and q are given
242!************************************
243
244
245  uu(:,:,1,n)=uuh(:,:,1)
246  vv(:,:,1,n)=vvh(:,:,1)
247  tt(:,:,1,n)=tth(:,:,1,n)
248  qv(:,:,1,n)=qvh(:,:,1,n)
249!hg adding the cloud water
250  if (readclouds) then
251    clwc(:,:,1,n)=clwch(:,:,1,n)
252    if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n)
253  end if
254!hg
255  pv(:,:,1,n)=pvh(:,:,1)
256  rho(:,:,1,n)=rhoh(:,:,1)
257
258  uu(:,:,nz,n)=uuh(:,:,nuvz)
259  vv(:,:,nz,n)=vvh(:,:,nuvz)
260  tt(:,:,nz,n)=tth(:,:,nuvz,n)
261  qv(:,:,nz,n)=qvh(:,:,nuvz,n)
262!hg adding the cloud water
263  if (readclouds) then
264    clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
265    if (.not.sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
266  end if
267!hg
268  pv(:,:,nz,n)=pvh(:,:,nuvz)
269  rho(:,:,nz,n)=rhoh(:,:,nuvz)
270
271
272  kmin=2
273  idx=kmin
274  do iz=2,nz-1
275    do jy=0,nymin1
276      do ix=0,nxmin1
277        if(height(iz).gt.uvzlev(ix,jy,nuvz)) then
278          uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
279          vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
280          tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
281          qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
282!hg adding the cloud water
283          if (readclouds) then
284            clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
285            if (.not.sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
286          end if
287!hg
288          pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
289          rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
290        else
291          innuvz: do kz=idx(ix,jy),nuvz
292            if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
293                 (height(iz).le.uvzlev(ix,jy,kz))) then
294              idx(ix,jy)=kz
295              exit innuvz
296            endif
297          enddo innuvz
298        endif
299      enddo
300    enddo
301    do jy=0,nymin1
302      do ix=0,nxmin1
303        if(height(iz).le.uvzlev(ix,jy,nuvz)) then
304          kz=idx(ix,jy)
305          dz1=height(iz)-uvzlev(ix,jy,kz-1)
306          dz2=uvzlev(ix,jy,kz)-height(iz)
307          dz=dz1+dz2
308          uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
309          vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
310          tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
311               +tth(ix,jy,kz,n)*dz1)/dz
312          qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
313               +qvh(ix,jy,kz,n)*dz1)/dz
314!hg adding the cloud water
315          if (readclouds) then
316            clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
317            if (.not.sumclouds) &
318                 &ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
319          end if
320!hg
321          pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
322          rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)*dz2+rhoh(ix,jy,kz)*dz1)/dz
323        endif
324      enddo
325    enddo
326  enddo
327
328
329! Levels, where w is given
330!*************************
331
332  ww(:,:,1,n)=wwh(:,:,1)*pinmconv(:,:,1)
333  ww(:,:,nz,n)=wwh(:,:,nwz)*pinmconv(:,:,nz)
334  kmin=2
335  idx=kmin
336  do iz=2,nz
337    do jy=0,nymin1
338      do ix=0,nxmin1
339        inn:         do kz=idx(ix,jy),nwz
340          if(idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1).and. &
341               height(iz).le.wzlev(ix,jy,kz)) then
342            idx(ix,jy)=kz
343            exit inn
344          endif
345        enddo inn
346      enddo
347    enddo
348    do jy=0,nymin1
349      do ix=0,nxmin1
350        kz=idx(ix,jy)
351        dz1=height(iz)-wzlev(ix,jy,kz-1)
352        dz2=wzlev(ix,jy,kz)-height(iz)
353        dz=dz1+dz2
354        ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(ix,jy,kz-1)*dz2 &
355             +wwh(ix,jy,kz)*pinmconv(ix,jy,kz)*dz1)/dz
356      enddo
357    enddo
358  enddo
359
360! Compute density gradients at intermediate levels
361!*************************************************
362
363  drhodz(:,:,1,n)=(rho(:,:,2,n)-rho(:,:,1,n))/ &
364       (height(2)-height(1))
365  do kz=2,nz-1
366    drhodz(:,:,kz,n)=(rho(:,:,kz+1,n)-rho(:,:,kz-1,n))/ &
367         (height(kz+1)-height(kz-1))
368  end do
369  drhodz(:,:,nz,n)=drhodz(:,:,nz-1,n)
370
371!    end do
372!  end do
373
374
375!****************************************************************
376! Compute slope of eta levels in windward direction and resulting
377! vertical wind correction
378!****************************************************************
379
380  do jy=1,ny-2
381    cosf(jy)=1./cos((real(jy)*dy+ylat0)*pi180)
382  enddo
383
384  kmin=2
385  idx=kmin
386  do iz=2,nz-1
387    do jy=1,ny-2
388      do ix=1,nx-2
389
390        inneta: do kz=idx(ix,jy),nz
391          if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
392               (height(iz).le.uvzlev(ix,jy,kz))) then
393            idx(ix,jy)=kz
394            exit inneta
395          endif
396        enddo inneta
397      enddo
398    enddo
399
400    do jy=1,ny-2
401      do ix=1,nx-2
402        kz=idx(ix,jy)
403        dz1=height(iz)-uvzlev(ix,jy,kz-1)
404        dz2=uvzlev(ix,jy,kz)-height(iz)
405        dz=dz1+dz2
406        ix1=ix-1
407        jy1=jy-1
408        ixp=ix+1
409        jyp=jy+1
410
411        dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2.
412        dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2.
413        dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
414
415        dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2.
416        dzdy2=(uvzlev(ix,jyp,kz)-uvzlev(ix,jy1,kz))/2.
417        dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
418
419        ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*uu(ix,jy,iz,n)*dxconst*cosf(jy)+dzdy*vv(ix,jy,iz,n)*dyconst)
420
421      end do
422
423    end do
424  end do
425
426! If north pole is in the domain, calculate wind velocities in polar
427! stereographic coordinates
428!*******************************************************************
429
430  if (nglobal) then
431    do iz=1,nz
432      do jy=int(switchnorthg)-2,nymin1
433        ylat=ylat0+real(jy)*dy
434        do ix=0,nxmin1
435          xlon=xlon0+real(ix)*dx
436          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
437               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
438               vvpol(ix,jy,iz,n))
439        end do
440      end do
441    end do
442
443
444    do iz=1,nz
445
446! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
447!
448!   AMSnauffer Nov 18 2004 Added check for case vv=0
449!
450      xlon=xlon0+real(nx/2-1)*dx
451      xlonr=xlon*pi/180.
452      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
453           vv(nx/2-1,nymin1,iz,n)**2)
454      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
455        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
456             vv(nx/2-1,nymin1,iz,n))-xlonr
457      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
458        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
459             vv(nx/2-1,nymin1,iz,n))-xlonr
460      else
461        ddpol=pi/2-xlonr
462      endif
463      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
464      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
465
466! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
467      xlon=180.0
468      xlonr=xlon*pi/180.
469      ylat=90.0
470      uuaux=-ffpol*sin(xlonr+ddpol)
471      vvaux=-ffpol*cos(xlonr+ddpol)
472      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
473           vvpolaux)
474
475      jy=nymin1
476      do ix=0,nxmin1
477        uupol(ix,jy,iz,n)=uupolaux
478        vvpol(ix,jy,iz,n)=vvpolaux
479      end do
480    end do
481
482
483! Fix: Set W at pole to the zonally averaged W of the next equator-
484! ward parallel of latitude
485
486    do iz=1,nz
487      wdummy=0.
488      jy=ny-2
489      do ix=0,nxmin1
490        wdummy=wdummy+ww(ix,jy,iz,n)
491      end do
492      wdummy=wdummy/real(nx)
493      jy=nymin1
494      do ix=0,nxmin1
495        ww(ix,jy,iz,n)=wdummy
496      end do
497    end do
498
499  endif
500
501
502! If south pole is in the domain, calculate wind velocities in polar
503! stereographic coordinates
504!*******************************************************************
505
506  if (sglobal) then
507    do iz=1,nz
508      do jy=0,int(switchsouthg)+3
509        ylat=ylat0+real(jy)*dy
510        do ix=0,nxmin1
511          xlon=xlon0+real(ix)*dx
512          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
513               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
514               vvpol(ix,jy,iz,n))
515        end do
516      end do
517    end do
518
519    do iz=1,nz
520
521! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
522!
523!   AMSnauffer Nov 18 2004 Added check for case vv=0
524!
525      xlon=xlon0+real(nx/2-1)*dx
526      xlonr=xlon*pi/180.
527      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
528           vv(nx/2-1,0,iz,n)**2)
529      if (vv(nx/2-1,0,iz,n).lt.0.) then
530        ddpol=atan(uu(nx/2-1,0,iz,n)/ &
531             vv(nx/2-1,0,iz,n))+xlonr
532      else if (vv(nx/2-1,0,iz,n).gt.0.) then
533        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
534             vv(nx/2-1,0,iz,n))+xlonr
535      else
536        ddpol=pi/2-xlonr
537      endif
538      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
539      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
540
541! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
542      xlon=180.0
543      xlonr=xlon*pi/180.
544      ylat=-90.0
545      uuaux=+ffpol*sin(xlonr-ddpol)
546      vvaux=-ffpol*cos(xlonr-ddpol)
547      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
548           vvpolaux)
549
550      jy=0
551      do ix=0,nxmin1
552        uupol(ix,jy,iz,n)=uupolaux
553        vvpol(ix,jy,iz,n)=vvpolaux
554      end do
555    end do
556
557
558! Fix: Set W at pole to the zonally averaged W of the next equator-
559! ward parallel of latitude
560
561    do iz=1,nz
562      wdummy=0.
563      jy=1
564      do ix=0,nxmin1
565        wdummy=wdummy+ww(ix,jy,iz,n)
566      end do
567      wdummy=wdummy/real(nx)
568      jy=0
569      do ix=0,nxmin1
570        ww(ix,jy,iz,n)=wdummy
571      end do
572    end do
573  endif
574
575
576!*********************************************************************************** 
577  if (readclouds) then !HG METHOD
578! The method is loops all grids vertically and constructs the 3D matrix for clouds
579! Cloud top and cloud bottom gid cells are assigned as well as the total column
580! cloud water. For precipitating grids, the type and whether it is in or below
581! cloud scavenging are assigned with numbers 2-5 (following the old metod).
582! Distinction is done for lsp and convp though they are treated the same in regards
583! to scavenging. Also clouds that are not precipitating are defined which may be
584! to include future cloud processing by non-precipitating-clouds.
585!***********************************************************************************
586    write(*,*) 'Global ECMWF fields: using cloud water'
587    clw(:,:,:,n)=0.0
588    icloud_stats(:,:,:,n)=0.0
589    clouds(:,:,:,n)=0
590! If water/ice are read separately into clwc and ciwc, store sum in clwc
591    if (.not.sumclouds) then
592      clwc(:,:,:,n) = clwc(:,:,:,n) + ciwc(:,:,:,n)
593    end if
594    do jy=0,nymin1
595      do ix=0,nxmin1
596        lsp=lsprec(ix,jy,1,n)
597        convp=convprec(ix,jy,1,n)
598        prec=lsp+convp
599        tot_cloud_h=0
600! Find clouds in the vertical
601        do kz=1, nz-1 !go from top to bottom
602          if (clwc(ix,jy,kz,n).gt.0) then     
603! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
604            clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
605            tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
606            icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]
607            icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz))                     ! Cloud BOT height stats      [m]
608!ZHG 2015 extra for testing
609!            clh(ix,jy,kz,n)=height(kz+1)-height(kz)
610            icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
611            icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz))          ! Cloud TOP height            [m]
612!ZHG
613          endif
614        end do
615
616! If Precipitation. Define removal type in the vertical
617        if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
618
619          do kz=nz,1,-1 !go Bottom up!
620            if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
621              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
622              clouds(ix,jy,kz,n)=1                               ! is a cloud
623              if (lsp.ge.convp) then
624                clouds(ix,jy,kz,n)=3                            ! lsp in-cloud
625              else
626                clouds(ix,jy,kz,n)=2                             ! convp in-cloud
627              endif                                              ! convective or large scale
628            elseif((clw(ix,jy,kz,n).le.0) .and. (icloud_stats(ix,jy,3,n).ge.height(kz)) ) then ! is below cloud
629              if (lsp.ge.convp) then
630                clouds(ix,jy,kz,n)=5                             ! lsp dominated washout
631              else
632                clouds(ix,jy,kz,n)=4                             ! convp dominated washout
633              endif                                              ! convective or large scale
634            endif
635
636            if (height(kz).ge. 19000) then                        ! set a max height for removal
637              clouds(ix,jy,kz,n)=0
638            endif !clw>0
639          end do !nz
640        endif ! precipitation
641      end do
642    end do
643
644! eso: copy the relevant data to clw4 to reduce amount of communicated data for MPI
645    clw4(:,:,n) = icloud_stats(:,:,4,n)
646
647!**************************************************************************
648  else       ! use old definitions
649!**************************************************************************
650!   create a cloud and rainout/washout field, clouds occur where rh>80%
651!   total cloudheight is stored at level 0
652    write(*,*) 'Global fields: using cloud water from Parameterization'
653    do jy=0,nymin1
654      do ix=0,nxmin1
655! OLD METHOD
656        rain_cloud_above(ix,jy)=0
657        lsp=lsprec(ix,jy,1,n)
658        convp=convprec(ix,jy,1,n)
659        cloudsh(ix,jy,n)=0
660        do kz_inv=1,nz-1
661          kz=nz-kz_inv+1
662          pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
663          rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
664          clouds(ix,jy,kz,n)=0
665          if (rh.gt.0.8) then ! in cloud
666            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
667              rain_cloud_above(ix,jy)=1
668              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
669                   height(kz)-height(kz-1)
670              if (lsp.ge.convp) then
671                clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
672              else
673                clouds(ix,jy,kz,n)=2 ! convp dominated rainout
674              endif
675            else ! no precipitation
676              clouds(ix,jy,kz,n)=1 ! cloud
677            endif
678          else ! no cloud
679            if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
680              if (lsp.ge.convp) then
681                clouds(ix,jy,kz,n)=5 ! lsp dominated washout
682              else
683                clouds(ix,jy,kz,n)=4 ! convp dominated washout
684              endif
685            endif
686          endif
687        end do
688!END OLD METHOD
689      end do
690    end do
691  endif !readclouds
692
693
694     !********* TEST ***************
695     ! WRITE OUT SOME TEST VARIABLES
696     !********* TEST ************'**
697!teller(:)=0
698!virr=virr+1
699!WRITE(aspec, '(i3.3)'), virr
700
701!if (readclouds) then
702!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement.txt'
703!else
704!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement_old.txt'
705!endif
706!
707!OPEN(UNIT=118, FILE=fnameH,FORM='FORMATTED',STATUS = 'UNKNOWN')
708!do kz_inv=1,nz-1
709!  kz=nz-kz_inv+1
710!  !kz=91
711!  do jy=0,nymin1
712!     do ix=0,nxmin1
713!          if (clouds(ix,jy,kz,n).eq.1) teller(1)=teller(1)+1 ! no precipitation cloud
714!          if (clouds(ix,jy,kz,n).eq.2) teller(2)=teller(2)+1 ! convp dominated rainout
715!          if (clouds(ix,jy,kz,n).eq.3) teller(3)=teller(3)+1 ! lsp dominated rainout
716!          if (clouds(ix,jy,kz,n).eq.4) teller(4)=teller(4)+1 ! convp dominated washout
717!          if (clouds(ix,jy,kz,n).eq.5) teller(5)=teller(5)+1 ! lsp dominated washout
718!         
719!        !  write(*,*) height(kz),teller
720!     end do
721!  end do
722!  write(118,*) height(kz),teller
723!  teller(:)=0
724!end do
725!teller(:)=0
726!write(*,*) teller
727!write(*,*) aspec
728!
729!fnameA=trim(zhgpath)//trim(aspec)//'cloudV.txt'
730!fnameB=trim(zhgpath)//trim(aspec)//'cloudT.txt'
731!fnameC=trim(zhgpath)//trim(aspec)//'cloudB.txt'
732!fnameD=trim(zhgpath)//trim(aspec)//'cloudW.txt'
733!fnameE=trim(zhgpath)//trim(aspec)//'old_cloudV.txt'
734!fnameF=trim(zhgpath)//trim(aspec)//'lsp.txt'
735!fnameG=trim(zhgpath)//trim(aspec)//'convp.txt'
736!if (readclouds) then
737!OPEN(UNIT=111, FILE=fnameA,FORM='FORMATTED',STATUS = 'UNKNOWN')
738!OPEN(UNIT=112, FILE=fnameB,FORM='FORMATTED',STATUS = 'UNKNOWN')
739!OPEN(UNIT=113, FILE=fnameC,FORM='FORMATTED',STATUS = 'UNKNOWN')
740!OPEN(UNIT=114, FILE=fnameD,FORM='FORMATTED',STATUS = 'UNKNOWN')
741!else
742!OPEN(UNIT=115, FILE=fnameE,FORM='FORMATTED',STATUS = 'UNKNOWN')
743!OPEN(UNIT=116, FILE=fnameF,FORM='FORMATTED',STATUS = 'UNKNOWN')
744!OPEN(UNIT=117, FILE=fnameG,FORM='FORMATTED',STATUS = 'UNKNOWN')
745!endif
746!
747!do ix=0,nxmin1
748!if (readclouds) then
749!write(111,*) (icloud_stats(ix,jy,1,n),jy=0,nymin1)
750!write(112,*) (icloud_stats(ix,jy,2,n),jy=0,nymin1)
751!write(113,*) (icloud_stats(ix,jy,3,n),jy=0,nymin1)
752!write(114,*) (icloud_stats(ix,jy,4,n),jy=0,nymin1)
753!else
754!write(115,*) (cloudsh(ix,jy,n),jy=0,nymin1)    !integer
755!write(116,*) (lsprec(ix,jy,1,n),jy=0,nymin1)   !7.83691406E-02
756!write(117,*) (convprec(ix,jy,1,n),jy=0,nymin1) !5.38330078E-02
757!endif
758!end do
759!
760!if (readclouds) then
761!CLOSE(111)
762!CLOSE(112)
763!CLOSE(113)
764!CLOSE(114)
765!else
766!CLOSE(115)
767!CLOSE(116)
768!CLOSE(117)
769!endif
770!
771!END ********* TEST *************** END
772! WRITE OUT SOME TEST VARIABLES
773!END ********* TEST *************** END
774
775
776! PS 2012
777!      lsp=lsprec(ix,jy,1,n)
778!      convp=convprec(ix,jy,1,n)
779!          prec=lsp+convp
780!          if (lsp.gt.convp) then !  prectype='lsp'
781!            lconvectprec = .false.
782!          else ! prectype='cp'
783!            lconvectprec = .true.
784!           endif
785!      else ! windfields does not contain cloud data
786!          rhmin = 0.90 ! standard condition for presence of clouds
787!PS       note that original by Sabine Eckhart was 80%
788!PS       however, for T<-20 C we consider saturation over ice
789!PS       so I think 90% should be enough         
790!          icloudbot(ix,jy,n)=icmv
791!          icloudtop=icmv ! this is just a local variable
792!98        do kz=1,nz
793!            pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
794!            rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
795!ps            if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
796!            if (rh .gt. rhmin) then
797!              if (icloudbot(ix,jy,n) .eq. icmv) then
798!                icloudbot(ix,jy,n)=nint(height(kz))
799!              endif
800!              icloudtop=nint(height(kz)) ! use int to save memory
801!            endif
802    ! end do
803!PS try to get a cloud thicker than 50 m
804!PS if there is at least .01 mm/h  - changed to 0.002 and put into
805!PS parameter precpmin       
806!          if ((icloudbot(ix,jy,n) .eq. icmv .or. &
807!              icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
808!              prec .gt. precmin) then
809!            rhmin = rhmin - 0.05
810!            if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
811!   end if
812
813!PS is based on looking at a limited set of comparison data
814!          if (lconvectprec .and. icloudtop .lt. 6000 .and. &
815!             prec .gt. precmin) then
816!
817!            if (convp .lt. 0.1) then
818!              icloudbot(ix,jy,n) = 500
819!              icloudtop =         8000
820!            else
821!              icloudbot(ix,jy,n) = 0
822!              icloudtop =      10000
823!            endif
824!          endif
825!          if (icloudtop .ne. icmv) then
826!            icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
827!          else
828!            icloudthck(ix,jy,n) = icmv
829!          endif
830!PS  get rid of too thin clouds     
831!          if (icloudthck(ix,jy,n) .lt. 50) then
832!            icloudbot(ix,jy,n)=icmv
833!            icloudthck(ix,jy,n)=icmv
834!          endif
835!hg__________________________________
836!           rcw(ix,jy)=2E-7*prec**0.36
837!           rpc(ix,jy)=prec
838!hg end______________________________
839
840!      endif !hg read clouds
841
842
843
844
845!eso measure CPU time
846!  call mpif_mtime('verttransform',1)
847
848!eso print out the same measure as in Leo's routine
849    ! write(*,*) 'verttransform: ', &
850    !      sum(tt(:,:,:,n)*tt(:,:,:,n)), &
851    !      sum(uu(:,:,:,n)*uu(:,:,:,n)),sum(vv(:,:,:,n)*vv(:,:,:,n)),&
852    !      sum(qv(:,:,:,n)*qv(:,:,:,n)),sum(pv(:,:,:,n)*pv(:,:,:,n)),&
853    !      sum(rho(:,:,:,n)*rho(:,:,:,n)),sum(drhodz(:,:,:,n)*drhodz(:,:,:,n)),&
854    !      sum(ww(:,:,:,n)*ww(:,:,:,n)), &
855    !      sum(clouds(:,:,:,n)), sum(cloudsh(:,:,n)),sum(idx),sum(pinmconv)
856end subroutine verttransform
857
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG