Merge pull request NOAA-EMC#6 from ClaraDraper-NOAA/feature/2mDA_4atm…

…osupdate

Feature/2m da 4atmosupdate

src/enkf/controlvec.f90

@@ -47,13 +47,12 @@ module controlvec
                 mpi_integer,mpi_wtime,mpi_status,mpi_real8
 
 use gridio,    only: readgriddata, readgriddata_pnc, writegriddata, writegriddata_pnc, &
-                     writeincrement, writeincrement_pnc, readgriddata_2mDA, &
-                     writegriddata_2mDA, writeincrement_2mDA
+                     writeincrement, writeincrement_pnc 
 use gridinfo,  only: getgridinfo, gridinfo_cleanup,                    &
                      npts, vars3d_supported, vars2d_supported
 use params,    only: nlevs, nbackgrounds, fgfileprefixes, reducedgrid, &
                      nanals, pseudo_rh, use_qsatensmean, nlons, nlats,&
-                     nanals_per_iotask, ntasks_io, nanal1, nanal2, global_2mDA, &
+                     nanals_per_iotask, ntasks_io, nanal1, nanal2,&
                      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
@@ -132,7 +131,7 @@ subroutine init_controlvec()
       cvars3d(nc3d) = trim(adjustl(var))
       clevels(nc3d) = ilev + clevels(nc3d-1)
    else 
-      if (nproc .eq. 0) print *,'Error: only ', nlevs, ' and ', nlevs+1,' number of levels is supported in current version, got ',ilev
+      if (nproc .eq. 0) print *,'Error: controlvec only ', nlevs, ' and ', nlevs+1,' number of levels is supported in current version, got ',ilev
       call stop2(503)
    endif
 enddo
@@ -213,30 +212,22 @@ subroutine read_control()
 ! read in whole control vector on i/o procs - keep in memory 
 ! (needed in write_ensemble)
 allocate(grdin(npts,ncdim,nbackgrounds,nanals_per_iotask))
-if (.not. global_2mDA) allocate(qsat(npts,nlevs,nbackgrounds,nanals_per_iotask))
+q_ind = getindex(cvars3d, 'q')
+if (q_ind > 0)  allocate(qsat(npts,nlevs,nbackgrounds,nanals_per_iotask))
 if (paranc) then
-   if (global_2mDA) then
-      print *,'paranc not supported for global_2mDA'
-      call mpi_barrier(mpi_comm_world,ierr)
-      call mpi_finalize(ierr)
-   endif
    if (nproc == 0) t1 = mpi_wtime()
    call readgriddata_pnc(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,nbackgrounds, &
            fgfileprefixes,fgsfcfileprefixes,reducedgrid,grdin,qsat)
 end if
 if (nproc <= ntasks_io-1) then
    if (.not. paranc) then
       if (nproc == 0) t1 = mpi_wtime()
-      if (global_2mDA) then
-         call readgriddata_2mDA(nanal1(nproc),nanal2(nproc),cvars2d,nc2d,ncdim,nbackgrounds, &
-              fgsfcfileprefixes,reducedgrid,grdin)
-      else
-         call readgriddata(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,nbackgrounds, &
-              fgfileprefixes,fgsfcfileprefixes,reducedgrid,grdin,qsat)
-      endif
+      call readgriddata(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,nbackgrounds, &
+           fgfileprefixes,fgsfcfileprefixes,reducedgrid,grdin,qsat)
    end if
    !print *,'min/max qsat',nanal,'=',minval(qsat),maxval(qsat)
-   if (use_qsatensmean .and. .not. global_2mDA) then
+   q_ind = getindex(cvars3d, 'q')
+   if (use_qsatensmean .and. q_ind>0 ) then ! flag for if q is in control vector
        allocate(qsatmean(npts,nlevs,nbackgrounds))
        allocate(qsat_tmp(npts))
        ! compute ensemble mean qsat
@@ -268,8 +259,6 @@ subroutine read_control()
    !   print *,'min/max qsatmean proc',nproc,'=',&
    !            minval(qsatmean(:,:,nbackgrounds/2+1)),maxval(qsatmean(:,:,nbackgrounds/2+1))
    !endif
-   if (.not. global_2mDA) then
-   q_ind = getindex(cvars3d, 'q')
    if (pseudo_rh .and. q_ind > 0) then
       if (use_qsatensmean) then
          do ne=1,nanals_per_iotask
@@ -289,7 +278,6 @@ subroutine read_control()
          enddo
       endif
    end if
-   end if
 
 endif
 
@@ -345,7 +333,6 @@ subroutine write_control(no_inflate_flag)
 100 format('ens. mean anal. increment min/max  ',a,2x,g19.12,2x,g19.12)
    deallocate(grdin_mean_tmp)
 
-   if (.not. global_2mDA) then
    q_ind = getindex(cvars3d, 'q')
    if (pseudo_rh .and. q_ind > 0) then
       if (use_qsatensmean) then
@@ -374,35 +361,20 @@ subroutine write_control(no_inflate_flag)
          enddo
       endif
    end if
+
    if (.not. paranc) then
       if (write_fv3_incr) then
-         if (global_2mDA) then
-             call writeincrement_2mDA(nanal1(nproc),nanal2(nproc),cvars2d,nc2d,ncdim,grdin,no_inflate_flag)
-         else
              call writeincrement(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag)
-         endif
       else
-         if (global_2mDA) then
-             call writegriddata_2mDA(nanal1(nproc),nanal2(nproc),cvars2d,nc2d,ncdim,grdin,no_inflate_flag)
-         else
              call writegriddata(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag)
-         endif
       end if
       if (nproc == 0) then
         if (write_ensmean) then
            ! also write out ens mean on root task.
            if (write_fv3_incr) then
-              if (global_2mDA) then
-                 call writeincrement_2mDA(0,0,cvars2d,nc2d,ncdim,grdin_mean,no_inflate_flag)
-              else
                  call writeincrement(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag)
-              endif
            else
-              if (global_2mDA) then
-                 call writegriddata_2mDA(0,0,cvars2d,nc2d,ncdim,grdin_mean,no_inflate_flag)
-              else
                  call writegriddata(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag)
-              endif
            end if
         endif
         deallocate(grdin_mean)
@@ -414,11 +386,6 @@ subroutine write_control(no_inflate_flag)
 end if ! io task
 
 if (paranc) then
-   if (global_2mDA) then
-      print *,'paranc not supported for global_2mDA'
-      call mpi_barrier(mpi_comm_world,ierr)
-      call mpi_finalize(ierr)
-   endif
    if (write_fv3_incr) then
       call writeincrement_pnc(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag)
    else
@@ -454,3 +421,4 @@ subroutine controlvec_cleanup()
 end subroutine controlvec_cleanup
 
 end module controlvec
+

src/enkf/gridinfo_gfs.f90

@@ -45,7 +45,7 @@ module gridinfo
 
 use mpisetup, only: nproc, mpi_integer, mpi_real4
 use mpimod, only: mpi_comm_world
-use params, only: datapath,nlevs,nlons,nlats,use_gfs_nemsio,use_gfs_ncio,fgfileprefixes, global_2mDA
+use params, only: datapath,nlevs,nlons,nlats,use_gfs_nemsio,use_gfs_ncio,fgfileprefixes
 use kinds, only: r_kind, i_kind, r_double, r_single
 use constants, only: one,zero,pi,cp,rd,grav,rearth,max_varname_length
 use specmod, only: sptezv_s, sptez_s, init_spec_vars, isinitialized, asin_gaulats, &
@@ -138,11 +138,7 @@ subroutine getgridinfo(fileprefix, reducedgrid)
      dset = open_dataset(filename)
      londim = get_dim(dset,'grid_xt'); nlonsin = londim%len
      latdim = get_dim(dset,'grid_yt'); nlatsin = latdim%len
-     if (.not. global_2mDA) then 
-        levdim = get_dim(dset,'pfull');   nlevsin = levdim%len
-     else
-        nlevsin = 1
-     endif
+     levdim = get_dim(dset,'pfull');   nlevsin = levdim%len
      idvc = 2; ntrunc = nlatsin-2
      if (nlons /= nlonsin .or. nlats /= nlatsin .or. nlevs /= nlevsin) then
        print *,'incorrect dims in netcdf file'
@@ -233,13 +229,8 @@ subroutine getgridinfo(fileprefix, reducedgrid)
       endif
       ! convert to 1d array, units to millibars, flip so lats go N to S.
       spressmn = 0.01_r_kind*reshape(values_2d,(/nlons*nlats/))
-      if (.not. global_2mDA) then
-         call read_attribute(dset, 'ak', ak)
-         call read_attribute(dset, 'bk', bk)
-      else
-         allocate(ak(nlevs+1),bk(nlevs+1))
-         ak=0; bk=1
-      endif
+      call read_attribute(dset, 'ak', ak)
+      call read_attribute(dset, 'bk', bk)
       call close_dataset(dset)
       ! pressure at interfaces
       do k=1,nlevs+1
@@ -311,15 +302,11 @@ subroutine getgridinfo(fileprefix, reducedgrid)
          enddo
       enddo
    endif
-   if (global_2mDA) then
-      presslmn(:,nlevs) = spressmn
-   else
-      do k=1,nlevs
-        ! layer pressure from Phillips vertical interpolation.
-        presslmn(:,k) = ((pressimn(:,k)**kap1-pressimn(:,k+1)**kap1)/&
+   do k=1,nlevs
+      ! layer pressure from Phillips vertical interpolation.
+      presslmn(:,k) = ((pressimn(:,k)**kap1-pressimn(:,k+1)**kap1)/&
                          (kap1*(pressimn(:,k)-pressimn(:,k+1))))**kapr
-      end do
-   endif
+   end do
    print *,'ensemble mean first guess surface pressure:'
    print *,minval(spressmn),maxval(spressmn)
    !do k=1,nlevs