Optimize GSI

src/gsi/adjtest.f90

@@ -33,10 +33,12 @@ module adjtest
 use control_vectors, only: control_vector,allocate_cv,random_cv, &
     deallocate_cv,dot_product,assignment(=)
 use state_vectors, only: allocate_state,deallocate_state,dot_product
+use gridmod, only : minmype
 use gsi_bundlemod, only: gsi_bundle
 use gsi_bundlemod, only: assignment(=)
 use bias_predictors, only: predictors,allocate_preds,deallocate_preds, &
     assignment(=)
+use control2state_mod, only: control2state,c2sset,control2state_ad
 
 implicit none
 private
@@ -81,7 +83,7 @@ subroutine adtest(xhat)
 integer(i_kind) :: ii,idig
 real(r_kind) :: zz1,zz2,zz3
 
-if (mype==0) write(6,*)'ADTEST starting'
+if (mype==minmype) write(6,*)'ADTEST starting'
 
 ! ----------------------------------------------------------------------
 ! Allocate local variables
@@ -97,10 +99,10 @@ subroutine adtest(xhat)
 ! Initialize control space vectors
 if (present(xhat)) then
    xtest1=xhat
-   if (mype==0) write(6,*)'ADTEST use input xhat'
+   if (mype==minmype) write(6,*)'ADTEST use input xhat'
 else
    call random_cv(xtest1)
-   if (mype==0) write(6,*)'ADTEST use random_cv(xhat)'
+   if (mype==minmype) write(6,*)'ADTEST use random_cv(xhat)'
 endif
 xtest2=zero
 
@@ -135,18 +137,20 @@ subroutine adtest(xhat)
 do ii=1,nsubwin
    zz2=zz2+dot_product(stest1(ii),stest1(ii))
 enddo
-DO ii=1,nrclen
+do ii=1,nrclen
    zz2=zz2+sbias1%values(ii)*sbias1%values(ii)
-ENDDO
+enddo
 
-if ( abs(zz1+zz2) > sqrt(tiny(zz3)) ) then
-   zz3=two*abs(zz1-zz2)/(zz1+zz2)
-else
-   zz3=abs(zz1-zz2)
-endif
-idig= int(-log(zz3+tiny(zz3))/log(10.0_r_kind))
+if (mype==minmype) then
+   if ( abs(zz1+zz2) > sqrt(tiny(zz3)) ) then
+      zz3=two*abs(zz1-zz2)/(zz1+zz2)
+   else
+      zz3=abs(zz1-zz2)
+   end if
+   idig= int(-log(zz3+tiny(zz3))/log(10.0_r_kind))
 
-if (mype==0) then
+! Note that this result is not completely correct especially on processors
+! other than minmype.  See issue 548.
    write(6,'(A)')' ADTEST            0.123456789012345678'
    write(6,'(A,ES25.18)')' ADTEST <F*F.Y,X>= ',zz1
    write(6,'(A,ES25.18)')' ADTEST <F.Y,F.Y>= ',zz2
@@ -166,7 +170,7 @@ subroutine adtest(xhat)
 call deallocate_preds(sbias2)
 ! ----------------------------------------------------------------------
 
-if (mype==0) write(6,*)'ADTEST finished'
+if (mype==minmype) write(6,*)'ADTEST finished'
 
 return
 end subroutine adtest

src/gsi/adjtest_obs.f90

@@ -78,6 +78,7 @@ subroutine adtest_obs
 
   use m_obsdiags, only: obsLLists
   use m_obsLList, only: obsLList_getTLDdotprod
+  use control2state_mod, only: control2state
 
   implicit none
 

src/gsi/aeroinfo.f90

@@ -313,12 +313,10 @@ subroutine aeroinfo_read
 
 !      Successful read, return to calling routine
      else
-!       File does not exist, write warning message to alert users
+!       File does not exist, write warning message to unit 6 to alert users
         if (mype==mype_aero) then
-           open(iout_aero)
-           write(iout_aero,*)'AEROINFO_READ:  ***WARNING*** FILE ',trim(fname),' does not exist'
-           write(iout_aero,*)'AEROINFO_READ:  jpch_aero=',jpch_aero
-           close(iout_aero)
+           write(6,*)'AEROINFO_READ:  ***WARNING*** FILE ',trim(fname),' does not exist'
+           write(6,*)'AEROINFO_READ:  jpch_aero=',jpch_aero
         endif
      end if
 

src/gsi/atms_spatial_average_mod.f90

@@ -153,7 +153,7 @@ SUBROUTINE ATMS_Spatial_Average(Num_Obs, NChanl, FOV, Time, BT_InOut, &
        Scanline_Back(FOV(I),Scanline(I))=I
     END DO
 
-!$omp parallel do schedule(dynamic,1) private(ichan,iscan,ios,ifov)
+!$omp parallel do schedule(dynamic,1) private(i,ichan,iscan,ios,ifov)
     DO IChan=1,nchanl
 
        err(ichan)=0

src/gsi/bicg.f90

@@ -30,7 +30,7 @@ subroutine bicg()
 
 use kinds,     only: r_kind,i_kind,r_quad
 use gsi_4dvar, only: l4dvar, &
-                     ladtest, lgrtest, lanczosave, ltcost, nwrvecs
+                     ladtest, lgrtest, lanczosave, ltcost, nwrvecs, lsqrtb
 use jfunc,     only: jiter,miter,niter,xhatsave,yhatsave,jiterstart
 use constants, only: zero,tiny_r_kind
 use mpimod,    only: mype
@@ -39,6 +39,7 @@ subroutine bicg()
 use obsmod,    only: lsaveobsens,l_do_adjoint
 use adjtest,   only: adtest
 use grdtest,   only: grtest
+use gsi_bundlemod, only : gsi_bundlegetpointer
 use control_vectors, only: control_vector
 use control_vectors, only: allocate_cv,deallocate_cv,write_cv,inquire_cv
 use control_vectors, only: dot_product,assignment(=)
@@ -89,6 +90,13 @@ subroutine bicg()
 call allocate_cv(gradf)
 call allocate_cv(grads)
 
+if(l_hyb_ens .and. .not. aniso_a_en) then
+   if (lsqrtb) then
+      write(6,*)'l_hyb_ens: not for use with lsqrtb'
+      call stop2(317)
+   end if
+end if
+
 
 ! Get initial cost function and gradient
 nprt=2

src/gsi/bicglanczos.F90

@@ -57,13 +57,14 @@ module bicglanczos
 use constants, only : zero, one, half,two, zero_quad,tiny_r_kind
 use timermod , only : timer_ini, timer_fnl
 use lanczos  , only : save_precond
-use gsi_4dvar, only : iorthomax
+use gsi_4dvar, only : iorthomax,lsqrtb
 use control_vectors, only: control_vector
 use control_vectors, only: allocate_cv,deallocate_cv,inquire_cv
 use control_vectors, only: read_cv,write_cv
 use control_vectors, only: dot_product,assignment(=)
 use gsi_bundlemod, only: gsi_bundle
 use gsi_bundlemod, only: assignment(=)
+use gsi_bundlemod, only : gsi_bundlegetpointer
 use mpimod  ,  only : mpi_comm_world
 use mpimod,    only: mype
 use jfunc   ,  only : iter, jiter
@@ -248,7 +249,13 @@ subroutine pcglanczos(xhat,yhat,pcost,gradx,grady,preduc,kmaxit,lsavevecs)
 if(nprt>=1.and.ltcost_) call allocate_cv(gradf)
 call allocate_cv(dirw)
 
-!--- 'zeta' is an upper bound on the relative error of the gradient.
+if(l_hyb_ens .and. .not. aniso_a_en) then
+   if (lsqrtb) then
+      write(6,*)'l_hyb_ens: not for use with lsqrtb'
+      call stop2(317)
+   end if
+end if
+ !--- 'zeta' is an upper bound on the relative error of the gradient.
 
 zeta  = 1.0e-4_r_kind
 zreqrd = preduc

src/gsi/bkerror.f90

@@ -71,7 +71,6 @@ subroutine bkerror(grady)
 
 ! Declare local variables
   integer(i_kind) i,ii
-  integer(i_kind) i_t,i_p,i_st,i_vp
   integer(i_kind) ipnts(4),istatus
 ! integer(i_kind) nval_lenz,ndim2d
   real(r_kind),pointer,dimension(:,:,:):: p_t  =>NULL()
@@ -97,11 +96,7 @@ subroutine bkerror(grady)
 ! Only need to get pointer for ii=1 - all other are the same
   call gsi_bundlegetpointer ( grady%step(1), (/'t ','sf','vp','ps'/), &
                               ipnts, istatus )
-  i_t  = ipnts(1)
-  i_st = ipnts(2)
-  i_vp = ipnts(3)
-  i_p  = ipnts(4)
-  dobal = i_t>0.and.i_p>0.and.i_st>0.and.i_vp>0
+  dobal = ipnts(1)>0 .and. ipnts(2)>0 .and. ipnts(3)>0 .and. ipnts(4)>0
 
 ! if ensemble run, multiply by sqrt_beta_s
   if(l_hyb_ens) call sqrt_beta_s_mult(grady)

src/gsi/calctends_no_tl.f90

@@ -244,28 +244,21 @@ subroutine calctends_no_tl(st,vp,t,p,mype,u_t,v_t,t_t,p_t,uvflag)
       end if
     end if
 
-!   top/bottom boundary condition:
 
     do j=jtstart(kk),jtstop(kk)
       do i=1,lat2
+
+!   top/bottom boundary condition:
+
         what(i,j,1)=zero
         what(i,j,nsig+1)=zero
-      enddo
-    enddo
-
 
 !   load actual dp/dt
 
-    do j=jtstart(kk),jtstop(kk)
-      do i=1,lat2
         p_t(i,j)=prsth(i,j,1)
-      end do
-    end do
 
 !   before big k loop, zero out the km1 summation arrays
 
-    do j=jtstart(kk),jtstop(kk)
-      do i=1,lat2
         sumkm1  (i,j)=zero
         sum2km1 (i,j)=zero
         sumvkm1 (i,j)=zero

src/gsi/compact_diffs.f90

@@ -268,7 +268,6 @@ subroutine stvp2uv(work,idim)
   integer(i_kind) ix,iy
   integer(i_kind) ny,i,j
   real(r_kind) polsu,polnu,polnv,polsv
-  real(r_kind),dimension(nlon):: grid3n,grid3s,grid1n,grid1s
   real(r_kind),dimension(nlat-2,nlon):: a,b,grid1,grid2,grid3,grid4
 
   if(idim <=1) write(6,*) ' error in call to stvp2uv ',idim
@@ -318,27 +317,17 @@ subroutine stvp2uv(work,idim)
   polnv=polnv/float(nlon)
   polsu=polsu/float(nlon)
   polsv=polsv/float(nlon)
-  do ix=1,nlon
-     grid3n(ix)= polnu*coslon(ix)+polnv*sinlon(ix)
-     grid1n(ix)=-polnu*sinlon(ix)+polnv*coslon(ix)
-     grid3s(ix)= polsu*coslon(ix)+polsv*sinlon(ix)
-     grid1s(ix)= polsu*sinlon(ix)-polsv*coslon(ix)
-  end do
 ! work(1 is u, work(2 is v
   do j=1,nlon
-     do i=1,nlat
-        if(i /= 1 .and. i /= nlat)then
-           work(1,i,j)=grid3(i-1,j)
-           work(2,i,j)=grid1(i-1,j)
-        else if(i == 1)then
-           work(1,i,j)=grid3s(j)
-           work(2,i,j)=grid1s(j)
-        else
-           work(1,i,j)=grid3n(j)
-           work(2,i,j)=grid1n(j)
-        end if
+     do i=2,nlat-1
+        work(1,i,j)=grid3(i-1,j)
+        work(2,i,j)=grid1(i-1,j)
      end do
-  enddo
+     work(1,1,j)= polsu*coslon(j)+polsv*sinlon(j)
+     work(2,1,j)= polsu*sinlon(j)-polsv*coslon(j)
+     work(1,nlat,j)= polnu*coslon(j)+polnv*sinlon(j)
+     work(2,nlat,j)= -polnu*sinlon(j)+polnv*coslon(j)
+  end do
 
   return
   end subroutine stvp2uv
@@ -749,18 +738,14 @@ subroutine tstvp2uv(work,idim)
   ny=nlat-2
 
   do j=1,nlon
-     do i=1,nlat
-        if(i /= 1 .and. i /= nlat)then
-           grid3(i-1,j)=work(1,i,j)
-           grid1(i-1,j)=work(2,i,j)
-        else if(i == 1)then
-           grid3s(j)=work(1,i,j)
-           grid1s(j)=work(2,i,j)
-        else
-           grid3n(j)=work(1,i,j)
-           grid1n(j)=work(2,i,j)
-        end if
+     do i=2,nlat-1
+        grid3(i-1,j)=work(1,i,j)
+        grid1(i-1,j)=work(2,i,j)
      end do
+     grid3s(j)=work(1,1,j)
+     grid1s(j)=work(2,1,j)
+     grid3n(j)=work(1,nlat,j)
+     grid1n(j)=work(2,nlat,j)
   end do
 
   polnu=zero
@@ -815,16 +800,15 @@ subroutine tstvp2uv(work,idim)
        nlon,ny,noq)
 !$omp end parallel sections
   do j=1,nlon
-     do i=1,nlat
-        if(i /= 1 .and. i /= nlat)then
-!          NOTE:  Adjoint of first derivative is its negative
-           work(1,i,j)=-(a(i-1,j)+d(i-1,j))
-           work(2,i,j)=-(b(i-1,j)+c(i-1,j))
-        else
-           work(1,i,j)=zero
-           work(2,i,j)=zero
-        end if
+     do i=2,nlat-1
+!       NOTE:  Adjoint of first derivative is its negative
+        work(1,i,j)=-(a(i-1,j)+d(i-1,j))
+        work(2,i,j)=-(b(i-1,j)+c(i-1,j))
      end do
+     work(1,1,j)=zero
+     work(2,1,j)=zero
+     work(1,nlat,j)=zero
+     work(2,nlat,j)=zero
   end do
 
   return

src/gsi/constants.f90

@@ -90,7 +90,7 @@ module constants
 
 ! Declare derived constants
   integer(i_kind):: huge_i_kind
-  integer(i_kind), parameter :: max_varname_length=64
+  integer(i_kind), parameter :: max_varname_length=20
   real(r_single):: tiny_single, huge_single
   real(r_kind):: xai, xa, xbi, xb, dldt, rozcon,ozcon,fv, tpwcon,eps, rd_over_g
   real(r_kind):: el2orc, g_over_rd, rd_over_cp, cpr, omeps, epsm1, factor2