update NMC_Bkerror utility to handle netcdf Gaussian files (#1)

src/NMC_Bkerror/py-bkerror/avg_bkerror.py

@@ -0,0 +1,42 @@
+from bkerror import bkerror
+import numpy as np
+
+nmon = 0
+for mon in ['jan','feb','mar','apr','may','june','july','aug','sept','oct','nov','dec']:
+    filename = '../../../../../staticB/24h/global_berror.l127y770.f77_%ssmooth0p5' % (mon,)
+    print nmon,filename
+    nsig,nlat,nlon = bkerror.get_header(filename)
+    ivar,agvin,bgvin,wgvin,corzin,hscalesin,vscalesin,corq2in,corsstin,hsstin,corpin,hscalespin = bkerror.get_bkerror(filename,nsig,nlat,nlon)
+    if not nmon:
+        print 'initalize arrays'
+        agvout = np.zeros(agvin.shape, agvin.dtype)
+        bgvout = np.zeros(bgvin.shape, bgvin.dtype)
+        wgvout = np.zeros(wgvin.shape, wgvin.dtype)
+        corzout = np.zeros(corzin.shape, corzin.dtype)
+        hscalesout = np.zeros(hscalesin.shape, hscalesin.dtype)
+        vscalesout = np.zeros(vscalesin.shape, vscalesin.dtype)
+        corq2out = np.zeros(corq2in.shape, corq2in.dtype)
+        corsstout = np.zeros(corsstin.shape, corsstin.dtype)
+        hsstout = np.zeros(hsstin.shape, hsstin.dtype)
+        corpout = np.zeros(corpin.shape, corpin.dtype)
+        hscalespout = np.zeros(hscalespin.shape, hscalespin.dtype)
+    agvout += agvin/12.
+    bgvout += bgvin/12.
+    wgvout += wgvin/12.
+    corzout += corzin/12.
+    hscalesout += hscalesin/12.
+    vscalesout += vscalesin/12.
+    corq2out += corq2in/12.
+    corsstout += corsstin/12.
+    hsstout += hsstin/12.
+    corpout += corpin/12.
+    hscalespout += hscalespin/12.
+    nmon += 1
+
+filename = '../../../../../staticB/24h/global_berror.l127y770.f77_annmeansmooth0p5'
+print filename
+bkerror.put_bkerror(filename,ivar,\
+        agvout,bgvout,wgvout,\
+        corzout,hscalesout,vscalesout,\
+        corq2out,corsstout,hsstout,corpout,hscalespout)
+

src/NMC_Bkerror/sorc/CMakeLists.txt

@@ -40,5 +40,6 @@ target_link_libraries(calcstats.x PRIVATE w3emc::w3emc_d)
 target_link_libraries(calcstats.x PRIVATE sp::sp_4)
 target_link_libraries(calcstats.x PRIVATE sigio::sigio)
 target_link_libraries(calcstats.x PRIVATE nemsio::nemsio)
+target_link_libraries(calcstats.x PRIVATE ncio::ncio)
 
 install(TARGETS calcstats.x RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})

src/NMC_Bkerror/sorc/getcases.f90

@@ -5,17 +5,20 @@ subroutine getcases(numcases,mype)
   use kinds, only: r_kind,r_single
   use variables, only: ak5,bk5,ck5,maxcases,nsig,dimbig,hybrid,&
       filename,na,nb,zero,idpsfc5,idvm5,idthrm5,idvc5,ntrac5,cp5,&
-      use_enkf,use_gfs_nemsio,ncepgfs_head,nlonin,nlatin
+      nlon,nlat,use_gfs_ncio,use_enkf,use_gfs_nemsio,ncepgfs_head,nlonin,nlatin
   use sigio_module, only: sigio_head,sigio_srhead,sigio_sclose,&
        sigio_sropen
   use nemsio_module, only:  nemsio_init,nemsio_open,nemsio_close
   use nemsio_module, only:  nemsio_gfile,nemsio_getfilehead
+  use module_ncio, only: Dataset, Variable, Dimension, open_dataset,&
+                         get_idate_from_time_units,&
+                         read_attribute, close_dataset, get_dim, read_vardata
   implicit none
 
   integer,parameter:: i_missing=-9999 
   integer,dimension(4):: idate4
-  integer nmin24(dimbig),nmin48(dimbig),idate5(5)
-  integer nmina,nminb
+  integer nmin24(dimbig),nmin48(dimbig),idate5(5),idate6(6)
+  integer nmina,nminb,nsigin
 
   real*4 fhour4
 
@@ -31,10 +34,13 @@ subroutine getcases(numcases,mype)
   integer,dimension(7):: idate
   integer :: nfhour, nfminute, nfsecondn, nfsecondd, ncount24, ncount48
   real(r_single),allocatable,dimension(:,:,:) :: nems_vcoord
+  real(r_single), allocatable, dimension(:) :: values_1d,ak,bk
 
   type(sigio_head):: sighead
   type(ncepgfs_head):: gfshead
   type(nemsio_gfile) :: gfile
+  type(Dataset) :: dset
+  type(Dimension) :: londim,latdim,levdim
   logical :: isfile
 
   if (mype==0) write(6,*) 'BEGIN TESTCASES'
@@ -92,6 +98,40 @@ subroutine getcases(numcases,mype)
         nlonin=gfshead%lonb
         nlatin=gfshead%latb
 
+     else if (use_gfs_ncio) then
+        dset = open_dataset(filename(loop),errcode=iret2)
+        if ( iret2 /= 0 ) then
+           write(6,*)' getcases:  ***ERROR*** problem open_dataset file = ', &
+              trim(filename(loop)),', Status = ',iret2
+           stop
+        end if
+        idate6 = get_idate_from_time_units(dset)
+        idate5 = 0; idate5(1:4) = idate6(1:4)
+        call read_vardata(dset,'time',values_1d,errcode=iret2)
+        if (iret2 /= 0) then
+           print *,'error reading time from',trim(filename(loop))
+           stop
+        endif
+        fhour5 = values_1d(1)
+        deallocate(values_1d)
+        londim = get_dim(dset,'grid_xt'); nlonin = londim%len
+        latdim = get_dim(dset,'grid_yt'); nlatin = latdim%len
+        levdim = get_dim(dset,'pfull');   nsigin = levdim%len
+        if (nsig .ne. nsigin) then
+           print *,'number of vert levels in namelist does not match file'
+           stop
+        endif
+        if (nlat-2 .ne. nlatin) then
+           print *,'number of lats in namelist does not match file',&
+           nlat,nlatin
+           stop
+        endif
+        if (nlon .ne. nlonin) then
+           print *,'number of lons in namelist does not match file'
+           stop
+        endif
+        if (loop == 1 .and. mype == 0) print *,'ncio idate5,fhour5',idate5,fhour5
+        call close_dataset(dset)
      else ! not use_gfs_nemsio
 
        call sigio_sropen(inges,filename(loop),iret)
@@ -226,7 +266,31 @@ subroutine getcases(numcases,mype)
      vcoord5(:,1:nvcoord5)=nems_vcoord(:,1:nvcoord5,1)
 
      deallocate(nems_vcoord)
-
+  else if ( use_gfs_ncio ) then
+     ! hardwired for now
+     idpsfc5 = 2
+     idvc5   = 2
+     idthrm5 = 2
+     ntrac5  = 3 ! (spfh,o3mr,total condensate)
+     dset = open_dataset(filename(1),errcode=iret2)
+     call read_attribute(dset, 'ak', values_1d)
+     nvcoord5 = 2
+     allocate(vcoord5(nsig+1,nvcoord5))
+     do k=1,nsig+1
+        vcoord5(nsig-k+2,1) = values_1d(k)
+     enddo
+     call read_attribute(dset, 'bk', values_1d)
+     do k=1,nsig+1
+        vcoord5(nsig-k+2,2) = values_1d(k)
+     enddo
+     deallocate(values_1d)
+     if (mype == 0) then
+        print *,'ncio ak,bk:'
+        print *,vcoord5(:,1)
+        print *,vcoord5(:,2)
+     endif
+     call close_dataset(dset)
+
   else ! not use_gfs_nemsio
 
 ! DTK NEW  EXTRACT SOME STUFF FROM THE FIRST LISTED FILE
@@ -261,7 +325,7 @@ subroutine getcases(numcases,mype)
      end if
 
      call sigio_sclose(inges,iret)
-  endif !no use_gfs_nemsio
+  endif 
 
   do k=1,nsig+1
      ak5(k)=zero