Merge pull request NOAA-EMC#125 from MichaelLueken-NOAA/master

GitHub Issue NOAA-EMC#108. Fix restart and clean up pcgsoi.

src/gsi/anberror.f90

@@ -284,15 +284,15 @@ subroutine create_anberror_vars(mype)
 !$$$ end documentation block
 
     use fgrid2agrid_mod, only: create_fgrid2agrid
-    use jfunc, only: nrclen,nclen,diag_precon
+    use jfunc, only: nrclen,nclen
     use berror, only: varprd,vprecond,bnf=>nf,bnr=>nr
     use gridmod, only: nlat,nlon
     implicit none
 
     integer(i_kind),intent(in   ) :: mype
 
     allocate(varprd(max(1,nrclen)))
-    if(diag_precon)allocate(vprecond(nclen))
+    allocate(vprecond(nclen))
     allocate(an_amp(max_ngauss,nvars))
     an_amp=one/three
 
@@ -388,14 +388,13 @@ subroutine destroy_anberror_vars
 !
 !$$$
     use fgrid2agrid_mod, only: destroy_fgrid2agrid
-    use jfunc, only: diag_precon
     use berror, only: vprecond
     implicit none
 
     deallocate(an_amp)
     deallocate(afact0)
     deallocate(qvar3d)
-    if(diag_precon)deallocate(vprecond)
+    deallocate(vprecond)
 
     call destroy_fgrid2agrid(pf2aP1)
     call destroy_fgrid2agrid(pf2aP2)
@@ -431,7 +430,7 @@ subroutine create_anberror_vars_reg(mype)
 !
 !$$$
     use fgrid2agrid_mod, only: create_fgrid2agrid
-    use jfunc, only: nrclen,nclen,diag_precon
+    use jfunc, only: nrclen,nclen
     use berror, only: varprd,vprecond
     use gridmod, only: nlat,nlon,istart,jstart
     use general_commvars_mod, only: s2g_raf
@@ -442,7 +441,7 @@ subroutine create_anberror_vars_reg(mype)
     logical regional
 
     allocate(varprd(max(1,nrclen)))
-    if(diag_precon)allocate(vprecond(nclen))
+    allocate(vprecond(nclen))
     allocate(an_amp(max_ngauss,nvars))
     an_amp=one/three
 
@@ -654,15 +653,14 @@ subroutine destroy_anberror_vars_reg
 !$$$ end documentation block
 
     use fgrid2agrid_mod, only: destroy_fgrid2agrid
-    use jfunc, only: diag_precon
     use berror, only: vprecond
     use general_sub2grid_mod, only: general_sub2grid_destroy_info
     implicit none
 
     deallocate(an_amp)
     deallocate(afact0)
     deallocate(qvar3d)
-    if(diag_precon)deallocate(vprecond)
+    deallocate(vprecond)
     call destroy_fgrid2agrid(pf2aP1)
     call general_sub2grid_destroy_info(s2g_rff)
 !write(6,'(" FOR TEST ONLY--REMOVE THIS MESSAGE BEFORE FINAL COMMIT--SUCCESSFUL CALL TO general_sub2grid_destroy_info to remove s2g_rff in destroy_anberror_vars_reg")')

src/gsi/anbkerror.f90

@@ -1,4 +1,4 @@
-subroutine anbkerror(gradx,grady)
+subroutine anbkerror(grady)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    anbkerror apply anisotropic background error covariance  
@@ -33,10 +33,10 @@ subroutine anbkerror(gradx,grady)
 !   2015-07-02  pondeca - update slab mode option to work with any number of control variables
 !
 !   input argument list:
-!     gradx    - input field  
+!     grady    - input field  
 !
 !   output
-!     grady    - background structure * gradx 
+!     grady    - background structure * grady 
 !
 ! attributes:
 !   language: f90
@@ -57,7 +57,6 @@ subroutine anbkerror(gradx,grady)
   implicit none
 
 ! Declare passed variables
-  type(control_vector),intent(inout) :: gradx
   type(control_vector),intent(inout) :: grady
 
 ! Declare local variables
@@ -82,9 +81,6 @@ subroutine anbkerror(gradx,grady)
 ! Initialize timer
   call timer_ini('anbkerror')
 
-! Put things in grady first since operations change input variables
-  grady=gradx
-
 ! Since each internal vector [step(jj)] of grad has the same structure, pointers
 ! are the same independent of the subwindow jj
 call gsi_bundlegetpointer (grady%step(1),myvnames,ipnts,istatus)

src/gsi/berror.f90

@@ -249,7 +249,7 @@ subroutine create_berror_vars
 !$$$
   use balmod, only: llmin,llmax
   use gridmod, only: nlat,nlon,lat2,lon2,nsig,nnnn1o
-  use jfunc, only: nrclen,nclen,diag_precon
+  use jfunc, only: nrclen,nclen
   use constants, only: zero,one
   implicit none
 
@@ -293,7 +293,7 @@ subroutine create_berror_vars
      dssvs = zero
   endif
   allocate(varprd(nrclen))
-  if(diag_precon)allocate(vprecond(nclen))
+  allocate(vprecond(nclen))
   allocate(inaxs(nf,nlon/8),inxrs(nlon/8,mr:nr) )
 
   allocate(slw(ny*nx,nnnn1o),&
@@ -330,7 +330,6 @@ subroutine destroy_berror_vars
 !   machine:  ibm RS/6000 SP
 !
 !$$$
-    use jfunc, only: diag_precon
     implicit none
     if(allocated(table)) deallocate(table)
     deallocate(wtaxs)
@@ -342,7 +341,7 @@ subroutine destroy_berror_vars
     if(allocated(alv))   deallocate(alv)
     if(allocated(dssv))  deallocate(dssv)
     if(allocated(dssvs)) deallocate(dssvs)
-    if(diag_precon)deallocate(vprecond)
+    deallocate(vprecond)
     deallocate(slw,slw1,slw2)
     deallocate(ii,jj,ii1,jj1,ii2,jj2)
 
@@ -566,7 +565,7 @@ subroutine pcinfo
     use kinds, only: r_kind,i_kind
     use radinfo, only: ostats,rstats,varA,jpch_rad,npred,newpc4pred
     use aircraftinfo, only: aircraft_t_bc_pof,aircraft_t_bc,ntail,npredt,ostats_t,rstats_t,varA_t
-    use jfunc, only: nclen,nrclen,diag_precon,step_start,ntclen
+    use jfunc, only: nclen,nrclen,step_start,ntclen,diag_precon
     use constants, only:  zero,one
     implicit none
 
@@ -580,17 +579,18 @@ subroutine pcinfo
 !   Only diagonal elements are considered
 
 !   set a coeff. factor for variances of control variables
-    if(diag_precon)then
-      lfact=step_start
-      vprecond=lfact
+    lfact=step_start
+    vprecond=lfact
 
+    if(diag_precon)then
       if(newpc4pred)then
 !       for radiance bias predictor coeff.
         nclen1=nclen-nrclen
         ii=0
         do i=1,jpch_rad
            do j=1,npred
               ii=ii+1
+!             if (ostats(i)>zero) vprecond(nclen1+ii)=vprecond(nclen1+ii)/(one+rstats(j,i)*varprd(ii))
               if (ostats(i)>zero) vprecond(nclen1+ii)=one/(one+rstats(j,i)*varprd(ii))
               if (ostats(i)>20.0_r_kind) then
                  if (rstats(j,i)>zero) then
@@ -615,6 +615,7 @@ subroutine pcinfo
                 if (aircraft_t_bc_pof) obs_count = ostats_t(j,i)
                 if (aircraft_t_bc) obs_count = ostats_t(1,i)
 
+!               if (obs_count>zero) vprecond(nclen1+ii)=vprecond(nclen1+ii)/(one+rstats_t(j,i)*varprd(jj))
                 if (obs_count>zero) vprecond(nclen1+ii)=one/(one+rstats_t(j,i)*varprd(jj))
                 if (obs_count>3.0_r_kind) then
                    varA_t(j,i)=one/(one/varprd(jj)+rstats_t(j,i))
@@ -917,7 +918,7 @@ subroutine create_berror_vars_reg
     use constants, only: zero
     use balmod, only: llmin,llmax
     use gridmod, only: nlat,nlon,nsig,nnnn1o,lat2,lon2
-    use jfunc, only: nrclen,nclen,diag_precon
+    use jfunc, only: nrclen,nclen
     implicit none
 
     nx=nlon
@@ -940,7 +941,7 @@ subroutine create_berror_vars_reg
     endif
 
     allocate(varprd(max(1,nrclen) ) )     
-    if(diag_precon)allocate(vprecond(nclen))
+    allocate(vprecond(nclen))
 
     allocate(slw(ny*nx,nnnn1o) )
     allocate(ii(ny,nx,3,nnnn1o),jj(ny,nx,3,nnnn1o) )
@@ -973,7 +974,6 @@ subroutine destroy_berror_vars_reg
 !   machine:  ibm RS/6000 SP
 !
 !$$$
-    use jfunc, only:diag_precon
     implicit none
 
     deallocate(be,qvar3d)
@@ -984,7 +984,7 @@ subroutine destroy_berror_vars_reg
     deallocate(ii,jj)
     deallocate(slw)
     deallocate(varprd)
-    if(diag_precon)deallocate(vprecond)
+    deallocate(vprecond)
 
     return
   end subroutine destroy_berror_vars_reg

src/gsi/bicg.f90

@@ -31,8 +31,7 @@ subroutine bicg()
 use kinds,     only: r_kind,i_kind,r_quad
 use gsi_4dvar, only: l4dvar, &
                      ladtest, lgrtest, lanczosave, ltcost, nwrvecs
-use jfunc,     only: jiter,miter,niter,xhatsave,yhatsave,jiterstart, &
-                     diag_precon
+use jfunc,     only: jiter,miter,niter,xhatsave,yhatsave,jiterstart
 use constants, only: zero,tiny_r_kind
 use mpimod,    only: mype
 use obs_sensitivity, only: lobsensmin, lobsensfc, lobsensincr, &
@@ -99,23 +98,24 @@ subroutine bicg()
 if(LMPCGL) then 
    call pcgprecond(gradx,grady)
 else 
-   call bkerror(gradx,grady)
+   grady=gradx
+   call bkerror(grady)
 
    ! If hybrid ensemble run, then multiply ensemble control variable a_en 
    !                                 by its localization correlation
    if(l_hyb_ens) then
 
      if(aniso_a_en) then
-   !   call anbkerror_a_en(gradx,grady)    !  not available yet
+   !   call anbkerror_a_en(grady)    !  not available yet
        write(6,*)' ANBKERROR_A_EN not written yet, program stops'
        call stop2 (999)
      else
-       call bkerror_a_en(gradx,grady)
+       call bkerror_a_en(grady)
      end if
 
    end if
    ! Add potential additional preconditioner
-   if(diag_precon) call precond(grady)
+   call precond(grady)
 endif
 
 zg0=dot_product(gradx,grady,r_quad)
@@ -188,24 +188,25 @@ subroutine bicg()
 
 else ! not sensitivity run
 
-   call bkerror(gradf,grads)
+   grads=gradf
+   call bkerror(grads)
 
 ! If hybrid ensemble run, then multiply ensemble control variable a_en 
 !                                 by its localization correlation
    if(l_hyb_ens) then
 
      if(aniso_a_en) then
-!      call anbkerror_a_en(gradf,grads)    !  not available yet
+!      call anbkerror_a_en(grads)    !  not available yet
        write(6,*)' ANBKERROR_A_EN not written yet, program stops'
        stop
      else
-       call bkerror_a_en(gradf,grads)
+       call bkerror_a_en(grads)
      end if
 
    end if
 
 !  Add potential additional preconditioner 
-   if(diag_precon) call precond(grads)
+   call precond(grads)
 
 !  Update xhatsave
    do ii=1,xhat%lencv