source: flexpart.git/src/verttransform_ecmwf.f90 @ 79e0349

devrelease-10
Last change on this file since 79e0349 was 79e0349, checked in by Sabine <sabine.eckhardt@…>, 9 months ago

Doing the fit for scavenging when using metfields with no clouds

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