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

devrelease-10univie
Last change on this file since 6ecb30a was 6ecb30a, checked in by Espen Sollum ATMOS <eso@…>, 2 years ago

Merged changes from CTBTO project

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