Index: src/readwind_gfs.f90
===================================================================
--- src/readwind_gfs.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
+++ src/readwind_gfs.f90	(revision db91eb78afeb1bb8bafe66068877db8841859062)
@@ -1,3 +1,3 @@
-!**********************************************************************
+
 ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
@@ -108,4 +108,5 @@
   character(len=24) :: gribErrorMsg = 'Error reading grib file'
   character(len=20) :: gribFunction = 'readwind_gfs'
+  character(len=20) :: shortname
 
 
@@ -120,5 +121,5 @@
 
   !HSO
-5   call grib_open_file(ifile,path(3)(1:length(3)) &
+  call grib_open_file(ifile,path(3)(1:length(3)) &
          //trim(wfname(indj)),'r',iret)
   if (iret.ne.GRIB_SUCCESS) then
@@ -162,4 +163,6 @@
 
   !read the grib2 identifiers
+  call grib_get_string(igrib,'shortName',shortname,iret)
+
   call grib_get_int(igrib,'discipline',discipl,iret)
 !  call grib_check(iret,gribFunction,gribErrorMsg)
@@ -173,5 +176,6 @@
        valSurf,iret)
 !  call grib_check(iret,gribFunction,gribErrorMsg)
-
+  
+!  write(*,*) 'Field: ',ifield,parCat,parNum,typSurf,shortname
   !convert to grib1 identifiers
   isec1(6)=-1
@@ -214,4 +218,8 @@
     isec1(7)=105         ! indicatorOfTypeOfLevel
     isec1(8)=10
+  elseif ((parCat.eq.1).and.(parNum.eq.22).and.(typSurf.eq.100)) then ! CLWMR Cloud Mixing Ratio [kg/kg]:
+    isec1(6)=153         ! indicatorOfParameter
+    isec1(7)=100         ! indicatorOfTypeOfLevel
+    isec1(8)=valSurf/100 ! level, convert to hPa
   elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
     isec1(6)=2           ! indicatorOfParameter
@@ -548,4 +556,11 @@
         endif
       endif
+! SEC & IP 12/2018 read GFS clouds
+      if(isec1(6).eq.153) then  !! CLWCR  Cloud liquid water content [kg/kg]
+        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
+        readclouds=.true.
+        sumclouds=.true.
+      endif
+
 
     end do
@@ -675,4 +690,6 @@
     call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
     call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
+! IP & SEC adding GFS Clouds 20181205
+    call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,2,n)
   endif
 
Index: src/verttransform_gfs.f90
===================================================================
--- src/verttransform_gfs.f90	(revision 6ecb30a41fb8a07e153abadc712699a09c15f90c)
+++ src/verttransform_gfs.f90	(revision db91eb78afeb1bb8bafe66068877db8841859062)
@@ -75,5 +75,5 @@
   integer :: rain_cloud_above,kz_inv
   real :: f_qvsat,pressure
-  real :: rh,lsp,convp
+  real :: rh,lsp,cloudh_min,convp,prec
   real :: rhoh(nuvzmax),pinmconv(nzmax)
   real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
@@ -224,4 +224,8 @@
       tt(ix,jy,1,n)=tth(ix,jy,llev,n)
       qv(ix,jy,1,n)=qvh(ix,jy,llev,n)
+! IP & SEC, 201812 add clouds
+      if (readclouds) then
+         clwc(ix,jy,1,n)=clwch(ix,jy,llev,n)
+      endif 
       pv(ix,jy,1,n)=pvh(ix,jy,llev)
       rho(ix,jy,1,n)=rhoh(llev)
@@ -231,4 +235,8 @@
       tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
       qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
+! IP & SEC, 201812 add clouds
+      if (readclouds) then
+         clwc(ix,jy,nz,n)=clwch(ix,jy,nuvz,n)
+      endif
       pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
       rho(ix,jy,nz,n)=rhoh(nuvz)
@@ -242,4 +250,8 @@
             tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
             qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
+! IP & SEC, 201812 add clouds
+            if (readclouds) then
+               clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
+            endif
             pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
             rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
@@ -258,4 +270,9 @@
             qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
             +qvh(ix,jy,kz,n)*dz1)/dz
+! IP & SEC, 201812 add clouds
+            if (readclouds) then
+               clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2 &
+               +clwch(ix,jy,kz,n)*dz1)/dz
+            endif
             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
@@ -496,4 +513,65 @@
 
 
+
+!***********************************************************************************
+! IP & SEC, 201812 GFS clouds read
+  if (readclouds) then
+! The method is loops all grids vertically and constructs the 3D matrix for clouds
+! Cloud top and cloud bottom gid cells are assigned as well as the total column
+! cloud water. For precipitating grids, the type and whether it is in or below
+! cloud scavenging are assigned with numbers 2-5 (following the old metod).
+! Distinction is done for lsp and convp though they are treated the same in regards
+! to scavenging. Also clouds that are not precipitating are defined which may be
+! to include future cloud processing by non-precipitating-clouds.
+!***********************************************************************************
+    write(*,*) 'Global NCEP fields: using cloud water'
+    clw(:,:,:,n)=0.0
+    ctwc(:,:,n)=0.0
+    clouds(:,:,:,n)=0
+! If water/ice are read separately into clwc and ciwc, store sum in clwc
+    do jy=0,nymin1
+      do ix=0,nxmin1
+        lsp=lsprec(ix,jy,1,n)
+        convp=convprec(ix,jy,1,n)
+        prec=lsp+convp
+! Find clouds in the vertical
+        do kz=1, nz-1 !go from top to bottom
+          if (clwc(ix,jy,kz,n).gt.0) then
+! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
+            clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
+            ctwc(ix,jy,n) = ctwc(ix,jy,n)+clw(ix,jy,kz,n)
+            cloudh_min=min(height(kz+1),height(kz))
+          endif
+        end do
+
+! If Precipitation. Define removal type in the vertical
+        if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+
+          do kz=nz,1,-1 !go Bottom up!
+            if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
+              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
+              clouds(ix,jy,kz,n)=1                               ! is a cloud
+              if (lsp.ge.convp) then
+                clouds(ix,jy,kz,n)=3                            ! lsp in-cloud
+              else
+                clouds(ix,jy,kz,n)=2                             ! convp in-cloud
+              endif                                              ! convective or large scale
+            elseif((clw(ix,jy,kz,n).le.0) .and. (cloudh_min.ge.height(kz))) then ! is below cloud
+              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                                              ! convective or large scale
+            endif
+
+            if (height(kz).ge. 19000) then                        ! set a max height for removal
+              clouds(ix,jy,kz,n)=0
+            endif !clw>0
+          end do !nz
+        endif ! precipitation
+      end do
+    end do
+  else
+  write(*,*) 'Global NCEP fields: using cloud water from Parameterization'
   !   write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
   !   create a cloud and rainout/washout field, clouds occur where rh>80%
@@ -534,4 +612,5 @@
     end do
   end do
+  endif  ! IP & SEC 201812, GFS clouds read
 
 
