Add Grib2 writing capabilities to the ocnicepost module

.vscode/settings.json

@@ -0,0 +1,3 @@
+{
+    "fortran.fortls.disabled": true
+}

parm/ocnicepost/ice.csv

@@ -1,12 +1,11 @@
-! variable name, dimension, grid location, remapping method, vector pair, grid location of vector pair
-       'hi_h', 2 , 'Ct', 'bilinear', '', ''
-       'hs_h', 2 , 'Ct', 'bilinear', '', ''
-     'aice_h', 2 , 'Ct', 'bilinear', '', ''
-     'Tsfc_h', 2 , 'Ct', 'bilinear', '', ''
-     'uvel_h', 2 , 'Bu_x', 'bilinear', 'vvel_h', 'Bu_y'
-     'vvel_h', 2 , 'Bu_y', 'bilinear', 'uvel_h', 'Bu_x'
-   'frzmlt_h', 2 , 'Ct', 'conserve', '', ''
-   'albsni_h', 2 , 'Ct', 'bilinear', '', ''
-'mlt_onset_h', 2 , 'Ct', 'bilinear', '', ''
-'frz_onset_h', 2 , 'Ct', 'bilinear', '', ''
-      'sst_h', 2 , 'Ct', 'bilinear', '', ''
+! variable name, dimension, grid location, remapping method, vector pair, grid location of vector pair,gb2 varname,discription,unit,gc1,gc2,gc3,gc4,gc5,gc6,gc7,gc8
+       'hi_h',2, 'Ct', 'bilinear', '', '', 'ICETK', 'Ice Thickness', 'm',10,2,7,0,2,1,1,0
+       'hs_h',2, 'Ct', 'bilinear', '', '', 'SNVOLSI', 'Snow Volume Over Sea Ice per Unit Area', 'm3m-2',10,2,7,0,2,16,1,0
+     'aice_h',2, 'Ct', 'bilinear', '', '', 'ICEC', 'Ice Cover', 'Proportion',10,2,7,0,2,0,1,0
+     'Tsfc_h',2, 'Ct', 'bilinear', '', '', 'ICETMP', 'Ice Temperature', 'K',10,2,7,0,2,8,1,0
+     'uvel_h',2, 'Bu_x', 'bilinear', 'vvel_h', 'Bu_y', 'UICE', 'U-Component of Ice Drift', 'm/s',10,2,7,0,2,4,1,0
+     'vvel_h',2, 'Bu_y', 'bilinear', 'uvel_h', 'Bu_x', 'VICE', 'V-Component of Ice Drift', 'm/s',10,2,7,0,2,5,1,0
+   'frzmlt_h',2, 'Ct', 'conserve', '', '', 'FRZMLTPOT', 'freeze/melt potential', 'W m-2',10,2,7,0,2,27,1,0
+   'albsni_h',2, 'Ct', 'bilinear', '', '', 'ALBDOICE', 'Albedo', 'Numeric',10,2,7,0,2,14,1,0
+'mlt_onset_h',2, 'Ct', 'bilinear', '', '', 'MLTDATE', 'melt onset date', 'Numeric',10,2,7,0,2,28,1,0
+'frz_onset_h',2, 'Ct', 'bilinear', '', '', 'FRZDATE', 'freeze onset date', 'Numeric',10,2,7,0,2,29,1,0

parm/ocnicepost/ice_mx025_to_0p25.nml

@@ -7,5 +7,6 @@ maskvar='tmask'
 sinvar=''
 cosvar=''
 angvar='ANGLET'
+grib2=.true.
 debug=.false.
 /

parm/ocnicepost/ocean.csv

@@ -1,26 +1,20 @@
-! variable name, dimension, grid location, remapping method, vector pair, grid location of vector pair
-      'SSH', 2 ,'Ct' , 'bilinear' , '', ''
-      'SST', 2 ,'Ct' , 'bilinear' , '', ''
-      'SSS', 2 ,'Ct' , 'bilinear' , '', ''
-    'speed', 2 ,'Ct' , 'bilinear' , '', ''
-     'ePBL', 2 ,'Ct' , 'bilinear' , '', ''
-  'MLD_003', 2 ,'Ct' , 'bilinear' , '', ''
- 'MLD_0125', 2 ,'Ct' , 'bilinear' , '', ''
-   'latent', 2 ,'Ct' , 'conserve' , '', ''
- 'sensible', 2 ,'Ct' , 'conserve' , '', ''
-       'SW', 2 ,'Ct' , 'conserve' , '', ''
-       'LW', 2 ,'Ct' , 'conserve' , '', ''
-     'evap', 2 ,'Ct' , 'conserve' , '', ''
-    'lprec', 2 ,'Ct' , 'conserve' , '', ''
-    'fprec', 2 ,'Ct' , 'conserve' , '', ''
-'LwLatSens', 2 ,'Ct' , 'conserve' , '', ''
- 'Heat_PmE', 2 ,'Ct' , 'conserve' , '', ''
-      'SSU', 2 ,'Cu' , 'bilinear' , 'SSV'  , 'Cv'
-      'SSV', 2 ,'Cv' , 'bilinear' , 'SSU'  , 'Cu'
-     'taux', 2 ,'Cu' , 'conserve' , 'tauy' , 'Cv'
-     'tauy', 2 ,'Cv' , 'conserve' , 'taux' , 'Cu'
-     'temp', 3 ,'Ct' , 'bilinear' , '', ''
-      'tob', 2 ,'Ct' , 'bilinear' , '', ''
-       'so', 3 ,'Ct' , 'bilinear' , '', ''
-       'uo', 3 ,'Cu' , 'bilinear' , 'vo' , 'Cv'
-       'vo', 3 ,'Cv' , 'bilinear' , 'uo' , 'Cu'
+! variable name, dimension, grid location, remapping method, vector pair, grid location of vector pair,gb2 varname,discription,unit,gc1,gc2,gc3,gc4,gc5,gc6,gc7,gc8
+      'SSH',2,'Ct' , 'bilinear' , '', '','SSHG','Sea Surface Height Relative to Geoid','m',10,0,7,1,3,195,101,0
+      'SST',2,'Ct' , 'bilinear' , '', '','WTMP','Water Temperature','K',10,1,7,0,3,0,101,0
+      'SSS',2,'Ct' , 'bilinear' , '', '',SALIN','Salinity','psu',10,0,7,1,4,193,101,0
+    'speed',2,'Ct' , 'bilinear' , '', '','SPC','Current Speed','m/s',10,1,7,0,1,1,101,0
+  'MLD_003',2,'Ct' , 'bilinear' , '', '','WDEPTH','Water Depth','m',10,1,7,0,4,14,101,0
+   'latent',2,'Ct' , 'conserve' , '', '','LHTFL','Latent Heat Net Flux','W/m^2',0,1,7,0,0,10,101,0
+ 'sensible',2,'Ct' , 'conserve' , '', '','SHTFL','Sensible Heat Net Flux','W/m^2',0,1,7,0,0,11,101,0
+       'SW',2,'Ct' , 'conserve' , '', '','NSWRF','Net Short Wave Radiation Flux','W/m^2',0,1,7,0,4,9,101,0
+       'LW',2,'Ct' , 'conserve' , '', '','NLWRF','Net Long-Wave Radiation Flux','W/m^2',0,1,7,0,5,5,101,0
+'LwLatSens',2,'Ct' , 'conserve' , '', '','THFLX','Total Downward Heat Flux at Surface','W/m^2',0,0,7,1,0,197,101,0
+ 'Heat_PmE',2,'Ct' , 'conserve' , '', '','DWHFLUX','Downward Heat Flux','W/m^2',10,1,7,0,3,4,101,0
+      'SSU',2,'Cu' , 'bilinear' , 'SSV'  , 'Cv','UOGRD','U-Component of Current','m/s',10,1,7,0,1,2,101,0
+      'SSV',2,'Cv' , 'bilinear' , 'SSU'  , 'Cu','VOGRD','V-Component of Current','m/s',10,1,7,0,1,3,101,0
+     'taux',2,'Cu' , 'conserve' , 'tauy' , 'Cv','XCOMPSS','x-component Surface Stress','N/m^2',10,1,7,0,3,7,101,0
+     'tauy',2,'Cv' , 'conserve' , 'taux' , 'Cu','YCOMPSS','y-component Surface Stress','N/m^2',10,1,7,0,3,7,101,0
+     'temp',3,'Ct' , 'bilinear' , '', '','WTMP','Water Temperature','K',10,1,7,0,3,0,160,0
+       'so',3,'Ct' , 'bilinear' , '', '','SALIN','3-D Salinity','psu',10,0,7,1,4,193,160,0
+       'uo',3,'Cu' , 'bilinear' , 'vo' , 'Cv','UOGRD','U-Component of Current','m/s',10,1,7,0,1,2,160,0
+       'vo',3,'Cv' , 'bilinear' , 'uo' , 'Cu','VOGRD','V-Component of Current','m/s',10,1,7,0,1,3,160,0

parm/ocnicepost/ocean_mx025_to_0p25.nml

@@ -7,5 +7,6 @@ maskvar='temp'
 sinvar='sin_rot'
 cosvar='cos_rot'
 angvar=''
+grib2=.true.
 debug=.false.
 /

parm/ocnicepost/ocnicepost.nml.jinja2

@@ -7,5 +7,6 @@
   sinvar = "{{ sinvar }}",
   cosvar = "{{ cosvar }}",
   angvar = "{{ angvar }}",
+  grib2 = {{ grib2 | default('.true.') }}
   debug = {{ debug | default('.false.') }}
 /

src/ocnicepost.fd/CMakeLists.txt

@@ -1,13 +1,16 @@
 list(APPEND fortran_src
-  arrays_mod.F90
   init_mod.F90
+  arrays_mod.F90
+  utils_mod.F90
   masking_mod.F90
   ocnicepost.F90
-  utils_mod.F90
+
 )
 
 set(exe_name ocnicepost.x)
 add_executable(${exe_name} ${fortran_src})
-target_link_libraries(${exe_name} PRIVATE NetCDF::NetCDF_Fortran)
+target_link_libraries(${exe_name} PRIVATE NetCDF::NetCDF_Fortran
+                                          g2::g2_d
+                                          bacio::bacio_4)
 
 install(TARGETS ${exe_name} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})

src/ocnicepost.fd/init_mod.F90

@@ -7,15 +7,27 @@ module init_mod
   real, parameter :: maxvars = 50              !< The maximum number of fields expected in a source file
 
   type :: vardefs
-     character(len= 20)   :: var_name          !< A variable's variable name
-     character(len=120)   :: long_name         !< A variable's long name
-     character(len= 20)   :: units             !< A variable's unit
-     character(len= 20)   :: var_remapmethod   !< A variable's mapping method
-     integer              :: var_dimen         !< A variable's dimensionality
-     character(len=  4)   :: var_grid          !< A variable's input grid location; all output locations are on cell centers
-     character(len= 20)   :: var_pair          !< A variable's pair
-     character(len=  4)   :: var_pair_grid     !< A pair variable grid
-     real                 :: var_fillvalue     !< A variable's fillvalue
+       character(len= 20)   :: var_name          !< A variable's variable name
+       character(len= 20)   :: var_remapmethod   !< A variable's mapping method
+       character(len=120)   :: long_name         !< A variable's long name
+       character(len= 20)   :: units             !< A variable's unit
+       integer              :: var_dimen         !< A variable's dimensionality
+       character(len=  4)   :: var_grid          !< A variable's input grid location; all output locations are on cell centers
+       character(len= 20)   :: var_pair          !< A variable's pair
+       character(len=  4)   :: var_pair_grid     !< A pair variable grid
+       real                 :: var_fillvalue     !< A variable's fillvalue
+       character(len= 20)   :: name_gb2          !< A variable's grib2 variable name
+       character(len=120)   :: discription_gb2   !< A variable's discription
+       character(len= 20)   :: unit_gb2          !< A variable's unit       
+  !!!may need to add a fillvalue for grib2 file
+       integer              ::var_g1             !< Variables' grib2 coefficients g1-Dissipline
+       integer              ::var_g2             !< Variables' grib2 coefficients g2-Master Tables Version Number
+       integer              ::var_g3             !< Variables' grib2 coefficients g3-Section 1 originating center, used for local tables
+       integer              ::var_g4             !< Variables' grib2 coefficients g4-Section 1 Local Tables Version Number
+       integer              ::var_g5             !< Variables' grib2 coefficients g5-Section 4 Template 4.0 Parameter category
+       integer              ::var_g6             !< Variables' grib2 coefficients g6-Section 4 Template 4.0 Parameter number
+       integer              ::var_g7             !< Variables' grib2 coefficients g7-Level ID
+       integer              ::var_g8             !< Variables' grib2 coefficients g8-
   end type vardefs
 
   type(vardefs) :: outvars(maxvars)            !< An empty structure filled by reading a csv file describing the fields
@@ -40,6 +52,7 @@ module init_mod
   integer :: nyr        !< The y-dimension of the destination rectilinear grid
 
   integer :: logunit    !< The log unit
+  logical :: grib2      !< If true, write grib2 message
   logical :: debug      !< If true, print debug messages and intermediate files
   logical :: do_ocnpost !< If true, the source file is ocean, otherwise ice
 
@@ -53,7 +66,7 @@ subroutine readnml
     integer :: srcdims(2), dstdims(2)
 
     namelist /ocnicepost_nml/ ftype, srcdims, wgtsdir, dstdims, maskvar, sinvar, cosvar, &
-         angvar, debug
+         angvar, grib2, debug
 
     ! --------------------------------------------------------
     ! read the name list
@@ -120,8 +133,8 @@ subroutine readcsv(nvalid)
 
     character(len= 40) :: fname
     character(len=100) :: chead
-    character(len= 20) :: c1,c3,c4,c5,c6
-    integer :: i2
+    character(len= 20) :: c1,c3,c4,c5,c6,c7,c8,c9
+    integer :: i2,i10,i11,i12,i13,i14,i15,i16,i17
     integer :: nn,n,ierr,iounit
 
     ! --------------------------------------------------------
@@ -138,7 +151,7 @@ subroutine readcsv(nvalid)
     read(iounit,*)chead
     nn=0
     do n = 1,maxvars
-       read(iounit,*,iostat=ierr)c1,i2,c3,c4,c5,c6
+       read(iounit,*,iostat=ierr)c1,i2,c3,c4,c5,c6,c7,c8,c9,i10,i11,i12,i13,i14,i15,i16,i17
        if (ierr .ne. 0) exit
        if (len_trim(c1) > 0) then
           nn = nn+1
@@ -148,6 +161,17 @@ subroutine readcsv(nvalid)
           outvars(nn)%var_remapmethod = trim(c4)
           outvars(nn)%var_pair = trim(c5)
           outvars(nn)%var_pair_grid = trim(c6)
+          outvars(nn)%name_gb2 = trim(c7)
+          outvars(nn)%discription_gb2 = trim(c8)
+          outvars(nn)%unit_gb2 = trim(c9)
+          outvars(nn)%var_g1 = i10
+          outvars(nn)%var_g2 = i11
+          outvars(nn)%var_g3 = i12
+          outvars(nn)%var_g4 = i13
+          outvars(nn)%var_g5 = i14
+          outvars(nn)%var_g6 = i15
+          outvars(nn)%var_g7 = i16
+          outvars(nn)%var_g8 = i17
        end if
     end do
     close(iounit)

src/ocnicepost.fd/masking_mod.F90

@@ -27,7 +27,7 @@ subroutine remap_masks(vfill)
     ! local variables
     integer            :: rc, ncid, varid, n
     character(len=240) :: wgtsfile
-    real               :: minlat = -78.63       
+    real               :: minlat = -78.00    ! This is masking Antarctica which is not masked in the weighting file
 
     real, allocatable, dimension(:) :: out1d
 
@@ -80,7 +80,7 @@ subroutine remap_masks(vfill)
 
        if (debug) then
           write(logunit,'(a,2g14.4)')'mask min/max on destination grid ',minval(rgmask3d),maxval(rgmask3d)
-          call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgmask3d.nc', 'rgmask3d', dims=(/nxr,nyr,nlevs/), &
+          call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgmask3d.nc', 'rgmask3d', dims=(/nxr*nyr,nlevs/), &  
                field=rgmask3d)
        end if
     else
@@ -93,7 +93,7 @@ subroutine remap_masks(vfill)
 
        if (debug) then
           write(logunit,'(a,2g14.4)')'mask min/max on destination grid ',minval(rgmask2d),maxval(rgmask2d)
-          call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgmask2d.nc', 'rgmask2d', dims=(/nxr,nyr/),       &
+          call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgmask2d.nc', 'rgmask2d', dims=(/nxr*nyr/),       &    
                field=rgmask2d)
        end if
     end if

src/ocnicepost.fd/ocnicepost.F90

@@ -27,16 +27,17 @@ program ocnicepost
   use netcdf
   use init_mod   , only : nxt, nyt, nlevs, nxr, nyr, outvars, readnml, readcsv
   use init_mod   , only : wgtsdir, ftype, fsrc, fdst, input_file, cosvar, sinvar, angvar
-  use init_mod   , only : do_ocnpost, debug, logunit
+  use init_mod   , only : do_ocnpost, debug, logunit, grib2
+  use init_mod   , only : vardefs
   use arrays_mod , only : b2d, c2d, b3d, rgb2d, rgc2d, rgb3d, dstlon, dstlat, setup_packing
   use arrays_mod , only : nbilin2d, nbilin3d, nconsd2d, bilin2d, bilin3d, consd2d
   use masking_mod, only : mask2d, mask3d, rgmask2d, rgmask3d, remap_masks
-  use utils_mod  , only : getfield, packarrays, remap, dumpnc, nf90_err
+  use utils_mod  , only : getfield, packarrays, remap, dumpnc, nf90_err, write_grib2_2d, write_grib2_3d
 
   implicit none
 
   character(len=120) :: wgtsfile
-  character(len=120) :: fout
+  character(len=120) :: fout, gout
 
   ! dimensions, units and variables from source file used in creation of
   ! output netcdf
@@ -49,6 +50,10 @@ program ocnicepost
   real, allocatable, dimension(:,:)   :: out2d !< 2D destination grid output array
   real, allocatable, dimension(:,:,:) :: out3d !< 3D destination grid output array
 
+  ! arrays for output grib2
+  real, allocatable, dimension(:,:)   :: grib2d    !< 2D destination grib2 concat array
+  type(vardefs), allocatable, dimension(:) :: g2d  !< concatinated variable metadata for 2D source fields remap
+
   real(kind=8)       :: timestamp
   character(len= 40) :: timeunit, timecal
   character(len= 20) :: vname, vunit
@@ -132,6 +137,8 @@ program ocnicepost
 
   call remap_masks(vfill)
 
+
+
   ! --------------------------------------------------------
   ! create packed arrays for mapping and remap packed arrays
   ! to the destination grid
@@ -250,6 +257,43 @@ program ocnicepost
           sqrt(rgb2d(:,idx2)**2 + rgb2d(:,idx3)**2)
   end if
 
+!--------------------------------------------------------
+! Write the grib2 files
+!--------------------------------------------------------
+
+  if(grib2) then
+   gout = trim(ftype)//'.'//trim(fdst)//'.gb2'
+   if (debug) write(logunit, '(a)')'GRIB2 2D output file: '//trim(gout)
+
+   if (allocated(rgb2d) .and. allocated(rgc2d)) then
+      allocate(grib2d(nxr*nyr,nconsd2d+nbilin2d), source=0.0)
+      allocate(g2d(1:nconsd2d+nbilin2d)) 
+      grib2d(:, 1:nbilin2d) = rgb2d
+      grib2d(:, nbilin2d+1:nconsd2d+nbilin2d) = rgc2d
+      g2d(1:nbilin2d) = b2d
+      g2d(nbilin2d+1:nconsd2d+nbilin2d) = c2d
+   else if (allocated(rgb2d)) then
+      allocate(grib2d(nxr*nyr,nconsd2d), source=0.0)
+      allocate(g2d(1:nconsd2d))
+      grib2d(:, 1:nconsd2d) = rgc2d
+      g2d(1:nconsd2d) = c2d
+   else if (allocated(rgc2d)) then
+      allocate(grib2d(nxr*nyr,nbilin2d), source=0.0)
+      allocate(g2d(1:nbilin2d))
+      grib2d(:, 1:nbilin2d) = rgb2d
+      g2d(1:nbilin2d) = b2d
+   end if
+
+   call write_grib2_2d(gout, g2d, (/nxr,nyr/), nconsd2d+nbilin2d, grib2d, vfill)
+
+   if (allocated(rgb3d)) then
+      gout = trim(ftype)//'.'//trim(fdst)//'_3D.gb2'
+      call write_grib2_3d(gout, b3d, (/nxr,nyr,nlevs/), nbilin3d, rgb3d, vfill)
+      if (debug) write(logunit, '(a)')'GRIB2 3D output file: '//trim(gout)
+   end if
+  end if
+
+
   ! --------------------------------------------------------
   ! write the mapped fields
   ! --------------------------------------------------------
@@ -357,6 +401,7 @@ program ocnicepost
         call nf90_err(nf90_put_var(ncid,   varid, out2d), 'put variable: '//vname)
      end do
   end if
+
   if (allocated(rgb3d)) then
      do n = 1,nbilin3d
         out3d(:,:,:) = reshape(rgb3d(:,:,n), (/nxr,nyr,nlevs/))
@@ -371,4 +416,4 @@ program ocnicepost
 
   stop
 
-end program ocnicepost
+end program ocnicepost