source: flexpart.git/src/verttransform_ecmwf.f90 @ 47f96e5

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

init_r was mistakenly set to .true. in commit de4c5e9885c62eece05be19e1c0165ade923e201

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