Changes based on shoyokota review NOAA-EMC#2.

src/gsi/hdraobmod.f90

@@ -1272,7 +1272,7 @@ subroutine read_hdraob(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
 
 ! Write header record and data to output file for further processing
 
-  call count_obs(ndata,nreal,ilat,ilon,cdata_all,nobs)
+  call count_obs(ndata,nreal,ilat,ilon,cdata_all(1,1:ndata),nobs)
   write(lunout) obstype,sis,nreal,nchanl,ilat,ilon,ndata
   write(lunout) ((cdata_all(k,i),k=1,nreal),i=1,ndata)
 

src/gsi/obs_para.f90

@@ -211,7 +211,6 @@ subroutine disobs(ndata,nobs,mm1,lunout,obsfile,obstypeall)
 !$$$
   use kinds, only: r_kind,i_kind
   use gridmod, only: periodic_s,nlon,nlat,jlon1,ilat1,istart,jstart
-  use mpimod, only: mype
   implicit none
 
 ! Declare passed variables

src/gsi/read_avhrr.f90

@@ -570,7 +570,6 @@ subroutine read_avhrr(mype,val_avhrr,ithin,rmesh,jsatid,&
      if ( nread > 0 ) then
         write(*,'(a,a10,I3,F6.1,3I10)') 'read_avhrr,satid,imesh,amesh,itxmax,nread,ndata_mesh : ',jsatid,imesh,amesh(imesh),itxmax,nread,ndata_mesh
      endif
-     deallocate(amesh,hsst_thd)
 !
 !    get data_all by combining data from all thinning box sizes
 !
@@ -586,6 +585,7 @@ subroutine read_avhrr(mype,val_avhrr,ithin,rmesh,jsatid,&
      deallocate(data_mesh,nrec)
 
   enddo     ! do imesh = 1, nmesh
+  deallocate(amesh,hsst_thd)
 
 ! Allow single task to check for bad obs, update superobs sum,
 ! and write out data to scratch file for further processing.

src/gsi/stpcalc.f90

@@ -860,6 +860,10 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, &
 ! Check for final stepsize negative (probable error)
   if(stpinout <= zero)then
      if(mype == minmype)then
+        do i=1,ipen
+            pen_est(i)=pbc(1,i)-(stpinout-stp(0))*(2.0_r_quad*bsum(i)- &
+                       (stpinout-stp(0))*csum(i))
+        end do
         write(iout_iter,130) istp_use,bx,cx,stp(istp_use)
         write(iout_iter,105) (bsum(i),i=1,ipen)
         write(iout_iter,110) (csum(i),i=1,ipen)
@@ -873,37 +877,30 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, &
 202 format(' penalties          = ',(10(e13.6,1x)))
 
 ! If convergence or failure of stepsize calculation return
-  if (.not. end_iter) then
 
 !    Estimate terms in penalty
-     if(mype == minmype)then
-        if(print_verbose)then
-           do i=1,ipen
-               pen_est(i)=pbc(1,i)-(stpinout-stp(0))*(2.0_r_quad*bsum(i)- &
-                          (stpinout-stp(0))*csum(i))
-           end do
-           write(iout_iter,101) (pbc(1,i)-pen_est(i),i=1,ipen)
-        end if
-        pjcostnew(1) = pbc(1,1)                                  ! Jb
-        pjcostnew(3) = pbc(1,2)+pbc(1,3)                         ! Jc
-        pjcostnew(4)=zero
-        do i=4,n0
-           pjcostnew(4) =  pjcostnew(4) + pbc(1,i) ! Jl
-        end do
-        pjcostnew(2) = zero 
-        do i=1,nobs_type
-           pjcostnew(2) = pjcostnew(2)+pbc(1,n0+i)               ! Jo
-        end do
-        penaltynew=pjcostnew(1)+pjcostnew(2)+pjcostnew(3)+pjcostnew(4)
+  if(mype == minmype)then
+     pjcostnew(1) = pbc(1,1)                                  ! Jb
+     pjcostnew(3) = pbc(1,2)+pbc(1,3)                         ! Jc
+     pjcostnew(4)=zero
+     do i=4,n0
+        pjcostnew(4) =  pjcostnew(4) + pbc(1,i) ! Jl
+     end do
+     pjcostnew(2) = zero 
+     do i=1,nobs_type
+        pjcostnew(2) = pjcostnew(2)+pbc(1,n0+i)               ! Jo
+     end do
+     penaltynew=pjcostnew(1)+pjcostnew(2)+pjcostnew(3)+pjcostnew(4)
 
-        if(print_verbose)then
-           write(iout_iter,200) (stp(i),i=0,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
+     if(print_verbose)then
+        write(iout_iter,200) (stp(i),i=0,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
 
+  if (.not. end_iter) then
 ! Update solution
 !$omp parallel do schedule(dynamic,1) private(i,ii)
      do ii=1,nobs_bins+2

src/gsi/turbl_tl.f90

@@ -58,9 +58,7 @@ subroutine turbl_tl(pges,tges,oges,u,v,prs,t,termu,termv,termt,jstart,jstop)
   real(r_kind),dimension(nsig_hlf):: dudz_bck,dvdz_bck,dodz_bck,ri_bck,rf_bck
   real(r_kind),dimension(nsig_hlf):: rdudz_bck,rdvdz_bck,sm_bck,sh_bck,rdzl_bck,rdzi_bck
   real(r_kind),dimension(nsig_hlf):: u_tl,v_tl,o_tl,zl_tl,t_tl,rssq,rofbck,rshbck
-  real(r_kind),dimension(nsig_hlf):: km_bck,kh_bck,zi_bck
-  real(r_kind),dimension(nsig_hlf+1):: p_bck,zi_tl,p_tl
-  real(r_kind),dimension(2:nsig_hlf):: km_tl,kh_tl
+  real(r_kind),dimension(nsig_hlf+1):: km_bck,kh_bck,p_bck,zi_bck,km_tl,kh_tl,zi_tl,p_tl
   real(r_kind),dimension(2:nsig_hlf):: dzl_tl,dodz_tl,dudz_tl,dvdz_tl,ri_tl
   real(r_kind),dimension(2:nsig_hlf):: rf_tl,sh_tl,sm_tl,lmix_tl
   real(r_kind):: a1,a2,ax,bx,px,rpx,zx,ssq
@@ -71,6 +69,7 @@ subroutine turbl_tl(pges,tges,oges,u,v,prs,t,termu,termv,termt,jstart,jstop)
   integer(i_kind) i,j,k
   integer(i_kind),dimension(nsig):: lssq
 
+
   do j=jstart,jstop
      do i=1,lat2
 
@@ -84,6 +83,7 @@ subroutine turbl_tl(pges,tges,oges,u,v,prs,t,termu,termv,termt,jstart,jstop)
            zi_bck(k)=zi(i,j,k)
         end do
         p_bck(nsig_hlf+1) =pges(i,j,nsig_hlf+1)
+        zi_bck(nsig_hlf+1)=zi  (i,j,nsig_hlf+1)
 
         do k=1,nsig_hlf
            dodz_bck(k)=dodz(i,j,k)
@@ -96,6 +96,8 @@ subroutine turbl_tl(pges,tges,oges,u,v,prs,t,termu,termv,termt,jstart,jstop)
            km_bck(k)=km(i,j,k)
            kh_bck(k)=kh(i,j,k)
         end do
+        km_bck(nsig_hlf+1)=zero
+        kh_bck(nsig_hlf+1)=zero
 
         do k=2,nsig_hlf
            ssq=dudz_bck(k)**2+dvdz_bck(k)**2
@@ -247,6 +249,8 @@ subroutine turbl_tl(pges,tges,oges,u,v,prs,t,termu,termv,termt,jstart,jstop)
            end if
         end do
 
+        km_tl(1)=zero; kh_tl(1)=zero
+
 ! update perturbation tendencies
 
 
@@ -255,26 +259,23 @@ subroutine turbl_tl(pges,tges,oges,u,v,prs,t,termu,termv,termt,jstart,jstop)
            ax=t_bck(k)/o_bck(k)
            hrdzbk=half*rdzibk
            ardzbk=ax*hrdzbk
+           zx= dzi_tl(k)*rdzibk
+           kmaz_bck=(km_bck(k)+km_bck(k+1))*hrdzbk
+           khaz_bck=(kh_bck(k)+kh_bck(k+1))*ardzbk
+
+           termu(i,j,k)=termu(i,j,k)-zx*dudtm(i,j,k)
+           termv(i,j,k)=termv(i,j,k)-zx*dvdtm(i,j,k)
+           termt(i,j,k)=termt(i,j,k)-zx*dtdtm(i,j,k)
            kmaz_tl=zero
            khaz_tl=zero
-           kmaz_bck=km_bck(k)*hrdzbk
-           khaz_bck=kh_bck(k)*ardzbk
            if(k<nsig_hlf) then
               kmaz_tl =kmaz_tl+km_tl (k+1)*hrdzbk
               khaz_tl =khaz_tl+kh_tl (k+1)*ardzbk
-              kmaz_bck=kmaz_bck+km_bck(k)*hrdzbk
-              khaz_bck=khaz_bck+kh_bck(k)*ardzbk
            end if
            if(k>1) then
               kmaz_tl= kmaz_tl+km_tl (k)*hrdzbk
               khaz_tl= khaz_tl+kh_tl (k)*ardzbk
            end if
-
-           zx= dzi_tl(k)*rdzibk
-
-           termu(i,j,k)=termu(i,j,k)-zx*dudtm(i,j,k)
-           termv(i,j,k)=termv(i,j,k)-zx*dvdtm(i,j,k)
-           termt(i,j,k)=termt(i,j,k)-zx*dtdtm(i,j,k)
            if(k<nsig_hlf) then
               termu(i,j,k)=termu(i,j,k)+&
                  kmaz_bck*dudz_tl(k+1) +&

src/gsi/windht.f90

@@ -149,7 +149,7 @@ subroutine find_wind_height(cprov,csubprov,finalheight,kcount)
 
     character(len=8),intent(in)::cprov,csubprov
     real(r_kind),intent(out)::finalheight
-    integer,dimension(3)::kcount
+    integer,dimension(3),intent(inout)::kcount
 
     !local vars
     integer(i_kind)::i