Fixed out of bounds and vector assignment bugs NOAA-EMC#471

hu5970 · May 19, 2023 · 8c2ead4 · 8c2ead4
1 parent f5d200d
commit 8c2ead4
Show file tree

Hide file tree

Showing 3 changed files with 31 additions and 9 deletions.
diff --git a/src/gsi/ensctl2state.f90 b/src/gsi/ensctl2state.f90
@@ -240,7 +240,6 @@ subroutine ensctl2state(xhat,mval,eval)
 !$omp section
 
 !  Get pointers to required state variables
-   call gsi_bundlegetpointer (eval(jj),'oz'  ,sv_oz , istatus)
    call gsi_bundlegetpointer (eval(jj),'sst' ,sv_sst, istatus)
    if(ls_w)then
      call gsi_bundlegetpointer (eval(jj),'w' ,sv_w, istatus)
@@ -249,7 +248,6 @@ subroutine ensctl2state(xhat,mval,eval)
      end if
    end if
 !  Copy variables
-   call gsi_bundlegetvar ( wbundle_c, 'oz' , sv_oz,  istatus )
    call gsi_bundlegetvar ( wbundle_c, 'sst', sv_sst, istatus )
    if(lc_w)then
       call gsi_bundlegetvar ( wbundle_c, 'w' , sv_w,  istatus )
@@ -258,6 +256,13 @@ subroutine ensctl2state(xhat,mval,eval)
       end if
    end if
 
+!  Get the ozone vector if it is defined
+   id=getindex(cvars3d,"oz")
+   if(id > 0) then
+      call gsi_bundlegetpointer (eval(jj),'oz'  ,sv_oz , istatus)
+      call gsi_bundlegetvar ( wbundle_c, 'oz' , sv_oz,  istatus )
+   endif
+
 !$omp end parallel sections
 
 ! Add contribution from static B, if necessary

diff --git a/src/gsi/ensctl2state_ad.f90 b/src/gsi/ensctl2state_ad.f90
@@ -206,9 +206,7 @@ subroutine ensctl2state_ad(eval,mval,grad)
 
 !$omp section
 
-   call gsi_bundlegetpointer (eval(jj),'oz'  ,rv_oz , istatus)
    call gsi_bundlegetpointer (eval(jj),'sst' ,rv_sst, istatus)
-   call gsi_bundleputvar ( wbundle_c, 'oz',  rv_oz,  istatus )
    call gsi_bundleputvar ( wbundle_c, 'sst', rv_sst, istatus )
    if(wdw_exist)then
      call gsi_bundlegetpointer (eval(jj),'w' ,rv_w, istatus)
@@ -219,6 +217,13 @@ subroutine ensctl2state_ad(eval,mval,grad)
      end if
    end if
 
+!  Get the ozone vector if it is defined
+   id=getindex(cvars3d,"oz")
+   if(id > 0) then
+      call gsi_bundlegetpointer (eval(jj),'oz'  ,rv_oz , istatus)
+      call gsi_bundleputvar ( wbundle_c, 'oz',  rv_oz,  istatus )
+   endif
+
 !$omp section
 
    if (do_cw_to_hydro_ad .and. .not.do_cw_to_hydro_ad_hwrf) then

diff --git a/src/gsi/stpcalc.f90 b/src/gsi/stpcalc.f90
@@ -263,7 +263,7 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, &
   real(r_quad),parameter:: one_tenth_quad = 0.1_r_quad 
 
 ! Declare local variables
-  integer(i_kind) i,j,mm1,ii,iis,ibin,ipenloc,it
+  integer(i_kind) i,j,mm1,ii,iis,final_ii,ibin,ipenloc,it
   integer(i_kind) istp_use,nstep,nsteptot,kprt
   real(r_quad),dimension(4,ipen):: pbc
   real(r_quad),dimension(4,nobs_type):: pbcjo 
@@ -299,6 +299,7 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, &
   kprt=3
   pjcalc=.false.
   pj=zero_quad
+  final_ii=1
 
 !   Begin calculating contributions to penalty and stepsize for various terms
 !
@@ -779,11 +780,13 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, &
               write(iout_iter,*) ' early termination due to cx or stp  <=0 ',cx,stp(ii)
               write(iout_iter,*) ' better stepsize found',cx,stp(ii)
            end if
+           final_ii=ii
            exit stepsize
         else if(ii == istp_iter)then
            write(iout_iter,*) ' early termination due to no decrease in penalty ',cx,stp(ii)
            stp(istp_use)=zero
            end_iter = .true.
+           final_ii=ii
            exit stepsize
         else
 !       Try different (better?) stepsize
@@ -808,12 +811,16 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, &
            end_iter = .true.
 !          Finalize timer
            call timer_fnl('stpcalc')
+           final_ii=ii
            exit stepsize
         end if
 !       Check for convergence in stepsize estimation
         stprat(ii)=zero
         if(stp(ii) > zero_quad)stprat(ii)=abs((stp(ii)-stp(ii-1))/stp(ii))
-        if(stprat(ii) < 1.e-4_r_kind) exit stepsize
+        if(stprat(ii) < 1.e-4_r_kind) then
+           final_ii=ii
+           exit stepsize
+        end if
         dels = one_tenth_quad*dels
      end if
 
@@ -840,16 +847,21 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, &
               istp_use=i
            end if
         end do
-        if(istp_use /= istp_iter)exit stepsize
+        if(istp_use /= istp_iter) then
+           final_ii=ii
+           exit stepsize
+        end if
 !       If no best stepsize set to zero and end minimization
         if(mype == minmype)then
            write(iout_iter,141)(outpen(i),i=1,nsteptot)
         end if
         end_iter = .true.
         stp(ii)=zero_quad
         istp_use=ii
+        final_ii=ii
         exit stepsize
      end if
+     final_ii=ii
   end do stepsize
   if(kprt >= 2 .and. iter == 0)then
      call mpl_allreduce(ipen,nobs_bins,pj)
@@ -880,15 +892,15 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, &
 
      if(print_verbose)then
         write(iout_iter,200) (stp(i),i=0,istp_use)
-        write(iout_iter,199) (stprat(ii),ii=1,istp_use)
+        write(iout_iter,199) (stprat(i),i=1,istp_use)
         write(iout_iter,201) (outstp(i),i=1,nsteptot)
         write(iout_iter,202) (outpen(i)-outpen(4),i=1,nsteptot)
      end if
   end if
 ! Check for final stepsize negative (probable error)
   if(stpinout <= zero)then
      if(mype == minmype)then
-        write(iout_iter,130) ii,bx,cx,stp(ii)
+        write(iout_iter,130) final_ii,bx,cx,stp(final_ii)
         write(iout_iter,105) (bsum(i),i=1,ipen)
         write(iout_iter,110) (csum(i),i=1,ipen)
         write(iout_iter,101) (pbc(1,i)-pen_est(i),i=1,ipen)