Merge pull request #281 from jswhit2/feature/enkf-64bit

GitHub Issue NOAA-EMC/GSI#277: Merge recent updates and bug fixes in EnKF

src/enkf/controlvec.f90

@@ -53,7 +53,7 @@ module controlvec
 use params,    only: nlevs, nbackgrounds, fgfileprefixes, reducedgrid, &
                      nanals, pseudo_rh, use_qsatensmean, nlons, nlats,&
                      nanals_per_iotask, ntasks_io, nanal1, nanal2, &
-                     fgsfcfileprefixes, paranc, write_fv3_incr
+                     fgsfcfileprefixes, paranc, write_fv3_incr, write_ensmean
 use kinds,     only: r_kind, i_kind, r_double, r_single
 use mpeu_util, only: gettablesize, gettable, getindex
 use constants, only: max_varname_length
@@ -291,13 +291,14 @@ subroutine write_control(no_inflate_flag)
 real(r_double)  :: t1,t2
 integer(i_kind) :: nb, nvar, ne
 integer(i_kind) :: q_ind, ierr
-real(r_single), allocatable, dimension(:,:) :: grdin_mean, grdin_mean_tmp
+real(r_single), allocatable, dimension(:,:) :: grdin_mean_tmp
+real(r_single), allocatable, dimension(:,:,:,:) :: grdin_mean
 
 if (nproc <= ntasks_io-1) then
 
    allocate(grdin_mean_tmp(npts,ncdim))
    if (nproc == 0) then
-     allocate(grdin_mean(npts,ncdim))
+     allocate(grdin_mean(npts,ncdim,nbackgrounds,1))
      grdin_mean = 0_r_single
      t1 = mpi_wtime()
    endif
@@ -311,27 +312,24 @@ subroutine write_control(no_inflate_flag)
       do ne=1,nanals_per_iotask
          call mpi_reduce(grdin(:,:,nb,ne), grdin_mean_tmp, npts*ncdim, mpi_real4,&
                          mpi_sum,0,mpi_comm_io,ierr)
-         if (nproc == 0) grdin_mean = grdin_mean + grdin_mean_tmp
+         if (nproc == 0) grdin_mean(:,:,nb,1) = grdin_mean(:,:,nb,1) + grdin_mean_tmp
       enddo
       ! print out ens mean increment info
       if (nproc == 0) then
-         grdin_mean = grdin_mean/real(nanals)
+         grdin_mean(:,:,nb,1) = grdin_mean(:,:,nb,1)/real(nanals)
          do nvar=1,nc3d
             write(6,100) trim(cvars3d(nvar)),   &
-                minval(grdin_mean(:,clevels(nvar-1)+1:clevels(nvar))),     &
-                maxval(grdin_mean(:,clevels(nvar-1)+1:clevels(nvar)))
+                minval(grdin_mean(:,clevels(nvar-1)+1:clevels(nvar),nb,1)),     &
+                maxval(grdin_mean(:,clevels(nvar-1)+1:clevels(nvar),nb,1))
          enddo
          do nvar=1,nc2d
             write(6,100) trim(cvars2d(nvar)),   &
-                minval(grdin_mean(:,clevels(nc3d) + nvar)),                &
-                maxval(grdin_mean(:,clevels(nc3d) + nvar))
+                minval(grdin_mean(:,clevels(nc3d) + nvar,nb,1)),                &
+                maxval(grdin_mean(:,clevels(nc3d) + nvar,nb,1))
          enddo
       endif
    enddo
 100 format('ens. mean anal. increment min/max  ',a,2x,g19.12,2x,g19.12)
-   if (nproc == 0) then
-      deallocate(grdin_mean)
-   endif
    deallocate(grdin_mean_tmp)
 
    q_ind = getindex(cvars3d, 'q')
@@ -353,6 +351,14 @@ subroutine write_control(no_inflate_flag)
          enddo
          enddo
       endif
+      if (nproc == 0 .and. write_ensmean) then
+         ! write_ensmean implies use_qsatensmean
+         do nb=1,nbackgrounds
+            ! re-scale normalized spfh with sat. sphf of ensmean first guess
+            grdin_mean(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,1) = &
+            grdin_mean(:,(q_ind-1)*nlevs+1:q_ind*nlevs,nb,1)*qsatmean(:,:,nb)
+         enddo
+      endif
    end if
    if (.not. paranc) then
       if (write_fv3_incr) then
@@ -361,6 +367,15 @@ subroutine write_control(no_inflate_flag)
          call writegriddata(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag)
       end if
       if (nproc == 0) then
+        if (write_ensmean) then
+           ! also write out ens mean on root task.
+           if (write_fv3_incr) then
+              call writeincrement(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag)
+           else
+              call writegriddata(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag)
+           end if
+        endif
+        deallocate(grdin_mean)
         t2 = mpi_wtime()
         print *,'time in write_control on root',t2-t1,'secs'
       endif 
@@ -375,6 +390,15 @@ subroutine write_control(no_inflate_flag)
       call writegriddata_pnc(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag)
    end if
    if (nproc == 0) then
+     ! also write out ens mean on root task
+     if (write_ensmean) then ! FIXME use parallel IO to write ensmean
+        if (write_fv3_incr) then
+           call writeincrement(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag)
+        else
+           call writegriddata(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag)
+        end if
+     endif
+     deallocate(grdin_mean)
      t2 = mpi_wtime()
      print *,'time in write_control on root',t2-t1,'secs'
    endif 

src/enkf/enkf_obsmod.f90

@@ -105,7 +105,7 @@ module enkf_obsmod
 use kinds, only : r_kind, r_double, i_kind, r_single
 use constants, only: zero, one, deg2rad, rad2deg, rd, cp, pi
 use params, only: & 
-      datestring,datapath,sprd_tol,nanals,saterrfact, &
+      letkf_flag,nobsl_max,datestring,datapath,sprd_tol,nanals,saterrfact, &
       lnsigcutoffnh, lnsigcutoffsh, lnsigcutofftr, corrlengthnh,&
       corrlengthtr, corrlengthsh, obtimelnh, obtimeltr, obtimelsh,&
       lnsigcutoffsatnh, lnsigcutoffsatsh, lnsigcutoffsattr,&
@@ -204,6 +204,13 @@ subroutine readobs()
 if (nproc == 0) then
  print *,'max time in mpireadobs  = ',tdiffmax
  print *,'total number of obs ',nobstot
+ print *,'min/max obtime ',minval(obtime),maxval(obtime)
+endif
+! if nobsl_max set for LETKF, and the total number of obs < nobsl_max,
+! reset nobsl_max to -1
+if (letkf_flag .and. nobsl_max > 0 .and. nobstot < nobsl_max) then
+   if (nproc == 0) print *,'resetting nobsl_max to -1'
+   nobsl_max=-1
 endif
 allocate(obfit_prior(nobstot))
 ! screen out some obs by setting ob error to a very large number

src/enkf/gridio_gfs.f90

@@ -333,14 +333,19 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, &
 
   if (iope==0) then
      do k=1,nlevs
-        krev = nlevs-k+1
-        ! layer pressure from phillips vertical interolation
+        ! pressure at bottom of layer interface (for gps jacobian, see prsltmp in setupbend.f90)
+        if (prse_ind > 0) then
+           ug(:) = pressi(:,k)
+           call copytogrdin(ug,pslg(:,k))
+           ! Jacobian for gps in pressure is saved in different units in GSI; need to
+           ! multiply pressure by 0.1
+           grdin(:,levels(prse_ind-1)+k,nb,ne) = 0.1*pslg(:,k)
+        endif
+        ! layer pressure from phillips vertical interolation (used for qsat
+        ! calculation)
         ug(:) = ((pressi(:,k)**kap1-pressi(:,k+1)**kap1)/&
                 (kap1*(pressi(:,k)-pressi(:,k+1))))**kapr
         call copytogrdin(ug,pslg(:,k))
-        ! Jacobian for gps in pressure is saved in different units in GSI; need to
-        ! multiply pressure by 0.1
-        if (prse_ind > 0)     grdin(:,levels(prse_ind-1)+k,nb,ne) = 0.1*pslg(:,k)
      end do
      if (pseudo_rh) then
         call genqsat1(q,qsat(:,:,nb,ne),pslg,tv,ice,npts,nlevs)
@@ -952,15 +957,19 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,
 
   ! compute saturation q.
   do k=1,nlevs
-    ! layer pressure from phillips vertical interolation
+    ! pressure at bottom of layer interface (for gps jacobian, see prsltmp in setupbend.f90)
+    if (prse_ind > 0) then
+       ug(:) = pressi(:,k)
+       call copytogrdin(ug,pslg(:,k))
+       ! Jacobian for gps in pressure is saved in different units in GSI; need to
+       ! multiply pressure by 0.1
+       grdin(:,levels(prse_ind-1)+k,nb,ne) = 0.1*pslg(:,k)
+    endif
+    ! layer pressure from phillips vertical interolation (used for qsat
+    ! calculation)
     ug(:) = ((pressi(:,k)**kap1-pressi(:,k+1)**kap1)/&
             (kap1*(pressi(:,k)-pressi(:,k+1))))**kapr
-
     call copytogrdin(ug,pslg(:,k))
-    ! Jacobian for gps in pressure is saved in different units in GSI; need to
-    ! multiply pressure by 0.1
-    if (prse_ind > 0)     grdin(:,levels(prse_ind-1)+k,nb,ne) = 0.1*pslg(:,k)
-
   end do
   if (pseudo_rh) then
      call genqsat1(q,qsat(:,:,nb,ne),pslg,tv,ice,npts,nlevs)
@@ -1115,7 +1124,9 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_
   ! need to distribute grdin to all PEs in this subcommunicator
   ! bring all the subdomains back to the main PE
   call mpi_barrier(iocomms(mem_pe(nproc)), iret)
-  call mpi_bcast(grdin,npts*ndim*nbackgrounds, mpi_real4, 0, iocomms(mem_pe(nproc)), iret)
+  do nb=1,nbackgrounds
+     call mpi_bcast(grdin(1,1,nb,1),npts*ndim, mpi_real4, 0, iocomms(mem_pe(nproc)), iret)
+  enddo
 
   ! loop through times and do the read
   ne = 1
@@ -1838,7 +1849,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n
                           write_attribute, quantize_data, has_var, has_attr
   use constants, only: grav
   use params, only: nbackgrounds,anlfileprefixes,fgfileprefixes,reducedgrid,&
-                    nccompress
+                    nccompress,write_ensmean
   implicit none
 
   integer, intent(in) :: nanal1,nanal2
@@ -1906,12 +1917,17 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n
   write(charnanal,'(i3.3)') nanal
   backgroundloop: do nb=1,nbackgrounds
 
-  if(no_inflate_flag) then
-    filenameout = trim(adjustl(datapath))//trim(adjustl(anlfileprefixes(nb)))//"nimem"//charnanal
+  if (nanal == 0 .and. write_ensmean) then
+     filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"ensmean"
+     filenameout = trim(adjustl(datapath))//trim(adjustl(anlfileprefixes(nb)))//"ensmean"
   else
-    filenameout = trim(adjustl(datapath))//trim(adjustl(anlfileprefixes(nb)))//"mem"//charnanal
-  end if
-  filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"mem"//charnanal
+     if(no_inflate_flag) then
+       filenameout = trim(adjustl(datapath))//trim(adjustl(anlfileprefixes(nb)))//"nimem"//charnanal
+     else
+       filenameout = trim(adjustl(datapath))//trim(adjustl(anlfileprefixes(nb)))//"mem"//charnanal
+     end if
+     filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"mem"//charnanal
+  endif
 
   if (use_gfs_nemsio) then
      clip = tiny(vg(1))
@@ -3289,7 +3305,7 @@ end subroutine writegriddata
  subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_flag)
   use netcdf
   use params, only: nbackgrounds,incfileprefixes,fgfileprefixes,reducedgrid,&
-                    datestring,nhr_anal
+                    datestring,nhr_anal,write_ensmean
   use constants, only: grav
   use mpi
   use module_fv3gfs_ncio, only: Dataset, Variable, Dimension, open_dataset,&
@@ -3354,12 +3370,17 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,
   write(charnanal,'(i3.3)') nanal
   backgroundloop: do nb=1,nbackgrounds
 
-  if(no_inflate_flag) then
-    filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"nimem"//charnanal
+  if (nanal == 0 .and. write_ensmean) then
+     filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"ensmean"
+     filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"ensmean"
   else
-    filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"mem"//charnanal
-  end if
-  filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"mem"//charnanal
+     if(no_inflate_flag) then
+       filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"nimem"//charnanal
+     else
+       filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"mem"//charnanal
+     end if
+     filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"mem"//charnanal
+  endif
 
   ! create the output netCDF increment file
   call nccheck_incr(nf90_create(path=trim(filenameout), cmode=nf90_netcdf4, ncid=ncid_out))
@@ -3773,7 +3794,9 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate
   ! need to distribute grdin to all PEs in this subcommunicator
   ! bring all the subdomains back to the main PE
   call mpi_barrier(iocomms(mem_pe(nproc)), iret)
-  call mpi_bcast(grdin,npts*ndim*nbackgrounds, mpi_real4, 0, iocomms(mem_pe(nproc)), iret)
+  do nb=1,nbackgrounds
+     call mpi_bcast(grdin(1,1,nb,1),npts*ndim, mpi_real4, 0, iocomms(mem_pe(nproc)), iret)
+  enddo
 
   ! loop through times and do the read
   ne = 1