source: branches/FLEXPART_9.1.3/src/verttransform_gfs.f90 @ 13

Last change on this file since 13 was 13, checked in by saeck, 9 years ago

update to wetdepo.f90

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