Merge branch 'NOAA-EMC:master' into feature/rm_sfcanl

src/GSD/gsdcloud/cloudCover_Surface.f90

@@ -1,4 +1,4 @@
-SUBROUTINE cloudCover_Surface(mype,nlat,nlon,nsig,r_radius,thunderRadius,&
+SUBROUTINE cloudCover_Surface(mype,nlat,nlon,nsig,thunderRadius,&
                         cld_bld_hgt,t_bk,p_bk,q,h_bk,zh,  &
                         mxst_p,NVARCLD_P,numsao,OI,OJ,OCLD,OWX,Oelvtn,Odist,&
                         cld_cover_3d,cld_type_3d,wthr_type,pcp_type_3d,     &
@@ -23,7 +23,6 @@ SUBROUTINE cloudCover_Surface(mype,nlat,nlon,nsig,r_radius,thunderRadius,&
 !     nlon        - no. of lons on subdomain (buffer points on ends)
 !     nlat        - no. of lats on subdomain (buffer points on ends)
 !     nsig        - no. of levels
-!     r_radius    - influence radius of the cloud observation
 !     thunderRadius -
 !     cld_bld_hgt - Height below which cloud building is done 
 !
@@ -71,7 +70,6 @@ SUBROUTINE cloudCover_Surface(mype,nlat,nlon,nsig,r_radius,thunderRadius,&
   implicit none
 
   integer(i_kind),intent(in) :: mype
-  REAL(r_single), intent(in) :: r_radius
   integer(i_kind),intent(in) :: nlat,nlon,nsig
   real(r_single), intent(in) :: thunderRadius
   real(r_kind),   intent(in) :: cld_bld_hgt

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.f90

@@ -123,7 +123,7 @@ module enkf
                   obtype, oberrvarmean, numobspersat, deltapredx, biaspreds,&
                   oberrvar_orig, probgrosserr, prpgerr,&
                   corrlengthsq,lnsigl,obtimel,obloclat,obloclon,obpress,stattype,&
-                  anal_ob,anal_ob_post
+                  anal_ob,anal_ob_post,assimltd_flag
 use constants, only: pi, one, zero
 use params, only: sprd_tol, paoverpb_thresh, datapath, nanals,&
                   iassim_order,sortinc,deterministic,numiter,nlevs,&
@@ -153,7 +153,8 @@ subroutine enkf_update()
 
 ! local variables.
 integer(i_kind) nob,nob1,nob2,nob3,npob,nf,nf2,ii,nobx,nskip,&
-                niter,i,nrej,npt,nuse,ncount,ncount_check,nb,np
+                niter,i,nrej,npt,nuse,ncount,ncount_check,nb,np,&
+                nuseconvoz,nusesat,nobs_convoz
 integer(i_kind) indxens1(nanals),indxens2(nanals)
 integer(i_kind) indxens1_modens(nanals*neigv),indxens2_modens(nanals*neigv)
 real(r_single) hxpost(nanals),hxprior(nanals),hxinc(nanals),&
@@ -757,17 +758,31 @@ subroutine enkf_update()
   tend = mpi_wtime()
   if (nproc .eq. 0) then
       write(6,8003) niter,'timing on proc',nproc,' = ',tend-tbegin,t2,t3,t4,t5,t6,nrej
+      allocate(assimltd_flag(nobstot))
+      assimltd_flag = 99999
       if (iassim_order == 2) then
           ncount_check = ncount
       else
           ncount_check = nobstot
       endif
       nuse = 0; covl_fact = 0.
+      nuseconvoz=0; nusesat = 0
+      nobs_convoz = nobs_conv + nobs_oz
       do nob1=1,ncount_check
          nob = indxassim(nob1)
          if (iskip(nob) .ne. 1) then
             covl_fact = covl_fact + sqrt(corrlengthsq(nob)/corrlengthsq_orig(nob))
             nuse = nuse + 1
+            assimltd_flag(nob) = 1
+            if (nob .le. nobs_convoz) then
+               nuseconvoz = nuseconvoz +1
+            else if (nob .gt. nobs_convoz) then
+               nusesat = nusesat + 1
+            else
+               print *,'nob ', nob ,' falling through'
+            endif
+         else
+            assimltd_flag(nob) = 0
          endif
       enddo
       nskip = nobstot-nuse
@@ -776,6 +791,8 @@ subroutine enkf_update()
       if (covl_fact < 0.99) print *,'mean covl_fact = ',covl_fact
       if (nskip > 0) print *,nskip,' out of',nobstot,'obs skipped,',nuse,' used'
       if (nsame > 0) print *,nsame,' out of', nobstot-nskip,' same lat/long'
+      if (nuseconvoz > 0 ) print *,nuseconvoz,' out of',nobs_conv + nobs_oz ,'convobs used'
+      if (nusesat > 0 ) print *,nusesat ,' out of',nobs_sat ,'satobs used'
       if (nrej >  0) print *,nrej,' obs rejected by varqc'
   endif
   8003  format(i2,1x,a14,1x,i5,1x,a3,6(f7.2,1x),i4)

src/enkf/enkf_obs_sensitivity.f90

@@ -45,7 +45,7 @@ module enkf_obs_sensitivity
                        obloclon,obpress,indxsat,oberrvar,stattype,obtime,ob, &
                        ensmean_ob,ensmean_obnobc,obsprd_prior,obfit_prior, &
                        oberrvar_orig,biaspreds,anal_ob_post,nobstot,lnsigl, &
-                       corrlengthsq,obtimel,oblnp,obloc
+                       corrlengthsq,obtimel,oblnp,obloc ,assimltd_flag
 use convinfo, only: convinfo_read,init_convinfo
 use ozinfo, only: ozinfo_read,init_oz
 use radinfo, only: radinfo_read,jpch_rad,nusis,nuchan,npred
@@ -91,6 +91,7 @@ module enkf_obs_sensitivity
   integer(i_kind) :: stattype           ! Observation type
   character(len=20) :: obtype           ! Observation element / Satellite name
   integer(i_kind) :: indxsat            ! Satellite index (channel)
+  integer(i_kind) :: assimltd_flag      ! is assimilated? flag
   real(r_single)  :: osense_kin         ! Observation sensitivity (kinetic energy) [J/kg]
   real(r_single)  :: osense_dry         ! Observation sensitivity (Dry total energy) [J/kg]
   real(r_single)  :: osense_moist       ! Observation sensitivity (Moist total energy) [J/kg]
@@ -215,6 +216,7 @@ subroutine read_ob_sens
   allocate(stattype(nobstot))
   allocate(obtype(nobstot))
   allocate(indxsat(nobs_sat))
+  allocate(assimltd_flag(nobstot))
   allocate(biaspreds(npred+1,nobs_sat))
   allocate(tmpanal_ob(nanals))
   allocate(tmpbiaspreds(npred+1))
@@ -235,6 +237,7 @@ subroutine read_ob_sens
      oberrvar_orig(nob) = real(indata%oberrvar_orig,r_kind)
      stattype(nob) = indata%stattype
      obtype(nob) = indata%obtype
+     assimltd_flag(nob) = indata%assimltd_flag
      if(nproc == 0) anal_ob_post(1:nanals,nob) = real(tmpanal_ob(1:nanals),r_kind)
   end do
   ! Read loop over satellite radiance observations
@@ -256,6 +259,7 @@ subroutine read_ob_sens
      stattype(nob) = indata%stattype
      obtype(nob) = indata%obtype
      indxsat(nn) = indata%indxsat
+     assimltd_flag(nob) = indata%assimltd_flag
      if(nproc == 0) anal_ob_post(1:nanals,nob) = real(tmpanal_ob(1:nanals),r_kind)
      biaspreds(1:npred+1,nn) = real(tmpbiaspreds(1:npred+1),r_kind)
   end do
@@ -430,6 +434,7 @@ subroutine print_ob_sens
         outdata%stattype = stattype(nob)
         outdata%obtype = obtype(nob)
         outdata%indxsat = 0
+        outdata%assimltd_flag = assimltd_flag(nob)
         if(efsoi_flag) then
            outdata%osense_kin = real(obsense_kin(nob),r_single)
            outdata%osense_dry = real(obsense_dry(nob),r_single)
@@ -500,6 +505,7 @@ subroutine print_ob_sens
         outdata%stattype = stattype(nob)
         outdata%obtype = obtype(nob)
         outdata%indxsat = nchan
+        outdata%assimltd_flag = assimltd_flag(nob)
         tmpbiaspreds(1:npred+1) = real(biaspreds(1:npred+1,nn),r_single)
         if(efsoi_flag) then
            outdata%osense_kin = real(obsense_kin(nob),r_single)

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,&
@@ -152,6 +152,9 @@ module enkf_obsmod
 
 ! ob-space posterior ensemble, needed for EFSOI
 real(r_single),public,allocatable, dimension(:,:) :: anal_ob_post   ! Fortran pointer
+! is the observation assimilated? logical would be preferable, but that confuses
+! Python
+integer(i_kind),public,allocatable, dimension(:) :: assimltd_flag ! Fortran pointer
 
 contains
 
@@ -204,6 +207,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
@@ -466,6 +476,7 @@ subroutine obsmod_cleanup()
 if (allocated(prpgerr)) deallocate(prpgerr)
 if (allocated(diagused)) deallocate(diagused)
 if (allocated(anal_ob_post)) deallocate(anal_ob_post)
+if (allocated(assimltd_flag)) deallocate(assimltd_flag)
 ! free shared memory segement, fortran pointer to that memory.
 nullify(anal_ob)
 call MPI_Barrier(mpi_comm_world,ierr)