diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 00000000..63e4335f --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,3 @@ +{ + "fortran.fortls.disabled": true +} \ No newline at end of file diff --git a/parm/ocnicepost/ice.csv b/parm/ocnicepost/ice.csv index c42ff4dc..c06c65fe 100644 --- a/parm/ocnicepost/ice.csv +++ b/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 diff --git a/parm/ocnicepost/ice_mx025_to_0p25.nml b/parm/ocnicepost/ice_mx025_to_0p25.nml index 89a7d016..e8cd3d99 100644 --- a/parm/ocnicepost/ice_mx025_to_0p25.nml +++ b/parm/ocnicepost/ice_mx025_to_0p25.nml @@ -7,5 +7,6 @@ maskvar='tmask' sinvar='' cosvar='' angvar='ANGLET' +grib2=.true. debug=.false. / diff --git a/parm/ocnicepost/ocean.csv b/parm/ocnicepost/ocean.csv index 1d853e98..e54ff2b5 100644 --- a/parm/ocnicepost/ocean.csv +++ b/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 diff --git a/parm/ocnicepost/ocean_mx025_to_0p25.nml b/parm/ocnicepost/ocean_mx025_to_0p25.nml index 0b800724..567bcd80 100644 --- a/parm/ocnicepost/ocean_mx025_to_0p25.nml +++ b/parm/ocnicepost/ocean_mx025_to_0p25.nml @@ -7,5 +7,6 @@ maskvar='temp' sinvar='sin_rot' cosvar='cos_rot' angvar='' +grib2=.true. debug=.false. / diff --git a/parm/ocnicepost/ocnicepost.nml.jinja2 b/parm/ocnicepost/ocnicepost.nml.jinja2 index 4cc6c596..8dd1e8bf 100644 --- a/parm/ocnicepost/ocnicepost.nml.jinja2 +++ b/parm/ocnicepost/ocnicepost.nml.jinja2 @@ -7,5 +7,6 @@ sinvar = "{{ sinvar }}", cosvar = "{{ cosvar }}", angvar = "{{ angvar }}", + grib2 = {{ grib2 | default('.true.') }} debug = {{ debug | default('.false.') }} / diff --git a/src/ocnicepost.fd/CMakeLists.txt b/src/ocnicepost.fd/CMakeLists.txt index 529454fa..c5d213a5 100644 --- a/src/ocnicepost.fd/CMakeLists.txt +++ b/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}) diff --git a/src/ocnicepost.fd/init_mod.F90 b/src/ocnicepost.fd/init_mod.F90 index 2949a088..3aa8824c 100644 --- a/src/ocnicepost.fd/init_mod.F90 +++ b/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) diff --git a/src/ocnicepost.fd/masking_mod.F90 b/src/ocnicepost.fd/masking_mod.F90 index 3613d6f5..b3ac9c99 100644 --- a/src/ocnicepost.fd/masking_mod.F90 +++ b/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 diff --git a/src/ocnicepost.fd/ocnicepost.F90 b/src/ocnicepost.fd/ocnicepost.F90 index 975219ac..c15c60cb 100644 --- a/src/ocnicepost.fd/ocnicepost.F90 +++ b/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 \ No newline at end of file diff --git a/src/ocnicepost.fd/utils_mod.F90 b/src/ocnicepost.fd/utils_mod.F90 index 590eec09..ee0d4f72 100644 --- a/src/ocnicepost.fd/utils_mod.F90 +++ b/src/ocnicepost.fd/utils_mod.F90 @@ -1,7 +1,9 @@ module utils_mod use netcdf - use init_mod, only : debug, logunit, vardefs, fsrc + use init_mod, only : debug, logunit, vardefs, fsrc, input_file, ftype + + implicit none @@ -39,6 +41,8 @@ module utils_mod public packarrays public remap public dumpnc + public write_grib2_2d + public write_grib2_3d public nf90_err contains @@ -578,6 +582,566 @@ subroutine dumpnc1d(fname, vname, dims, field) end subroutine dumpnc1d + + + !----------------------------------------------------------------------------------- + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Write Grib2 2D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!This subroutine write Grib2 file modified messages!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !----------------------------------------------------------------------------------- + + subroutine write_grib2_2d(fname, gcf, dims, nflds, field, vfill) + + implicit none + + character(len=*), intent(in) :: fname + type(vardefs), intent(in) :: gcf(:) + integer(4), intent(in) :: dims(2) + integer(4), intent(in) :: nflds + real(4), intent(inout) :: field(dims(1)*dims(2),nflds) + real(4), intent(in) :: vfill + + ! internal variables + integer(4) :: max_bytes, lengrib + integer(4) :: ref_time(6) + integer(4) :: lunout, ierr + integer(4) :: fortime, dij, npt + CHARACTER(len=1),allocatable,dimension(:) :: cgrib + real(8) :: tmpfld(size(field,1)) + + ! GRIB2 metadata arrays + integer(4) :: listsec0(2), listsec1(13) + integer(4) :: igdtnum, ipdtnum, idrtnum + integer(4) :: igdtlen, ipdtlen, idrtlen + integer(4) :: jgdt(19), jpdt(15), idrtmpl(16) + integer(4) :: igds(5) + integer(4) :: numcoord, ibmap + real(4) :: coordlist + integer(4) :: n, lon0, lon1, lat0, lat1 + integer(4) :: ideflist, idefnum + logical*1 :: bmp(dims(1)*dims(2)) + + real :: max_val, min_val, mean_val, count_val + + npt = dims(1) * dims(2) + + max_bytes = npt * 4 ! Estimated max bytes + bmp=.true. + + call getlun(lunout) + call baopenw(lunout, trim(fname), ierr) + if (ierr /= 0) then + write(0, *) 'Error opening grib2 file ', trim(fname) + return + end if + + call retrieve_time( fortime , ref_time ) + + ! Initialize GRIB2 message sections + listsec0(1) = gcf(1)%var_g1 ! Discipline - GRIB Master Table Number (Code Table 0.0) + listsec0(2) = 2 ! GRIB Edition Number (currently 2) + + listsec1(1) = gcf(1)%var_g3 ! Originating Centre (Common Code Table C-1) + listsec1(2) = 0 ! Originating Sub-centre (local table) EMC=4 +! listsec1(3) = 32 ! GRIB Master Tables Version Number (Code Table 1.0)-last one currently 32 + listsec1(3) = gcf(1)%var_g2 ! GRIB Master Tables Version Number (Code Table 1.0) + listsec1(4) = 1 ! GRIB Local Tables Version Number (Code Table 1.1) + listsec1(5) = 1 ! Significance of Reference Time (Code Table 1.2) + listsec1(6) = ref_time(1) ! Reference Time - Year -4digits + listsec1(7) = ref_time(2) ! Reference Time - Month + listsec1(8) = ref_time(3) ! Reference Time - Day + listsec1(9) = ref_time(4) ! Reference Time - Hour + listsec1(10) = ref_time(5) ! Reference Time - Minute + listsec1(11) = ref_time(6) ! Reference Time - Second + listsec1(12) = 0 ! Production status of data (Code Table 1.3) + listsec1(13) = 1 ! Type of processed data (Code Table 1.4) + + ! set grid res + if (dims(1) == 1440 .and. dims(2) == 721) dij= 250000 ! 1/4deg rectilinear + if (dims(1) == 720 .and. dims(2) == 361) dij= 500000 ! 1/2deg rectilinear + if (dims(1) == 360 .and. dims(2) == 181) dij= 1000000 ! 1deg rectilinear + if (dims(1) == 72 .and. dims(2) == 36) dij= 5000000 ! 5deg rectilinear + + lon0 = 0 + lon1 = 360000000 - dij + lat0 = -90000000 + lat1 = 90000000 + + ! Populate the jgdt array for Template 3.0 + jgdt(1) = 6 + jgdt(2) = 0 + jgdt(3) = 0 + jgdt(4) = 0 + jgdt(5) = 0 + jgdt(6) = 0 + jgdt(7) = 0 + jgdt(8) = dims(1) + jgdt(9) = dims(2) + jgdt(10) = 0 + jgdt(11) = -1 + jgdt(12) = lat0 + jgdt(13) = lon0 + jgdt(14) = 48 !0 + jgdt(15) = lat1 + jgdt(16) = lon1 + jgdt(17) = dij + jgdt(18) = dij + jgdt(19) = 64 + + igdtnum=0 + ! Define igds GRIB2 - SECTION 3 + igds(1) = 0 ! Source of grid definition + igds(2) = npt ! Number of grid points + igds(3) = 0 ! Number of octets for each additional grid points definition + igds(4) = 0 ! Interpretation of list for optional points definition + igds(5) = igdtnum ! GRIB2 - CODE TABLE 3.1 + + igdtlen=size(jgdt) + + if (debug) then + write(logunit, *) 'listsec0, listsec1: ', listsec0, listsec1 + write(logunit, *) 'igdtnum, igdtlen: ', igdtnum, igdtlen + write(logunit, *) 'jgdt: ', jgdt + write(logunit, *) 'igds: ', igds + write(logunit, *) 'dij: ', dij + write(logunit, *) 'max_bytes: ', max_bytes + write(logunit, *) 'forcast time: ', fortime + write(logunit, *) 'refference time: ', ref_time + end if + + ideflist=0 + idefnum=0 + + do n=1,nflds + + allocate(cgrib(max_bytes)) + + listsec0(1) = gcf(n)%var_g1 + + call gribcreate(cgrib, max_bytes, listsec0, listsec1, ierr) + if (ierr /= 0) then + write(0, *) 'Error initializing GRIB2 message', ierr + return + end if + + ! Compute max, min, and mean + max_val = maxval(field(:,n), mask = field(:,n) .ne. vfill) + min_val = minval(field(:,n), mask = field(:,n) .ne. vfill) + mean_val = sum(field(:,n), mask = field(:,n) .ne. vfill) / count(field(:,n) .ne. vfill) + + if (debug) then + write(logunit, *) 'Variable_name, max, min, mean: ', gcf(n)%var_name, max_val, min_val, mean_val + end if + + call addgrid(cgrib, max_bytes, igds, jgdt, igdtlen, ideflist, idefnum, ierr) ! there is an internal error here + if (ierr /= 0) then + write(0, *) 'Error adding grid to GRIB2 message', ierr + return + end if + +! Create Section 4 parametrs + ipdtnum=0 + + jpdt(1)=gcf(n)%var_g5 ! parm number catagory + jpdt(2)=gcf(n)%var_g6 ! parm number + jpdt(3)=2 ! (0-analysis, 1-initialazation, 2-forecast, .. GRIB2 - CODE TABLE 4.3 ) + jpdt(4)=0 ! + jpdt(5)=96 ! Code ON388 Table A- GFS + jpdt(6)=0 ! + jpdt(7)=0 ! + jpdt(8)=1 ! unit (Hour=1) 6hour=11 (ask later) Table 4.4 + jpdt(9)=fortime ! forecast time + jpdt(10)=gcf(n)%var_g7 ! level ID (1-Ground or Water Surface, 101 mean sea level, 160 depth bellow mean sea level , 168-Ocean Model Layer,...) + jpdt(11)=0 ! + jpdt(12)=0 ! + jpdt(13)=0 + jpdt(14)=0 + jpdt(15)=0 + + + if (debug) write(logunit, *) 'ipdtnum=', ipdtnum, ', jpdt= ', jpdt(1:16) + + ipdtlen=size(jpdt) + + numcoord=0 + coordlist=0. ! needed for hybrid vertical coordinate + + ibmap = 255 ! Bitmap indicator ( see Code Table 6.0 ) -255 no bitmap + bmp=.true. + + if (trim(gcf(n)%name_gb2) .eq. 'WTMP') then + where ( field(:,n) .ne. vfill ) field(:,n) = field(:,n) + 273.15 + endif + + where ( field(:,n) .eq. vfill ) field(:,n)= 9999.0 + + ! Create Section 5 parametrs + idrtnum = 0 ! Template 5.2 (Grid Point Data - complex Packing) + + idrtmpl(:)=0 + ! Populate idrtmpl + idrtmpl(1) = 0 ! Reference value (scaled value of the minimum data point) + idrtmpl(2) = 0 ! Binary scale factor (scale by 2^E) + idrtmpl(3) = 2 ! Decimal scale factor (scale by 10^D) + idrtmpl(4) = 0 ! + idrtmpl(5) = 0 ! + idrtmpl(6) = 0 ! + ! Reserved fields + idrtmpl(7:16) = 0 ! Reserved for future use + + idrtlen=size(idrtmpl) + + if (debug) write(logunit, *) 'idrtmpl: ', idrtmpl + + tmpfld=0 + tmpfld=real(field(:,n), 8) + + + write(logunit, *) 'Variable_name, max, min, mean, count: ', gcf(n)%var_name, max_val, min_val, mean_val, count_val + + call addfield(cgrib, max_bytes, ipdtnum, jpdt, ipdtlen, coordlist, numcoord, & + idrtnum, idrtmpl, idrtlen, tmpfld, npt, ibmap, bmp, ierr) + if (ierr /= 0) then + write(0, *) 'Error adding field to GRIB2 message', ierr + return + end if + + call gribend(cgrib, max_bytes, lengrib, ierr) + if (debug) write(logunit, *) 'gribend status=', ierr + if (debug) write(logunit, *) 'length of the final GRIB2 message in octets =', lengrib + call wryte(lunout, lengrib, cgrib) + + deallocate(cgrib) + + end do + + call baclose(lunout, ierr) + + return + + end subroutine write_grib2_2d + + + + !----------------------------------------------------------------------------------- + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Write Grib2 3D!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!This subroutine write Grib2 file modified messages!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !----------------------------------------------------------------------------------- + + subroutine write_grib2_3d(fname, gcf, dims, nflds, field, vfill) + + implicit none + + character(len=*), intent(in) :: fname + type(vardefs), intent(in) :: gcf(:) + integer, intent(in) :: dims(3) + integer, intent(in) :: nflds + real, intent(inout) :: field( dims(1) * dims(2) , dims(3) , nflds ) + real, intent(in) :: vfill + + ! internal variables + integer :: max_bytes, lengrib + integer :: ref_time(6) + integer :: lunout, ierr + integer :: fortime, dij, npt + CHARACTER(len=1),allocatable,dimension(:) :: cgrib + real(8) :: tmpfld(size(field,1)) + + ! GRIB2 metadata arrays + integer :: listsec0(2), listsec1(13) + integer :: igdtnum, ipdtnum, idrtnum + integer :: igdtlen, ipdtlen, idrtlen + integer :: jgdt(19), jpdt(15), idrtmpl(16) + integer :: igds(5) + integer :: numcoord, ibmap + real :: coordlist + integer :: ideflist, idefnum + logical*1 :: bmp( dims(1) * dims(2) ) + + integer :: n, lon0, lon1, lat0, lat1, nlay, lyr + integer, dimension(40) :: dep1 + integer, dimension(28) :: dep2 + integer, dimension(:), allocatable :: dep + + npt = dims(1) * dims(2) + + max_bytes = npt * 4 + bmp=.true. + + call getlun(lunout) + call baopenw(lunout, trim(fname), ierr) + if (ierr /= 0) then + write(0, *) 'Error opening grib2 file ', trim(fname) + return + end if + + call retrieve_time( fortime , ref_time ) + + ! Initialize GRIB2 message sections + listsec0(1) = gcf(1)%var_g1 ! Discipline - GRIB Master Table Number (Code Table 0.0) + listsec0(2) = 2 ! GRIB Edition Number (currently 2) + + listsec1(1) = gcf(1)%var_g3 ! Originating Centre (Common Code Table C-1) + listsec1(2) = 0 ! Originating Sub-centre (local table) EMC=4 + listsec1(3) = gcf(1)%var_g2 ! GRIB Master Tables Version Number (Code Table 1.0) + listsec1(4) = 1 ! GRIB Local Tables Version Number (Code Table 1.1) + listsec1(5) = 1 ! Significance of Reference Time (Code Table 1.2) + listsec1(6) = ref_time(1) ! Reference Time - Year -4digits + listsec1(7) = ref_time(2) ! Reference Time - Month + listsec1(8) = ref_time(3) ! Reference Time - Day + listsec1(9) = ref_time(4) ! Reference Time - Hour + listsec1(10) = ref_time(5) ! Reference Time - Minute + listsec1(11) = ref_time(6) ! Reference Time - Second + listsec1(12) = 0 ! Production status of data (Code Table 1.3) + listsec1(13) = 1 ! Type of processed data (Code Table 1.4) + + dep1=(/ 5, 15, 25, 35, 45, 55, 65, 75, 85, 95, 105, 115, 125, 135, 145, 155,& + 165, 175, 185, 195, 205, 215, 226, 241, 267, 309, 374, 467, 594, 757, 960,& + 1204, 1490, 1817, 2184, 2587, 3024, 3489, 3977, 4481 /) + + dep2=(/ 5, 15, 25, 35, 45, 55, 65, 75, 85, 95, 105, 115, 125, 135, 145, 155,& + 165, 175, 185, 195, 205, 215, 226, 241, 267, 309, 374, 467 /) + + if (dims(1) == 1440 .and. dims(2) == 721) then ! 1/4deg rectilinear + dij = 250000 + nlay = 40 + dep = dep1 + end if + + if (dims(1) == 720 .and. dims(2) == 361) then ! 1/2deg rectilinear + dij = 500000 + nlay = 40 + dep = dep1 + end if + + if (dims(1) == 360 .and. dims(2) == 181) then ! 1deg rectilinear + dij = 1000000 + nlay = 40 + dep = dep1 + end if + + if (dims(1) == 72 .and. dims(2) == 36 ) then ! 5deg rectilinear + dij = 5000000 + nlay = 25 + dep = dep2 + end if + + lon0 = 0 + lon1 = 360000000 - dij + lat0 = -90000000 + lat1 = 90000000 + + ! Populate the jgdt array for Template 3.0 (changed parameters to current grib2 files) + jgdt(1) = 6 + jgdt(2) = 0 + jgdt(3) = 0 + jgdt(4) = 0 + jgdt(5) = 0 + jgdt(6) = 0 + jgdt(7) = 0 + jgdt(8) = dims(1) + jgdt(9) = dims(2) + jgdt(10) = 0 + jgdt(11) = -1 + jgdt(12) = lat0 + jgdt(13) = lon0 + jgdt(14) = 48 + jgdt(15) = lat1 + jgdt(16) = lon1 + jgdt(17) = dij + jgdt(18) = dij + jgdt(19) = 64 + + igdtnum=0 + ! Define igds GRIB2 - SECTION 3 + igds(1) = 0 ! Source of grid definition + igds(2) = npt ! Number of grid points + igds(3) = 0 ! Number of octets for each additional grid points definition + igds(4) = 0 ! Interpretation of list for optional points definition + igds(5) = igdtnum ! GRIB2 - CODE TABLE 3.1 + + igdtlen=size(jgdt) + + if (debug) then + write(logunit, *) 'listsec0, listsec1: ', listsec0, listsec1 + write(logunit, *) 'igdtnum, igdtlen: ', igdtnum, igdtlen + write(logunit, *) 'jgdt: ', jgdt + write(logunit, *) 'igds: ', igds + write(logunit, *) 'dij: ', dij + write(logunit, *) 'max_bytes: ', max_bytes + write(logunit, *) 'forcast time: ', fortime + write(logunit, *) 'refference time: ', ref_time + end if + + ideflist=0 + idefnum=0 + + do lyr=1,nlay + + do n=1,nflds + + allocate(cgrib(max_bytes)) + + listsec0(1) = gcf(n)%var_g1 + + call gribcreate(cgrib, max_bytes, listsec0, listsec1, ierr) + if (ierr /= 0) then + write(0, *) 'Error initializing GRIB2 message', ierr + return + end if + + if (debug) write(logunit, *) 'n, nflds, npt, lay: ', n, nflds, npt, lyr, gcf(n)%discription_gb2, gcf(n)%var_fillvalue + + call addgrid(cgrib, max_bytes, igds, jgdt, igdtlen, ideflist, idefnum, ierr) ! there is an internal error here + if (ierr /= 0) then + write(0, *) 'Error adding grid to GRIB2 message', ierr + return + end if + +! Create Section 4 parametrs + ipdtnum=0 + + + jpdt(1)=gcf(n)%var_g5 ! parm number catagory + jpdt(2)=gcf(n)%var_g6 ! parm number + jpdt(3)=2 ! (0-analysis, 1-initialazation, 2-forecast, .. GRIB2 - CODE TABLE 4.3 ) + jpdt(4)=0 ! + jpdt(5)=96 ! Code ON388 Table A- GFS + jpdt(6)=0 ! + jpdt(7)=0 ! + jpdt(8)=1 ! unit (Hour=1) 6hour=11 (ask later) Table 4.4 + jpdt(9)=fortime ! forecast hour + jpdt(10)=gcf(n)%var_g7 ! level ID (1-Ground or Water Surface, 101 mean sea level, 160 depth bellow mean sea level , 168-Ocean Model Layer,...) + jpdt(11)=0 ! scale factor + jpdt(12)=dep(lyr) ! scale value + jpdt(13)=255 + jpdt(14)=0 + jpdt(15)=0 + + if (debug) write(logunit, *) 'ipdtnum=', ipdtnum, ', jpdt= ', jpdt(1:15) + + ipdtlen=size(jpdt) + + numcoord=0 + coordlist=0. ! needed for hybrid vertical coordinate + + ibmap=255 ! Bitmap indicator ( see Code Table 6.0 ) -255 no bitmap + bmp=.true. + + if (trim(gcf(n)%name_gb2) .eq. 'WTMP') then + where ( field(:,lyr,n) .ne. vfill ) field(:,lyr,n) = field(:,lyr,n) + 273.15 + endif + + ! Assign Template 5 + + idrtnum = 0 ! Template 5.2 (Grid Point Data - complex Packing) + + idrtmpl(:)=0 + ! Populate idrtmpl + idrtmpl(1) = 0 ! Reference value (scaled value of the minimum data point) + idrtmpl(2) = 0 ! Binary scale factor (scale by 2^E) + idrtmpl(3) = 3 ! Decimal scale factor (scale by 10^D) + idrtmpl(4) = 0 ! + idrtmpl(5) = 0 ! + idrtmpl(6) = 0 ! + ! Reserved fields + idrtmpl(7:16) = 0 ! Reserved for future use + + idrtlen=size(idrtmpl) + + tmpfld=0 + tmpfld=real(field(:,lyr,n), 8) + + call addfield(cgrib, max_bytes, ipdtnum, jpdt, ipdtlen, coordlist, numcoord, & + idrtnum, idrtmpl, idrtlen, tmpfld, npt, ibmap, bmp, ierr) + if (ierr /=- 0) then + write(0, *) 'Error adding field to GRIB2 message', ierr + return + end if + + call gribend(cgrib, max_bytes, lengrib, ierr) + if (debug) write(logunit, *) 'gribend status=', ierr + if (debug) write(logunit, *) 'length of the final GRIB2 message in octets =', lengrib + call wryte(lunout, lengrib, cgrib) + + deallocate(cgrib) + + end do + end do + + call baclose(lunout, ierr) + + return + +end subroutine write_grib2_3d + + +!-------------------------------------------------------------------------------------- + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!To get a lun used for bacio!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !-------------------------------------------------------------------------------------- + subroutine getlun(lun) + integer, intent(out) :: lun + logical :: is_open + lun = 50 + do + inquire(unit=lun, opened=is_open) + if (.not. is_open) then + return + else + lun = lun + 1 + end if + end do + end subroutine getlun + + + !-------------------------------------------------------------------------------------- + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!Retrieve Time From Input File!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !-------------------------------------------------------------------------------------- + + subroutine retrieve_time(forecast_hour, ref_time) + + implicit none + + integer, intent(out) :: forecast_hour ! Forecast hour as an integer + integer, dimension(6), intent(out) :: ref_time ! Array for GRIB2 reference time: [year, month, day, hour, minute, second] + + integer :: ncid, time_varid, T1_varid, T2_varid + character(len=30) :: units_str + double precision :: T1, T2 + + integer :: ref_year, ref_month, ref_day, ref_hour, ref_min, ref_sec + integer :: year, month, day, hour + double precision :: hours_offset + + call nf90_err(nf90_open(trim(input_file), nf90_nowrite, ncid), 'opening '//input_file) + call nf90_err(nf90_inq_varid(ncid, 'time', time_varid), 'get variable ID: time') + call nf90_err(nf90_get_var(ncid, time_varid, forecast_hour), 'get variable time') + call nf90_err(nf90_get_att(ncid, time_varid, 'units', units_str), 'get attribute: units') + + if (trim(ftype) == 'ocean') then + read(units_str(13:30), '(I4,1X,I2,1X,I2,1X,I2,1X,I2,1X,I2)') & + ref_year, ref_month, ref_day, ref_hour, ref_min, ref_sec + else + read(units_str(12:29), '(I4,1X,I2,1X,I2,1X,I2,1X,I2,1X,I2)') & ! remove it once ice time unit changed + ref_year, ref_month, ref_day, ref_hour, ref_min, ref_sec + forecast_hour=24*forecast_hour ! remove it once ice time unit changed) + end if + + ref_time(1) = ref_year + ref_time(2) = ref_month + ref_time(3) = ref_day + ref_time(4) = ref_hour + ref_time(5) = ref_min + ref_time(6) = ref_sec + + end subroutine retrieve_time + + !---------------------------------------------------------- ! handle netcdf errors !---------------------------------------------------------- @@ -594,4 +1158,5 @@ subroutine nf90_err(ierr, string) stop 99 end if end subroutine nf90_err + end module utils_mod