source: flexpart.git/src/verttransform.f90 @ e200b7a

10.4.1_peseiFPv9.3.1FPv9.3.1b_testingFPv9.3.2GFS_025NetCDFbugfixes+enhancementsdepositiondevflexpart-noresmfp9.3.1-20161214-nc4grib2nc4_repairinputlistlaptoprelease-10release-10.4.1scaling-bugsvn-petrasvn-trunkunivie
Last change on this file since e200b7a was e200b7a, checked in by Matthias Langer <matthias.langer@…>, 11 years ago

git-svn-id: http://flexpart.flexpart.eu:8088/svn/FlexPart90/trunk@4 ef8cc7e1-21b7-489e-abab-c1baa636049d

  • Property mode set to 100644
File size: 20.2 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  !*****************************************************************************
52  !                                                                            *
53  ! Variables:                                                                 *
54  ! nx,ny,nz                        field dimensions in x,y and z direction    *
55  ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition           *
56  ! uu(0:nxmax,0:nymax,nzmax,2)     wind components in x-direction [m/s]       *
57  ! vv(0:nxmax,0:nymax,nzmax,2)     wind components in y-direction [m/s]       *
58  ! ww(0:nxmax,0:nymax,nzmax,2)     wind components in z-direction [deltaeta/s]*
59  ! tt(0:nxmax,0:nymax,nzmax,2)     temperature [K]                            *
60  ! pv(0:nxmax,0:nymax,nzmax,2)     potential voriticity (pvu)                 *
61  ! ps(0:nxmax,0:nymax,2)           surface pressure [Pa]                      *
62  !                                                                            *
63  !*****************************************************************************
64
65  use par_mod
66  use com_mod
67  use cmapf_mod, only: cc2gll
68
69  implicit none
70
71  integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
72  integer :: rain_cloud_above,kz_inv
73  real :: f_qvsat,pressure
74  real :: rh,lsp,convp
75  real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
76  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
77  real :: xlon,ylat,xlonr,dzdx,dzdy
78  real :: dzdx1,dzdx2,dzdy1,dzdy2
79  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
80  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
81  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
82  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
83  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
84  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
85  real,parameter :: const=r_air/ga
86
87  logical :: init = .true.
88
89
90  !*************************************************************************
91  ! If verttransform is called the first time, initialize heights of the   *
92  ! z levels in meter. The heights are the heights of model levels, where  *
93  ! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
94  ! the vertical resolution in the z system is doubled. As reference point,*
95  ! the lower left corner of the grid is used.                             *
96  ! Unlike in the eta system, no difference between heights for u,v and    *
97  ! heights for w exists.                                                  *
98  !*************************************************************************
99
100
101  !      do 897 kz=1,nuvz
102  !       write (*,*) 'akz: ',akz(kz),'bkz',bkz(kz)
103  !897     continue
104
105  if (init) then
106
107  ! Search for a point with high surface pressure (i.e. not above significant topography)
108  ! Then, use this point to construct a reference z profile, to be used at all times
109  !*****************************************************************************
110
111    do jy=0,nymin1
112      do ix=0,nxmin1
113        if (ps(ix,jy,1,n).gt.100000.) then
114          ixm=ix
115          jym=jy
116          goto 3
117        endif
118      end do
119    end do
1203   continue
121
122
123    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
124         ps(ixm,jym,1,n))
125    pold=ps(ixm,jym,1,n)
126    height(1)=0.
127
128    do kz=2,nuvz
129      pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
130      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
131
132
133  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
134  ! upon the transformation to z levels. In order to save computer memory, this is
135  ! not done anymore in the standard version. However, this option can still be
136  ! switched on by replacing the following lines with those below, that are
137  ! currently commented out.
138  ! Note that two more changes are necessary in this subroutine below.
139  ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
140  !*****************************************************************************
141
142      if (abs(tv-tvold).gt.0.2) then
143        height(kz)= &
144             height(kz-1)+const*log(pold/pint)* &
145             (tv-tvold)/log(tv/tvold)
146      else
147        height(kz)=height(kz-1)+ &
148             const*log(pold/pint)*tv
149      endif
150
151  ! Switch on following lines to use doubled vertical resolution
152  !*************************************************************
153  !    if (abs(tv-tvold).gt.0.2) then
154  !      height((kz-1)*2)=
155  !    +      height(max((kz-2)*2,1))+const*log(pold/pint)*
156  !    +      (tv-tvold)/log(tv/tvold)
157  !    else
158  !      height((kz-1)*2)=height(max((kz-2)*2,1))+
159  !    +      const*log(pold/pint)*tv
160  !    endif
161  ! End doubled vertical resolution
162
163      tvold=tv
164      pold=pint
165    end do
166
167
168  ! Switch on following lines to use doubled vertical resolution
169  !*************************************************************
170  !  do 7 kz=3,nz-1,2
171  !    height(kz)=0.5*(height(kz-1)+height(kz+1))
172  !  height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
173  ! End doubled vertical resolution
174
175
176  ! Determine highest levels that can be within PBL
177  !************************************************
178
179    do kz=1,nz
180      if (height(kz).gt.hmixmax) then
181        nmixz=kz
182        goto 9
183      endif
184    end do
1859   continue
186
187  ! Do not repeat initialization of the Cartesian z grid
188  !*****************************************************
189
190    init=.false.
191
192  endif
193
194
195  ! Loop over the whole grid
196  !*************************
197
198  do jy=0,nymin1
199    do ix=0,nxmin1
200      tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
201           ps(ix,jy,1,n))
202      pold=ps(ix,jy,1,n)
203      uvzlev(1)=0.
204      wzlev(1)=0.
205      rhoh(1)=pold/(r_air*tvold)
206
207
208  ! Compute heights of eta levels
209  !******************************
210
211      do kz=2,nuvz
212        pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
213        tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
214        rhoh(kz)=pint/(r_air*tv)
215
216        if (abs(tv-tvold).gt.0.2) then
217          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
218               (tv-tvold)/log(tv/tvold)
219        else
220          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
221        endif
222
223        tvold=tv
224        pold=pint
225      end do
226
227
228      do kz=2,nwz-1
229        wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
230      end do
231      wzlev(nwz)=wzlev(nwz-1)+ &
232           uvzlev(nuvz)-uvzlev(nuvz-1)
233
234      uvwzlev(ix,jy,1)=0.0
235      do kz=2,nuvz
236        uvwzlev(ix,jy,kz)=uvzlev(kz)
237      end do
238
239  ! Switch on following lines to use doubled vertical resolution
240  ! Switch off the three lines above.
241  !*************************************************************
242  !22          uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
243  !     do 23 kz=2,nwz
244  !23          uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
245  ! End doubled vertical resolution
246
247  ! pinmconv=(h2-h1)/(p2-p1)
248
249      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
250           ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
251           (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
252      do kz=2,nz-1
253        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
254             ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
255             (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
256      end do
257      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
258           ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
259           (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
260
261  ! Levels, where u,v,t and q are given
262  !************************************
263
264      uu(ix,jy,1,n)=uuh(ix,jy,1)
265      vv(ix,jy,1,n)=vvh(ix,jy,1)
266      tt(ix,jy,1,n)=tth(ix,jy,1,n)
267      qv(ix,jy,1,n)=qvh(ix,jy,1,n)
268      pv(ix,jy,1,n)=pvh(ix,jy,1)
269      rho(ix,jy,1,n)=rhoh(1)
270      uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
271      vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
272      tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
273      qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
274      pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
275      rho(ix,jy,nz,n)=rhoh(nuvz)
276      kmin=2
277      do iz=2,nz-1
278        do kz=kmin,nuvz
279          if(height(iz).gt.uvzlev(nuvz)) then
280            uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
281            vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
282            tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
283            qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
284            pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
285            rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
286            goto 30
287          endif
288          if ((height(iz).gt.uvzlev(kz-1)).and. &
289               (height(iz).le.uvzlev(kz))) then
290           dz1=height(iz)-uvzlev(kz-1)
291           dz2=uvzlev(kz)-height(iz)
292           dz=dz1+dz2
293           uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
294           vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
295           tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
296                +tth(ix,jy,kz,n)*dz1)/dz
297           qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
298                +qvh(ix,jy,kz,n)*dz1)/dz
299           pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
300           rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
301           kmin=kz
302           goto 30
303          endif
304        end do
30530      continue
306      end do
307
308
309  ! Levels, where w is given
310  !*************************
311
312      ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
313      ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
314      kmin=2
315      do iz=2,nz
316        do kz=kmin,nwz
317          if ((height(iz).gt.wzlev(kz-1)).and. &
318               (height(iz).le.wzlev(kz))) then
319           dz1=height(iz)-wzlev(kz-1)
320           dz2=wzlev(kz)-height(iz)
321           dz=dz1+dz2
322           ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
323                +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
324           kmin=kz
325           goto 40
326          endif
327        end do
32840      continue
329      end do
330
331  ! Compute density gradients at intermediate levels
332  !*************************************************
333
334      drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
335           (height(2)-height(1))
336      do kz=2,nz-1
337        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
338             (height(kz+1)-height(kz-1))
339      end do
340      drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
341
342    end do
343  end do
344
345
346  !****************************************************************
347  ! Compute slope of eta levels in windward direction and resulting
348  ! vertical wind correction
349  !****************************************************************
350
351  do jy=1,ny-2
352    do ix=1,nx-2
353
354      kmin=2
355      do iz=2,nz-1
356
357        ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
358        vi=vv(ix,jy,iz,n)*dyconst
359
360        do kz=kmin,nz
361          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
362               (height(iz).le.uvwzlev(ix,jy,kz))) then
363            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
364            dz2=uvwzlev(ix,jy,kz)-height(iz)
365            dz=dz1+dz2
366            kl=kz-1
367            klp=kz
368            kmin=kz
369            goto 47
370          endif
371        end do
372
37347      ix1=ix-1
374        jy1=jy-1
375        ixp=ix+1
376        jyp=jy+1
377
378        dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
379        dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
380        dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
381
382        dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
383        dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
384        dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
385
386        ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*ui+dzdy*vi)
387
388      end do
389
390    end do
391  end do
392
393
394  ! If north pole is in the domain, calculate wind velocities in polar
395  ! stereographic coordinates
396  !*******************************************************************
397
398  if (nglobal) then
399    do jy=int(switchnorthg)-2,nymin1
400      ylat=ylat0+real(jy)*dy
401      do ix=0,nxmin1
402        xlon=xlon0+real(ix)*dx
403        do iz=1,nz
404          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
405               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
406               vvpol(ix,jy,iz,n))
407        end do
408      end do
409    end do
410
411
412    do iz=1,nz
413
414  ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
415  !
416  !   AMSnauffer Nov 18 2004 Added check for case vv=0
417  !
418      xlon=xlon0+real(nx/2-1)*dx
419      xlonr=xlon*pi/180.
420      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
421           vv(nx/2-1,nymin1,iz,n)**2)
422      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
423        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
424             vv(nx/2-1,nymin1,iz,n))-xlonr
425      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
426        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
427             vv(nx/2-1,nymin1,iz,n))-xlonr
428      else
429        ddpol=pi/2-xlonr
430      endif
431      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
432      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
433
434  ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
435      xlon=180.0
436      xlonr=xlon*pi/180.
437      ylat=90.0
438      uuaux=-ffpol*sin(xlonr+ddpol)
439      vvaux=-ffpol*cos(xlonr+ddpol)
440      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
441           vvpolaux)
442
443      jy=nymin1
444      do ix=0,nxmin1
445        uupol(ix,jy,iz,n)=uupolaux
446        vvpol(ix,jy,iz,n)=vvpolaux
447      end do
448    end do
449
450
451  ! Fix: Set W at pole to the zonally averaged W of the next equator-
452  ! ward parallel of latitude
453
454  do iz=1,nz
455      wdummy=0.
456      jy=ny-2
457      do ix=0,nxmin1
458        wdummy=wdummy+ww(ix,jy,iz,n)
459      end do
460      wdummy=wdummy/real(nx)
461      jy=nymin1
462      do ix=0,nxmin1
463        ww(ix,jy,iz,n)=wdummy
464      end do
465  end do
466
467  endif
468
469
470  ! If south pole is in the domain, calculate wind velocities in polar
471  ! stereographic coordinates
472  !*******************************************************************
473
474  if (sglobal) then
475    do jy=0,int(switchsouthg)+3
476      ylat=ylat0+real(jy)*dy
477      do ix=0,nxmin1
478        xlon=xlon0+real(ix)*dx
479        do iz=1,nz
480          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
481               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
482               vvpol(ix,jy,iz,n))
483        end do
484      end do
485    end do
486
487    do iz=1,nz
488
489  ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
490  !
491  !   AMSnauffer Nov 18 2004 Added check for case vv=0
492  !
493      xlon=xlon0+real(nx/2-1)*dx
494      xlonr=xlon*pi/180.
495      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
496           vv(nx/2-1,0,iz,n)**2)
497      if (vv(nx/2-1,0,iz,n).lt.0.) then
498        ddpol=atan(uu(nx/2-1,0,iz,n)/ &
499             vv(nx/2-1,0,iz,n))+xlonr
500      else if (vv(nx/2-1,0,iz,n).gt.0.) then
501        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
502             vv(nx/2-1,0,iz,n))+xlonr
503      else
504        ddpol=pi/2-xlonr
505      endif
506      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
507      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
508
509  ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
510      xlon=180.0
511      xlonr=xlon*pi/180.
512      ylat=-90.0
513      uuaux=+ffpol*sin(xlonr-ddpol)
514      vvaux=-ffpol*cos(xlonr-ddpol)
515      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
516           vvpolaux)
517
518      jy=0
519      do ix=0,nxmin1
520        uupol(ix,jy,iz,n)=uupolaux
521        vvpol(ix,jy,iz,n)=vvpolaux
522      end do
523    end do
524
525
526  ! Fix: Set W at pole to the zonally averaged W of the next equator-
527  ! ward parallel of latitude
528
529    do iz=1,nz
530      wdummy=0.
531      jy=1
532      do ix=0,nxmin1
533        wdummy=wdummy+ww(ix,jy,iz,n)
534      end do
535      wdummy=wdummy/real(nx)
536      jy=0
537      do ix=0,nxmin1
538        ww(ix,jy,iz,n)=wdummy
539      end do
540    end do
541  endif
542
543
544  !write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
545  !   create a cloud and rainout/washout field, clouds occur where rh>80%
546  !   total cloudheight is stored at level 0
547  do jy=0,nymin1
548    do ix=0,nxmin1
549      rain_cloud_above=0
550      lsp=lsprec(ix,jy,1,n)
551      convp=convprec(ix,jy,1,n)
552      cloudsh(ix,jy,n)=0
553      do kz_inv=1,nz-1
554         kz=nz-kz_inv+1
555         pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
556         rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
557         clouds(ix,jy,kz,n)=0
558         if (rh.gt.0.8) then ! in cloud
559            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
560               rain_cloud_above=1
561               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
562                    height(kz)-height(kz-1)
563               if (lsp.ge.convp) then
564                  clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
565               else
566                  clouds(ix,jy,kz,n)=2 ! convp dominated rainout
567               endif
568            else ! no precipitation
569                  clouds(ix,jy,kz,n)=1 ! cloud
570            endif
571         else ! no cloud
572            if (rain_cloud_above.eq.1) then ! scavenging
573               if (lsp.ge.convp) then
574                  clouds(ix,jy,kz,n)=5 ! lsp dominated washout
575               else
576                  clouds(ix,jy,kz,n)=4 ! convp dominated washout
577               endif
578            endif
579         endif
580      end do
581    end do
582  end do
583
584  !do 102 kz=1,nuvz
585  !write(an,'(i02)') kz+10
586  !write(*,*) nuvz,nymin1,nxmin1,'--',an,'--'
587  !open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')
588  !do 101 jy=0,nymin1
589  !    write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)
590  !101   continue
591  ! close(4)
592  !102   continue
593
594  ! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')
595  ! do 103 jy=0,nymin1
596  !     write (4,*)
597  !+       (height(kz),kz=1,nuvz)
598  !103   continue
599  ! close(4)
600
601  !open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')
602  ! do 104 jy=0,nymin1
603  !     write (4,*)
604  !+       (r_air*tt(ix,jy,1,n)*rho(ix,jy,1,n),ix=1,nxmin1)
605  !104   continue
606  ! close(4)
607end subroutine verttransform
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG