Index: trunk/src/verttransform_gfs.f90
===================================================================
--- trunk/src/verttransform_gfs.f90	(revision 4)
+++ trunk/src/verttransform_gfs.f90	(revision 24)
@@ -21,5 +21,5 @@
 
 subroutine verttransform(n,uuh,vvh,wwh,pvh)
-  !                         i  i   i   i   i
+  !                      i  i   i   i   i
   !*****************************************************************************
   !                                                                            *
@@ -71,8 +71,8 @@
   real :: f_qvsat,pressure
   real :: rh,lsp,convp
-  real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
+  real :: rhoh(nuvzmax),pinmconv(nzmax)
   real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
   real :: xlon,ylat,xlonr,dzdx,dzdy
-  real :: dzdx1,dzdx2,dzdy1,dzdy2
+  real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
   real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
   real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
@@ -92,9 +92,5 @@
   ! If verttransform is called the first time, initialize heights of the   *
   ! z levels in meter. The heights are the heights of model levels, where  *
-  ! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
-  ! the vertical resolution in the z system is doubled. As reference point,*
-  ! the lower left corner of the grid is used.                             *
-  ! Unlike in the eta system, no difference between heights for u,v and    *
-  ! heights for w exists.                                                  *
+  ! u,v,T and qv are given.                                                *
   !*************************************************************************
 
@@ -118,5 +114,5 @@
 
     tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
-         ps(ixm,jym,1,n))
+    ps(ixm,jym,1,n))
     pold=ps(ixm,jym,1,n)
     height(1)=0.
@@ -126,46 +122,14 @@
       tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
 
-
-  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
-  ! upon the transformation to z levels. In order to save computer memory, this is
-  ! not done anymore in the standard version. However, this option can still be
-  ! switched on by replacing the following lines with those below, that are
-  ! currently commented out.
-  ! Note that two more changes are necessary in this subroutine below.
-  ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
-  !*****************************************************************************
-
       if (abs(tv-tvold).gt.0.2) then
-        height(kz)= &
-             height(kz-1)+const*log(pold/pint)* &
-             (tv-tvold)/log(tv/tvold)
+        height(kz)=height(kz-1)+const*log(pold/pint)* &
+        (tv-tvold)/log(tv/tvold)
       else
-        height(kz)=height(kz-1)+ &
-             const*log(pold/pint)*tv
+        height(kz)=height(kz-1)+const*log(pold/pint)*tv
       endif
-
-  ! Switch on following lines to use doubled vertical resolution
-  !*************************************************************
-  !    if (abs(tv-tvold).gt.0.2) then
-  !      height((kz-1)*2)=
-  !    +      height(max((kz-2)*2,1))+const*log(pold/pint)*
-  !    +      (tv-tvold)/log(tv/tvold)
-  !    else
-  !      height((kz-1)*2)=height(max((kz-2)*2,1))+
-  !    +      const*log(pold/pint)*tv
-  !    endif
-  ! End doubled vertical resolution
 
       tvold=tv
       pold=pint
     end do
-
-
-  ! Switch on following lines to use doubled vertical resolution
-  !*************************************************************
-  !  do 7 kz=3,nz-1,2
-  !    height(kz)=0.5*(height(kz-1)+height(kz+1))
-  !  height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
-  ! End doubled vertical resolution
 
 
@@ -198,5 +162,5 @@
       llev = 0
       do i=1,nuvz
-       if (ps(ix,jy,1,n).lt.akz(i)) llev=i
+        if (ps(ix,jy,1,n).lt.akz(i)) llev=i
       end do
        llev = llev+1
@@ -212,5 +176,4 @@
       tvold=tth(ix,jy,llev,n)*(1.+0.608*qvh(ix,jy,llev,n))
       pold=akz(llev)
-      uvzlev(llev)=0.
       wzlev(llev)=0.
       uvwzlev(ix,jy,llev)=0.
@@ -223,24 +186,14 @@
 
         if (abs(tv-tvold).gt.0.2) then
-          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
-               (tv-tvold)/log(tv/tvold)
+          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
+          (tv-tvold)/log(tv/tvold)
         else
-          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
+          uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
         endif
-        wzlev(kz)=uvzlev(kz)
-        uvwzlev(ix,jy,kz)=uvzlev(kz)
+        wzlev(kz)=uvwzlev(ix,jy,kz)
 
         tvold=tv
         pold=pint
       end do
-
-
-  ! Switch on following lines to use doubled vertical resolution
-  ! Switch off the three lines above.
-  !*************************************************************
-  !22          uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
-  !     do 23 kz=2,nwz
-  !23          uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
-  ! End doubled vertical resolution
 
   ! pinmconv=(h2-h1)/(p2-p1)
@@ -279,5 +232,5 @@
       do iz=2,nz-1
         do kz=kmin,nuvz
-          if(height(iz).gt.uvzlev(nuvz)) then
+          if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
             uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
             vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
@@ -289,18 +242,18 @@
             goto 30
           endif
-          if ((height(iz).gt.uvzlev(kz-1)).and. &
-               (height(iz).le.uvzlev(kz))) then
-           dz1=height(iz)-uvzlev(kz-1)
-           dz2=uvzlev(kz)-height(iz)
-           dz=dz1+dz2
-           uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
-           vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
-           tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
-                +tth(ix,jy,kz,n)*dz1)/dz
-           qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
-                +qvh(ix,jy,kz,n)*dz1)/dz
-           pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
-           rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
-           pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz
+          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
+          (height(iz).le.uvwzlev(ix,jy,kz))) then
+            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
+            dz2=uvwzlev(ix,jy,kz)-height(iz)
+            dz=dz1+dz2
+            uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
+            vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
+            tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
+            +tth(ix,jy,kz,n)*dz1)/dz
+            qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
+            +qvh(ix,jy,kz,n)*dz1)/dz
+            pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
+            rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
+            pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz
           endif
         end do
@@ -318,11 +271,10 @@
         do kz=kmin,nwz
           if ((height(iz).gt.wzlev(kz-1)).and. &
-               (height(iz).le.wzlev(kz))) then
-           dz1=height(iz)-wzlev(kz-1)
-           dz2=wzlev(kz)-height(iz)
-           dz=dz1+dz2
-           ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
-                +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
-
+          (height(iz).le.wzlev(kz))) then
+            dz1=height(iz)-wzlev(kz-1)
+            dz2=wzlev(kz)-height(iz)
+            dz=dz1+dz2
+            ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
+            +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
           endif
         end do
@@ -337,5 +289,5 @@
       do kz=2,nz-1
         drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
-             (height(kz+1)-height(kz-1))
+        (height(kz+1)-height(kz-1))
       end do
       drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
@@ -351,4 +303,5 @@
 
   do jy=1,ny-2
+    cosf=cos((real(jy)*dy+ylat0)*pi180)
     do ix=1,nx-2
 
@@ -367,10 +320,10 @@
       do iz=2,nz-1
 
-        ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
+        ui=uu(ix,jy,iz,n)*dxconst/cosf
         vi=vv(ix,jy,iz,n)*dyconst
 
         do kz=kmin,nz
           if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
-               (height(iz).le.uvwzlev(ix,jy,kz))) then
+          (height(iz).le.uvwzlev(ix,jy,kz))) then
             dz1=height(iz)-uvwzlev(ix,jy,kz-1)
             dz2=uvwzlev(ix,jy,kz)-height(iz)
@@ -426,12 +379,10 @@
       xlon=xlon0+real(nx/2-1)*dx
       xlonr=xlon*pi/180.
-      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
-           vv(nx/2-1,nymin1,iz,n)**2)
-      if(vv(nx/2-1,nymin1,iz,n).lt.0.) then
-        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
-             vv(nx/2-1,nymin1,iz,n))-xlonr
+      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2)
+      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
+        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
       elseif (vv(nx/2-1,nymin1,iz,n).gt.0.) then
         ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
-             vv(nx/2-1,nymin1,iz,n))-xlonr
+        vv(nx/2-1,nymin1,iz,n))-xlonr
       else
         ddpol=pi/2-xlonr
@@ -446,7 +397,5 @@
       uuaux=-ffpol*sin(xlonr+ddpol)
       vvaux=-ffpol*cos(xlonr+ddpol)
-      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
-           vvpolaux)
-
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
       jy=nymin1
       do ix=0,nxmin1
@@ -460,5 +409,5 @@
   ! ward parallel of latitude
 
-  do iz=1,nz
+    do iz=1,nz
       wdummy=0.
       jy=ny-2
@@ -471,5 +420,5 @@
         ww(ix,jy,iz,n)=wdummy
       end do
-  end do
+    end do
 
   endif
@@ -487,6 +436,5 @@
         do iz=1,nz
           call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
-               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
-               vvpol(ix,jy,iz,n))
+          vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
         end do
       end do
@@ -498,12 +446,9 @@
       xlon=xlon0+real(nx/2-1)*dx
       xlonr=xlon*pi/180.
-      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
-           vv(nx/2-1,0,iz,n)**2)
+      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2)
       if(vv(nx/2-1,0,iz,n).lt.0.) then
-        ddpol=atan(uu(nx/2-1,0,iz,n)/ &
-             vv(nx/2-1,0,iz,n))+xlonr
+        ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
       elseif (vv(nx/2-1,0,iz,n).gt.0.) then
-        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
-             vv(nx/2-1,0,iz,n))-xlonr
+        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))-xlonr
       else
         ddpol=pi/2-xlonr
@@ -518,6 +463,5 @@
       uuaux=+ffpol*sin(xlonr-ddpol)
       vvaux=-ffpol*cos(xlonr-ddpol)
-      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
-           vvpolaux)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
 
       jy=0
@@ -562,24 +506,23 @@
          clouds(ix,jy,kz,n)=0
          if (rh.gt.0.8) then ! in cloud
-            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
-               rain_cloud_above=1
-               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
-                    height(kz)-height(kz-1)
-               if (lsp.ge.convp) then
-                  clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
-               else
-                  clouds(ix,jy,kz,n)=2 ! convp dominated rainout
-               endif
-            else ! no precipitation
-                  clouds(ix,jy,kz,n)=1 ! cloud
-            endif
+           if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+              rain_cloud_above=1
+              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
+              if (lsp.ge.convp) then
+                 clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
+              else
+                 clouds(ix,jy,kz,n)=2 ! convp dominated rainout
+              endif
+           else ! no precipitation
+             clouds(ix,jy,kz,n)=1 ! cloud
+           endif
          else ! no cloud
-            if (rain_cloud_above.eq.1) then ! scavenging
-               if (lsp.ge.convp) then
-                  clouds(ix,jy,kz,n)=5 ! lsp dominated washout
-               else
-                  clouds(ix,jy,kz,n)=4 ! convp dominated washout
-               endif
-            endif
+           if (rain_cloud_above.eq.1) then ! scavenging
+             if (lsp.ge.convp) then
+               clouds(ix,jy,kz,n)=5 ! lsp dominated washout
+             else
+               clouds(ix,jy,kz,n)=4 ! convp dominated washout
+             endif
+           endif
          endif
       end do
