Issue #2: Merge branch 'master' at 85f38c0 into release/gfsda.v16.0.0

cmake/Modules/setIntelFlags.cmake

@@ -64,7 +64,7 @@ function(set_GSI_ENKF_Intel)
 endfunction(set_GSI_ENKF_Intel)
 
 function (set_GSI_ENKF_Debug_Intel)
-    set(GSI_Fortran_FLAGS "-DPOUND_FOR_STRINGIFY -O0 -fp-model source -convert big_endian -assume byterecl -implicitnone -mcmodel medium -shared-intel -g -traceback -debug -ftrapuv -check all,noarg_temp_created -fp-stack-check -fstack-protector -warn all,nointerfaces -convert big_endian -implicitnone -D_REAL8_ ${OpenMP_Fortran_FLAGS} ${MPI_Fortran_COMPILE_FLAGS} ${GSDCLOUDOPT}" CACHE INTERNAL "")
+    set(GSI_Fortran_FLAGS "-DPOUND_FOR_STRINGIFY -O0 -fp-model source -convert big_endian -assume byterecl -implicitnone -g -traceback -debug -ftrapuv -check all,noarg_temp_created -fp-stack-check -fstack-protector -warn all,nointerfaces -convert big_endian -implicitnone -D_REAL8_ ${OpenMP_Fortran_FLAGS} ${MPI_Fortran_COMPILE_FLAGS} ${GSDCLOUDOPT}" CACHE INTERNAL "")
     set(ENKF_Fortran_FLAGS "-O0 ${HOST_FLAG} -warn all -implicitnone -traceback -g -debug all -check all,noarg_temp_created -fp-model strict -convert big_endian -assume byterecl -D_REAL8_  ${OpenMP_Fortran_FLAGS} ${GSDCLOUDOPT}" CACHE INTERNAL "")
     set(GSDCLOUD_Fortran_FLAGS "-DPOUND_FOR_STRINGIFY -O3 -convert big_endian" CACHE INTERNAL "")
     #Common debug flags

modulefiles/modulefile.ProdGSI.jet

@@ -11,21 +11,19 @@
 #module-whatis " GDAS_ENKF whatis description"
 #
 #
-module load newdefaults
-module load intel/2016.2.181 
-module load impi/5.1.1.109
+
+module load intel/18.0.5.274
+module load impi
 #
-module load szip/2.1   
+module load szip/2.1
 module load hdf5/1.8.9
-module load netcdf4/4.2.1.1  
 
-export COREPATH=/lfs3/projects/hfv3gfs/nwprod/lib
-export COREPATH=/mnt/lfs3/projects/hfv3gfs/gwv/ljtjet/lib
+export COREPATH=/lfs4/HFIP/hfv3gfs/gwv/l0530/lib
 export NCEPLIBS=$COREPATH
 module use $NCEPLIBS/modulefiles
 
-export FFLAGS="-openmp"
-export LDFLAGS="-openmp"
+export FFLAGS="-qopenmp"
+export LDFLAGS="-qopenmp"
 module load bacio-intel-sandybridge/2.0.2
 module load crtm-intel-sandybridge/2.2.5
 module load bufr-intel-sandybridge/11.1.0
@@ -35,10 +33,7 @@ module load w3emc-intel-sandybridge/2.2.0
 module load w3nco-intel-sandybridge/2.0.6
 module load sigio-intel-sandybridge/2.0.1
 module load sfcio-intel-sandybridge/1.0.0
+module load ip-intel-sandybridge/3.0.0
 
-module use /home/George.Vandenberghe/t1/l1025/lib/modulefiles
-module load hdf5
-module load netcdf
 
-export NETCDF=${NETCDF4}
-export NETCDF_DIR=${NETCDF4}
+module load netcdfp/4.7.4

modulefiles/modulefile.ProdGSI.wcoss_c

@@ -40,7 +40,7 @@ module unload intel/15.0.3.187
 module load intel/18.1.163
 
 # Loading pe environment
-module load craype-haswell
+module load craype-sandybridge
 module load cray-mpich/7.2.0
 
 module use /usrx/local/dev/modulefiles

regression/regression_var.sh

@@ -138,6 +138,8 @@ case $machine in
       export noscrub="/gpfs/hps3/emc/global/noscrub/$LOGNAME"
    elif [ -d /gpfs/hps3/emc/da/noscrub/$LOGNAME ]; then
       export noscrub="/gpfs/hps3/emc/da/noscrub/$LOGNAME"
+   elif [ -d /gpfs/hps3/emc/hwrf/noscrub/$LOGNAME ]; then
+       export noscrub="/gpfs/hps3/emc/hwrf/noscrub/$LOGNAME"
    fi
    if [[ "$cmaketest" = "false" ]]; then
      export basedir="/gpfs/hps3/emc/global/noscrub/$LOGNAME/svn/gsi"

src/enkf/gridio_wrf.f90

@@ -46,7 +46,8 @@ module gridio
   ! Define all public subroutines within this module
   private
   public :: readgriddata, readgriddata_pnc
-  public :: writegriddata, writegriddata_pnc, WRITEINCREMENT, WRITEINCREMENT_PNC
+  public :: writegriddata, writegriddata_pnc
+  public :: writeincrement, writeincrement_pnc
 
   !-------------------------------------------------------------------------
 
@@ -74,6 +75,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes,f
 
   end subroutine readgriddata
 
+
   !========================================================================
   ! readgriddata_arw.f90: read WRF-ARW state or control vector
   !-------------------------------------------------------------------------
@@ -1036,6 +1038,7 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,vargrid
 
   end subroutine writegriddata
 
+
   !======================================================================
   ! readwrfvar.f90: This subroutine reads a varname variable from WRF
   ! ARW or NMM netcdf file and returns the variable interpolated to 

src/gsi/get_gefs_for_regional.f90

@@ -22,6 +22,7 @@ subroutine get_gefs_for_regional
 !   2016-05-19  Carley/s.liu   - prevent the GSI from printing out erroneous error  
 !                               when using ensembles from different time
 !   2016-12-12  tong    - add code to get nemsio meta data, if use_gfs_nemsio=True
+!   2020-07-01  Bi   - add code to get netCDF data, if use_gfs_ncio=.true.
 !
 !   input argument list:
 !
@@ -33,7 +34,8 @@ subroutine get_gefs_for_regional
 !
 !$$$ end documentation block
 
-  use gridmod, only: idsl5,regional,use_gfs_nemsio
+  use gridmod, only: idsl5,regional,use_gfs_nemsio,use_gfs_ncio,&
+                     ncepgfs_head,ncepgfs_headv
   use gridmod, only: nlon,nlat,lat2,lon2,nsig,rotate_wind_ll2xy
   use hybrid_ensemble_parameters, only: region_lat_ens,region_lon_ens
   use hybrid_ensemble_parameters, only: en_perts,ps_bar,nelen
@@ -74,6 +76,9 @@ subroutine get_gefs_for_regional
   use nemsio_module, only: nemsio_init,nemsio_open,nemsio_close
   use ncepnems_io, only: error_msg
   use nemsio_module, only: nemsio_gfile,nemsio_getfilehead
+  use module_fv3gfs_ncio, only: Dimension, Dataset, open_dataset, get_dim, &
+                                read_vardata, get_idate_from_time_units,&
+                                read_attribute, close_dataset
   use get_wrf_mass_ensperts_mod, only: get_wrf_mass_ensperts_class
   use gsi_io, only: verbose
   use obsmod, only: l_wcp_cwm
@@ -109,9 +114,12 @@ subroutine get_gefs_for_regional
   character(len=*),parameter::myname='get_gefs_for_regional'
   real(r_kind) bar_norm,sig_norm,kapr,kap1,trk
   integer(i_kind) iret,i,j,k,k2,n,mm1,iderivative
+  integer(i_kind) mype_out
   integer(i_kind) ic2,ic3,it
   integer(i_kind) ku,kv,kt,kq,koz,kcw,kz,kps
   character(255) filename,filelists(ntlevs_ens)
+  character(6) sfilename
+
   logical ice
   integer(sigio_intkind):: lunges = 11
   type(sigio_head):: sighead
@@ -121,10 +129,11 @@ subroutine get_gefs_for_regional
   logical,allocatable :: vector(:)
   real(r_kind),parameter::  zero_001=0.001_r_kind
   real(r_kind),allocatable,dimension(:) :: xspli,yspli,xsplo,ysplo
-  integer(i_kind) iyr,ihourg
+  integer(i_kind) iyr,ihourg,kr
   integer(i_kind),dimension(4):: idate4
   integer(i_kind),dimension(8) :: ida,jda 
   integer(i_kind),dimension(5) :: iadate_gfs
+  integer(i_kind),dimension(6):: idate6
   real(r_kind) hourg
   real(r_kind),dimension(5):: fha
   integer(i_kind) istatus
@@ -136,6 +145,13 @@ subroutine get_gefs_for_regional
   integer(i_kind),dimension(7):: idate
   real(r_kind) :: fhour
   type(nemsio_gfile) :: gfile
+  type(Dataset) :: atmges      
+  type(ncepgfs_head):: gfshead
+  type(ncepgfs_headv):: gfsheadv
+  real(r_single),allocatable,dimension(:) :: aknc, bknc, fhour1
+
+
+  type(Dimension) :: ncdim
   integer(i_kind) :: nvcoord
   real(r_single),allocatable:: nems_vcoord(:,:,:)
   real(r_single),allocatable:: vcoord(:,:)
@@ -223,7 +239,7 @@ subroutine get_gefs_for_regional
 
   rewind (10) 
   read(10,'(a)',err=20,end=20)filename 
-  if(.not. use_gfs_nemsio)then
+  if ((.not. use_gfs_nemsio) .and. (.not. use_gfs_ncio))then
      open(lunges,file=trim(filename),form='unformatted')
      call sigio_srhead(lunges,sighead,iret)
      close(lunges)
@@ -264,7 +280,7 @@ subroutine get_gefs_for_regional
      idate4(2)= sighead%idate(2)
      idate4(3)= sighead%idate(3)
      idate4(4)= sighead%idate(4)
-  else
+  else if ( use_gfs_nemsio ) then
      call nemsio_init(iret=iret)
      if (iret /= 0) call error_msg(trim(my_name),trim(filename),' ','init',istop,iret)
 
@@ -336,12 +352,77 @@ subroutine get_gefs_for_regional
      if(mype == 0) then
         write(6,*) ' nemsio: fhour,idate=',fhour,idate
         write(6,*) ' iadate(y,m,d,hr,min)=',iadate
-        write(6,*) ' nemsio: jcap,levs=',njcap,levs
+        write(6,*) ' nemsio: jcap,levs=',levs
         write(6,*) ' nemsio: latb,lonb=',latb,lonb
         write(6,*) ' nemsio: idvc,nvcoord=',idvc,nvcoord
         write(6,*) ' nemsio: idsl=',idsl
      end if
 
+!  add netCDF header information: 
+     else ! use_gfs_ncio and get this information
+        write(sfilename,'("sfcf",i2.2)')nhr_assimilation
+        ! open the netCDF file
+        atmges = open_dataset(filename)
+        ! get dimension sizes
+        ncdim = get_dim(atmges, 'grid_xt'); gfshead%lonb = ncdim%len
+        ncdim = get_dim(atmges, 'grid_yt'); gfshead%latb = ncdim%len
+        ncdim = get_dim(atmges, 'pfull') ; gfshead%levs = ncdim%len
+        nsig_gfs = gfshead%levs
+        ! hard code jcap,idsl,idvc
+        gfshead%jcap = -9999
+        gfshead%idsl= 1
+        gfshead%idvc = 2
+
+        ! FV3GFS write component does not include JCAP, infer from DIMY-2
+        njcap=latb-2
+
+        nlat_gfs=gfshead%latb+2
+        nlon_gfs=gfshead%lonb
+        nsig_gfs=gfshead%levs
+
+        jcap_gfs=gfshead%latb-2
+
+        if (mype==mype_out) write(6,*)'GESINFO:  Read NCEP FV3GFS netCDF ', &
+           'format file, ',trim(filename)
+        ! hard code nvcoord to be 2
+        gfshead%nvcoord=2 ! ak and bk
+        if (allocated(gfsheadv%vcoord)) deallocate(gfsheadv%vcoord)
+        allocate(gfsheadv%vcoord(gfshead%levs+1,gfshead%nvcoord))
+        call read_attribute(atmges, 'ak', aknc)
+        call read_attribute(atmges, 'bk', bknc)
+        do k=1,gfshead%levs+1
+           kr = gfshead%levs+2-k
+           gfsheadv%vcoord(k,1) = aknc(kr)
+           gfsheadv%vcoord(k,2) = bknc(kr)
+        end do
+        deallocate(aknc,bknc) 
+        ! get time information
+        idate6 = get_idate_from_time_units(atmges)
+        gfshead%idate(1) = idate6(4)  !hour
+        gfshead%idate(2) = idate6(2)  !month
+        gfshead%idate(3) = idate6(3)  !day
+        gfshead%idate(4) = idate6(1)  !year
+        call read_vardata(atmges, 'time', fhour1) ! might need to change this to attribute later
+                                               ! depends on model changes from Jeff Whitaker
+        gfshead%fhour = fhour1(1)
+
+        call close_dataset(atmges)
+
+        if(mype == 0) then
+          write(6,*) ' netCDF:fhour,idate=',fhour1,idate6
+          write(6,*) ' netCDF:iadate(y,m,d,hr,min)=',iadate
+          write(6,*) ' netCDF: jcap,levs=',gfshead%levs
+          write(6,*) ' netCDF: latb,lonb=',gfshead%latb,gfshead%lonb
+          write(6,*) ' netCDF: nvcoord=',gfshead%nvcoord
+          write(6,*) ' netCDF: idvc,idsl=',gfshead%idvc,gfshead%idsl
+        endif
+
+        hourg = fhour1(1)
+        idate4(1) = idate6(4)
+        idate4(2) = idate6(2)
+        idate4(3) = idate6(3)
+        idate4(4) = idate6(1)
+
   end if
 
 ! Compute valid time from ensemble date and forecast length and compare to iadate, the analysis time
@@ -354,6 +435,7 @@ subroutine get_gefs_for_regional
         iyr=iyr+2000
      end if
   end if
+
   fha=zero ; ida=0; jda=0
   fha(2)=ihourg    ! relative time interval in hours
   ida(1)=iyr       ! year
@@ -371,30 +453,36 @@ subroutine get_gefs_for_regional
      iadate_gfs(4)=jda(5) ! hour
   endif
   iadate_gfs(5)=0      ! minute
+
+
   if(mype == 0) then
      write(6,*)' in get_gefs_for_regional, iadate_gefs=',iadate_gfs
      write(6,*)' in get_gefs_for_regional, iadate    =',iadate
   end if
-           call w3fs21(iadate,nming1)
-           call w3fs21(iadate_gfs,nming2)
+
+     call w3fs21(iadate,nming1)
+     call w3fs21(iadate_gfs,nming2)
+
   if( (nming1/=nming2) .and. (.not.l_ens_in_diff_time) ) then
      if(mype == 0) write(6,*)' GEFS ENSEMBLE MEMBER DATE NOT EQUAL TO ANALYSIS DATE, PROGRAM STOPS'
      call stop2(85)
   end if
 
-
 !         set up ak5,bk5,ck5 for use in computing 3d pressure field (needed for vertical interp to regional)
 !                            following is code segment from gesinfo.F90
   allocate(ak5(nsig_gfs+1))
   allocate(bk5(nsig_gfs+1))
   allocate(ck5(nsig_gfs+1))
   allocate(tref5(nsig_gfs))
+
+  idvc=gfshead%idvc
+  idsl=gfshead%idsl
   do k=1,nsig_gfs+1
      ak5(k)=zero
      bk5(k)=zero
      ck5(k)=zero
   end do
-  if (.not. use_gfs_nemsio) then
+  if ((.not. use_gfs_nemsio) .and. (.not. use_gfs_ncio))then
      if (sighead%nvcoord == 1) then
         do k=1,sighead%levs+1
            bk5(k) = sighead%vcoord(k,1)
@@ -414,6 +502,26 @@ subroutine get_gefs_for_regional
         write(6,*)'GET_GEFS_FOR_REGIONAL:  ***ERROR*** INVALID value for nvcoord=',sighead%nvcoord
         call stop2(85)
      endif
+  else if ( use_gfs_ncio ) then 
+     if (gfshead%nvcoord == 1) then
+        do k=1,nsig_gfs+1
+           bk5(k) = gfsheadv%vcoord(k,1)
+        end do
+     elseif (gfshead%nvcoord == 2) then
+        do k = 1,nsig_gfs+1
+           ak5(k) = gfsheadv%vcoord(k,1)*zero_001
+           bk5(k) = gfsheadv%vcoord(k,2)
+        end do
+     elseif (gfshead%nvcoord == 3) then
+        do k = 1,nsig_gfs+1
+           ak5(k) = gfsheadv%vcoord(k,1)*zero_001
+           bk5(k) = gfsheadv%vcoord(k,2)
+           ck5(k) = gfsheadv%vcoord(k,3)*zero_001
+        end do
+     else
+        write(6,*)'GET_GEFS_FOR_REGIONAL netCDF:  ***ERROR*** INVALID value for nvcoord=',gfshead%nvcoord
+        call stop2(85)
+     endif
   else
      if (nvcoord == 1) then
         do k=1,nsig_gfs+1
@@ -436,6 +544,7 @@ subroutine get_gefs_for_regional
      endif
   end if
 
+
   if(mype == 0 .and. print_verbose)then
      do k=1,nsig_gfs+1
         write(6,*)' ak5,bk5,ck5=',ak5(k),bk5(k),ck5(k)
@@ -509,6 +618,9 @@ subroutine get_gefs_for_regional
      if(use_gfs_nemsio)then
         call general_read_gfsatm_nems(grd_gfst,sp_gfs,filename,uv_hyb_ens,.false.,.true., &
                atm_bundle,.true.,iret)
+     else if (use_gfs_ncio) then
+        call general_read_gfsatm_nc(grd_gfst,sp_gfs,filename,uv_hyb_ens,.false.,.true., &
+               atm_bundle,.true.,iret)
      else
         call general_read_gfsatm(grd_gfst,sp_gfs,sp_gfs,filename,uv_hyb_ens,.false.,.true., &
                atm_bundle,inithead,iret)
@@ -543,6 +655,8 @@ subroutine get_gefs_for_regional
            end do
         end do
      end do
+
+
      kz=num_fields ; kps=kz-1
      do j=1,grd_gfs%lon2
         do i=1,grd_gfs%lat2
@@ -805,6 +919,7 @@ subroutine get_gefs_for_regional
 
   end do   !  end loop over ensemble members.
 
+
 !   next, compute mean of ensembles.
 
   allocate(stbar(grd_mix%lat2,grd_mix%lon2,grd_mix%nsig))

src/gsi/gfs_stratosphere.f90

@@ -109,7 +109,7 @@ subroutine init_gfs_stratosphere
 
    use_gfs_stratosphere=.false.
    good_o3mr=.false.
-   nsig_max=120
+   nsig_max=127
    k0m=0
    k1m=0
    k0r=0
@@ -1153,7 +1153,12 @@ subroutine add_gfs_stratosphere
          ncdim = get_dim(atmges, 'grid_xt'); lonb = ncdim%len
          ncdim = get_dim(atmges, 'grid_yt'); latb = ncdim%len
          ncdim = get_dim(atmges, 'pfull'); levs = ncdim%len
+         if ( mype == 0 ) then
+            write(6,*) 'ncio lonb latb lev =',lonb, latb, levs
+         endif
          njcap = -9999
+         ! FV3GFS write component does not include JCAP, infer from DIMY-2
+         if (njcap<0) njcap=latb-2
          jcap_org = njcap
          idate6 = get_idate_from_time_units(atmges)
          call read_vardata(atmges, 'time', fhour2)