diff --git a/src/gsi/obs_para.f90 b/src/gsi/obs_para.f90
index e560fafe1f..f9deb83d0b 100644
--- a/src/gsi/obs_para.f90
+++ b/src/gsi/obs_para.f90
@@ -15,7 +15,7 @@ subroutine obs_para(ndata,mype)
 !   2004-07-28  treadon - add only to module use, add intent in/out
 !   2004-11-19  derber - modify to eliminate file and change to use logical 
 !                        rather than weights
-!   2005-06-14  wu      - add OMI total ozone
+!   2005-06-14  wu      - add omi total ozone
 !   2005-09-08  derber - simplify data set handling
 !   2005-10-17  treadon - fix bug in defnition of mype_diaghdr
 !   2006-02-03  derber  - modify for new obs control and obs count
@@ -25,7 +25,7 @@ subroutine obs_para(ndata,mype)
 !   2008-09-08  lueken  - merged ed's changes into q1fy09 code
 !   2009-04-21  derber  - reformulate to remove communication
 !   2008-05-10  meunier - handle for lagrangian data
-!   2013-01-26  parrish - attempt fix for bug flagged by WCOSS debug compiler.
+!   2013-01-26  parrish - attempt fix for bug flagged by wcoss debug compiler.
 !                            Replace 
 !                              "call dislag(.....,nobs_s)"
 !                            with
@@ -33,10 +33,10 @@ subroutine obs_para(ndata,mype)
 !                           nobs_s is an array in current subroutine, but is a
 !                           scalar inside subroutine dislag.
 !   2014-10-03  carley  - add creation mpi subcommunicator needed for 
-!                          buddy check QC to distinguish among pe subdomains 
+!                          buddy check qc to distinguish among pe subdomains 
 !                          with and without obs (only for t obs and twodvar_regional 
 !                          at the moment) 
-!   2015-10-28  guo     - added ioid_s(:) to support split-observer GSI with a
+!   2015-10-28  guo     - added ioid_s(:) to support split-observer gsi with a
 !                         control vector grid different from the guess state
 !                         grid.
 !
@@ -104,13 +104,13 @@ subroutine obs_para(ndata,mype)
         if (dtype(is)=='lag') then    ! lagrangian data
            call dislag(ndata(is,1),mm1,lunout,obsfile_all(is),dtype(is),&
                 nobs_s) 
-        nsat1(is)= nobs_sub(mm1,is)
+           nsat1(is)= nobs_sub(mm1,is)
         else 
            obproc:do ii=1,npe
-             if(nobs_sub(ii,is) > 0)then
-               mype_diaghdr(is) = ii-1
-               exit obproc
-             end if
+              if(nobs_sub(ii,is) > 0)then
+                 mype_diaghdr(is) = ii-1
+                 exit obproc
+              end if
            end do obproc
                 
            if(nobs_sub(mm1,is) > zero)then    ! classical observations
@@ -125,7 +125,7 @@ subroutine obs_para(ndata,mype)
         end if
 
         ! Simple logic to organize which tasks do and do not have obs. 
-        !  Needed for buddy check QC.   
+        !  Needed for buddy check qc.   
         if (twodvar_regional .and. dtype(is) == 't' .and.  buddycheck_t) then 
            ikey_yes=0 
            ikey_no=0 
@@ -142,13 +142,13 @@ subroutine obs_para(ndata,mype)
               end if 
            end do 
  
-           ! With organized colors and keys, now create the new MPI communicator 
+           ! With organized colors and keys, now create the new mpi communicator 
            !   which only talks to pe's who have obs on their subdomains.  This is 
-           !   needed for MPI communication within the setup* routines (e.g. a buddy check). 
+           !   needed for mpi communication within the setup* routines (e.g. a buddy check). 
                
            call mpi_comm_split(mpi_comm_world,icolor(mm1),ikey(mm1),obs_sub_comm(is),ierror)   
-           CALL MPI_COMM_SIZE(obs_sub_comm(is), newprocs, ierror) 
-           CALL MPI_COMM_RANK(obs_sub_comm(is), newrank, ierror) 
+           call mpi_comm_size(obs_sub_comm(is), newprocs, ierror) 
+           call mpi_comm_rank(obs_sub_comm(is), newrank, ierror) 
            if (buddydiag_save) write(6,'(A,I3,I10,A,I20,A,I3,A,I3)') 'obs_para: mype/myobs=',& 
                               mype,nobs_sub(mm1,is),'newcomm=',obs_sub_comm(is),'newprocs=', & 
                               newprocs,'newrank=',newrank            
@@ -226,9 +226,9 @@ subroutine disobs(ndata,nobs,mm1,lunout,obsfile,obstypeall)
   integer(i_kind),dimension(mm1):: ibe,ibw,ibn,ibs
   logical,allocatable,dimension(:):: luse_s
   integer(i_kind),allocatable,dimension(:):: ioid_s
-        ! IOID represents Initial-Obs-ID, which is the serial index of given
+        ! ioid represents initial-obs-id, which is the serial index of given
         ! observation record of a given observation source, before they are
-        ! distributed according to the guess-grid partition.  This ID is saved
+        ! distributed according to the guess-grid partition.  This id is saved
         ! for all future needs of obs. redistribution, such that the sorted
         ! sequences can be reproduced.
   real(r_kind),allocatable,dimension(:,:):: obs_data,data1_s
@@ -256,7 +256,7 @@ subroutine disobs(ndata,nobs,mm1,lunout,obsfile,obstypeall)
   nn_obs = nreal + nchanl
 
   allocate(obs_data(nn_obs,ndata))
-!  Read in all observations of a given type along with subdomain flags
+! Read in all observations of a given type along with subdomain flags
   read(lunin) obs_data
   close(lunin)
 
diff --git a/src/gsi/obs_sensitivity.f90 b/src/gsi/obs_sensitivity.f90
index 880ee6384a..3c38f6341f 100644
--- a/src/gsi/obs_sensitivity.f90
+++ b/src/gsi/obs_sensitivity.f90
@@ -12,7 +12,7 @@ module obs_sensitivity
 !   2007-07-19 tremolet - increment sensitivity to observations
 !   2009-08-07 lueken   - updated documentation
 !   2010-04-30 tangborn - add pointer to carbon monoxide
-!   2010-05-27 todling  - remove all user-specific TL-related references
+!   2010-05-27 todling  - remove all user-specific tl-related references
 !   2010-07-16 todling  - add reference to aero and aerol
 !   2010-08-19 lueken   - add only to module use;no machine code, so use .f90
 !   2011-03-29 todling  - add reference to pm2_5
@@ -21,10 +21,10 @@ module obs_sensitivity
 !                         howv ,tcamt, lcbas, cldch
 !   2016-02-20 pagowski - add pm10
 !   2016-05-05 pondeca  - add reference to uwnd10m, vwnd10m
-!   2017-05-12 Y. Wang and X. Wang - add reflectivity (dBZ), POC: xuguang.wang@ou.edu
+!   2017-05-12 Y. Wang and X. Wang - add reflectivity (dbz)
 !   2017-01-16 Apodaca  - add reference to lightning
-!   2017-05-06 todling  - reload ensemble when FSOI calc doing EFSOI-like
-!   2017-05-21 todling  - implement ability to do 2nd EFSOI-like calc
+!   2017-05-06 todling  - reload ensemble when fsoi calc doing efsoi-like
+!   2017-05-21 todling  - implement ability to do 2nd efsoi-like calc
 !
 ! Subroutines Included:
 !   init_fc_sens  - Initialize computations
@@ -44,8 +44,8 @@ module obs_sensitivity
 use gsi_4dvar, only: tau_fcst,efsoi_order,efsoi_afcst,efsoi_ana
 use jfunc, only: jiter, miter, niter, iter
 
-use gsi_obOperTypeManager, only: nobs_type => obOper_count
-use gsi_obOperTypeManager, only: obOper_typeinfo
+use gsi_obopertypemanager, only: nobs_type => oboper_count
+use gsi_obopertypemanager, only: oboper_typeinfo
 use mpimod, only: mype
 use control_vectors, only: control_vector,allocate_cv,read_cv,deallocate_cv, &
     dot_product,assignment(=)
@@ -71,9 +71,9 @@ module obs_sensitivity
        init_obsens, init_fc_sens, save_fc_sens, dot_prod_obs
 public efsoi_o2_update
 
-public:: obsensCounts_realloc
-public:: obsensCounts_set
-public:: obsensCounts_dealloc
+public:: obsenscounts_realloc
+public:: obsenscounts_set
+public:: obsenscounts_dealloc
 
 logical lobsensfc,lobsensjb,lobsensincr, &
         lobsensadj,lobsensmin,llancdone,lsensrecompute
@@ -89,47 +89,47 @@ module obs_sensitivity
 character(len=*),parameter:: myname="obs_sensitivity"
 ! ------------------------------------------------------------------------------
 contains
-!>> object obsensCounts_:
+!>> object obsenscounts_:
 !>> this object was public obsmod::obscounts(:,:), but now private module
 !>> variable in this module.  It is accessed through following module procedures
 !>> []_alloc(), []_set() and []_dealloc().
 
-subroutine obsensCounts_realloc(ntype,nbin)
+subroutine obsenscounts_realloc(ntype,nbin)
 !>> was implemented in setuprhsall()
   implicit none
   integer(i_kind),intent(in):: ntype
   integer(i_kind),intent(in):: nbin
-  character(len=*),parameter:: myname_=myname//"::obsensCounts_realloc"
+  character(len=*),parameter:: myname_=myname//"::obsenscounts_realloc"
   if(allocated(obscounts)) deallocate(obscounts)
   allocate(obscounts(ntype,nbin))
-end subroutine obsensCounts_realloc
+end subroutine obsenscounts_realloc
 
-subroutine obsensCounts_set(iobsglb)
+subroutine obsenscounts_set(iobsglb)
 !>> was implemented in evaljo()
   implicit none
   integer(i_kind),dimension(:,:),intent(in):: iobsglb
-  character(len=*),parameter:: myname_=myname//"::obsensCounts_set"
+  character(len=*),parameter:: myname_=myname//"::obsenscounts_set"
   if(.not.allocated(obscounts)) then
-    call perr(myname_,'not allocated, obscounts')
-    call perr(myname_,'was evaljo() exception 125')
-    call  die(myname_)
+     call perr(myname_,'not allocated, obscounts')
+     call perr(myname_,'was evaljo() exception 125')
+     call  die(myname_)
   endif
   if(any(shape(obscounts)/=shape(iobsglb))) then
-    call perr(myname_,'mismatched, storage size(obscounts,1) =',size(obscounts,1))
-    call perr(myname_,'             argument size(iobsglb,1) =',size(iobsglb,1))
-    call perr(myname_,'            storage size(obscounts,2) =',size(obscounts,2))
-    call perr(myname_,'             argument size(iobsglb,2) =',size(iobsglb,2))
-    call  die(myname_)
+     call perr(myname_,'mismatched, storage size(obscounts,1) =',size(obscounts,1))
+     call perr(myname_,'             argument size(iobsglb,1) =',size(iobsglb,1))
+     call perr(myname_,'            storage size(obscounts,2) =',size(obscounts,2))
+     call perr(myname_,'             argument size(iobsglb,2) =',size(iobsglb,2))
+     call  die(myname_)
   endif
   obscounts(:,:)=iobsglb(:,:)
-end subroutine obsensCounts_set
+end subroutine obsenscounts_set
 
-subroutine obsensCounts_dealloc()
+subroutine obsenscounts_dealloc()
 !>> was a part of obsmod::destroyobs_().
   implicit none
-  character(len=*),parameter:: myname_=myname//"::obsensCounts_dealloc"
+  character(len=*),parameter:: myname_=myname//"::obsenscounts_dealloc"
   if(allocated(obscounts)) deallocate(obscounts)
-end subroutine obsensCounts_dealloc
+end subroutine obsenscounts_dealloc
 
 ! ------------------------------------------------------------------------------
 subroutine init_obsens
@@ -176,9 +176,9 @@ subroutine init_fc_sens
 ! program history log:
 !   2007-06-26  tremolet - initial code
 !   2009-08-07  lueken - added subprogram doc block
-!   2010-05-27  todling - gsi_4dcoupler; remove dependence on GMAO specifics
+!   2010-05-27  todling - gsi_4dcoupler; remove dependence on gmao specifics
 !   2012-05-22  todling - update interface to getpert
-!   2015-12-01  todling - add several obs-types that Pondeca forget to add here
+!   2015-12-01  todling - add several obs-types that pondeca forget to add here
 !
 !   input argument list:
 !
@@ -247,7 +247,7 @@ subroutine init_fc_sens
             write(clfile(10:12),'(I3.3)') miter
             call read_cv(fcsens,clfile)
          else
-!           read and convert output of GCM adjoint
+!           read and convert output of gcm adjoint
             allocate(fname(nsubwin))
             fname='NULL'
             if (tau_fcst>0) then
@@ -297,7 +297,7 @@ subroutine init_fc_sens
       else
 !        read gradient from outer loop jiter+1
          clfile='fgsens.ZZZ'
-         WRITE(clfile(8:10),'(I3.3)') jiter+1
+         write(clfile(8:10),'(I3.3)') jiter+1
          call read_cv(fcsens,clfile)
       endif
    endif
@@ -320,9 +320,9 @@ subroutine efsoi_o2_update(sval)
 ! program history log:
 !   2007-06-26  tremolet - initial code
 !   2009-08-07  lueken - added subprogram doc block
-!   2010-05-27  todling - gsi_4dcoupler; remove dependence on GMAO specifics
+!   2010-05-27  todling - gsi_4dcoupler; remove dependence on gmao specifics
 !   2012-05-22  todling - update interface to getpert
-!   2015-12-01  todling - add several obs-types that Pondeca forget to add here
+!   2015-12-01  todling - add several obs-types that pondeca forget to add here
 !
 !   input argument list:
 !
@@ -360,75 +360,75 @@ subroutine efsoi_o2_update(sval)
 end if
 
 ! Zero out whatever is in fcsens
- fcsens=zero
-
-! read and convert output of GCM adjoint
- allocate(fname(nsubwin))
- fname='A'
- do ii=1,nsubwin
-    call allocate_state(fcgrad(ii))
- end do
- call allocate_preds(zbias)
- zbias=zero
- do ii=1,ntlevs_ens
-    call allocate_state(eval(ii))
- end do
- do ii=1,nsubwin
-    fcgrad(ii)=zero
- end do
- call gsi_4dcoupler_getpert(fcgrad,nsubwin,'adm',fname)
+fcsens=zero
+
+! read and convert output of gcm adjoint
+allocate(fname(nsubwin))
+fname='A'
+do ii=1,nsubwin
+   call allocate_state(fcgrad(ii))
+end do
+call allocate_preds(zbias)
+zbias=zero
+do ii=1,ntlevs_ens
+   call allocate_state(eval(ii))
+end do
+do ii=1,nsubwin
+   fcgrad(ii)=zero
+end do
+call gsi_4dcoupler_getpert(fcgrad,nsubwin,'adm',fname)
 ! Wipe out loaded ensemble members to allow reloading ensemble of analysis
- call destroy_ensemble
+call destroy_ensemble
 ! Read in ens forecasts issues from "analysis" 
- efsoi_afcst=.true.
- call create_ensemble
- call load_ensemble
- do ii=1,ntlevs_ens
-    call allocate_state(eval(ii))
- end do
- eval(1)=fcgrad(1)
- fcgrad(1)=zero
- call ensctl2state_ad(eval,fcgrad(1),fcsens)
- call control2state_ad(fcgrad,zbias,fcsens)
- do ii=1,ntlevs_ens
-    call deallocate_state(eval(ii))
- end do
- do ii=1,nsubwin
-    call deallocate_state(fcgrad(ii))
- end do
- call deallocate_preds(zbias)
- deallocate(fname)
-! Wipe  outensemble from memory
- call destroy_ensemble
+efsoi_afcst=.true.
+call create_ensemble
+call load_ensemble
+do ii=1,ntlevs_ens
+   call allocate_state(eval(ii))
+end do
+eval(1)=fcgrad(1)
+fcgrad(1)=zero
+call ensctl2state_ad(eval,fcgrad(1),fcsens)
+call control2state_ad(fcgrad,zbias,fcsens)
+do ii=1,ntlevs_ens
+   call deallocate_state(eval(ii))
+end do
+do ii=1,nsubwin
+   call deallocate_state(fcgrad(ii))
+end do
+call deallocate_preds(zbias)
+deallocate(fname)
+! Wipe outensemble from memory
+call destroy_ensemble
 
 ! Report magnitude of input vector
- zjx=dot_product(fcsens,fcsens)
- if (mype==0) write(6,888)'order2_fc_sens: Norm fcsens=',sqrt(zjx)
+zjx=dot_product(fcsens,fcsens)
+if (mype==0) write(6,888)'order2_fc_sens: Norm fcsens=',sqrt(zjx)
 888 format(A,3(1X,ES25.18))
 
-! Read in ensemble of analysis (these are EnKF, not GSI, analyses obviously)
- efsoi_ana=.true.
- tau_fcst=-1
- call create_ensemble
- call load_ensemble
+! Read in ensemble of analysis (these are enkf, not gsi, analyses obviously)
+efsoi_ana=.true.
+tau_fcst=-1
+call create_ensemble
+call load_ensemble
 
 ! Allocate local variables
- do ii=1,nsubwin
-    call allocate_state(mval(ii))
- end do
- do ii=1,ntlevs_ens
-    call allocate_state(eval(ii))
- end do
- call allocate_preds(zbias)
+do ii=1,nsubwin
+   call allocate_state(mval(ii))
+end do
+do ii=1,ntlevs_ens
+   call allocate_state(eval(ii))
+end do
+call allocate_preds(zbias)
 
 ! Convert from control variable to state space
- call control2state(fcsens,mval,zbias)
+call control2state(fcsens,mval,zbias)
 
 ! Convert ensemble control variable to state space and update from previous value
- call ensctl2state(fcsens,mval(1),eval)
- do ii=1,ntlevs_ens
-    call self_add(sval(ii),eval(ii))
- end do
+call ensctl2state(fcsens,mval(1),eval)
+do ii=1,ntlevs_ens
+   call self_add(sval(ii),eval(ii))
+end do
 
 ! Release memory
 call deallocate_preds(zbias)
@@ -472,7 +472,7 @@ subroutine save_fc_sens
 
 !  Full stats
    do jj=1,nobs_type
-      write(6,'(A,2X,I3,2X,A)')'Obs types:',jj,obOper_typeinfo(jj)
+      write(6,'(A,2X,I3,2X,A)')'Obs types:',jj,oboper_typeinfo(jj)
    enddo
    write(6,'(A,2X,I4)')'Obs bins:',nobs_bins
    write(6,*)'Obs Count Begin'
@@ -486,8 +486,8 @@ subroutine save_fc_sens
    write(6,*)'Obs Count End'
 
    write(6,*)'Obs Impact Begin'
-   do kk=1,SIZE(sensincr,3)
-      if (SIZE(sensincr,3)==1) then
+   do kk=1,size(sensincr,3)
+      if (size(sensincr,3)==1) then
          write(6,'(A,I4)')'Obs Impact iteration= ',niter(jiter)
       else
          write(6,'(A,I4)')'Obs Impact iteration= ',kk
@@ -498,14 +498,14 @@ subroutine save_fc_sens
    enddo
    write(6,*)'Obs Impact End'
 
-   kk=SIZE(sensincr,3)
+   kk=size(sensincr,3)
 !  Summary by obs type
    do jj=1,nobs_type
       zz=zero
       do ii=1,nobs_bins
          zz=zz+sensincr(ii,jj,kk)
       enddo
-      if (zz/=zero) write(6,'(A,2X,A3,2X,ES12.5)')'Obs Impact type',obOper_typeinfo(jj),zz
+      if (zz/=zero) write(6,'(A,2X,A3,2X,ES12.5)')'Obs Impact type',oboper_typeinfo(jj),zz
    enddo
 
 !  Summary by obs bins
@@ -548,7 +548,7 @@ real(r_kind) function dot_prod_obs()
 !   machine:
 !
 !$$$ end documentation block
-use m_obsdiagNode, only: obs_diag
+use m_obsdiagnode, only: obs_diag
 use m_obsdiags, only: obsdiags
 implicit none
 
diff --git a/src/gsi/omegas_ad.f90 b/src/gsi/omegas_ad.f90
index 7a7274540c..f8103da83a 100644
--- a/src/gsi/omegas_ad.f90
+++ b/src/gsi/omegas_ad.f90
@@ -6,7 +6,7 @@ module omegas_ad_mod
 !
 ! abstract: This module has been added as a wrapper around subroutine omegas_ad
 !            to eliminate type mismatch compile errors when using the debug
-!            compile option on WCOSS.
+!            compile option on wcoss.
 !
 ! program history log:
 !   2012-01-26  parrish
@@ -73,7 +73,7 @@ subroutine omegas_ad_1_1_( im, ix, km, dphi_i_, dlam_i_, u_i_, v_i_, div_i_, &
   integer(i_kind) k
 
   if( im /= 1 .or. ix /= 1 ) then
-     write(6,*)' GSCOND_AD_1_1_, IM,IX=',IM,IX,' -- BOTH MUST BE 1.  PROGRAM FAILS'
+     write(6,*)' GSCOND_AD_1_1_, IM,IX=',im,ix,' -- BOTH MUST BE 1.  PROGRAM FAILS'
      stop
   end if
 
@@ -117,7 +117,7 @@ subroutine omegas_ad_im_ix_( im, ix, km, dphi_i, dlam_i, u_i, v_i, div_i, &
 !     prgmmr:    treadon     org: np23                date: 2003-12-18
 !
 ! abstract:  This subroutine contains the forward and ajoint models for the
-!            GFS vertical velocity calculation
+!            gfs vertical velocity calculation
 !
 ! program history log:
 !   2003-12-18  treadon - initial routine
@@ -127,7 +127,7 @@ subroutine omegas_ad_im_ix_( im, ix, km, dphi_i, dlam_i, u_i, v_i, div_i, &
 !   2013-01-26  parrish - module added as a wrapper around subroutine omegas_ad
 !                            to eliminate type mismatch compile errors when
 !                            using the debug
-!                            compile option on WCOSS.
+!                            compile option on wcoss.
 !
 !   input argument list:
 !     im      - actual number of profiles to be processed
@@ -155,7 +155,7 @@ subroutine omegas_ad_im_ix_( im, ix, km, dphi_i, dlam_i, u_i, v_i, div_i, &
 !
 ! remarks:
 !    The adjoint portion of this routine was generated by the
-!    Tangent linear and Adjoint Model Compiler,  TAMC 5.3.0
+!    tangent linear and adjoint model compiler,  tamc 5.3.0
 !
 ! attributes:
 !   language: f90
@@ -210,7 +210,7 @@ subroutine omegas_ad_im_ix_( im, ix, km, dphi_i, dlam_i, u_i, v_i, div_i, &
   integer(i_kind) k
 
 !----------------------------------------------
-! RESET LOCAL ADJOINT VARIABLES
+! Reset local adjoint variables
 !----------------------------------------------
   do ip2 = 1, ix
      do ip1 = 1, km
@@ -234,10 +234,10 @@ subroutine omegas_ad_im_ix_( im, ix, km, dphi_i, dlam_i, u_i, v_i, div_i, &
   end do
 
 !----------------------------------------------
-! ROUTINE BODY
+! Routine body
 !----------------------------------------------
 !----------------------------------------------
-! FUNCTION AND TAPE COMPUTATIONS
+! Function and tape computations
 !----------------------------------------------
   do i = 1, im
      do k = 1, km+1
@@ -280,7 +280,7 @@ subroutine omegas_ad_im_ix_( im, ix, km, dphi_i, dlam_i, u_i, v_i, div_i, &
   if (.not.adjoint) return
 
 !----------------------------------------------
-! ADJOINT COMPUTATIONS
+! Adjoint computations
 !----------------------------------------------
   do i = 1, im
      do k = 1, km
diff --git a/src/gsi/oneobmod.F90 b/src/gsi/oneobmod.F90
index 571b175972..b6146d5f97 100644
--- a/src/gsi/oneobmod.F90
+++ b/src/gsi/oneobmod.F90
@@ -39,7 +39,7 @@ module oneobmod
 !   def oblon      - observation longitude for one ob exp
 !   def obhourset  - observation delta time from analysis time for 
 !                    one ob exp
-!   def obpres     - observation pressure (hPa) or one ob exp
+!   def obpres     - observation pressure (hpa) or one ob exp
 !   def obdattim   - observation date for one ob exp
 !   def oneob_type - observation type for one ob exp
 !   def oneobtest  - single observation test flag (true=on)
@@ -146,7 +146,7 @@ subroutine oneobmakebufr
 !   2014-08-04  carley - modify for tcamt and howv obs
 !   2014-08-18  carley - added td2m, mxtm, mitm, pmsl, and wsdp10m
 !   2016-01-18  pondeca - added cldch
-!   2016-06-17  Sienkiewicz- virtmp switch - if true write VIRTMP program code to indicate
+!   2016-06-17  Sienkiewicz- virtmp switch - if true write virtmp program code to indicate
 !                            that the observation is virtual temperature
 !
 !   input argument list:
@@ -186,7 +186,7 @@ subroutine oneobmakebufr
     character(80):: cld2seqstr='TOCC HBLCS'      ! total cloud cover and height above surface of base of lowest cloud seen
     character(80):: cldseqstr='VSSO CLAM HOCB'   ! vertical significance, cloud amount and cloud base height
     character(80):: owavestr='HOWV'              ! significant wave height
-    character(80):: maxtmintstr='MXTM MITM'      ! Max T and Min T
+    character(80):: maxtmintstr='MXTM MITM'      ! max t and min t
 !   character(80):: cldceilhstr='CEILING'        ! cldch
     if (oneobmade) return
 
@@ -194,14 +194,14 @@ subroutine oneobmakebufr
        call oneobo3lv
        return
     end if
-! set values from parameter list
+!   set values from parameter list
     xob=oblon
     yob=oblat
     dhr=obhourset
     idate=iadate(1)*1000000+iadate(2)*10000+iadate(3)*100+iadate(4)
     write(6,*)'OneObMake: ', idate
     pob=obpres
-! set default values for this routine
+!   set default values for this routine
     ludx=22
     nobs=1
     nlev=1
@@ -230,7 +230,7 @@ subroutine oneobmakebufr
        subset='ADPSFC'
        typ(1)=87._r_kind
        cat(1,1)=zero
-       cld2seq(1,1)=25._r_kind !TOCC - total cloud amount (%)       
+       cld2seq(1,1)=25._r_kind !tocc - total cloud amount (%)       
     else if (oneob_type=='howv') then
        subset='SFCSHP'
        typ(1)=80._r_kind
@@ -240,24 +240,24 @@ subroutine oneobmakebufr
        subset='ADPSFC'
        typ(1)=87._r_kind
        cat(1,1)=zero
-       obs(12,1)=280.0_r_kind !Dewpoint in Kelvin
+       obs(12,1)=280.0_r_kind !Dewpoint in kelvin
     else if (oneob_type=='mxtm') then
        subset='ADPSFC'
        typ(1)=87._r_kind
        cat(1,1)=zero
-       maxtmint(1,1)=280.0_r_kind !Max T in Kelvin
+       maxtmint(1,1)=280.0_r_kind !Max t in kelvin
     else if (oneob_type=='mitm') then
        subset='ADPSFC'
        typ(1)=87._r_kind
        cat(1,1)=zero
-       maxtmint(2,1)=273.15_r_kind !Min T in Kelvin
+       maxtmint(2,1)=273.15_r_kind !Min t in kelvin
     else if (oneob_type=='pmsl') then
        subset='ADPSFC'
        typ(1)=87._r_kind
        cat(1,1)=zero
-       obs(13,1)=1008.10_r_kind !PMSL in MB
+       obs(13,1)=1008.10_r_kind !pmsl in mb
        qms(7,1)=one
-    else if (oneob_type=='wspd10m') then !10m AGL wind speed - need to store both u and v (already by this point done)
+    else if (oneob_type=='wspd10m') then !10m agl wind speed - need to store both u and v (already by this point done)
        subset='ADPSFC'
        typ(1)=87._r_kind
        cat(1,1)=zero
@@ -270,7 +270,7 @@ subroutine oneobmakebufr
        typ(1)=20._r_kind
        cat(1,1)=one
     endif
-! keep errs small so the single ob passes the QC check
+!   keep errs small so the single ob passes the qc check
     poe=r0_01
     qoe=one_tenth
     toe=one_tenth
@@ -322,12 +322,12 @@ subroutine oneobmakebufr
        call ufbint(lendian_in,err,10,nlev,iret,errstr)
        call ufbint(lendian_in,pcd,10,nlev,iret,pcdstr)
        if (oneob_type=='tcamt') then
-         call ufbint(lendian_in,cldseq,3,10,iret,cldseqstr)
-         call ufbint(lendian_in,cld2seq,2,1,iret,cld2seqstr)
+          call ufbint(lendian_in,cldseq,3,10,iret,cldseqstr)
+          call ufbint(lendian_in,cld2seq,2,1,iret,cld2seqstr)
        else if (oneob_type=='howv') then
-         call ufbint(lendian_in,owave,1,nlev,iret,owavestr)
+          call ufbint(lendian_in,owave,1,nlev,iret,owavestr)
        else if ( oneob_type=='mxtm' .or. oneob_type=='mitm') then
-         call ufbint(lendian_in,maxtmint,2,nlev,iret,maxtmintstr)
+          call ufbint(lendian_in,maxtmint,2,nlev,iret,maxtmintstr)
        end if                              
        call writsb(lendian_in)
        hdr(1)=transfer(sid(n),hdr(1))
@@ -440,7 +440,7 @@ subroutine oneobmakerwsupob
     h=ha-epsh
     thishgt=this_stahgt+h
 
-!     Get corrected tilt angle
+!   Get corrected tilt angle
     celev=celev0
     selev=selev0
     if(thisrange>=one) then
@@ -455,7 +455,7 @@ subroutine oneobmakerwsupob
     deldistmin=min(gamma-thisrange,deldistmin)
 
 
-!     Get earth lat lon of superob
+!   Get earth lat lon of superob
     thisazimuthr=thisazimuth*deg2rad
     rlonloc=rad_per_meter*gamma*cos(thisazimuthr)
     rlatloc=rad_per_meter*gamma*sin(thisazimuthr)
@@ -464,7 +464,7 @@ subroutine oneobmakerwsupob
     thislon=rlonglob*rad2deg
 
 
-!     Get corrected azimuth
+!   Get corrected azimuth
     clat1=cos(rlatglob)
     caz0=cos(thisazimuthr)
     saz0=sin(thisazimuthr)
@@ -496,7 +496,6 @@ subroutine oneobmakerwsupob
     write(6,*) 'thislon     = ',thislon
     write(6,*) 'thishgt     = ',thishgt
     write(6,*) 'maginnov    = ',maginnov
-    !write(6,*) 'thisvr      = ',thisvr ! This quantity is not used for SOT.
     write(6,*) 'corrected_az= ',corrected_azimuth
     write(6,*) 'corrected_tl= ',corrected_tilt
     write(6,*) 'thiserr     = ',thiserr
@@ -520,7 +519,7 @@ subroutine oneobo3lv
 ! abstract: create ozone level text file for single ob experiment
 !
 ! program history log:
-!   2007-09-11  Sienkiewicz - extend to create MLS text file on request
+!   2007-09-11  Sienkiewicz - extend to create mls text file on request
 !
 !   input argument list:
 !
@@ -547,8 +546,6 @@ subroutine oneobo3lv
     isnd = 1
     ppmv = one                ! dummy value 
 
-!    obdattim=2000010100
-
     rlnc = zero
     rlnc(2) = obhourset
     ildat(1) = obdattim  / 1000000            ! year
@@ -575,27 +572,27 @@ subroutine oneobo3lv
 
   end subroutine oneobo3lv
 
-      SUBROUTINE invtllv(ALM,APH,TLMO,CTPH0,STPH0,TLM,TPH)
-!     borrowed from read_radar
-      use kinds, only:  r_kind
-      implicit none
-
-      real(r_kind),intent(in   ) :: alm,aph,tlmo,ctph0,stph0
-      real(r_kind),intent(  out) :: tlm,tph
-
-      real(r_kind):: relm,srlm,crlm,sph,cph,cc,anum,denom
-
-      RELM=ALM
-      SRLM=SIN(RELM)
-      CRLM=COS(RELM)
-      SPH=SIN(APH)
-      CPH=COS(APH)
-      CC=CPH*CRLM
-      ANUM=CPH*SRLM
-      DENOM=CTPH0*CC-STPH0*SPH
-      TLM=tlmo+ATAN2(ANUM,DENOM)
-      TPH=ASIN(CTPH0*SPH+STPH0*CC)
-      END SUBROUTINE invtllv
+  subroutine invtllv(alm,aph,tlmo,ctph0,stph0,tlm,tph)
+!   borrowed from read_radar
+    use kinds, only:  r_kind
+    implicit none
+
+    real(r_kind),intent(in   ) :: alm,aph,tlmo,ctph0,stph0
+    real(r_kind),intent(  out) :: tlm,tph
+
+    real(r_kind):: relm,srlm,crlm,sph,cph,cc,anum,denom
+
+    relm=alm
+    srlm=sin(relm)
+    crlm=cos(relm)
+    sph=sin(aph)
+    cph=cos(aph)
+    cc=cph*crlm
+    anum=cph*srlm
+    denom=ctph0*cc-stph0*sph
+    tlm=tlmo+atan2(anum,denom)
+    tph=asin(ctph0*sph+stph0*cc)
+  end subroutine invtllv
 
 
 end module oneobmod
diff --git a/src/gsi/ozinfo.f90 b/src/gsi/ozinfo.f90
index b42bc29e5c..052d5fbd3f 100644
--- a/src/gsi/ozinfo.f90
+++ b/src/gsi/ozinfo.f90
@@ -16,7 +16,7 @@ module ozinfo
 !   2005-09-28  derber  - Modify for new ozinfo file, add var qc parameters
 !   2006-02-03  derber  - modify for new obs control and obs count
 !   2007-06-29  Zhou    - change total number of ozone enteries (jpch_oz) from
-!                         53 (version 6 SBUV/2) to 67 (version 8 SBUV/2)
+!                         53 (version 6 sbuv/2) to 67 (version 8 sbuv/2)
 !
 ! Subroutines Included:
 !   sub init_oz       - set ozone related variables to defaults
@@ -28,10 +28,10 @@ module ozinfo
 !   def diag_ozone     - logical to turn off or on the diagnostic ozone file (true=on)
 !   def jpch_oz        - number of (levels+1) * number of satellites
 !   def mype_oz        - task id for writing out radiance diagnostics
-!   def pob_oz         - pressure level of observation (hPa)
+!   def pob_oz         - pressure level of observation (hpa)
 !   def gross_oz       - gross error limit
 !   def error_oz       - observation error
-!   def nusis_oz       - sensor/intrument/satellite id (14=NOAA-14, 15=NOAA-15, 16=NOAA-16, etc)
+!   def nusis_oz       - sensor/intrument/satellite id (14=noaa-14, 15=noaa-15, 16=noaa-16, etc)
 !   def nulev          - integer level of ozone observation
 !   def iuse_oz        - integer flag to control usage of ozone data (-1=don't use, 1=use)
 !
@@ -96,7 +96,7 @@ subroutine init_oz
     jpch_oz = 0                        ! number of enteries read from ozinfo
     diag_ozone = .true.                ! default is to generate ozone diagnostic file
     mype_oz     = max(0,npe-6)         ! mpi task to write ozone summary report
-    ihave_oz=(getindex(svars3d,'oz')>0)! .t. when OZ present in state-vector
+    ihave_oz=(getindex(svars3d,'oz')>0)! .t. when oz present in state-vector
 
   end subroutine init_oz
   
@@ -155,8 +155,8 @@ subroutine ozinfo_read
     endif
     jpch_oz = j
     if(jpch_oz == 0)then
-      close(lunin)
-      return
+       close(lunin)
+       return
     end if
 
 !   Allocate arrays to hold ozone information
diff --git a/src/gsi/patch2grid_mod.f90 b/src/gsi/patch2grid_mod.f90
index cccff75772..b6faa0d7c3 100644
--- a/src/gsi/patch2grid_mod.f90
+++ b/src/gsi/patch2grid_mod.f90
@@ -29,7 +29,7 @@ module patch2grid_mod
 !$$$ end documentation block
   use kinds, only: r_kind,i_kind
   use constants,only: zero,one
-  use anberror, only: pf2aP1,pf2aP2,pf2aP3,&
+  use anberror, only: pf2ap1,pf2ap2,pf2ap3,&
                       nx,ny,mr,nr,nf
   use gridmod,  only: nlat,nlon
 
@@ -54,9 +54,9 @@ module patch2grid_mod
 
   real(r_kind),allocatable,dimension(:):: bl,bl2
 
-  contains
+contains
 
-subroutine setup_patch2grid
+  subroutine setup_patch2grid
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    setup_patch2grid
@@ -76,37 +76,37 @@ subroutine setup_patch2grid
 !   machine:
 !
 !$$$ end documentation block
-  implicit none
-
-  integer(i_kind):: ier
-
-  ndx=(nx-nlon)/2
-  ndy=(nlat-ny)/2
-  ndx2=2*ndx
-  nmix=nr+1+(ny-nlat)/2
-  nrmxb=ndy-1
-  nmixp=nmix+1
-  nymx=ny-nmix
-  nlatxb=nlat-nrmxb
-
-  allocate(p1all(ny,nx))
-  allocate(p2all(nlon+1,mr:nr))
-  allocate(p3all(nlon+1,mr:nr))
-  allocate(p2pol(nf*2+1,nf*2+1))
-  allocate(p3pol(nf*2+1,nf*2+1),stat=ier)
-
-  if( ier /= 0 ) then
-     write(6,*) 'setup_patch2grid: could not allocate memories'
-     call stop2(99)
-  end if
+    implicit none
+
+    integer(i_kind):: ier
+
+    ndx=(nx-nlon)/2
+    ndy=(nlat-ny)/2
+    ndx2=2*ndx
+    nmix=nr+1+(ny-nlat)/2
+    nrmxb=ndy-1
+    nmixp=nmix+1
+    nymx=ny-nmix
+    nlatxb=nlat-nrmxb
+
+    allocate(p1all(ny,nx))
+    allocate(p2all(nlon+1,mr:nr))
+    allocate(p3all(nlon+1,mr:nr))
+    allocate(p2pol(nf*2+1,nf*2+1))
+    allocate(p3pol(nf*2+1,nf*2+1),stat=ier)
+
+    if( ier /= 0 ) then
+       write(6,*) 'setup_patch2grid: could not allocate memories'
+       call stop2(99)
+    end if
 
 !   initialize blend array (bl,bl2)
-  allocate( bl(nx-nlon), bl2(nr+1+(ny-nlat)/2) )
-  call setup_blend
+    allocate( bl(nx-nlon), bl2(nr+1+(ny-nlat)/2) )
+    call setup_blend
 
-end subroutine
+  end subroutine
 
-subroutine setup_blend
+  subroutine setup_blend
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    setup_blend
@@ -127,61 +127,61 @@ subroutine setup_blend
 !
 !$$$ end documentation block
 
-  use blendmod, only: blend
-  implicit none
-
-  integer(i_kind),dimension(0:40):: iblend
-  integer(i_kind):: mm,nolp,nbuf,nmixbl,ndxbl
-  real(r_kind):: dxx,x,y
-
-  integer(i_kind):: i,j,k
-
-! Setup blending
-  mm=4
-  call blend(mm,iblend)
-
-  nolp=nr+1+(ny-nlat)/2
-  nbuf=0
-  nmixbl=nolp-nbuf*2
-  dxx=one/(nmixbl+1)
-  bl2=zero
-  k=0
-  do i=1,nmixbl
-     k=k+1
-     x=i*dxx
-     y=iblend(mm)
-     do j=mm-1,0,-1
-        y=x*y+iblend(j)
-     end do
-     y=y*x**(mm+1)
-     bl2(k)=one-y
-  end do
-  do k=1,nmixbl
-     bl2(k)=sqrt(bl2(k))
-  end do
-
-  nmixbl=(nx-nlon)
-  dxx=one/(nmixbl+1)
-  ndxbl=(nx-nlon)
-  bl=zero
-  k=ndxbl-nmixbl
-  do i=1,nmixbl
-     k=k+1
-     x=i*dxx
-     y=iblend(mm)
-     do j=mm-1,0,-1
-        y=x*y+iblend(j)
-     end do
-     y=y*x**(mm+1)
-     bl(k)=one-y
-  enddo
-  do k=1,nmixbl
-     bl(k)=sqrt(bl(k))
-  end do
-
-end subroutine setup_blend
-
-subroutine destroy_patch2grid
+    use blendmod, only: blend
+    implicit none
+
+    integer(i_kind),dimension(0:40):: iblend
+    integer(i_kind):: mm,nolp,nbuf,nmixbl,ndxbl
+    real(r_kind):: dxx,x,y
+
+    integer(i_kind):: i,j,k
+
+!   Setup blending
+    mm=4
+    call blend(mm,iblend)
+
+    nolp=nr+1+(ny-nlat)/2
+    nbuf=0
+    nmixbl=nolp-nbuf*2
+    dxx=one/(nmixbl+1)
+    bl2=zero
+    k=0
+    do i=1,nmixbl
+       k=k+1
+       x=i*dxx
+       y=iblend(mm)
+       do j=mm-1,0,-1
+          y=x*y+iblend(j)
+       end do
+       y=y*x**(mm+1)
+       bl2(k)=one-y
+    end do
+    do k=1,nmixbl
+       bl2(k)=sqrt(bl2(k))
+    end do
+
+    nmixbl=(nx-nlon)
+    dxx=one/(nmixbl+1)
+    ndxbl=(nx-nlon)
+    bl=zero
+    k=ndxbl-nmixbl
+    do i=1,nmixbl
+       k=k+1
+       x=i*dxx
+       y=iblend(mm)
+       do j=mm-1,0,-1
+          y=x*y+iblend(j)
+       end do
+       y=y*x**(mm+1)
+       bl(k)=one-y
+    enddo
+    do k=1,nmixbl
+       bl(k)=sqrt(bl(k))
+    end do
+
+  end subroutine setup_blend
+
+  subroutine destroy_patch2grid
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    destroy_patch2grid
@@ -201,18 +201,18 @@ subroutine destroy_patch2grid
 !   machine:
 !
 !$$$ end documentation block
-  implicit none
+    implicit none
 
-  deallocate(p1all)
-  deallocate(p2all,p3all)
-  deallocate(p2pol,p3pol)
-  deallocate(bl,bl2)
+    deallocate(p1all)
+    deallocate(p2all,p3all)
+    deallocate(p2pol,p3pol)
+    deallocate(bl,bl2)
 
-end subroutine
+  end subroutine
 
 !=======================================================================
 !=======================================================================
-subroutine grid2patch(grid_wrk,hflt_all,hflt_p2,hflt_p3)
+  subroutine grid2patch(grid_wrk,hflt_all,hflt_p2,hflt_p3)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    grid2patch
@@ -237,60 +237,60 @@ subroutine grid2patch(grid_wrk,hflt_all,hflt_p2,hflt_p3)
 !
 !$$$ end documentation block
 
-  use smooth_polcarf, only: smooth_caspol
-  use fgrid2agrid_mod, only: agrid2fgrid
-
-  implicit none
-
-  real(r_kind),dimension(nlat,nlon),intent(in   ) :: grid_wrk
-
-  real(r_kind)                     ,intent(  out) :: hflt_all(pf2aP1%nlatf,pf2aP1%nlonf)
-  real(r_kind)                     ,intent(  out) :: hflt_p2 (pf2aP2%nlatf ,pf2aP2%nlonf )
-  real(r_kind)                     ,intent(  out) :: hflt_p3 (pf2aP3%nlatf ,pf2aP3%nlonf )
-
-  integer(i_kind):: i,i1,i2,j,j1
-
-! Extract central patch (band) from full grid_wrk (work --> p1)
-! Blending zones
-  p1all=zero
-  do i=1,ndx
-     i1=i-ndx+nlon
-     i2=nx-ndx+i
-     do j=1,ny
-        j1=j+ndy
-        p1all(j,i) =grid_wrk(j1,i1)      ! left (west) blending zone
-        p1all(j,i2)=grid_wrk(j1,i)       ! right (east) blending zone
-     end do
-  end do
-
-! Middle zone (no blending)
-  do i=ndx+1,nx-ndx
-     i1=i-ndx
-     do j=1,ny
-        p1all(j,i)=grid_wrk(j+ndy,i1)
-     end do
-  end do
-
-! Handle polar patches
-  do i=1,nlon
-     do j=mr,nrmxb+nmix
-        p2all(i,j)=grid_wrk(nlat-j,i)
-        p3all(i,j)=grid_wrk(j+1,i)
-     end do
-  end do
-
-  call agrid2fgrid(pf2aP1,p1all,hflt_all) !analysis to filter grid_wrk
-
-  call smooth_caspol(p2pol,p2all)
-  call agrid2fgrid(pf2aP2 ,p2pol,hflt_p2)  !analysis to filter grid_wrk
-
-  call smooth_caspol(p3pol,p3all)
-  call agrid2fgrid(pf2aP3 ,p3pol,hflt_p3)  !analysis to filter grid_wrk
-
-end subroutine grid2patch
+    use smooth_polcarf, only: smooth_caspol
+    use fgrid2agrid_mod, only: agrid2fgrid
+
+    implicit none
+
+    real(r_kind),dimension(nlat,nlon),intent(in   ) :: grid_wrk
+
+    real(r_kind)                     ,intent(  out) :: hflt_all(pf2ap1%nlatf,pf2ap1%nlonf)
+    real(r_kind)                     ,intent(  out) :: hflt_p2 (pf2ap2%nlatf,pf2ap2%nlonf)
+    real(r_kind)                     ,intent(  out) :: hflt_p3 (pf2ap3%nlatf,pf2ap3%nlonf)
+
+    integer(i_kind):: i,i1,i2,j,j1
+
+!   Extract central patch (band) from full grid_wrk (work --> p1)
+!   Blending zones
+    p1all=zero
+    do i=1,ndx
+       i1=i-ndx+nlon
+       i2=nx-ndx+i
+       do j=1,ny
+          j1=j+ndy
+          p1all(j,i) =grid_wrk(j1,i1)      ! left (west) blending zone
+          p1all(j,i2)=grid_wrk(j1,i)       ! right (east) blending zone
+       end do
+    end do
+
+!   Middle zone (no blending)
+    do i=ndx+1,nx-ndx
+       i1=i-ndx
+       do j=1,ny
+          p1all(j,i)=grid_wrk(j+ndy,i1)
+       end do
+    end do
+
+!   Handle polar patches
+    do i=1,nlon
+       do j=mr,nrmxb+nmix
+          p2all(i,j)=grid_wrk(nlat-j,i)
+          p3all(i,j)=grid_wrk(j+1,i)
+       end do
+    end do
+
+    call agrid2fgrid(pf2ap1,p1all,hflt_all) !analysis to filter grid_wrk
+
+    call smooth_caspol(p2pol,p2all)
+    call agrid2fgrid(pf2ap2 ,p2pol,hflt_p2)  !analysis to filter grid_wrk
+
+    call smooth_caspol(p3pol,p3all)
+    call agrid2fgrid(pf2ap3 ,p3pol,hflt_p3)  !analysis to filter grid_wrk
+
+  end subroutine grid2patch
 !=======================================================================
 !=======================================================================
-subroutine tpatch2grid(grid_wrk,hflt_all,hflt_p2,hflt_p3)
+  subroutine tpatch2grid(grid_wrk,hflt_all,hflt_p2,hflt_p3)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    tpatch2grid
@@ -315,93 +315,93 @@ subroutine tpatch2grid(grid_wrk,hflt_all,hflt_p2,hflt_p3)
 !
 !$$$ end documentation block
 
-  use smooth_polcarf, only: smooth_polcasa
-  use fgrid2agrid_mod, only: tfgrid2agrid
-
-  implicit none
-
-  real(r_kind),dimension(nlat,nlon),intent(in   ) :: grid_wrk
-
-  real(r_kind)                     ,intent(  out) :: hflt_all(pf2aP1%nlatf,pf2aP1%nlonf)
-  real(r_kind)                     ,intent(  out) :: hflt_p2 (pf2aP2%nlatf ,pf2aP2%nlonf )
-  real(r_kind)                     ,intent(  out) :: hflt_p3 (pf2aP3%nlatf ,pf2aP3%nlonf )
-
-  integer(i_kind):: i,i1,i2,j,j1
-
-! Extract central patch (band) from full grid_wrk (work --> p1)
-! Blending zones
-  p1all=zero
-
-  do i=1,ndx
-     i1=i-ndx+nlon
-     i2=nx-ndx+i
-     do j=1,ny
-        j1=j+ndy
-        p1all(j,i) =grid_wrk(j1,i1)      ! left (west) blending zone
-        p1all(j,i2)=grid_wrk(j1,i)       ! right (east) blending zone
-     end do
-  end do
-
-! Middle zone (no blending)
-  do i=ndx+1,nx-ndx
-     i1=i-ndx
-     do j=1,ny
-        p1all(j,i)=grid_wrk(j+ndy,i1)
-     end do
-  end do
-
-! Apply blending coefficients to central patch
-  do i=1,ndx2
-     i1=ndx2+1-i
-     i2=nx-ndx2+i
-     do j=1,ny
-        p1all(j,i) =p1all(j,i) *bl(i1)  ! left (west) blending zone
-        p1all(j,i2)=p1all(j,i2)*bl(i)   ! right (east) blending zone
-     end do
-  end do
-
-! bl2 of p1
-  do i=1,nx
-     do j=1,nmix
-        p1all(j,i)=p1all(j,i)*bl2(nmixp-j)
-     end do
-     do j=nymx+1,ny
-       p1all(j,i)=p1all(j,i)*bl2(j-nymx)
-     end do
-  end do
-
-! Handle polar patches
-  p2all=zero
-  p3all=zero
-
-  do j=mr,nrmxb+nmix
-     do i=1,nlon
-        p2all(i,j)=grid_wrk(nlat-j,i)
-        p3all(i,j)=grid_wrk(j+1,i)
-     end do
-  end do
-
-! Apply blending coefficients
-  do j=nrmxb+1,nrmxb+nmix
-     j1=j-nrmxb
-     do i=1,nlon
-        p2all(i,j)=p2all(i,j)*bl2(j1)
-        p3all(i,j)=p3all(i,j)*bl2(j1)
-     end do
-  end do
-
-  call tfgrid2agrid(pf2aP1,p1all,hflt_all) !analysis to filter grid_wrk
-
-  call smooth_polcasa(p2pol,p2all)
-  call tfgrid2agrid(pf2aP2,p2pol,hflt_p2) !analysis to filter grid_wrk
-
-  call smooth_polcasa(p3pol,p3all)
-  call tfgrid2agrid(pf2aP3,p3pol,hflt_p3) !analysis to filter grid_wrk
-
-end subroutine tpatch2grid
+    use smooth_polcarf, only: smooth_polcasa
+    use fgrid2agrid_mod, only: tfgrid2agrid
+
+    implicit none
+
+    real(r_kind),dimension(nlat,nlon),intent(in   ) :: grid_wrk
+
+    real(r_kind)                     ,intent(  out) :: hflt_all(pf2ap1%nlatf,pf2ap1%nlonf)
+    real(r_kind)                     ,intent(  out) :: hflt_p2 (pf2ap2%nlatf,pf2ap2%nlonf)
+    real(r_kind)                     ,intent(  out) :: hflt_p3 (pf2ap3%nlatf,pf2ap3%nlonf)
+
+    integer(i_kind):: i,i1,i2,j,j1
+
+!   Extract central patch (band) from full grid_wrk (work --> p1)
+!   Blending zones
+    p1all=zero
+
+    do i=1,ndx
+       i1=i-ndx+nlon
+       i2=nx-ndx+i
+       do j=1,ny
+          j1=j+ndy
+          p1all(j,i) =grid_wrk(j1,i1)      ! left (west) blending zone
+          p1all(j,i2)=grid_wrk(j1,i)       ! right (east) blending zone
+       end do
+    end do
+
+!   Middle zone (no blending)
+    do i=ndx+1,nx-ndx
+       i1=i-ndx
+       do j=1,ny
+          p1all(j,i)=grid_wrk(j+ndy,i1)
+       end do
+    end do
+
+!   Apply blending coefficients to central patch
+    do i=1,ndx2
+       i1=ndx2+1-i
+       i2=nx-ndx2+i
+       do j=1,ny
+          p1all(j,i) =p1all(j,i) *bl(i1)  ! left (west) blending zone
+          p1all(j,i2)=p1all(j,i2)*bl(i)   ! right (east) blending zone
+       end do
+    end do
+
+!   bl2 of p1
+    do i=1,nx
+       do j=1,nmix
+          p1all(j,i)=p1all(j,i)*bl2(nmixp-j)
+       end do
+       do j=nymx+1,ny
+         p1all(j,i)=p1all(j,i)*bl2(j-nymx)
+       end do
+    end do
+
+!   Handle polar patches
+    p2all=zero
+    p3all=zero
+
+    do j=mr,nrmxb+nmix
+       do i=1,nlon
+          p2all(i,j)=grid_wrk(nlat-j,i)
+          p3all(i,j)=grid_wrk(j+1,i)
+       end do
+    end do
+
+!   Apply blending coefficients
+    do j=nrmxb+1,nrmxb+nmix
+       j1=j-nrmxb
+       do i=1,nlon
+          p2all(i,j)=p2all(i,j)*bl2(j1)
+          p3all(i,j)=p3all(i,j)*bl2(j1)
+       end do
+    end do
+
+    call tfgrid2agrid(pf2ap1,p1all,hflt_all) !analysis to filter grid_wrk
+
+    call smooth_polcasa(p2pol,p2all)
+    call tfgrid2agrid(pf2ap2,p2pol,hflt_p2) !analysis to filter grid_wrk
+
+    call smooth_polcasa(p3pol,p3all)
+    call tfgrid2agrid(pf2ap3,p3pol,hflt_p3) !analysis to filter grid_wrk
+
+  end subroutine tpatch2grid
 !=======================================================================
 !=======================================================================
-subroutine patch2grid(grid_wrk,hflt_all,hflt_p2,hflt_p3)
+  subroutine patch2grid(grid_wrk,hflt_all,hflt_p2,hflt_p3)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    patch2grid
@@ -426,103 +426,103 @@ subroutine patch2grid(grid_wrk,hflt_all,hflt_p2,hflt_p3)
 !
 !$$$ end documentation block
 
-  use smooth_polcarf, only: smooth_polcas
-  use fgrid2agrid_mod, only: fgrid2agrid
-
-  implicit none
-
-  real(r_kind),dimension(nlat,nlon),intent(  out) :: grid_wrk
-
-  real(r_kind)                     ,intent(in   ) :: hflt_all(pf2aP1%nlatf,pf2aP1%nlonf)
-  real(r_kind)                     ,intent(in   ) :: hflt_p2 (pf2aP2%nlatf ,pf2aP2%nlonf )
-  real(r_kind)                     ,intent(in   ) :: hflt_p3 (pf2aP3%nlatf ,pf2aP3%nlonf )
-
-  integer(i_kind):: i,i1,i2,j,j1
-
-  call fgrid2agrid(pf2aP1,hflt_all,p1all) !analysis to filter grid_wrk
-
-  call fgrid2agrid(pf2aP2,hflt_p2,p2pol) !analysis to filter grid_wrk
-  call smooth_polcas(p2pol,p2all)
-
-  call fgrid2agrid(pf2aP3,hflt_p3,p3pol) !analysis to filter grid_wrk
-  call smooth_polcas(p3pol,p3all)
-
-! zero output array
-  grid_wrk=zero
-
-! Equatorial patch
-! Adjoint of central patch blending on left/right sides of patch
-  do i=1,ndx2
-     i1=ndx2+1-i
-     i2=nx-ndx2+i
-     do j=1,ny
-        p1all(j,i) =p1all(j,i) *bl(i1)   ! left (west) blending zone
-        p1all(j,i2)=p1all(j,i2)*bl(i)    ! right (east) blending zone
-     end do
-  end do
-
-! bl2 of p1
-  do i=1,nx
-     do j=1,nmix
-        p1all(j,i)=p1all(j,i)*bl2(nmixp-j)
-     end do
-     do j=nymx+1,ny
-        p1all(j,i)=p1all(j,i)*bl2(j-nymx)
-     end do
-  end do
-
-! Adjoint of transfer between central band and full grid (p1 --> work)
-  do i=1,ndx
-     i1=i-ndx+nlon
-     i2=nx-ndx+i
-     do j=1,ny
-        j1=j+ndy
-        grid_wrk(j1,i1)=grid_wrk(j1,i1)+p1all(j,i)  ! left (west) blending zone
-        grid_wrk(j1,i )=grid_wrk(j1,i )+p1all(j,i2) ! right (east) blending zone
-     end do
-  end do
-
-! Middle zone (no blending)
-  do i=ndx+1,nx-ndx
-     i1=i-ndx
-     do j=1,ny
-        j1=j+ndy
-        grid_wrk(j1,i1)=grid_wrk(j1,i1)+p1all(j,i)
-     end do
-  end do
-
-! Adjoint of North pole patch(p2) -- blending and transfer to grid
-! Adjoint of South pole patch(p3) -- blending and transfer to grid
-
-  do j=nlatxb-nmix,nlatxb-1
-! Adjoint of blending
-     do i=1,nlon
-        p2all(i,nlat-j)=p2all(i,nlat-j)*bl2(nlatxb-j)
-     end do
-  end do
-
-  do j=nrmxb+1,nrmxb+nmix
-!    Adjoint of blending
-     do i=1,nlon
-        p3all(i,j)=p3all(i,j)*bl2(j-nrmxb)
-     end do
-  end do
-
-  do i=1,nlon
-!    Adjoint of transfer
-     do j=mr,nrmxb+nmix
-        grid_wrk(j+1,i)=grid_wrk(j+1,i)+p3all(i,j)
-     end do
-     do j=nlatxb-nmix,nlat-mr
-        grid_wrk(j  ,i)=grid_wrk(j  ,i)+p2all(i,nlat-j)
-     end do
-  end do
-
-end subroutine patch2grid
+    use smooth_polcarf, only: smooth_polcas
+    use fgrid2agrid_mod, only: fgrid2agrid
+
+    implicit none
+
+    real(r_kind),dimension(nlat,nlon),intent(  out) :: grid_wrk
+
+    real(r_kind)                     ,intent(in   ) :: hflt_all(pf2ap1%nlatf,pf2ap1%nlonf)
+    real(r_kind)                     ,intent(in   ) :: hflt_p2 (pf2ap2%nlatf,pf2ap2%nlonf)
+    real(r_kind)                     ,intent(in   ) :: hflt_p3 (pf2ap3%nlatf,pf2ap3%nlonf)
+
+    integer(i_kind):: i,i1,i2,j,j1
+
+    call fgrid2agrid(pf2ap1,hflt_all,p1all) !analysis to filter grid_wrk
+
+    call fgrid2agrid(pf2ap2,hflt_p2,p2pol) !analysis to filter grid_wrk
+    call smooth_polcas(p2pol,p2all)
+
+    call fgrid2agrid(pf2ap3,hflt_p3,p3pol) !analysis to filter grid_wrk
+    call smooth_polcas(p3pol,p3all)
+
+!   zero output array
+    grid_wrk=zero
+
+!   Equatorial patch
+!   Adjoint of central patch blending on left/right sides of patch
+    do i=1,ndx2
+       i1=ndx2+1-i
+       i2=nx-ndx2+i
+       do j=1,ny
+          p1all(j,i) =p1all(j,i) *bl(i1)   ! left (west) blending zone
+          p1all(j,i2)=p1all(j,i2)*bl(i)    ! right (east) blending zone
+       end do
+    end do
+
+!   bl2 of p1
+    do i=1,nx
+       do j=1,nmix
+          p1all(j,i)=p1all(j,i)*bl2(nmixp-j)
+       end do
+       do j=nymx+1,ny
+          p1all(j,i)=p1all(j,i)*bl2(j-nymx)
+       end do
+    end do
+
+!   Adjoint of transfer between central band and full grid (p1 --> work)
+    do i=1,ndx
+       i1=i-ndx+nlon
+       i2=nx-ndx+i
+       do j=1,ny
+          j1=j+ndy
+          grid_wrk(j1,i1)=grid_wrk(j1,i1)+p1all(j,i)  ! left (west) blending zone
+          grid_wrk(j1,i )=grid_wrk(j1,i )+p1all(j,i2) ! right (east) blending zone
+       end do
+    end do
+
+!   Middle zone (no blending)
+    do i=ndx+1,nx-ndx
+       i1=i-ndx
+       do j=1,ny
+          j1=j+ndy
+          grid_wrk(j1,i1)=grid_wrk(j1,i1)+p1all(j,i)
+       end do
+    end do
+
+!   Adjoint of north pole patch(p2) -- blending and transfer to grid
+!   Adjoint of south pole patch(p3) -- blending and transfer to grid
+
+    do j=nlatxb-nmix,nlatxb-1
+!      Adjoint of blending
+       do i=1,nlon
+          p2all(i,nlat-j)=p2all(i,nlat-j)*bl2(nlatxb-j)
+       end do
+    end do
+
+    do j=nrmxb+1,nrmxb+nmix
+!      Adjoint of blending
+       do i=1,nlon
+          p3all(i,j)=p3all(i,j)*bl2(j-nrmxb)
+       end do
+    end do
+
+    do i=1,nlon
+!      Adjoint of transfer
+       do j=mr,nrmxb+nmix
+          grid_wrk(j+1,i)=grid_wrk(j+1,i)+p3all(i,j)
+       end do
+       do j=nlatxb-nmix,nlat-mr
+          grid_wrk(j  ,i)=grid_wrk(j  ,i)+p2all(i,nlat-j)
+       end do
+    end do
+
+  end subroutine patch2grid
 
 !=======================================================================
 !=======================================================================
-subroutine vpatch2grid(grid_wrk,hflt_all,hflt_p2,hflt_p3)
+  subroutine vpatch2grid(grid_wrk,hflt_all,hflt_p2,hflt_p3)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    vpatch2grid
@@ -547,130 +547,130 @@ subroutine vpatch2grid(grid_wrk,hflt_all,hflt_p2,hflt_p3)
 !
 !$$$ end documentation block
 
-  use smooth_polcarf, only: smooth_polcasv
-  use fgrid2agrid_mod, only: fgrid2agrid
-
-  implicit none
-
-  real(r_kind),dimension(nlat,nlon),intent(  out) :: grid_wrk
-
-  real(r_kind)                     ,intent(in   ) :: hflt_all(pf2aP1%nlatf,pf2aP1%nlonf)
-  real(r_kind)                     ,intent(in   ) :: hflt_p2 (pf2aP2%nlatf ,pf2aP2%nlonf )
-  real(r_kind)                     ,intent(in   ) :: hflt_p3 (pf2aP3%nlatf ,pf2aP3%nlonf )
-
-  real(r_kind),allocatable,dimension(:,:) :: wgt_wrk
-  real(r_kind),allocatable,dimension(:,:) :: p1wgt
-  real(r_kind),allocatable,dimension(:,:) :: p2wgt,p3wgt
-
-  integer(i_kind):: i,i1,i2,j,j1
-
-  allocate(wgt_wrk(nlat,nlon))
-
-  call fgrid2agrid(pf2aP1,hflt_all,p1all) !analysis to filter grid_wrk
-
-  call fgrid2agrid(pf2aP2,hflt_p2,p2pol) !analysis to filter grid_wrk
-  call smooth_polcasv(p2pol,p2all)
-
-  call fgrid2agrid(pf2aP3,hflt_p3,p3pol) !analysis to filter grid_wrk
-  call smooth_polcasv(p3pol,p3all)
-
-  wgt_wrk=zero
-
-! Equatorial patch
-! Adjoint of central patch blending on left/right sides of patch
-  allocate(p1wgt(ny,nx)); p1wgt=one
-
-  do i=1,ndx2
-     i1=ndx2+1-i
-     i2=nx-ndx2+i
-     do j=1,ny
-        p1all(j,i) =p1all(j,i) *bl(i1)   ! left (west) blending zone
-        p1all(j,i2)=p1all(j,i2)*bl(i)    ! right (east) blending zone
-        p1wgt(j,i) =p1wgt(j,i) *bl(i1)
-        p1wgt(j,i2)=p1wgt(j,i2)*bl(i)
-     end do
-  end do
-
-! bl2 of p1
-  do i=1,nx
-     do j=1,nmix
-        p1all(j,i)=p1all(j,i)*bl2(nmixp-j)
-        p1wgt(j,i)=p1wgt(j,i)*bl2(nmixp-j)
-     end do
-     do j=nymx+1,ny
-        p1all(j,i)=p1all(j,i)*bl2(j-nymx)
-        p1wgt(j,i)=p1wgt(j,i)*bl2(j-nymx)
-     end do
-  end do
-
-! zero output array
-  grid_wrk=zero
-  wgt_wrk=zero
-
-! Adjoint of transfer between central band and full grid (p1 --> work)
-  do i=1,ndx
-     i1=i-ndx+nlon
-     i2=nx-ndx+i
-     do j=1,ny
-        j1=j+ndy
-        grid_wrk(j1,i1)=grid_wrk(j1,i1)+p1all(j,i)  ! left  (west) blending zone
-        grid_wrk(j1,i )=grid_wrk(j1,i )+p1all(j,i2) ! right (east) blending zone
-        wgt_wrk (j1,i1)=wgt_wrk (j1,i1)+p1wgt(j,i)  ! left  (west) blending zone
-        wgt_wrk (j1,i )=wgt_wrk (j1,i )+p1wgt(j,i2) ! right (east) blending zone
-     end do
-  end do
-
-! Middle zone (no blending)
-  do i=ndx+1,nx-ndx
-     i1=i-ndx
-     do j=1,ny
-        j1=j+ndy
-        grid_wrk(j1,i1)=grid_wrk(j1,i1)+p1all(j,i)
-        wgt_wrk (j1,i1)=wgt_wrk (j1,i1)+p1wgt(j,i)
-     end do
-  end do
-
-  deallocate(p1wgt)
-
-! Adjoint of North pole patch(p2) -- blending and transfer to grid
-! Adjoint of South pole patch(p3) -- blending and transfer to grid
-  allocate(p2wgt(nlon+1,mr:nr)); p2wgt=one
-  allocate(p3wgt(nlon+1,mr:nr)); p3wgt=one
-
-  do j=nlatxb-nmix,nlatxb-1
-! Adjoint of blending
-     do i=1,nlon
-        p2all(i,nlat-j)=p2all(i,nlat-j)*bl2(nlatxb-j)
-        p2wgt(i,nlat-j)=p2wgt(i,nlat-j)*bl2(nlatxb-j)
-     end do
-  end do
-
-  do j=nrmxb+1,nrmxb+nmix
-!    Adjoint of blending
-     do i=1,nlon
-        p3all(i,j)=p3all(i,j)*bl2(j-nrmxb)
-        p3wgt(i,j)=p3wgt(i,j)*bl2(j-nrmxb)
-     end do
-  end do
-
-  do i=1,nlon
-!    Adjoint of transfer
-     do j=mr,nrmxb+nmix
-        grid_wrk(j+1,i)=grid_wrk(j+1,i)+p3all(i,j)
-        wgt_wrk (j+1,i)=wgt_wrk (j+1,i)+p3wgt(i,j)
-     end do
-     do j=nlatxb-nmix,nlat-mr
-        grid_wrk(j  ,i)=grid_wrk(j  ,i)+p2all(i,nlat-j)
-        wgt_wrk (j  ,i)=wgt_wrk (j  ,i)+p2wgt(i,nlat-j)
-     end do
-  end do
-
-  deallocate(p2wgt,p3wgt)
-
-  where(wgt_wrk>zero) grid_wrk=grid_wrk/wgt_wrk
-
-  deallocate(wgt_wrk)
-
-end subroutine vpatch2grid
+    use smooth_polcarf, only: smooth_polcasv
+    use fgrid2agrid_mod, only: fgrid2agrid
+
+    implicit none
+
+    real(r_kind),dimension(nlat,nlon),intent(  out) :: grid_wrk
+
+    real(r_kind)                     ,intent(in   ) :: hflt_all(pf2ap1%nlatf,pf2ap1%nlonf)
+    real(r_kind)                     ,intent(in   ) :: hflt_p2 (pf2ap2%nlatf,pf2ap2%nlonf)
+    real(r_kind)                     ,intent(in   ) :: hflt_p3 (pf2ap3%nlatf,pf2ap3%nlonf)
+
+    real(r_kind),allocatable,dimension(:,:) :: wgt_wrk
+    real(r_kind),allocatable,dimension(:,:) :: p1wgt
+    real(r_kind),allocatable,dimension(:,:) :: p2wgt,p3wgt
+
+    integer(i_kind):: i,i1,i2,j,j1
+
+    allocate(wgt_wrk(nlat,nlon))
+
+    call fgrid2agrid(pf2ap1,hflt_all,p1all) !analysis to filter grid_wrk
+
+    call fgrid2agrid(pf2ap2,hflt_p2,p2pol) !analysis to filter grid_wrk
+    call smooth_polcasv(p2pol,p2all)
+
+    call fgrid2agrid(pf2ap3,hflt_p3,p3pol) !analysis to filter grid_wrk
+    call smooth_polcasv(p3pol,p3all)
+
+    wgt_wrk=zero
+
+!   Equatorial patch
+!   Adjoint of central patch blending on left/right sides of patch
+    allocate(p1wgt(ny,nx)); p1wgt=one
+
+    do i=1,ndx2
+       i1=ndx2+1-i
+       i2=nx-ndx2+i
+       do j=1,ny
+          p1all(j,i) =p1all(j,i) *bl(i1)   ! left (west) blending zone
+          p1all(j,i2)=p1all(j,i2)*bl(i)    ! right (east) blending zone
+          p1wgt(j,i) =p1wgt(j,i) *bl(i1)
+          p1wgt(j,i2)=p1wgt(j,i2)*bl(i)
+       end do
+    end do
+
+!   bl2 of p1
+    do i=1,nx
+       do j=1,nmix
+          p1all(j,i)=p1all(j,i)*bl2(nmixp-j)
+          p1wgt(j,i)=p1wgt(j,i)*bl2(nmixp-j)
+       end do
+       do j=nymx+1,ny
+          p1all(j,i)=p1all(j,i)*bl2(j-nymx)
+          p1wgt(j,i)=p1wgt(j,i)*bl2(j-nymx)
+       end do
+    end do
+
+!   zero output array
+    grid_wrk=zero
+    wgt_wrk=zero
+
+!   Adjoint of transfer between central band and full grid (p1 --> work)
+    do i=1,ndx
+       i1=i-ndx+nlon
+       i2=nx-ndx+i
+       do j=1,ny
+          j1=j+ndy
+          grid_wrk(j1,i1)=grid_wrk(j1,i1)+p1all(j,i)  ! left  (west) blending zone
+          grid_wrk(j1,i )=grid_wrk(j1,i )+p1all(j,i2) ! right (east) blending zone
+          wgt_wrk (j1,i1)=wgt_wrk (j1,i1)+p1wgt(j,i)  ! left  (west) blending zone
+          wgt_wrk (j1,i )=wgt_wrk (j1,i )+p1wgt(j,i2) ! right (east) blending zone
+       end do
+    end do
+
+!   Middle zone (no blending)
+    do i=ndx+1,nx-ndx
+       i1=i-ndx
+       do j=1,ny
+          j1=j+ndy
+          grid_wrk(j1,i1)=grid_wrk(j1,i1)+p1all(j,i)
+          wgt_wrk (j1,i1)=wgt_wrk (j1,i1)+p1wgt(j,i)
+       end do
+    end do
+
+    deallocate(p1wgt)
+
+!   Adjoint of north pole patch(p2) -- blending and transfer to grid
+!   Adjoint of south pole patch(p3) -- blending and transfer to grid
+    allocate(p2wgt(nlon+1,mr:nr)); p2wgt=one
+    allocate(p3wgt(nlon+1,mr:nr)); p3wgt=one
+
+    do j=nlatxb-nmix,nlatxb-1
+!      Adjoint of blending
+       do i=1,nlon
+          p2all(i,nlat-j)=p2all(i,nlat-j)*bl2(nlatxb-j)
+          p2wgt(i,nlat-j)=p2wgt(i,nlat-j)*bl2(nlatxb-j)
+       end do
+    end do
+
+    do j=nrmxb+1,nrmxb+nmix
+!      Adjoint of blending
+       do i=1,nlon
+          p3all(i,j)=p3all(i,j)*bl2(j-nrmxb)
+          p3wgt(i,j)=p3wgt(i,j)*bl2(j-nrmxb)
+       end do
+    end do
+
+    do i=1,nlon
+!      Adjoint of transfer
+       do j=mr,nrmxb+nmix
+          grid_wrk(j+1,i)=grid_wrk(j+1,i)+p3all(i,j)
+          wgt_wrk (j+1,i)=wgt_wrk (j+1,i)+p3wgt(i,j)
+       end do
+       do j=nlatxb-nmix,nlat-mr
+          grid_wrk(j  ,i)=grid_wrk(j  ,i)+p2all(i,nlat-j)
+          wgt_wrk (j  ,i)=wgt_wrk (j  ,i)+p2wgt(i,nlat-j)
+       end do
+    end do
+
+    deallocate(p2wgt,p3wgt)
+
+    where(wgt_wrk>zero) grid_wrk=grid_wrk/wgt_wrk
+
+    deallocate(wgt_wrk)
+
+  end subroutine vpatch2grid
 
 end module patch2grid_mod
diff --git a/src/gsi/pcgsoi.f90 b/src/gsi/pcgsoi.f90
index 31b1c6ac2d..69ba49bfd7 100644
--- a/src/gsi/pcgsoi.f90
+++ b/src/gsi/pcgsoi.f90
@@ -14,7 +14,7 @@ module pcgsoimod
 !   2009-08-12  lueken  - update documentation
 !   2009-09-17  parrish - add bkerror_a_en and anbkerror_reg_a_en for hybrid ensemble control variable a_en
 !   2014-12-03  derber - thread dot products and modify so obsdiag can be turned off
-!   2018-08-10  guo     - removed m_obsHeadBundle references
+!   2018-08-10  guo     - removed m_obsheadbundle references
 !                       - replaced stpjo_setup() with a new stpjomod::stpjo_setup()
 !
 ! subroutines included:
@@ -30,8 +30,8 @@ module pcgsoimod
 
 implicit none
 
-PRIVATE
-PUBLIC pcgsoi
+private
+public pcgsoi
 
 contains
 
@@ -69,13 +69,13 @@ subroutine pcgsoi()
 !   2005-09-29  kleist,parrish - include _t (time derivatives) array
 !   2006-04-06  treadon - move bias cor. and tskin update into update_guess
 !   2006-04-16  derber - change call to stpcalc - move stepsize print to stpcalc
-!   2006-04-21  kleist - add calls to update Jc terms
+!   2006-04-21  kleist - add calls to update jc terms
 !   2006-05-26  derber - modify to improve convergence testing
 !   2006-07-28  derber - remove calls to makeobs
 !   2006-08-04  parrish - add changes for strong constraint option
 !   2007-03-09  su      - add option for observation perturbation
 !   2007-04-13  tremolet - use control vectors and state vectors                         
-!   2007-05-15  todling - add AGCM/TLM/ADM tester; fix one-ob case NaN's
+!   2007-05-15  todling - add agcm/tlm/adm tester; fix one-ob case nan's
 !   2007-07-05  todling - allow 4dvar to write out increment
 !   2007-09-30  todling - add timer
 !   2008-03-24  wu      - oberror tuning
@@ -97,14 +97,14 @@ subroutine pcgsoi()
 !   2011-04-25  el akkraoui - add option for re-orthogonalization.
 !   2011-07-10  todling - minor fixes for general precision handling. 
 !   2011-11-17  kleist - add handling for separate state vector for ensemble bits (hybrid ens/var)
-!   2013-01-26  parrish - WCOSS debug compile flags type mismatch for calls to ensctl2state_ad
+!   2013-01-26  parrish - wcoss debug compile flags type mismatch for calls to ensctl2state_ad
 !                          and ensctl2state.  I put in temporary fix to allow debug compile
 !                          by replacing mval with mval(1).  This is likely not
 !                          correct for multiple obs bins.
 !   2014-10-25  todling - reposition final clean to allow proper complition of 4dvar
 !   2014-12-22  Hu      -  add option i_gsdcldanal_type to control cloud analysis  
 !   2016-03-02  s.liu/carley  - remove use_reflectivity and use i_gsdcldanal_type 
-!   2016-03-25  todling - beta-mult param now within cov (following Dave Parrish corrections)
+!   2016-03-25  todling - beta-mult param now within cov (following dave parrish corrections)
 !   2016-05-13  parrish -  remove beta12mult.  Replace with sqrt_beta_s_mult, sqrt_beta_e_mult, inside
 !                          bkerror and bkerror_a_en.
 !
@@ -130,7 +130,7 @@ subroutine pcgsoi()
        iguess,read_guess_solution,diag_precon,step_start, &
        niter_no_qc,print_diag_pcg,lgschmidt
   use gsi_4dvar, only: nobs_bins, nsubwin, l4dvar, iwrtinc, ladtest, &
-                       iorthomax
+       iorthomax
   use gridmod, only: twodvar_regional
   use constants, only: zero,one,tiny_r_kind
   use anberror, only: anisotropic
@@ -149,7 +149,7 @@ subroutine pcgsoi()
   use xhat_vordivmod, only : xhat_vordiv_init, xhat_vordiv_calc, xhat_vordiv_clean
   use timermod, only: timer_ini,timer_fnl
   use projmethod_support, only: init_mgram_schmidt, &
-                                mgram_schmidt,destroy_mgram_schmidt
+       mgram_schmidt,destroy_mgram_schmidt
   use hybrid_ensemble_parameters,only : l_hyb_ens,aniso_a_en,ntlevs_ens
   use hybrid_ensemble_isotropic, only: bkerror_a_en
   use gsi_bundlemod, only : gsi_bundle
@@ -240,21 +240,21 @@ subroutine pcgsoi()
 
   if(iorthomax>0) then 
      allocate(cglwork(iorthomax+1))
-     DO ii=1,iorthomax+1
-        CALL allocate_cv(cglwork(ii))
+     do ii=1,iorthomax+1
+        call allocate_cv(cglwork(ii))
         cglwork(ii)=zero
-     ENDDO
+     enddo
      allocate(cglworkhat(iorthomax+1))
-     DO ii=1,iorthomax+1
-        CALL allocate_cv(cglworkhat(ii))
+     do ii=1,iorthomax+1
+        call allocate_cv(cglworkhat(ii))
         cglworkhat(ii)=zero
-     END DO
+     end do
   end if
 
 ! Perform inner iteration
   inner_iteration: do iter=0,niter(jiter)
 
-! Gradually turn on variational qc to avoid possible convergence problems
+!    Gradually turn on variational qc to avoid possible convergence problems
      if(vqc) then
         nlnqc_iter = iter >= niter_no_qc(jiter)
         if(jiter == jiterstart) then
@@ -280,10 +280,10 @@ subroutine pcgsoi()
            mval(1)=eval(1)
         end if
 
-!       Perform test of AGCM TLM and ADM
+!       Perform test of agcm tlm and adm
         call gsi_4dcoupler_grtests(mval,sval,nsubwin,nobs_bins)
 
-!       Run TL model to fill sval
+!       Run tl model to fill sval
         call model_tl(mval,sval,llprt)
      else
         if (l_hyb_ens) then
@@ -358,9 +358,9 @@ subroutine pcgsoi()
 !    Re-orthonormalization if requested
      if(iorthomax>0) then 
         iortho=min(iorthomax,iter) 
-        if(iter .ne. 0) then 
+        if(iter /= 0) then 
            do ii=iortho,1,-1
-              zdla = DOT_PRODUCT(gradx,cglworkhat(ii))
+              zdla = dot_product(gradx,cglworkhat(ii))
               do i=1,nclen
                  gradx%values(i) = gradx%values(i) - zdla*cglwork(ii)%values(i)
               end do
@@ -377,7 +377,7 @@ subroutine pcgsoi()
      end if
 
 !    If hybrid ensemble run, then multiply ensemble control variable a_en 
-!                                    by its localization correlation
+!    by its localization correlation
      if(l_hyb_ens) then
         if(aniso_a_en) then
     !      call anbkerror_a_en(gradx,grady)    !  not available yet
@@ -409,7 +409,7 @@ subroutine pcgsoi()
      dprod(1) = qdot_prod_sub(gradx,grady)
      if(diag_precon)then
         if (lanlerr) then
-! xdiff used as a temporary array
+!          xdiff used as a temporary array
            do i=1,nclen
               xdiff%values(i)=vprecond(i)*gradx%values(i)
            end do
@@ -424,7 +424,7 @@ subroutine pcgsoi()
            end do
            dprod(2) = qdot_prod_sub(xdiff,grady)
            dprod(3) = qdot_prod_sub(ydiff,gradx)
-! xdiff used as a temporary array
+!          xdiff used as a temporary array
            do i=1,nclen
               xdiff%values(i)=vprecond(i)*gradx%values(i)
            end do
@@ -471,7 +471,7 @@ subroutine pcgsoi()
 
 !    Calculate new search direction
      if (.not. restart) then
-       if(.not. lanlerr)then
+        if(.not. lanlerr)then
            do i=1,nclen
               xdiff%values(i)=gradx%values(i)
               ydiff%values(i)=grady%values(i)
@@ -489,7 +489,7 @@ subroutine pcgsoi()
            end do
         end if
      else
-!    If previous solution available, transfer into local arrays.
+!       If previous solution available, transfer into local arrays.
         if( .not. lanlerr)then
            xdiff=zero
            ydiff=zero
@@ -612,9 +612,9 @@ subroutine pcgsoi()
      pennorm=penx
 
   end do inner_iteration
-!  End of inner iteration
+! End of inner iteration
 
-!  Deallocate space for renormalization
+! Deallocate space for renormalization
   if(iorthomax>0) then 
      do ii=1,iorthomax+1
         call deallocate_cv(cglwork(ii))
@@ -656,83 +656,83 @@ subroutine pcgsoi()
   llprt=(mype==0)
   call control2state(xhat,mval,sbias)
   if (l4dvar) then
-    if (l_hyb_ens) then
-       call ensctl2state(xhat,mval(1),eval)
-       mval(1)=eval(1)
-    end if
-    call model_tl(mval,sval,llprt)
+     if (l_hyb_ens) then
+        call ensctl2state(xhat,mval(1),eval)
+        mval(1)=eval(1)
+     end if
+     call model_tl(mval,sval,llprt)
   else
-    if (l_hyb_ens) then
-       call ensctl2state(xhat,mval(1),eval)
-       do ii=1,nobs_bins
-          sval(ii)=eval(ii)
-       end do
-    else
-       do ii=1,nobs_bins
-          sval(ii)=mval(1)
-       end do
-    end if
+     if (l_hyb_ens) then
+        call ensctl2state(xhat,mval(1),eval)
+        do ii=1,nobs_bins
+           sval(ii)=eval(ii)
+        end do
+     else
+        do ii=1,nobs_bins
+           sval(ii)=mval(1)
+        end do
+     end if
   end if
 
   if(print_diag_pcg)then
 
-! Evaluate final cost function and gradient
+!    Evaluate final cost function and gradient
      if (mype==0) write(6,*)'Minimization final diagnostics'
 
      do ii=1,nobs_bins
-       rval(ii)=zero
+        rval(ii)=zero
      end do
      call intall(sval,sbias,rval,rbias)
      gradx=zero
      if (l4dvar) then
-       call model_ad(mval,rval,llprt)
-       if (l_hyb_ens) then
-          eval(1)=mval(1)
-          call ensctl2state_ad(eval,mval(1),gradx)
-       end if
+        call model_ad(mval,rval,llprt)
+        if (l_hyb_ens) then
+           eval(1)=mval(1)
+           call ensctl2state_ad(eval,mval(1),gradx)
+        end if
      else
-       if (l_hyb_ens) then
-          do ii=1,nobs_bins
-            eval(ii)=rval(ii)
-          end do
-          call ensctl2state_ad(eval,mval(1),gradx)
-       else
-          mval(1)=rval(1)
-          if (nobs_bins > 1 ) then
-             do ii=2,nobs_bins
-                call self_add(mval(1),rval(ii))
-             enddo
-          end if
-       end if
+        if (l_hyb_ens) then
+           do ii=1,nobs_bins
+              eval(ii)=rval(ii)
+           end do
+           call ensctl2state_ad(eval,mval(1),gradx)
+        else
+           mval(1)=rval(1)
+           if (nobs_bins > 1 ) then
+              do ii=2,nobs_bins
+                 call self_add(mval(1),rval(ii))
+              enddo
+           end if
+        end if
      end if
      call control2state_ad(mval,rbias,gradx)
   
 !    Add contribution from background term
      do i=1,nclen
-       gradx%values(i)=gradx%values(i)+yhatsave%values(i)
+        gradx%values(i)=gradx%values(i)+yhatsave%values(i)
      end do
   
 !    Multiply by background error
      if(anisotropic) then
-       call anbkerror(gradx,grady)
+        call anbkerror(gradx,grady)
      else
-       call bkerror(gradx,grady)
+        call bkerror(gradx,grady)
      end if
   
 !    If hybrid ensemble run, then multiply ensemble control variable a_en 
-!                                    by its localization correlation
+!    by its localization correlation
      if(l_hyb_ens) then
-       if(aniso_a_en) then
-    !     call anbkerror_a_en(gradx,grady)    !  not available yet
-          write(6,*)' ANBKERROR_A_EN not written yet, program stops'
-          stop
-       else
-          call bkerror_a_en(gradx,grady)
-       end if
+        if(aniso_a_en) then
+    !      call anbkerror_a_en(gradx,grady)    !  not available yet
+           write(6,*)' ANBKERROR_A_EN not written yet, program stops'
+           stop
+        else
+           call bkerror_a_en(gradx,grady)
+        end if
 
      end if
 
-! Print final Jo table
+!    Print final Jo table
      zgend=dot_product(gradx,grady,r_quad)
 !    nprt=2
 !    call evaljo(zjo,iobs,nprt,llouter)
@@ -740,7 +740,7 @@ subroutine pcgsoi()
 
      if(l_hyb_ens) then
 
-!       If hybrid ensemble run, compute contribution to Jb and Je separately
+!       If hybrid ensemble run, compute contribution to jb and je separately
 
         fjcost_e=   dot_product(xhatsave,yhatsave,r_quad,'cost_e')
 !       fjcost(1) = dot_product(xhatsave,yhatsave,r_quad,'cost_b')
@@ -755,7 +755,7 @@ subroutine pcgsoi()
      if (mype==0) then
         if(l_hyb_ens) then
 
-!          If hybrid ensemble run, print out contribution to Jb and Je separately
+!          If hybrid ensemble run, print out contribution to jb and je separately
 
            write(iout_iter,999)'costterms Jb,Je,Jo,Jc,Jl =',jiter,iter,fjcostnew(1)- fjcost_e, &
                fjcost_e,fjcostnew(2:4)
@@ -775,7 +775,7 @@ subroutine pcgsoi()
         endif
 
 
-!    Print final diagnostics
+!       Print final diagnostics
         write(6,888)'Final   cost function=',zfend
         write(6,888)'Final   gradient norm=',sqrt(zgend)
         if(zfini>tiny_r_kind) write(6,888)'Final/Initial cost function=',zfend/zfini
@@ -795,24 +795,24 @@ subroutine pcgsoi()
 ! cloud analysis  after iteration
 ! if(jiter == miter .and. i_gsdcldanal_type==1) then
   if(jiter == miter) then
-    if(i_gsdcldanal_type==2) then
-       call gsdcloudanalysis4nmmb(mype)
-    else if(i_gsdcldanal_type==1) then
-       call gsdcloudanalysis(mype)
-    else if(i_gsdcldanal_type==30) then
-       call gsdcloudanalysis4gfs(mype)
-    endif
+     if(i_gsdcldanal_type==2) then
+        call gsdcloudanalysis4nmmb(mype)
+     else if(i_gsdcldanal_type==1) then
+        call gsdcloudanalysis(mype)
+     else if(i_gsdcldanal_type==30) then
+        call gsdcloudanalysis4gfs(mype)
+     endif
   endif
 
 ! Write output analysis files
   if(.not.l4dvar) call prt_guess('analysis')
   call prt_state_norms(sval(1),'increment')
   if (twodvar_regional) then
-      call write_all(-1)
-    else
-      if(jiter == miter) then
-         call write_all(-1)
-      endif
+     call write_all(-1)
+  else
+     if(jiter == miter) then
+        call write_all(-1)
+     endif
   endif
 
 ! Overwrite guess with increment (4d-var only, for now)
@@ -821,9 +821,9 @@ subroutine pcgsoi()
      call inc2guess(sval)
      call write_all(iwrtinc)
      call prt_guess('increment')
-     ! NOTE: presently in 4dvar, we handle the biases in a slightly inconsistent way
+     ! Note: presently in 4dvar, we handle the biases in a slightly inconsistent way
      ! as when in 3dvar - that is, the state is not updated, but the biases are.
-     ! This assumes GSI handles a single iteration of the outer loop at a time
+     ! This assumes gsi handles a single iteration of the outer loop at a time
      ! when doing 4dvar (that is, multiple iterations require stop-and-go).
      call update_bias_preds(twodvar_regional,sbias)
   endif
diff --git a/src/gsi/pcgsqrt.f90 b/src/gsi/pcgsqrt.f90
index 28336de44c..e754a57801 100644
--- a/src/gsi/pcgsqrt.f90
+++ b/src/gsi/pcgsqrt.f90
@@ -6,7 +6,7 @@ subroutine pcgsqrt(xhat,costf,gradx,itermax,nprt)
 !   prgmmr:      tremolet
 !
 ! abstract: solve inner loop of analysis equation using conjugate
-!           gradient preconditioned with sqrt(B).
+!           gradient preconditioned with sqrt(b).
 !
 ! program history log:
 !   2007-04-27  tremolet - initial code
@@ -107,7 +107,7 @@ subroutine pcgsqrt(xhat,costf,gradx,itermax,nprt)
 !  Evaluate cost and gradient
    call evaljgrad(xtry,zfk,grtry,lsavinc,0,myname)
 
-!  Get A q_k
+!  Get a q_k
    do ii=1,grtry%lencv
       grtry%values(ii)=grtry%values(ii)-grad0%values(ii)
    end do
@@ -127,11 +127,11 @@ subroutine pcgsqrt(xhat,costf,gradx,itermax,nprt)
    if(iorthomax>0) then 
       iortho=min(iter,iorthomax)
       do jj=iortho,1,-1
-        zdla = dot_product(gradx,cglwork(jj))
-        do ii=1,gradx%lencv
-           gradx%values(ii) = gradx%values(ii) - zdla*cglwork(jj)%values(ii)
-        enddo
-     enddo
+         zdla = dot_product(gradx,cglwork(jj))
+         do ii=1,gradx%lencv
+            gradx%values(ii) = gradx%values(ii) - zdla*cglwork(jj)%values(ii)
+         enddo
+      enddo
    end if 
 
 !  Save gradients for orthonormalization