source: trunk/src/verttransform.f90 @ 20

Last change on this file since 20 was 20, checked in by igpis, 9 years ago

move version 9.1.8 form branches to trunk. Contributions from HSO, saeck, pesei, NIK, RT, XKF, IP and others

File size: 22.8 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine verttransform(n,uuh,vvh,wwh,pvh)
23  !                         i  i   i   i   i
24  !*****************************************************************************
25  !                                                                            *
26  !     This subroutine transforms temperature, dew point temperature and      *
27  !     wind components from eta to meter coordinates.                         *
28  !     The vertical wind component is transformed from Pa/s to m/s using      *
29  !     the conversion factor pinmconv.                                        *
30  !     In addition, this routine calculates vertical density gradients        *
31  !     needed for the parameterization of the turbulent velocities.           *
32  !                                                                            *
33  !     Author: A. Stohl, G. Wotawa                                            *
34  !                                                                            *
35  !     12 August 1996                                                         *
36  !     Update: 16 January 1998                                                *
37  !                                                                            *
38  !     Major update: 17 February 1999                                         *
39  !     by G. Wotawa                                                           *
40  !                                                                            *
41  !     - Vertical levels for u, v and w are put together                      *
42  !     - Slope correction for vertical velocity: Modification of calculation  *
43  !       procedure                                                            *
44  !                                                                            *
45  !*****************************************************************************
46  !  Changes, Bernd C. Krueger, Feb. 2001:
47  !   Variables tth and qvh (on eta coordinates) from common block
48  !*****************************************************************************
49  ! Sabine Eckhardt, March 2007
50  ! added the variable cloud for use with scavenging - descr. in com_mod
51  ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
52  ! note that also other subroutines are affected by the fix
53  !*****************************************************************************
54  !                                                                            *
55  ! Variables:                                                                 *
56  ! nx,ny,nz                        field dimensions in x,y and z direction    *
57  ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition           *
58  ! uu(0:nxmax,0:nymax,nzmax,2)     wind components in x-direction [m/s]       *
59  ! vv(0:nxmax,0:nymax,nzmax,2)     wind components in y-direction [m/s]       *
60  ! ww(0:nxmax,0:nymax,nzmax,2)     wind components in z-direction [deltaeta/s]*
61  ! tt(0:nxmax,0:nymax,nzmax,2)     temperature [K]                            *
62  ! pv(0:nxmax,0:nymax,nzmax,2)     potential voriticity (pvu)                 *
63  ! ps(0:nxmax,0:nymax,2)           surface pressure [Pa]                      *
64  !                                                                            *
65  !*****************************************************************************
66
67  use par_mod
68  use com_mod
69  use cmapf_mod, only: cc2gll
70
71  implicit none
72
73  integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
74 integer :: rain_cloud_above,kz_inv !SE
75  integer icloudtop !PS
76  real :: f_qvsat,pressure
77 !real :: rh,lsp,convp
78  real :: rh,lsp,convp,prec,rhmin 
79  real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
80  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
81  real :: xlon,ylat,xlonr,dzdx,dzdy
82  real :: dzdx1,dzdx2,dzdy1,dzdy2, precmin
83  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
84  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
85  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
86  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
87  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
88  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
89  logical lconvectprec
90  real,parameter :: const=r_air/ga
91  parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
92
93  logical :: init = .true.
94
95
96  !*************************************************************************
97  ! If verttransform is called the first time, initialize heights of the   *
98  ! z levels in meter. The heights are the heights of model levels, where  *
99  ! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
100  ! the vertical resolution in the z system is doubled. As reference point,*
101  ! the lower left corner of the grid is used.                             *
102  ! Unlike in the eta system, no difference between heights for u,v and    *
103  ! heights for w exists.                                                  *
104  !*************************************************************************
105
106
107  !      do 897 kz=1,nuvz
108  !       write (*,*) 'akz: ',akz(kz),'bkz',bkz(kz)
109  !897     continue
110
111  if (init) then
112
113  ! Search for a point with high surface pressure (i.e. not above significant topography)
114  ! Then, use this point to construct a reference z profile, to be used at all times
115  !*****************************************************************************
116
117    do jy=0,nymin1
118      do ix=0,nxmin1
119        if (ps(ix,jy,1,n).gt.100000.) then
120          ixm=ix
121          jym=jy
122          goto 3
123        endif
124      end do
125    end do
1263   continue
127
128
129    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
130         ps(ixm,jym,1,n))
131    pold=ps(ixm,jym,1,n)
132    height(1)=0.
133
134    do kz=2,nuvz
135      pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
136      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
137
138
139  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
140  ! upon the transformation to z levels. In order to save computer memory, this is
141  ! not done anymore in the standard version. However, this option can still be
142  ! switched on by replacing the following lines with those below, that are
143  ! currently commented out.
144  ! Note that two more changes are necessary in this subroutine below.
145  ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
146  !*****************************************************************************
147
148      if (abs(tv-tvold).gt.0.2) then
149        height(kz)= &
150             height(kz-1)+const*log(pold/pint)* &
151             (tv-tvold)/log(tv/tvold)
152      else
153        height(kz)=height(kz-1)+ &
154             const*log(pold/pint)*tv
155      endif
156
157  ! Switch on following lines to use doubled vertical resolution
158  !*************************************************************
159  !    if (abs(tv-tvold).gt.0.2) then
160  !      height((kz-1)*2)=
161  !    +      height(max((kz-2)*2,1))+const*log(pold/pint)*
162  !    +      (tv-tvold)/log(tv/tvold)
163  !    else
164  !      height((kz-1)*2)=height(max((kz-2)*2,1))+
165  !    +      const*log(pold/pint)*tv
166  !    endif
167  ! End doubled vertical resolution
168
169      tvold=tv
170      pold=pint
171    end do
172
173
174  ! Switch on following lines to use doubled vertical resolution
175  !*************************************************************
176  !  do 7 kz=3,nz-1,2
177  !    height(kz)=0.5*(height(kz-1)+height(kz+1))
178  !  height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
179  ! End doubled vertical resolution
180
181
182  ! Determine highest levels that can be within PBL
183  !************************************************
184
185    do kz=1,nz
186      if (height(kz).gt.hmixmax) then
187        nmixz=kz
188        goto 9
189      endif
190    end do
1919   continue
192
193  ! Do not repeat initialization of the Cartesian z grid
194  !*****************************************************
195
196    init=.false.
197
198  endif
199
200
201  ! Loop over the whole grid
202  !*************************
203
204  do jy=0,nymin1
205    do ix=0,nxmin1
206      tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
207           ps(ix,jy,1,n))
208      pold=ps(ix,jy,1,n)
209      uvzlev(1)=0.
210      wzlev(1)=0.
211      rhoh(1)=pold/(r_air*tvold)
212
213
214  ! Compute heights of eta levels
215  !******************************
216
217      do kz=2,nuvz
218        pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
219        tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
220        rhoh(kz)=pint/(r_air*tv)
221
222        if (abs(tv-tvold).gt.0.2) then
223          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
224               (tv-tvold)/log(tv/tvold)
225        else
226          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
227        endif
228
229        tvold=tv
230        pold=pint
231      end do
232
233
234      do kz=2,nwz-1
235        wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
236      end do
237      wzlev(nwz)=wzlev(nwz-1)+ &
238           uvzlev(nuvz)-uvzlev(nuvz-1)
239
240      uvwzlev(ix,jy,1)=0.0
241      do kz=2,nuvz
242        uvwzlev(ix,jy,kz)=uvzlev(kz)
243      end do
244
245  ! Switch on following lines to use doubled vertical resolution
246  ! Switch off the three lines above.
247  !*************************************************************
248  !22          uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
249  !     do 23 kz=2,nwz
250  !23          uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
251  ! End doubled vertical resolution
252
253  ! pinmconv=(h2-h1)/(p2-p1)
254
255      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
256           ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
257           (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
258      do kz=2,nz-1
259        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
260             ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
261             (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
262      end do
263      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
264           ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
265           (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
266
267  ! Levels, where u,v,t and q are given
268  !************************************
269
270      uu(ix,jy,1,n)=uuh(ix,jy,1)
271      vv(ix,jy,1,n)=vvh(ix,jy,1)
272      tt(ix,jy,1,n)=tth(ix,jy,1,n)
273      qv(ix,jy,1,n)=qvh(ix,jy,1,n)
274      pv(ix,jy,1,n)=pvh(ix,jy,1)
275      rho(ix,jy,1,n)=rhoh(1)
276      uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
277      vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
278      tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
279      qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
280      pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
281      rho(ix,jy,nz,n)=rhoh(nuvz)
282      kmin=2
283      do iz=2,nz-1
284        do kz=kmin,nuvz
285          if(height(iz).gt.uvzlev(nuvz)) then
286            uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
287            vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
288            tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
289            qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
290            pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
291            rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
292            goto 30
293          endif
294          if ((height(iz).gt.uvzlev(kz-1)).and. &
295               (height(iz).le.uvzlev(kz))) then
296           dz1=height(iz)-uvzlev(kz-1)
297           dz2=uvzlev(kz)-height(iz)
298           dz=dz1+dz2
299           uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
300           vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
301           tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
302                +tth(ix,jy,kz,n)*dz1)/dz
303           qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
304                +qvh(ix,jy,kz,n)*dz1)/dz
305           pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
306           rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
307           kmin=kz
308           goto 30
309          endif
310        end do
31130      continue
312      end do
313
314
315  ! Levels, where w is given
316  !*************************
317
318      ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
319      ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
320      kmin=2
321      do iz=2,nz
322        do kz=kmin,nwz
323          if ((height(iz).gt.wzlev(kz-1)).and. &
324               (height(iz).le.wzlev(kz))) then
325           dz1=height(iz)-wzlev(kz-1)
326           dz2=wzlev(kz)-height(iz)
327           dz=dz1+dz2
328           ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
329                +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
330           kmin=kz
331           goto 40
332          endif
333        end do
33440      continue
335      end do
336
337  ! Compute density gradients at intermediate levels
338  !*************************************************
339
340      drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
341           (height(2)-height(1))
342      do kz=2,nz-1
343        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
344             (height(kz+1)-height(kz-1))
345      end do
346      drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
347
348    end do
349  end do
350
351
352  !****************************************************************
353  ! Compute slope of eta levels in windward direction and resulting
354  ! vertical wind correction
355  !****************************************************************
356
357  do jy=1,ny-2
358    do ix=1,nx-2
359
360      kmin=2
361      do iz=2,nz-1
362
363        ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
364        vi=vv(ix,jy,iz,n)*dyconst
365
366        do kz=kmin,nz
367          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
368               (height(iz).le.uvwzlev(ix,jy,kz))) then
369            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
370            dz2=uvwzlev(ix,jy,kz)-height(iz)
371            dz=dz1+dz2
372            kl=kz-1
373            klp=kz
374            kmin=kz
375            goto 47
376          endif
377        end do
378
37947      ix1=ix-1
380        jy1=jy-1
381        ixp=ix+1
382        jyp=jy+1
383
384        dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
385        dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
386        dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
387
388        dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
389        dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
390        dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
391
392        ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*ui+dzdy*vi)
393
394      end do
395
396    end do
397  end do
398
399
400  ! If north pole is in the domain, calculate wind velocities in polar
401  ! stereographic coordinates
402  !*******************************************************************
403
404  if (nglobal) then
405    do jy=int(switchnorthg)-2,nymin1
406      ylat=ylat0+real(jy)*dy
407      do ix=0,nxmin1
408        xlon=xlon0+real(ix)*dx
409        do iz=1,nz
410          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
411               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
412               vvpol(ix,jy,iz,n))
413        end do
414      end do
415    end do
416
417
418    do iz=1,nz
419
420  ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
421  !
422  !   AMSnauffer Nov 18 2004 Added check for case vv=0
423  !
424      xlon=xlon0+real(nx/2-1)*dx
425      xlonr=xlon*pi/180.
426      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
427           vv(nx/2-1,nymin1,iz,n)**2)
428      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
429        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
430             vv(nx/2-1,nymin1,iz,n))-xlonr
431      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
432        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
433             vv(nx/2-1,nymin1,iz,n))-xlonr
434      else
435        ddpol=pi/2-xlonr
436      endif
437      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
438      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
439
440  ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
441      xlon=180.0
442      xlonr=xlon*pi/180.
443      ylat=90.0
444      uuaux=-ffpol*sin(xlonr+ddpol)
445      vvaux=-ffpol*cos(xlonr+ddpol)
446      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
447           vvpolaux)
448
449      jy=nymin1
450      do ix=0,nxmin1
451        uupol(ix,jy,iz,n)=uupolaux
452        vvpol(ix,jy,iz,n)=vvpolaux
453      end do
454    end do
455
456
457  ! Fix: Set W at pole to the zonally averaged W of the next equator-
458  ! ward parallel of latitude
459
460  do iz=1,nz
461      wdummy=0.
462      jy=ny-2
463      do ix=0,nxmin1
464        wdummy=wdummy+ww(ix,jy,iz,n)
465      end do
466      wdummy=wdummy/real(nx)
467      jy=nymin1
468      do ix=0,nxmin1
469        ww(ix,jy,iz,n)=wdummy
470      end do
471  end do
472
473  endif
474
475
476  ! If south pole is in the domain, calculate wind velocities in polar
477  ! stereographic coordinates
478  !*******************************************************************
479
480  if (sglobal) then
481    do jy=0,int(switchsouthg)+3
482      ylat=ylat0+real(jy)*dy
483      do ix=0,nxmin1
484        xlon=xlon0+real(ix)*dx
485        do iz=1,nz
486          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
487               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
488               vvpol(ix,jy,iz,n))
489        end do
490      end do
491    end do
492
493    do iz=1,nz
494
495  ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
496  !
497  !   AMSnauffer Nov 18 2004 Added check for case vv=0
498  !
499      xlon=xlon0+real(nx/2-1)*dx
500      xlonr=xlon*pi/180.
501      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
502           vv(nx/2-1,0,iz,n)**2)
503      if (vv(nx/2-1,0,iz,n).lt.0.) then
504        ddpol=atan(uu(nx/2-1,0,iz,n)/ &
505             vv(nx/2-1,0,iz,n))+xlonr
506      else if (vv(nx/2-1,0,iz,n).gt.0.) then
507        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
508             vv(nx/2-1,0,iz,n))+xlonr
509      else
510        ddpol=pi/2-xlonr
511      endif
512      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
513      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
514
515  ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
516      xlon=180.0
517      xlonr=xlon*pi/180.
518      ylat=-90.0
519      uuaux=+ffpol*sin(xlonr-ddpol)
520      vvaux=-ffpol*cos(xlonr-ddpol)
521      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
522           vvpolaux)
523
524      jy=0
525      do ix=0,nxmin1
526        uupol(ix,jy,iz,n)=uupolaux
527        vvpol(ix,jy,iz,n)=vvpolaux
528      end do
529    end do
530
531
532  ! Fix: Set W at pole to the zonally averaged W of the next equator-
533  ! ward parallel of latitude
534
535    do iz=1,nz
536      wdummy=0.
537      jy=1
538      do ix=0,nxmin1
539        wdummy=wdummy+ww(ix,jy,iz,n)
540      end do
541      wdummy=wdummy/real(nx)
542      jy=0
543      do ix=0,nxmin1
544        ww(ix,jy,iz,n)=wdummy
545      end do
546    end do
547  endif
548
549
550  !write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
551  !   create a cloud and rainout/washout field, clouds occur where rh>80%
552  !   total cloudheight is stored at level 0
553
554
555
556  do jy=0,nymin1
557    do ix=0,nxmin1
558
559
560
561  !    rain_cloud_above=0
562  !    lsp=lsprec(ix,jy,1,n)
563  !    convp=convprec(ix,jy,1,n)
564  !    cloudsh(ix,jy,n)=0
565  !    do kz_inv=1,nz-1
566  !       kz=nz-kz_inv+1
567  !       pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
568  !       rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
569  !       clouds(ix,jy,kz,n)=0
570  !       if (rh.gt.0.8) then ! in cloud
571  !          if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
572  !             rain_cloud_above=1
573  !             cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
574  !                  height(kz)-height(kz-1)
575  !             if (lsp.ge.convp) then
576  !                clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
577  !             else
578  !                clouds(ix,jy,kz,n)=2 ! convp dominated rainout
579  !             endif
580  !          else ! no precipitation
581  !                clouds(ix,jy,kz,n)=1 ! cloud
582  !          endif
583  !       else ! no cloud
584  !          if (rain_cloud_above.eq.1) then ! scavenging
585  !             if (lsp.ge.convp) then
586  !                clouds(ix,jy,kz,n)=5 ! lsp dominated washout
587  !             else
588  !                clouds(ix,jy,kz,n)=4 ! convp dominated washout
589  !             endif
590  !          endif
591  !       endif
592  !    end do
593
594
595   ! PS 3012
596
597             lsp=lsprec(ix,jy,1,n)
598          convp=convprec(ix,jy,1,n)
599          prec=lsp+convp
600          if (lsp.gt.convp) then !  prectype='lsp'
601            lconvectprec = .false.
602          else ! prectype='cp '
603            lconvectprec = .true.
604          endif
605          rhmin = 0.90 ! standard condition for presence of clouds
606!PS       note that original by Sabine Eckhart was 80%
607!PS       however, for T<-20 C we consider saturation over ice
608!PS       so I think 90% should be enough         
609          icloudbot(ix,jy,n)=icmv
610          icloudtop=icmv ! this is just a local variable
61198        do kz=1,nz
612            pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
613            rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
614!ps            if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
615            if (rh .gt. rhmin) then
616              if (icloudbot(ix,jy,n) .eq. icmv) then
617                icloudbot(ix,jy,n)=nint(height(kz))
618              endif
619              icloudtop=nint(height(kz)) ! use int to save memory
620            endif
621          enddo
622
623!PS try to get a cloud thicker than 50 m
624!PS if there is at least .01 mm/h  - changed to 0.002 and put into
625!PS parameter precpmin       
626          if ((icloudbot(ix,jy,n) .eq. icmv .or. &
627              icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
628              prec .gt. precmin) then
629            rhmin = rhmin - 0.05
630            if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
631          endif
632!PS implement a rough fix for badly represented convection
633!PS is based on looking at a limited set of comparison data
634          if (lconvectprec .and. icloudtop .lt. 6000 .and. &
635             prec .gt. precmin) then
636
637            if (convp .lt. 0.1) then
638              icloudbot(ix,jy,n) = 500
639              icloudtop =         8000
640            else
641              icloudbot(ix,jy,n) = 0
642              icloudtop =      10000
643            endif
644          endif
645          if (icloudtop .ne. icmv) then
646            icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
647          else
648            icloudthck(ix,jy,n) = icmv
649          endif
650!PS  get rid of too thin clouds     
651          if (icloudthck(ix,jy,n) .lt. 50) then
652            icloudbot(ix,jy,n)=icmv
653            icloudthck(ix,jy,n)=icmv
654          endif
655
656
657    end do
658  end do
659
660  !do 102 kz=1,nuvz
661  !write(an,'(i02)') kz+10
662  !write(*,*) nuvz,nymin1,nxmin1,'--',an,'--'
663  !open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')
664  !do 101 jy=0,nymin1
665  !    write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)
666  !101   continue
667  ! close(4)
668  !102   continue
669
670  ! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')
671  ! do 103 jy=0,nymin1
672  !     write (4,*)
673  !+       (height(kz),kz=1,nuvz)
674  !103   continue
675  ! close(4)
676
677  !open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')
678  ! do 104 jy=0,nymin1
679  !     write (4,*)
680  !+       (r_air*tt(ix,jy,1,n)*rho(ix,jy,1,n),ix=1,nxmin1)
681  !104   continue
682  ! close(4)
683
684
685end subroutine verttransform
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG