source: flexpart.git/src/verttransform.f90 @ ceb7e4d

devrelease-10univie
Last change on this file since ceb7e4d was ceb7e4d, checked in by Sabine <sabine.eckhardt@…>, 3 years ago

commented testing vars in verttransf

  • Property mode set to 100644
File size: 29.2 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,cloudh_min,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    ctwc(:,:,n)=0.0
590    clouds(:,:,:,n)=0
591! If water/ice are read separately into clwc and ciwc, store sum in clwc
592    if (.not.sumclouds) then
593      clwc(:,:,:,n) = clwc(:,:,:,n) + ciwc(:,:,:,n)
594    end if
595    do jy=0,nymin1
596      do ix=0,nxmin1
597        lsp=lsprec(ix,jy,1,n)
598        convp=convprec(ix,jy,1,n)
599        prec=lsp+convp
600!        tot_cloud_h=0
601! Find clouds in the vertical
602        do kz=1, nz-1 !go from top to bottom
603          if (clwc(ix,jy,kz,n).gt.0) then     
604! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
605            clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
606!            tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
607           
608!            icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]
609            ctwc(ix,jy,n) = ctwc(ix,jy,n)+clw(ix,jy,kz,n)
610!            icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz))                     ! Cloud BOT height stats      [m]
611            cloudh_min=min(height(kz+1),height(kz))
612!ZHG 2015 extra for testing
613!            clh(ix,jy,kz,n)=height(kz+1)-height(kz)
614!            icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
615!            icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz))          ! Cloud TOP height            [m]
616!ZHG
617          endif
618        end do
619
620! If Precipitation. Define removal type in the vertical
621        if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
622
623          do kz=nz,1,-1 !go Bottom up!
624            if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
625              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
626              clouds(ix,jy,kz,n)=1                               ! is a cloud
627              if (lsp.ge.convp) then
628                clouds(ix,jy,kz,n)=3                            ! lsp in-cloud
629              else
630                clouds(ix,jy,kz,n)=2                             ! convp in-cloud
631              endif                                              ! convective or large scale
632            elseif((clw(ix,jy,kz,n).le.0) .and. (cloudh_min.ge.height(kz))) then ! is below cloud
633              if (lsp.ge.convp) then
634                clouds(ix,jy,kz,n)=5                             ! lsp dominated washout
635              else
636                clouds(ix,jy,kz,n)=4                             ! convp dominated washout
637              endif                                              ! convective or large scale
638            endif
639
640            if (height(kz).ge. 19000) then                        ! set a max height for removal
641              clouds(ix,jy,kz,n)=0
642            endif !clw>0
643          end do !nz
644        endif ! precipitation
645      end do
646    end do
647
648! eso: copy the relevant data to clw4 to reduce amount of communicated data for MPI
649!    ctwc(:,:,n) = icloud_stats(:,:,4,n)
650
651!**************************************************************************
652  else       ! use old definitions
653!**************************************************************************
654!   create a cloud and rainout/washout field, clouds occur where rh>80%
655!   total cloudheight is stored at level 0
656    write(*,*) 'Global fields: using cloud water from Parameterization'
657    do jy=0,nymin1
658      do ix=0,nxmin1
659! OLD METHOD
660        rain_cloud_above(ix,jy)=0
661        lsp=lsprec(ix,jy,1,n)
662        convp=convprec(ix,jy,1,n)
663        cloudsh(ix,jy,n)=0
664        do kz_inv=1,nz-1
665          kz=nz-kz_inv+1
666          pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
667          rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
668          clouds(ix,jy,kz,n)=0
669          if (rh.gt.0.8) then ! in cloud
670            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
671              rain_cloud_above(ix,jy)=1
672              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
673                   height(kz)-height(kz-1)
674              if (lsp.ge.convp) then
675                clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
676              else
677                clouds(ix,jy,kz,n)=2 ! convp dominated rainout
678              endif
679            else ! no precipitation
680              clouds(ix,jy,kz,n)=1 ! cloud
681            endif
682          else ! no cloud
683            if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
684              if (lsp.ge.convp) then
685                clouds(ix,jy,kz,n)=5 ! lsp dominated washout
686              else
687                clouds(ix,jy,kz,n)=4 ! convp dominated washout
688              endif
689            endif
690          endif
691        end do
692!END OLD METHOD
693      end do
694    end do
695  endif !readclouds
696
697
698     !********* TEST ***************
699     ! WRITE OUT SOME TEST VARIABLES
700     !********* TEST ************'**
701!teller(:)=0
702!virr=virr+1
703!WRITE(aspec, '(i3.3)'), virr
704
705!if (readclouds) then
706!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement.txt'
707!else
708!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement_old.txt'
709!endif
710!
711!OPEN(UNIT=118, FILE=fnameH,FORM='FORMATTED',STATUS = 'UNKNOWN')
712!do kz_inv=1,nz-1
713!  kz=nz-kz_inv+1
714!  !kz=91
715!  do jy=0,nymin1
716!     do ix=0,nxmin1
717!          if (clouds(ix,jy,kz,n).eq.1) teller(1)=teller(1)+1 ! no precipitation cloud
718!          if (clouds(ix,jy,kz,n).eq.2) teller(2)=teller(2)+1 ! convp dominated rainout
719!          if (clouds(ix,jy,kz,n).eq.3) teller(3)=teller(3)+1 ! lsp dominated rainout
720!          if (clouds(ix,jy,kz,n).eq.4) teller(4)=teller(4)+1 ! convp dominated washout
721!          if (clouds(ix,jy,kz,n).eq.5) teller(5)=teller(5)+1 ! lsp dominated washout
722!         
723!        !  write(*,*) height(kz),teller
724!     end do
725!  end do
726!  write(118,*) height(kz),teller
727!  teller(:)=0
728!end do
729!teller(:)=0
730!write(*,*) teller
731!write(*,*) aspec
732!
733!fnameA=trim(zhgpath)//trim(aspec)//'cloudV.txt'
734!fnameB=trim(zhgpath)//trim(aspec)//'cloudT.txt'
735!fnameC=trim(zhgpath)//trim(aspec)//'cloudB.txt'
736!fnameD=trim(zhgpath)//trim(aspec)//'cloudW.txt'
737!fnameE=trim(zhgpath)//trim(aspec)//'old_cloudV.txt'
738!fnameF=trim(zhgpath)//trim(aspec)//'lsp.txt'
739!fnameG=trim(zhgpath)//trim(aspec)//'convp.txt'
740!if (readclouds) then
741!OPEN(UNIT=111, FILE=fnameA,FORM='FORMATTED',STATUS = 'UNKNOWN')
742!OPEN(UNIT=112, FILE=fnameB,FORM='FORMATTED',STATUS = 'UNKNOWN')
743!OPEN(UNIT=113, FILE=fnameC,FORM='FORMATTED',STATUS = 'UNKNOWN')
744!OPEN(UNIT=114, FILE=fnameD,FORM='FORMATTED',STATUS = 'UNKNOWN')
745!else
746!OPEN(UNIT=115, FILE=fnameE,FORM='FORMATTED',STATUS = 'UNKNOWN')
747!OPEN(UNIT=116, FILE=fnameF,FORM='FORMATTED',STATUS = 'UNKNOWN')
748!OPEN(UNIT=117, FILE=fnameG,FORM='FORMATTED',STATUS = 'UNKNOWN')
749!endif
750!
751!do ix=0,nxmin1
752!if (readclouds) then
753!write(111,*) (icloud_stats(ix,jy,1,n),jy=0,nymin1)
754!write(112,*) (icloud_stats(ix,jy,2,n),jy=0,nymin1)
755!write(113,*) (icloud_stats(ix,jy,3,n),jy=0,nymin1)
756!write(114,*) (icloud_stats(ix,jy,4,n),jy=0,nymin1)
757!else
758!write(115,*) (cloudsh(ix,jy,n),jy=0,nymin1)    !integer
759!write(116,*) (lsprec(ix,jy,1,n),jy=0,nymin1)   !7.83691406E-02
760!write(117,*) (convprec(ix,jy,1,n),jy=0,nymin1) !5.38330078E-02
761!endif
762!end do
763!
764!if (readclouds) then
765!CLOSE(111)
766!CLOSE(112)
767!CLOSE(113)
768!CLOSE(114)
769!else
770!CLOSE(115)
771!CLOSE(116)
772!CLOSE(117)
773!endif
774!
775!END ********* TEST *************** END
776! WRITE OUT SOME TEST VARIABLES
777!END ********* TEST *************** END
778
779
780! PS 2012
781!      lsp=lsprec(ix,jy,1,n)
782!      convp=convprec(ix,jy,1,n)
783!          prec=lsp+convp
784!          if (lsp.gt.convp) then !  prectype='lsp'
785!            lconvectprec = .false.
786!          else ! prectype='cp'
787!            lconvectprec = .true.
788!           endif
789!      else ! windfields does not contain cloud data
790!          rhmin = 0.90 ! standard condition for presence of clouds
791!PS       note that original by Sabine Eckhart was 80%
792!PS       however, for T<-20 C we consider saturation over ice
793!PS       so I think 90% should be enough         
794!          icloudbot(ix,jy,n)=icmv
795!          icloudtop=icmv ! this is just a local variable
796!98        do kz=1,nz
797!            pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
798!            rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
799!ps            if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
800!            if (rh .gt. rhmin) then
801!              if (icloudbot(ix,jy,n) .eq. icmv) then
802!                icloudbot(ix,jy,n)=nint(height(kz))
803!              endif
804!              icloudtop=nint(height(kz)) ! use int to save memory
805!            endif
806    ! end do
807!PS try to get a cloud thicker than 50 m
808!PS if there is at least .01 mm/h  - changed to 0.002 and put into
809!PS parameter precpmin       
810!          if ((icloudbot(ix,jy,n) .eq. icmv .or. &
811!              icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
812!              prec .gt. precmin) then
813!            rhmin = rhmin - 0.05
814!            if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
815!   end if
816
817!PS is based on looking at a limited set of comparison data
818!          if (lconvectprec .and. icloudtop .lt. 6000 .and. &
819!             prec .gt. precmin) then
820!
821!            if (convp .lt. 0.1) then
822!              icloudbot(ix,jy,n) = 500
823!              icloudtop =         8000
824!            else
825!              icloudbot(ix,jy,n) = 0
826!              icloudtop =      10000
827!            endif
828!          endif
829!          if (icloudtop .ne. icmv) then
830!            icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
831!          else
832!            icloudthck(ix,jy,n) = icmv
833!          endif
834!PS  get rid of too thin clouds     
835!          if (icloudthck(ix,jy,n) .lt. 50) then
836!            icloudbot(ix,jy,n)=icmv
837!            icloudthck(ix,jy,n)=icmv
838!          endif
839!hg__________________________________
840!           rcw(ix,jy)=2E-7*prec**0.36
841!           rpc(ix,jy)=prec
842!hg end______________________________
843
844!      endif !hg read clouds
845
846
847
848
849!eso measure CPU time
850!  call mpif_mtime('verttransform',1)
851
852!eso print out the same measure as in Leo's routine
853    ! write(*,*) 'verttransform: ', &
854    !      sum(tt(:,:,:,n)*tt(:,:,:,n)), &
855    !      sum(uu(:,:,:,n)*uu(:,:,:,n)),sum(vv(:,:,:,n)*vv(:,:,:,n)),&
856    !      sum(qv(:,:,:,n)*qv(:,:,:,n)),sum(pv(:,:,:,n)*pv(:,:,:,n)),&
857    !      sum(rho(:,:,:,n)*rho(:,:,:,n)),sum(drhodz(:,:,:,n)*drhodz(:,:,:,n)),&
858    !      sum(ww(:,:,:,n)*ww(:,:,:,n)), &
859    !      sum(clouds(:,:,:,n)), sum(cloudsh(:,:,n)),sum(idx),sum(pinmconv)
860end subroutine verttransform
861
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG