source: trunk/src/verttransform.f90 @ 20

Last change on this file since 20 was 20, checked in by igpis, 10 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