!/ ------------------------------------------------------------------- / !/ !/ !/ version | date | description !/ 1.0 | 01-nov-2006 | netcdf support in wavewatch III 2.22 !/ 1.1 | 03-dec-2009 | netcdf support in wavewatch III 3.04 !/ 1.2 | 01-sep-2012 | netcdf support in SWAN 40.91 !/ 1.3 | 01-may-2013 | overhauled to simplify maintenance !/ 1.4 | 11-apr-2014 | bug fixes !/ | reduced memory consumption !/ !/ copyright 2009-2014 bmtargoss !/ no unauthorized use without permission. !/ ! 1. purpose : ! ! agioncmd is a module that provides methods to write cf1.5 compliant ! netcdf files containing integrated parameters and spectra. ! Furthermore, some methods are available to read data from cf-1.5 and ! coads compliant files. ! ! For more information: ! http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.5 ! ! 2. method : ! ! 3.a.types ! --name-------------type-------description----------------------- ! mapgrid data on a lon / lat grid ! %longitude r.a x-coordinates ! %latitude r.a y-coordinates ! %nx int number of columns ! %ny int number of rows ! %dx real horizontal cellsize ! %dy real vertical cellsize ! %lunit char likely unit (meter, degrees) ! %alpc real angle of the grid ! %mdc bool multi dimensional coordinates ! ! recordaxe derived type with fields ! %content int.a seconds since 1-jan-1970 / run ! %delta int timestep (s) / 1 ! %ncontent int number of times / runs ! %nstatm bool non-stationary mode (time /run switch) ! %dimid int netcdf dimension id ! %varid int netcdf variable id ! ! pointgrid data along a point dimension ! %points int 1:npoints ! %longitude r.a x coordinates points ! %latitude r.a y coordinates points ! %lunit char likely unit (meter, degrees) ! %npoints ! ! spcgrid spectral grid ! %frequency r.a frequency axe ! %direction r.a. directional axe ! %nfreq int number of frequencies ! %ndir int number of directions ! ! nctable ncttype struct array from nctablemd.f90 ! %standard_name char see remarks a ! %long_name char " ! %units char " ! %nctype integer agnc_short, agnc_byte etc. ! %datamin real minimum expected data value ! %datamax real maximum expected data value ! %varid int nc90_var_id (remark e) ! ! fill_values fvtype struct with fields ! %byte int*1 ! %short int*2 ! %float real ! %double real*8 ! ! 3.a.variables ! ncid int identifier of opened / created file ! pi real 3.1415927 ! r2d real 180. / pi ! d2r real pi / 180. ! ! 4.procedures ! ---- create new netcdf files ---- ! create_ncfile(ncfile, ncid) ! agnc_define_mapgrid ! agnc_define_pntgrid ! agnc_define_spcgrid ! agnc_define_mapvariable ! agnc_define_pntvariables ! agnc_set_mapgrid(ncid, mapgrid) ! agnc_set_pntgrid(ncid, pntgrid) ! agnc_set_spcgrid(ncid, spcgrid) ! agnc_set_recordaxe (ncid, recordaxe) ! agnc_add_variable_attributes ! agnc_add_coordinates_attribute ! ------------- or --------------- ! open_ncfile ! --------- get grids ------------- ! agnc_get_mapgrid ! agnc_get_pntgrid ! agnc_get_spcgrid ! agnc_get_recordaxe ! ------ insert/append data ------ ! agnc_add_mapdata ! agnc_add_pntdata ! agnc_add_spcdata ! agnc_add_density ! ------------- utils ------------- ! coordinate_names_units ! agnc_get_varid_by_name ! close_ncfile ! nccheck ! axe_is_regular ! compute_1d_spectra ! ! 5. used by : ! ! ww3_outnc.ftn (wavewatch) ! swoutp (swan) ! ! 6. exit codes ! ! 7. remarks : ! a. http://cf-pcmdi.llnl.gov/documents/cf-standard-names/ ! b. int16 = round((data - offset) / scale) ! scale = (datamax - datamin) / (2^n -2) ! offset = (datamax - datamin)/2 ! c. officially, netcdf apps should un-/pack data on request ! however, i didn't feel like writing the interface on ! http://www.lahey.com/lookat90.htm, interface to deal with ! the different data storage types. ! d. this module is fortran 95 due to the usage of allocatable ! arrays in derived data type (eg, structs) ! e. set when opening an existing file or creation of a variable in ! a new file ! module agioncmd use netcdf use nctablemd implicit none real, parameter :: AGIONCMD_VERSION = 1.4 real, parameter :: AGNC_DUMMY = NF90_FILL_FLOAT integer(kind=1), parameter :: AGNC_FILL_BYTE = -2**7 integer(kind=2), parameter :: AGNC_FILL_SHORT= -2**15 type mapgrid_type real, dimension(:,:), allocatable :: longitude, latitude character(len=10) :: lunit = 'degrees' integer :: nx = 0, ny = 0, & lon_dimid, lat_dimid, & lon_varid, lat_varid real :: dx, dy, alpc = 0. logical :: mdc = .false. end type mapgrid_type type pntgrid_type real, dimension(:), allocatable :: longitude, latitude integer, dimension(:), allocatable :: ips character(len=10) :: lunit = 'degrees' integer :: pnt_dimid, & lon_varid, lat_varid, & npoints = 0, ips_varid = 0, & xdimlen = 0, ydimlen = 0 end type pntgrid_type type spcgrid_type real, dimension(:), allocatable :: frequency, direction logical :: relative integer :: frq_dimid, frq_varid, & dir_dimid, dir_varid, & nfreq = 0, ndir = 0 end type spcgrid_type type recordaxe_type integer*8, dimension(:), allocatable :: content integer :: delta integer :: ncontent = 0, dimid, varid logical :: nstatm = .true. end type recordaxe_type ! Given the dimension size and/or type, select proper procedure/function interface agnc_get_griddef module procedure agnc_get_griddef_spc, agnc_get_griddef_map, & agnc_get_griddef_pnt end interface agnc_get_griddef interface axe_is_regular module procedure axe_is_regular_integer, & axe_is_regular_integer64, & axe_is_regular_float, & axe_is_regular_float_2d, & axe_is_regular_double end interface axe_is_regular interface agnc_add_mapdata module procedure agnc_add_mapdata_float, agnc_add_mapdata_double end interface agnc_add_mapdata interface agnc_add_spcdata module procedure agnc_add_spcdata_3d, agnc_add_spcdata_2d end interface agnc_add_spcdata interface agnc_add_pntdata module procedure agnc_add_pntdata1d_float, agnc_add_pntdata2d_float, & agnc_add_pntdata3d_float, & agnc_add_pntdata1d_double, agnc_add_pntdata2d_double, & agnc_add_pntdata3d_double end interface agnc_add_pntdata interface seconds_since_epoch module procedure seconds_since_epoch_datevec, seconds_since_epoch_twoint end interface seconds_since_epoch interface datevec module procedure datevec_from_twoint, datevec_from_epoch end interface datevec interface timeindex module procedure timeindex32, timeindex64 end interface timeindex contains ! ! ---- create new netcdf files ---- ! name ! agnc_define_mapgrid ! agnc_define_pntgrid ! agnc_define_spcgrid ! agnc_define_mapvariable ! agnc_define_pntvariables ! agnc_set_mapgrid(ncid, mapgrid) ! agnc_set_pntgrid(ncid, pntgrid) ! agnc_set_spcgrid(ncid, spcgrid) ! agnc_set_recordaxe (ncid, recordaxe) subroutine create_ncfile( ncfile, ncid, recordaxe, mapgrid, pntgrid, spcgrid, Escale, nautical ) character(len=*), intent(in) :: ncfile integer, intent(out) :: ncid type(recordaxe_type), intent(inout) :: recordaxe type(mapgrid_type), intent(inout), optional :: mapgrid type(spcgrid_type), intent(inout), optional :: spcgrid type(pntgrid_type), intent(inout), optional :: pntgrid logical ,intent( in), optional :: Escale logical ,intent( in), optional :: nautical character(len=100) :: history, dconv logical :: do_scale dconv = 'nautical' if ( present(nautical) ) then if ( .not.nautical ) then dconv = 'cartesian' ! nctablemd::nautical defaults to true call set_nctable_convention_nautical(.false.) end if end if do_scale = .true. if ( present(Escale) ) then do_scale = Escale end if write(history,'(a,f3.1)') 'Created with agioncmd version ', AGIONCMD_VERSION call nccheck ( nf90_create( ncfile, NF90_NOCLOBBER + NF90_64BIT_OFFSET, ncid) ) if ( recordaxe%nstatm ) then call nccheck ( nf90_def_dim( ncid, 'time', NF90_UNLIMITED, recordaxe%dimid ) ); call nccheck ( nf90_def_var( ncid, 'time', NF90_INT, recordaxe%dimid, recordaxe%varid ) ) call nccheck ( nf90_put_att( ncid, recordaxe%varid, 'units', 'seconds since 1970-01-01') ) call nccheck ( nf90_put_att( ncid, recordaxe%varid, 'calendar', 'gregorian') ) call nccheck ( nf90_put_att( ncid, recordaxe%varid, 'standard_name', 'time') ) call nccheck ( nf90_put_att( ncid, recordaxe%varid, 'long_name', 'time') ) else call nccheck ( nf90_def_dim( ncid, 'run', NF90_UNLIMITED, recordaxe%dimid ) ); call nccheck ( nf90_def_var( ncid, 'run', NF90_INT, recordaxe%dimid, recordaxe%varid ) ) call nccheck ( nf90_put_att( ncid, recordaxe%varid, 'long_name', 'run number') ) call nccheck ( nf90_put_att( ncid, recordaxe%varid, 'units', '1' ) ) end if ! global attributes call nccheck ( nf90_put_att( ncid, NF90_GLOBAL, 'Conventions', 'CF-1.5') ) call nccheck ( nf90_put_att( ncid, NF90_GLOBAL, 'History', history) ) call nccheck ( nf90_put_att( ncid, NF90_GLOBAL, 'Directional_convention', dconv) ) if ( present(mapgrid) ) then call agnc_define_mapgrid(ncid, mapgrid) if ( present(spcgrid) ) call agnc_define_spcgrid(ncid, spcgrid, do_scale, dconv, mapgrid=mapgrid) else if ( present(pntgrid) ) then call agnc_define_pntgrid(ncid, pntgrid) if ( present(spcgrid) ) call agnc_define_spcgrid(ncid, spcgrid, do_scale, dconv, pntgrid=pntgrid) else call agnc_set_error('Provide a point grid or a map grid to create_ncfile') end if end if ! reserve extra space in the header for adding attributes ! later without reorganizing the file structure. 4kb should be ! enough for most attributes and negligible compared to the ! total file size. call nccheck( nf90_enddef(ncid, h_minfree=4096) ) ! Set the file back to definition mode so variables can be defined later on call nccheck( nf90_redef(ncid) ) end subroutine create_ncfile subroutine coordinate_names_units( lunit, xname, xunit, yname, yunit ) character(len=10), intent( in):: lunit character(len=20), intent( out):: xname, xunit character(len=20), intent( out):: yname, yunit if ( trim(lunit) == 'degrees') then xname = 'longitude' xunit = 'degrees_east' yname = 'latitude' yunit = 'degrees_north' else xname = 'x' xunit = lunit yname = 'y' yunit = lunit end if end subroutine coordinate_names_units subroutine agnc_define_mapgrid( ncid, mapgrid ) integer, intent( in) :: ncid type (mapgrid_type),intent( inout) :: mapgrid character(len=20) :: xname, xunit character(len=20) :: yname, yunit call coordinate_names_units( mapgrid%lunit, xname, xunit, yname, yunit ) if ( mapgrid%mdc ) then call nccheck ( nf90_def_dim( ncid, 'xc', mapgrid%nx, mapgrid%lon_dimid ) ); call nccheck ( nf90_def_dim( ncid, 'yc', mapgrid%ny, mapgrid%lat_dimid ) ); call nccheck ( nf90_def_var( ncid, xname, NF90_FLOAT, & (/ mapgrid%lon_dimid, mapgrid%lat_dimid /), & mapgrid%lon_varid ) ) call nccheck ( nf90_def_var( ncid, yname, NF90_FLOAT, & (/ mapgrid%lon_dimid, mapgrid%lat_dimid /), & mapgrid%lat_varid ) ) call nccheck ( nf90_put_att( ncid, mapgrid%lon_varid, '_FillValue', NF90_FILL_FLOAT) ) call nccheck ( nf90_put_att( ncid, mapgrid%lat_varid, '_FillValue', NF90_FILL_FLOAT) ) else call nccheck ( nf90_def_dim( ncid, xname, mapgrid%nx, mapgrid%lon_dimid ) ); call nccheck ( nf90_def_var( ncid, xname, NF90_FLOAT, mapgrid%lon_dimid, mapgrid%lon_varid ) ); call nccheck ( nf90_def_dim( ncid, yname, mapgrid%ny, mapgrid%lat_dimid ) ); call nccheck ( nf90_def_var( ncid, yname, NF90_FLOAT, mapgrid%lat_dimid, mapgrid%lat_varid ) ); end if call nccheck ( nf90_put_att( ncid, mapgrid%lon_varid, 'units', xunit) ) call nccheck ( nf90_put_att( ncid, mapgrid%lon_varid, 'long_name', xname) ) if ( trim(xname) == 'longitude' ) then call nccheck ( nf90_put_att( ncid, mapgrid%lon_varid, 'standard_name', 'longitude') ) end if call nccheck ( nf90_put_att( ncid, mapgrid%lat_varid, 'units', yunit) ) call nccheck ( nf90_put_att( ncid, mapgrid%lat_varid, 'long_name', yname) ) if ( trim(yname) == 'latitude' ) then call nccheck ( nf90_put_att( ncid, mapgrid%lat_varid, 'standard_name', 'latitude') ) end if end subroutine agnc_define_mapgrid subroutine agnc_define_spcgrid( ncid, spcgrid, do_scale, dconv, mapgrid, pntgrid ) integer, intent(in ) :: ncid type(spcgrid_type), intent(inout) :: spcgrid logical, intent( in) :: do_scale character(len=100), intent( in) :: dconv type(mapgrid_type), intent(inout), optional :: mapgrid type(pntgrid_type), intent(inout), optional :: pntgrid integer :: evarid, ra_dimid, & nf90_var_type logical :: nautical nautical = .false. if ( dconv == 'nautical') nautical = .true. if ( do_scale ) then nf90_var_type = NF90_SHORT else nf90_var_type = NF90_FLOAT end if ! ! frequency ! call nccheck ( nf90_def_dim( ncid, 'frequency', spcgrid%nfreq, spcgrid%frq_dimid ) ); call nccheck ( nf90_def_var( ncid, 'frequency', NF90_FLOAT, spcgrid%frq_dimid, spcgrid%frq_varid ) ); call nccheck ( nf90_put_att( ncid, spcgrid%frq_varid, 'units', 's-1') ) call nccheck ( nf90_put_att( ncid, spcgrid%frq_varid, 'standard_name', 'wave_frequency') ) call nccheck ( nf90_put_att( ncid, spcgrid%frq_varid, 'flow', spcgrid%frequency(1) ) ) call nccheck ( nf90_put_att( ncid, spcgrid%frq_varid, 'fhigh', spcgrid%frequency(spcgrid%nfreq) ) ) call nccheck ( nf90_put_att( ncid, spcgrid%frq_varid, 'msc', spcgrid%nfreq-1 ) ) call agnc_get_recordaxe_ids( ncid, ra_dimid ) if ( spcgrid%ndir > 0 ) then ! ! 2D direction axe ! call nccheck ( nf90_def_dim( ncid, 'direction', spcgrid%ndir, spcgrid%dir_dimid ) ) call nccheck ( nf90_def_var( ncid, 'direction', NF90_FLOAT, spcgrid%dir_dimid, spcgrid%dir_varid ) ); call nccheck ( nf90_put_att( ncid, spcgrid%dir_varid, 'units', 'radians') ) call nccheck ( nf90_put_att( ncid, spcgrid%dir_varid, 'long_name', 'direction') ) if ( nautical ) then call nccheck ( nf90_put_att( ncid, spcgrid%dir_varid, 'standard_name', 'sea_surface_wave_from_direction') ) else call nccheck ( nf90_put_att( ncid, spcgrid%dir_varid, 'standard_name', 'sea_surface_wave_to_direction') ) end if call nccheck ( nf90_put_att( ncid, spcgrid%dir_varid, 'mdc', spcgrid%ndir ) ) ! ! 2D density spectrum ! if ( present(mapgrid) ) then call nccheck ( nf90_def_var( ncid, 'density', nf90_var_type, & (/ spcgrid%dir_dimid, spcgrid%frq_dimid, mapgrid%lon_dimid, mapgrid%lat_dimid, ra_dimid /), evarid ) ) else call nccheck ( nf90_def_var( ncid, 'density', nf90_var_type, & (/ spcgrid%dir_dimid, spcgrid%frq_dimid, pntgrid%pnt_dimid, ra_dimid /), evarid ) ) end if if ( do_scale ) then call nccheck ( nf90_put_att( ncid, evarid, '_FillValue',NF90_FILL_SHORT) ) else call nccheck ( nf90_put_att( ncid, evarid, '_FillValue',NF90_FILL_FLOAT) ) end if call nccheck ( nf90_put_att( ncid, evarid, 'long_name', 'density') ) call nccheck ( nf90_put_att( ncid, evarid, 'units', 'm2 s rad-1') ) call nccheck ( nf90_put_att( ncid, evarid, 'standard_name', & 'sea_surface_wave_directional_variance_spectral_density') ) if ( do_scale) then call nccheck ( nf90_put_att( ncid, evarid, 'note', & 'multiply with scale_density') ) end if if ( spcgrid%relative ) then call nccheck ( nf90_put_att( ncid, evarid, 'relative_to_current', 'true') ) else call nccheck ( nf90_put_att( ncid, evarid, 'relative_to_current', 'false') ) end if ! ! scale_density ! if ( do_scale ) then if ( present(mapgrid) ) then call nccheck ( nf90_def_var( ncid, 'scale_density', NF90_FLOAT, & (/ mapgrid%lon_dimid, mapgrid%lat_dimid, ra_dimid /), evarid ) ) else call nccheck ( nf90_def_var( ncid, 'scale_density', NF90_FLOAT, & (/ pntgrid%pnt_dimid, ra_dimid /), evarid ) ) end if call nccheck ( nf90_put_att( ncid, evarid, 'note', & 'multiply density values with this scale factor') ) call nccheck ( nf90_put_att( ncid, evarid, 'units', '1') ) end if else if ( present(mapgrid) ) then call nccheck ( nf90_def_var( ncid, 'energy_1d', nf90_var_type, & (/ spcgrid%frq_dimid, mapgrid%lon_dimid, mapgrid%lat_dimid, ra_dimid /), evarid ) ) else call nccheck ( nf90_def_var( ncid, 'energy_1d', nf90_var_type, & (/ spcgrid%frq_dimid, pntgrid%pnt_dimid, ra_dimid /), evarid ) ) end if if ( do_scale ) then call nccheck ( nf90_put_att( ncid, evarid, '_FillValue',NF90_FILL_SHORT) ) else call nccheck ( nf90_put_att( ncid, evarid, '_FillValue',NF90_FILL_FLOAT) ) end if call nccheck ( nf90_put_att( ncid, evarid, 'long_name', 'energy') ) call nccheck ( nf90_put_att( ncid, evarid, 'units', 'm2 s') ) call nccheck ( nf90_put_att( ncid, evarid, 'standard_name', 'sea_surface_wave_variance_spectral_density') ) if ( spcgrid%relative ) then call nccheck ( nf90_put_att( ncid, evarid, 'relative_to_current', 'true') ) else call nccheck ( nf90_put_att( ncid, evarid, 'relative_to_current', 'false') ) end if call agnc_add_coordinates_attribute(ncid, evarid) ! ! scale_density ! if ( do_scale ) then call nccheck ( nf90_put_att( ncid, evarid, 'note', & 'multiply with scale_energy_1d') ) if ( present(mapgrid) ) then call nccheck ( nf90_def_var( ncid, 'scale_energy_1d', NF90_FLOAT, & (/ mapgrid%lon_dimid, mapgrid%lat_dimid, ra_dimid /), evarid ) ) else call nccheck ( nf90_def_var( ncid, 'scale_energy_1d', NF90_FLOAT, & (/ pntgrid%pnt_dimid, ra_dimid /), evarid ) ) end if call nccheck ( nf90_put_att( ncid, evarid, 'long_name', & 'Multiply energy_1d values with this scale factor') ) call nccheck ( nf90_put_att( ncid, evarid, 'units', & '1') ) end if ! ! theta_1d ! if ( present(mapgrid) ) then call nccheck ( nf90_def_var( ncid, 'theta_1d', NF90_BYTE, & (/ spcgrid%frq_dimid, mapgrid%lon_dimid, mapgrid%lat_dimid, ra_dimid /), evarid ) ) else call nccheck ( nf90_def_var( ncid, 'theta_1d', NF90_BYTE, & (/ spcgrid%frq_dimid, pntgrid%pnt_dimid, ra_dimid /), evarid ) ) end if call nccheck ( nf90_put_att( ncid, evarid, 'units', 'degree') ) call nccheck ( nf90_put_att( ncid, evarid, 'long_name', 'principal wave direction') ) call nccheck ( nf90_put_att( ncid, evarid, '_FillValue', AGNC_FILL_BYTE) ) call nccheck ( nf90_put_att( ncid, evarid, 'scale_factor', 360. / (2**8 -2)) ) call nccheck ( nf90_put_att( ncid, evarid, 'add_offset', 360. / 2.) ) call agnc_add_coordinates_attribute(ncid, evarid) ! ! spread_1d ! if ( present(mapgrid) ) then call nccheck ( nf90_def_var( ncid, 'spread_1d', NF90_BYTE, & (/ spcgrid%frq_dimid, mapgrid%lon_dimid, mapgrid%lat_dimid, ra_dimid /), evarid ) ) else call nccheck ( nf90_def_var( ncid, 'spread_1d', NF90_BYTE, & (/ spcgrid%frq_dimid, pntgrid%pnt_dimid, ra_dimid /), evarid ) ) end if call nccheck ( nf90_put_att( ncid, evarid, 'units', 'degree') ) call nccheck ( nf90_put_att( ncid, evarid, 'long_name', & 'Longuet-Higgins short-crestedness parameter (s in cos(theta/2)^2s)') ) call nccheck ( nf90_put_att( ncid, evarid, '_FillValue', AGNC_FILL_BYTE) ) call agnc_add_coordinates_attribute(ncid, evarid) end if end subroutine agnc_define_spcgrid subroutine agnc_define_pntgrid( ncid, pntgrid ) integer, intent(in ) :: ncid type (pntgrid_type),intent(inout) :: pntgrid character(len=20) :: xname, xunit character(len=20) :: yname, yunit call coordinate_names_units( pntgrid%lunit, xname, xunit, yname, yunit ) call nccheck ( nf90_def_dim( ncid, 'points', pntgrid%npoints, pntgrid%pnt_dimid ) ) call nccheck ( nf90_def_var( ncid, xname, NF90_FLOAT, pntgrid%pnt_dimid, pntgrid%lon_varid ) ); call nccheck ( nf90_put_att( ncid, pntgrid%lon_varid, 'units', xunit) ) call nccheck ( nf90_put_att( ncid, pntgrid%lon_varid, 'long_name', xname) ) if ( trim(xname) == 'longitude' ) then call nccheck ( nf90_put_att( ncid, pntgrid%lon_varid, 'standard_name', 'longitude') ) end if call nccheck ( nf90_def_var( ncid, yname, NF90_FLOAT, pntgrid%pnt_dimid, pntgrid%lat_varid ) ); call nccheck ( nf90_put_att( ncid, pntgrid%lat_varid, 'units', yunit) ) call nccheck ( nf90_put_att( ncid, pntgrid%lat_varid, 'long_name', yname) ) if ( trim(yname) == 'latitude' ) then call nccheck ( nf90_put_att( ncid, pntgrid%lat_varid, 'standard_name', 'latitude') ) end if ! http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.5/ch08s02.html ! states that compression by gathering consists of having a pointlist with one dimensional ! indices in the two dimensional grid (lon, lat). In the example lon and lat are arrays sized ! nx and ny. This code delibertly breaks that rule because it is also used to write unstructured ! data and along a user defined pointlist. It makes more sense to have a longitude / latitude with the ! size of the pointlist. if ( allocated(pntgrid%ips) ) then call nccheck ( nf90_put_att( ncid, pntgrid%lon_varid, 'dimlen', pntgrid%xdimlen) ) call nccheck ( nf90_put_att( ncid, pntgrid%lat_varid, 'dimlen', pntgrid%ydimlen) ) call nccheck ( nf90_def_var( ncid, 'ips', NF90_INT, pntgrid%pnt_dimid, pntgrid%ips_varid ) ); call nccheck ( nf90_put_att( ncid, pntgrid%ips_varid, 'units', '1') ) call nccheck ( nf90_put_att( ncid, pntgrid%ips_varid, 'computed_as', '(i - 1) * dimlen1 + j') ) call nccheck ( nf90_put_att( ncid, pntgrid%ips_varid, 'description', 'i is index in first spatial dimension, dimlen1 its length') ) call nccheck ( nf90_put_att( ncid, pntgrid%ips_varid, 'long_name', 'one-dimensional index') ) end if end subroutine agnc_define_pntgrid subroutine agnc_define_mapvariable (ncid, varname) integer, intent( in) :: ncid character(len=*), intent( in) :: varname type(nctable_record) :: trecord integer :: ra_dimid, lon_dimid, lat_dimid,& varid call get_nctable_record(varname, trecord) call agnc_get_recordaxe_ids(ncid, ra_dimid) call agnc_get_ll_dimids(ncid, lon_dimid, lat_dimid) call nccheck ( nf90_def_var( ncid, trecord%name , trecord%nctype, & (/ lon_dimid, lat_dimid, ra_dimid /), varid ) ) call agnc_add_variable_attributes( ncid, varid, trecord ) end subroutine agnc_define_mapvariable subroutine agnc_define_pntvariable (ncid, varname) integer, intent( in) :: ncid character(len=*), intent( in) :: varname type(nctable_record) :: trecord integer :: pnt_dimid, ra_dimid, varid call get_nctable_record(varname, trecord) call agnc_get_recordaxe_ids(ncid, ra_dimid) ! both unstructured as spectral grid have the "point" dimension call nccheck ( nf90_inq_dimid(ncid, 'points', pnt_dimid) ) call nccheck ( nf90_def_var( ncid, trecord%name , trecord%nctype, & (/ pnt_dimid, ra_dimid /), varid ) ) call agnc_add_variable_attributes( ncid, varid, trecord ) end subroutine agnc_define_pntvariable subroutine agnc_add_variable_attributes( ncid, varid, trecord ) integer, intent(in) :: ncid, varid type(nctable_record), intent(in) :: trecord real(kind=4) :: add_offset, scale_factor, rng integer :: tmpid character(len=nf90_max_name) :: lonname call nccheck ( nf90_put_att( ncid, varid, 'units', trecord%units) ) if ( trim(trecord%standard_name) /= 'none' ) & call nccheck ( nf90_put_att( ncid, varid, 'standard_name', trecord%standard_name) ) call nccheck ( nf90_put_att( ncid, varid, 'long_name', trecord%long_name) ) call agnc_add_coordinates_attribute(ncid, varid) ! http://www.unidata.ucar.edu/software/netcdf/docs/BestPractices.html ! ! In either the signed or unsigned case, an alternate formula may be used for the add_offset ! and scale_factor packing parameters that reserves a packed value for a special value, such ! as an indicator of missing data. For example, to reserve the minimum packed value (-2n - 1) ! for use as a special value in the case of signed packed values: ! ! scale_factor =(dataMax - dataMin) / (2n - 2) ! add_offset = (dataMax + dataMin) / 2 ! ! Note that the NF90_FILL_ values are one off in comparison with the formulea used here ! ! if datamin equals datamax in an integer kind of way, do not apply scaling rng = (trecord%datamax - trecord%datamin) add_offset = (trecord%datamin + trecord%datamax) * .5 if ( rng /= 0 ) then select case ( trecord%nctype ) case (NF90_BYTE) scale_factor = rng/(2**8-2) case (NF90_SHORT) scale_factor = rng/(2**16-2) case default call agnc_set_error( 'The nctable contains an invalid agnc_TYPE') end select call nccheck ( nf90_put_att( ncid, varid, 'scale_factor', scale_factor) ) call nccheck ( nf90_put_att( ncid, varid, 'add_offset', add_offset) ) end if select case ( trecord%nctype ) case (NF90_BYTE) call nccheck ( nf90_put_att( ncid, varid, '_FillValue', AGNC_FILL_BYTE) ) case (NF90_SHORT) call nccheck ( nf90_put_att( ncid, varid, '_FillValue', AGNC_FILL_SHORT) ) case (NF90_FLOAT) call nccheck ( nf90_put_att( ncid, varid, '_FillValue', NF90_FILL_FLOAT) ) end select end subroutine agnc_add_variable_attributes subroutine agnc_add_coordinates_attribute(ncid, varid) integer, intent(in) :: ncid, varid integer :: tmpid ! CF-1.2+ prescibes a "coordinates" attribute on multi-dimension fields if ( nf90_inq_dimid(ncid, 'xc', tmpid) == nf90_noerr ) then ! multi-dimension field because of the hard coded dimension name. The coordinates ! are either in the longitude/latitude or x,y variables if ( nf90_inq_varid(ncid, 'x', tmpid) == nf90_noerr ) then call nccheck ( nf90_put_att( ncid, varid, 'coordinates', 'x y') ) else call nccheck ( nf90_put_att( ncid, varid, 'coordinates', 'longitude latitude') ) end if end if end subroutine agnc_add_coordinates_attribute ! ------------- or --------------- ! name status ! open_ncfile( ncfile, 'read/write', & ! recordaxe, monthly ) done ! subroutine open_ncfile( ncfile, permission, ncid, recordaxe, monthly ) ! ! recordaxe, is optional. When permission is "read", the information is read from the file. ! When permission is write, a check is made whether the recordaxe and grid is set in ! the dimension variables in the file. character(len=*), intent(in) :: ncfile, permission integer , intent(out) :: ncid type(recordaxe_type) , optional, intent(inout) :: recordaxe logical , optional, intent( in) :: monthly type(recordaxe_type) :: recordaxe_tmp integer :: oldmode if ( permission == "read") then call nccheck( nf90_open(ncfile, nf90_nowrite, ncid) ) else call nccheck( nf90_open(ncfile, nf90_write, ncid) ) end if if ( present(recordaxe)) then call agnc_get_recordaxe(ncid, recordaxe_tmp) ! Check if read time axe is regular and does not only contain _FillValues if ( axe_is_regular(recordaxe_tmp%content, recordaxe_tmp%ncontent, recordaxe_tmp%delta, .true.) ) then recordaxe = recordaxe_tmp else if ( axe_is_regular(recordaxe%content, recordaxe%ncontent, recordaxe%delta, .true.)) then if ( present(monthly) ) then call agnc_set_recordaxe(ncid, recordaxe, monthly) else call agnc_set_recordaxe(ncid, recordaxe) end if else ! If you're writing a program, you end up here if you didn't provide a ! regular grid. Check that dt and nt are set call agnc_set_error( 'Could not read valid record axe from file') end if end if end if end subroutine open_ncfile ! ! --------- set (quasi-) dimension variables ------------- ! name ! agnc_set_recordaxe ! agnc_set_mapgrid ! agnc_set_pntgrid ! agnc_set_spcgrid (ncid, spcgrid) done subroutine agnc_set_recordaxe (ncid, recordaxe, monthly) ! When the optional monthly is set to true, the time axis in the netcdf file ! is predeclared for the month of trecordaxe%content(1) integer, intent ( in) :: ncid type (recordaxe_type), intent (inout) :: recordaxe logical, optional , intent ( in) :: monthly integer :: t1(6), t2(6), i if ( recordaxe%varid < 1 ) call agnc_get_recordaxe_ids(ncid, recordaxe%dimid, recordaxe%varid) ! Some sanity checks if (recordaxe%ncontent == 0) then call agnc_set_error( 'Cannot set empty record axis') end if if ( present(monthly) ) then if (monthly) then t1 = datevec(recordaxe%content(1)) t1(3) = 1 t1(4:6) = (/0,0,0/) t2=t1 t2(2) = t2(2)+1 recordaxe%ncontent = (seconds_since_epoch(t2) - seconds_since_epoch(t1)) / recordaxe%delta deallocate(recordaxe%content) allocate(recordaxe%content(recordaxe%ncontent)) recordaxe%content = (/ (i, i=seconds_since_epoch(t1), seconds_since_epoch(t2) - recordaxe%delta, recordaxe%delta) /) end if end if if (.not. axe_is_regular(recordaxe%content, recordaxe%ncontent, recordaxe%delta) ) then call close_ncfile(ncid) call agnc_set_error( 'Failed to append data to this file. It would result in a irregular record axis'); end if call nccheck ( nf90_put_var(ncid, recordaxe%varid, recordaxe%content) ) end subroutine agnc_set_recordaxe subroutine agnc_set_pntgrid (ncid, pntgrid) integer, intent ( in) :: ncid type (pntgrid_type), intent (inout) :: pntgrid character(len=nf90_max_name) :: name if ( pntgrid%lon_varid < 1 ) then call agnc_get_griddef(ncid, pntgrid) end if call nccheck ( nf90_put_var(ncid, pntgrid%lon_varid, pntgrid%longitude) ) call nccheck ( nf90_put_var(ncid, pntgrid%lat_varid, pntgrid%latitude) ) if ( pntgrid%ips_varid > 0 ) then call nccheck ( nf90_put_var(ncid, pntgrid%ips_varid, pntgrid%ips) ) end if end subroutine agnc_set_pntgrid subroutine agnc_set_mapgrid (ncid, mapgrid) integer, intent ( in) :: ncid type (mapgrid_type), intent (inout) :: mapgrid character(len=nf90_max_name) :: name if ( mapgrid%lon_varid < 1 ) then call agnc_get_griddef(ncid, mapgrid) end if call nccheck ( nf90_put_var(ncid, mapgrid%lon_varid, mapgrid%longitude) ) if ( mapgrid%mdc ) then call nccheck ( nf90_put_var(ncid, mapgrid%lat_varid, mapgrid%latitude) ) else call nccheck ( nf90_put_var(ncid, mapgrid%lat_varid, & reshape(mapgrid%latitude, (/mapgrid%ny/)) )) end if end subroutine agnc_set_mapgrid subroutine agnc_set_spcgrid (ncid, spcgrid) integer, intent ( in) :: ncid type (spcgrid_type), intent (inout) :: spcgrid character(len=nf90_max_name) :: name if ( spcgrid%frq_varid < 1 ) then call agnc_get_griddef(ncid, spcgrid) end if call nccheck ( nf90_put_var(ncid, spcgrid%frq_varid, spcgrid%frequency) ) if (spcgrid%ndir > 0) & call nccheck ( nf90_put_var(ncid, spcgrid%dir_varid, spcgrid%direction) ) end subroutine agnc_set_spcgrid ! --------- get data ------------- ! agnc_get_mapgrid ! agnc_get_pntgrid ! agnc_get_spcgrid ! agnc_get_recordaxe subroutine agnc_get_mapgrid(ncid, mapgrid) ! Structure: ! find dimension id of dimenstion named longitude, lon or x ! get meta information on this dimension ! assume the name of the dimension corresponds with a variable ! read variable ! set results in mapgrid struct ! check that geographic axes are regular integer, intent ( in ) :: ncid type (mapgrid_type), intent (out) :: mapgrid integer :: ndims real :: xlen, ylen real, dimension(:), allocatable :: lattmp call agnc_get_griddef(ncid, mapgrid) call nccheck ( nf90_inquire_variable(ncid, mapgrid%lon_varid, ndims=ndims) ) if ( ndims == 1 ) then allocate(mapgrid%longitude(mapgrid%nx, 1)) allocate(mapgrid%latitude (1,mapgrid%ny)) allocate(lattmp(mapgrid%ny)) call nccheck ( nf90_get_var(ncid, mapgrid%lon_varid, mapgrid%longitude) ) call nccheck ( nf90_get_var(ncid, mapgrid%lat_varid, lattmp) ) mapgrid%latitude(1,:) = lattmp deallocate(lattmp) xlen = mapgrid%longitude(mapgrid%nx,1) - mapgrid%longitude(1,1) ylen = mapgrid%latitude(1,mapgrid%ny) - mapgrid%latitude(1,1) mapgrid%dx = xlen / ( mapgrid%nx - 1 ) mapgrid%dy = ylen / ( mapgrid%ny - 1 ) mapgrid%alpc = 0. else allocate(mapgrid%longitude(mapgrid%nx, mapgrid%ny)) allocate(mapgrid%latitude (mapgrid%nx, mapgrid%ny)) call nccheck ( nf90_get_var(ncid, mapgrid%lat_varid, mapgrid%latitude) ) call nccheck ( nf90_get_var(ncid, mapgrid%lon_varid, mapgrid%longitude) ) xlen = mapgrid%longitude(mapgrid%nx,1) - mapgrid%longitude(1,1) ylen = mapgrid%latitude(mapgrid%nx,1) - mapgrid%latitude(1,1) mapgrid%alpc = atan2(xlen, ylen) ! greatcircle?? mapgrid%dx = ( xlen**2 + ylen**2 )**0.5 / ( mapgrid%nx - 1 ) xlen = mapgrid%longitude(1,mapgrid%ny) - mapgrid%longitude(1,1) ylen = mapgrid%latitude(1,mapgrid%ny) - mapgrid%latitude(1,1) mapgrid%dy = ( xlen**2 + ylen**2)**0.5 / ( mapgrid%ny - 1 ) end if end subroutine agnc_get_mapgrid subroutine agnc_get_pntgrid(ncid, pntgrid) integer, intent ( in ) :: ncid type (pntgrid_type), intent (out) :: pntgrid call agnc_get_griddef(ncid, pntgrid) allocate(pntgrid%longitude(pntgrid%npoints)) call nccheck ( nf90_get_var(ncid, pntgrid%lon_varid, pntgrid%longitude) ); allocate(pntgrid%latitude(pntgrid%npoints)) call nccheck ( nf90_get_var(ncid, pntgrid%lat_varid, pntgrid%latitude) ) if ( pntgrid%ips_varid /= -1 ) then allocate(pntgrid%ips(pntgrid%npoints)) call nccheck ( nf90_get_var(ncid, pntgrid%ips_varid, pntgrid%ips) ) end if end subroutine agnc_get_pntgrid subroutine agnc_get_spcgrid(ncid, spcgrid) integer, intent ( in ) :: ncid type (spcgrid_type), intent (out) :: spcgrid call agnc_get_griddef(ncid, spcgrid) if ( spcgrid%frq_varid > 0 ) then allocate(spcgrid%frequency(spcgrid%nfreq)) call nccheck ( nf90_get_var(ncid, spcgrid%frq_varid, spcgrid%frequency) ) if ( spcgrid%dir_varid > 0 ) then allocate(spcgrid%direction(spcgrid%ndir)) call nccheck ( nf90_get_var(ncid, spcgrid%dir_varid, spcgrid%direction) ) end if end if end subroutine agnc_get_spcgrid subroutine agnc_get_recordaxe(ncid, recordaxe) integer, intent ( in ) :: ncid type (recordaxe_type), intent(out) :: recordaxe character(len=nf90_max_name) :: units, raname, emsg integer*8, allocatable, dimension(:) :: dvalues integer :: factor, dvec(6) ! find time dimid call agnc_get_recordaxe_ids(ncid, recordaxe%dimid, recordaxe%varid, recordaxe%ncontent) if ( recordaxe%varid < 0 ) then call agnc_set_error( 'Could not find record axe variable') end if call nccheck ( nf90_get_att(ncid, recordaxe%varid, 'units', units) ) call nccheck ( nf90_inquire_variable(ncid, recordaxe%varid, raname) ) allocate(recordaxe%content(recordaxe%ncontent)) if (raname /= 'run') then call agnc_scan_time(units, factor, dvec, emsg) if ( factor < 0 ) call agnc_set_error(emsg) allocate(dvalues(recordaxe%ncontent)) call nccheck ( nf90_get_var(ncid, recordaxe%varid, dvalues) ) recordaxe%content = nint(dvalues * dble(factor) + dble(seconds_since_epoch(dvec)),8) deallocate(dvalues) recordaxe%nstatm = .true. else call nccheck ( nf90_get_var(ncid, recordaxe%varid, recordaxe%content) ) recordaxe%nstatm = .false. end if if (recordaxe%ncontent > 1) then recordaxe%delta = recordaxe%content(2) - recordaxe%content(1) if (.not. axe_is_regular(recordaxe%content, recordaxe%ncontent, recordaxe%delta)) then call agnc_set_error( 'record axe does not appear to be regular') end if else ! You'll have set dt yourself when nt is 1 recordaxe%delta = 0 end if end subroutine agnc_get_recordaxe subroutine agnc_get_recordaxe_ids(ncid, dimid, varid, ncontent, raname) integer, intent( in) :: ncid integer, intent(out) :: dimid integer, optional, intent(out) :: varid, ncontent character(len=nf90_max_name), optional, intent(out) :: raname character*6 :: dnames(2), dname integer :: i logical :: dim_found dim_found = .false. dnames = (/ 'time ', 'run ' /) ! Note that the assumption is that the record axe dimension has a correspomding ! variable do i = 1, 2 if ( nf90_inq_dimid(ncid, dnames(i), dimid) == nf90_noerr ) then dname = dnames(i) dim_found = .true. exit end if end do if ( .not. dim_found ) then call agnc_set_error( 'record axe could not be found') end if if (present(ncontent)) then call nccheck ( nf90_inquire_dimension(ncid, dimid, dname, ncontent) ) end if if ( present(varid) ) call agnc_get_varid_by_name(ncid, dname, varid) if ( present(raname)) raname = dname end subroutine agnc_get_recordaxe_ids subroutine agnc_get_ll_dimids(ncid, lon_dimid, lat_dimid) integer , intent( in) :: ncid integer , intent(out) :: lon_dimid, lat_dimid if ( nf90_inq_dimid(ncid, "longitude", lon_dimid) /= nf90_noerr ) then if ( nf90_inq_dimid(ncid, "x", lon_dimid) /= nf90_noerr ) then if ( nf90_inq_dimid(ncid, "xc", lon_dimid) /= nf90_noerr ) then call nccheck ( nf90_inq_dimid(ncid, "lon", lon_dimid) ) end if end if end if if ( nf90_inq_dimid(ncid, "latitude", lat_dimid) /= nf90_noerr ) then if ( nf90_inq_dimid(ncid, "y", lat_dimid) /= nf90_noerr ) then if ( nf90_inq_dimid(ncid, "yc", lat_dimid) /= nf90_noerr ) then call nccheck ( nf90_inq_dimid(ncid, "lat", lat_dimid) ) end if end if end if end subroutine agnc_get_ll_dimids subroutine agnc_get_ll_varids(ncid, lon_varid, lat_varid) integer , intent(in ) :: ncid integer , intent( out) :: lon_varid, lat_varid character(len=nf90_max_name) :: lonnames(3), latnames(3) integer :: n lonnames = (/'longitude', & 'lon ', & 'x '/) latnames = (/'latitude ', & 'lat ', & 'y '/) do n=1,3 call agnc_get_varid_by_name(ncid, trim(lonnames(n)), lon_varid) if ( lon_varid /= -1 ) exit end do if ( lon_varid == -1 ) call agnc_set_error('Could not find horizontal coordinate axe') ! find latitude varid do n=1,3 call agnc_get_varid_by_name(ncid, trim(latnames(n)), lat_varid) if ( lat_varid /= -1 ) exit end do if ( lat_varid == -1 ) call agnc_set_error('Could not find vertical coordinate axe') end subroutine agnc_get_ll_varids subroutine agnc_get_griddef_map(ncid, mapgrid) integer , intent(in ) :: ncid type(mapgrid_type) , intent( out) :: mapgrid call agnc_get_ll_dimids(ncid, mapgrid%lon_dimid, mapgrid%lat_dimid) call nccheck ( nf90_inquire_dimension(ncid, mapgrid%lon_dimid, len=mapgrid%nx) ) call nccheck ( nf90_inquire_dimension(ncid, mapgrid%lat_dimid, len=mapgrid%ny) ) call agnc_get_ll_varids(ncid, mapgrid%lon_varid, mapgrid%lat_varid) end subroutine agnc_get_griddef_map subroutine agnc_get_griddef_pnt(ncid, pntgrid) integer , intent(in ) :: ncid type(pntgrid_type) , intent( out) :: pntgrid call nccheck ( nf90_inq_dimid(ncid, 'points', pntgrid%pnt_dimid) ) call nccheck ( nf90_inquire_dimension(ncid, pntgrid%pnt_dimid, len=pntgrid%npoints) ) call agnc_get_ll_varids(ncid, pntgrid%lon_varid, pntgrid%lat_varid) if ( nf90_inq_varid(ncid, 'ips', pntgrid%ips_varid) /= nf90_noerr ) then pntgrid%ips_varid = -1 else call nccheck ( nf90_get_att(ncid, pntgrid%lon_varid, 'dimlen', pntgrid%xdimlen) ) call nccheck ( nf90_get_att(ncid, pntgrid%lat_varid, 'dimlen', pntgrid%ydimlen) ) end if end subroutine agnc_get_griddef_pnt subroutine agnc_get_griddef_spc(ncid, spcgrid) integer , intent(in ) :: ncid type(spcgrid_type) , intent( out) :: spcgrid ! spectral grid call nccheck ( nf90_inq_dimid(ncid, "frequency", spcgrid%frq_dimid) ) call nccheck ( nf90_inquire_dimension(ncid, spcgrid%frq_dimid, len=spcgrid%nfreq) ) call nccheck ( nf90_inq_varid(ncid, "frequency", spcgrid%frq_varid) ) ! direction only available when frequency is set if ( nf90_inq_dimid(ncid, "direction", spcgrid%dir_dimid) == nf90_noerr ) then call nccheck ( nf90_inquire_dimension(ncid, spcgrid%dir_dimid, len=spcgrid%ndir) ) call nccheck ( nf90_inq_varid(ncid, "direction", spcgrid%dir_varid) ) else spcgrid%dir_dimid = -1 spcgrid%dir_varid = -1 spcgrid%ndir = 0 end if end subroutine agnc_get_griddef_spc subroutine get_scalies_float(ncid, varid, add_offset, scale_factor, fill_value, xtype) integer, intent( in) :: ncid, varid real(kind=4), intent(out) :: add_offset, scale_factor, fill_value integer, optional, intent(out) :: xtype call nccheck ( nf90_inquire_variable(ncid, varid, xtype=xtype)) if ( nf90_get_att(ncid, varid, "_FillValue", fill_value) /= NF90_NOERR ) then select case (xtype) case (NF90_BYTE) fill_value = NF90_FILL_BYTE case (NF90_SHORT) fill_value = NF90_FILL_SHORT case (NF90_INT) fill_value = NF90_FILL_INT case default fill_value = NF90_FILL_FLOAT end select end if if ( nf90_get_att(ncid, varid, "add_offset", add_offset) /= NF90_NOERR ) add_offset = 0. if ( nf90_get_att(ncid, varid, "scale_factor", scale_factor) /= NF90_NOERR ) scale_factor = 1. end subroutine get_scalies_float subroutine get_scalies_double(ncid, varid, add_offset, scale_factor, fill_value, xtype) integer, intent( in) :: ncid, varid real(kind=8), intent(out) :: add_offset, scale_factor, fill_value integer, optional, intent(out) :: xtype call nccheck ( nf90_inquire_variable(ncid, varid, xtype=xtype)) if ( nf90_get_att(ncid, varid, "_FillValue", fill_value) /= NF90_NOERR ) then select case (xtype) case (NF90_BYTE) fill_value = NF90_FILL_BYTE case (NF90_SHORT) fill_value = NF90_FILL_SHORT case (NF90_INT) fill_value = NF90_FILL_INT case default fill_value = NF90_FILL_DOUBLE end select end if if ( nf90_get_att(ncid, varid, "add_offset", add_offset) /= NF90_NOERR ) add_offset = 0. if ( nf90_get_att(ncid, varid, "scale_factor", scale_factor) /= NF90_NOERR ) scale_factor = 1. end subroutine get_scalies_double ! ! ------ insert/append data ------ ! ! agnc_add_mapdata(ncid, varname, time, values) done ! agnc_add_spcdata(ncid, e1, th1, spr1) done subroutine agnc_add_mapdata_float(ncid, varname, ti, values, dummyvalue, skip_range_error) ! if the optional dummyvalue is given, the values are searched for this value ! and replaced with the _FillValue ! ti is the index where the data write should be started integer, intent( in) :: ncid, ti character(len=*), intent( in) :: varname real(kind=4), optional, intent( in) :: dummyvalue logical, optional :: skip_range_error real(kind=4), dimension(:) :: values real(kind=4), dimension(:,:), allocatable :: vloc real(kind=4) :: add_offset, scale_factor, & fill_value type(mapgrid_type) :: mapgrid integer :: varid, xtype, chunksize, & sx, cx, sy, cy, & lxi, lyi, nf90_stat ! integrated parameters, chunks can be fairly large chunksize = 256 call agnc_get_varid_by_name(ncid, varname, varid) call get_scalies_float(ncid, varid, add_offset, scale_factor, fill_value, xtype) call agnc_get_griddef(ncid, mapgrid) ! ! _FillValue ! if ( present(dummyvalue) ) then where (abs(values - dummyvalue) < epsilon(1.)) values = fill_value end if ! ! packed_data_value = nint((unpacked_data_value - add_offset) / scale_factor) ! if ( xtype == NF90_SHORT .or. xtype == NF90_BYTE .or. xtype == NF90_INT ) then where ( abs(values - fill_value) > epsilon(1.) ) & values = nint( (values - add_offset) / scale_factor) end if ! ! Chunked write ! allocate(vloc(chunksize, chunksize)) do sx=1,mapgrid%nx,chunksize cx = minval( (/chunksize, mapgrid%nx-sx+1/) ) do sy=1,mapgrid%ny,chunksize cy = minval( (/chunksize, mapgrid%ny-sy+1/) ) vloc = 0. do lxi=1,cx do lyi=1,cy vloc(lxi,lyi) = values((sx + lxi - 1) + (sy + lyi - 2) * mapgrid%nx ) end do end do nf90_stat = nf90_put_var(ncid, varid, vloc(1:cx,1:cy), (/sx,sy,ti/), (/cx, cy,1/)) ! ! If the write failed, check the values against the limits configured in nctablemd.ftn90 ! if ( nf90_stat /= NF90_NOERR ) then if ( values_in_range_float(ncid, varid, values, & fill_value, scale_factor, add_offset, varname) ) then call agnc_set_error( nf90_strerror(nf90_stat)) else if (.not. present(skip_range_error)) call agnc_set_error( nf90_strerror(nf90_stat)) end if end if end do end do deallocate(vloc) end subroutine agnc_add_mapdata_float subroutine agnc_add_mapdata_double(ncid, varname, ti, values, dummyvalue, skip_range_error) integer, intent( in) :: ncid, ti character(len=*), intent( in) :: varname real(kind=8), optional, intent( in) :: dummyvalue logical, optional :: skip_range_error real(kind=8), dimension(:) :: values real(kind=8), dimension(:,:), allocatable :: vloc real(kind=8) :: add_offset, scale_factor, & fill_value type(mapgrid_type) :: mapgrid integer :: varid, xtype, chunksize, & sx, cx, sy, cy, & lxi, lyi, nf90_stat ! integrated parameters, chunks can be fairly large chunksize = 128 call agnc_get_varid_by_name(ncid, varname, varid) call get_scalies_double(ncid, varid, add_offset, scale_factor, fill_value, xtype) call agnc_get_griddef(ncid, mapgrid) ! ! _FillValue ! if ( present(dummyvalue) ) then where (abs(values - dummyvalue) < epsilon(1.)) values = fill_value end if ! ! packed_data_value = nint((unpacked_data_value - add_offset) / scale_factor) ! if ( xtype == NF90_SHORT .or. xtype == NF90_BYTE .or. xtype == NF90_INT ) then where ( abs(values - fill_value) > epsilon(1.) ) & values = nint( (values - add_offset) / scale_factor) end if ! ! Chunked write ! allocate(vloc(chunksize, chunksize)) do sx=1,mapgrid%nx,chunksize cx = minval( (/chunksize, mapgrid%nx-sx+1/) ) do sy=1,mapgrid%ny,chunksize cy = minval( (/chunksize, mapgrid%ny-sy+1/) ) vloc = 0. do lxi=1,cx do lyi=1,cy vloc(lxi,lyi) = values((sx + lxi - 1) + (sy + lyi - 2) * mapgrid%nx ) end do end do nf90_stat = nf90_put_var(ncid, varid, vloc(1:cx,1:cy), (/sx,sy,ti/), (/cx, cy,1/)) ! ! If the write failed, check the values against the limits configured in nctablemd.ftn90 ! if ( nf90_stat /= NF90_NOERR ) then if ( values_in_range_double(ncid, varid, values, & fill_value, scale_factor, add_offset, varname) ) then call agnc_set_error( nf90_strerror(nf90_stat)) else if (.not. present(skip_range_error)) call agnc_set_error( nf90_strerror(nf90_stat)) end if end if end do end do deallocate(vloc) end subroutine agnc_add_mapdata_double subroutine agnc_add_spcdata_density(ncid, ti, e2d, spc_as_map) integer, intent(in ) :: ncid, ti real(kind=4), dimension(:,:), intent(inout) :: e2d logical, intent(in ) :: spc_as_map real(kind=4), dimension(:,:,:), allocatable :: density real(kind=4), dimension(:,:,:,:), allocatable :: edloc real(kind=4), dimension(:), allocatable :: scale_density integer :: varid_scale_density, & varid_density real(kind=4) :: fv_density integer :: ip, ith, ik, isp, msc, & chunksize, sx, cx, sy, cy, & lxi, lyi type(spcgrid_type) :: spcgrid type(pntgrid_type) :: pntgrid type(mapgrid_type) :: mapgrid type(recordaxe_type) :: recordaxe logical :: do_scale chunksize = 64 do_scale = .false. call agnc_get_griddef(ncid, spcgrid) call agnc_get_recordaxe(ncid, recordaxe) call agnc_collect_spcmeta_2d(ncid, varid_density, varid_scale_density, fv_density) if ( spc_as_map ) then call agnc_get_griddef(ncid, mapgrid) msc = mapgrid%nx * mapgrid%ny else call agnc_get_griddef(ncid, pntgrid) msc = pntgrid%npoints end if ! should I scale and write in a big chuncked loop over msc? ! Allocating the density array creates a new copy for all 2D spectral points in the heap. allocate(density(spcgrid%ndir, spcgrid%nfreq, msc)) ! ! - scale and/or replace energy dummies with _FillValue ! - if scaling requested, write scale factors to netcdf ! if ( varid_scale_density > 0 ) then do_scale = .true. allocate(scale_density(msc)) ! Note the assumption that AGNC_FILL_SHORT is negative scale_density = maxval(e2d, 1) / (2**15-1) do ip=1, msc if (scale_density(ip) < epsilon(1.) ) scale_density(ip) = epsilon(1.) do ik=1, spcgrid%nfreq do ith=1, spcgrid%ndir isp = ith + (ik - 1) * spcgrid%ndir if ( abs(e2d(isp,ip) - AGNC_DUMMY) > epsilon(1.) ) then density(ith,ik,ip) = nint(e2d(isp,ip) / scale_density(ip)) else density(ith,ik,ip) = AGNC_FILL_SHORT end if end do end do end do if ( spc_as_map ) then call nccheck ( nf90_put_var(ncid, varid_scale_density, & reshape(scale_density, (/mapgrid%nx, mapgrid%ny, 1/)), (/1, 1, ti/)) ) else call nccheck ( nf90_put_var(ncid, varid_scale_density, scale_density, (/1, ti/)) ) end if else ! replace unscaled energy dummies with _FillValue do ip=1, msc do ik=1, spcgrid%nfreq do ith=1, spcgrid%ndir isp = ith + (ik - 1) * spcgrid%ndir if ( abs(e2d(isp,ip) - AGNC_DUMMY) > epsilon(1.) ) then density(ith,ik,ip) = e2d(isp,ip) else density(ith,ik,ip) = NF90_FILL_FLOAT end if end do end do end do end if ! ! chunked write to netcdf ! if ( spc_as_map ) then allocate(edloc(spcgrid%ndir, spcgrid%nfreq, chunksize, chunksize)) do sx=1,mapgrid%nx,chunksize cx = minval( (/chunksize, mapgrid%nx-sx+1/) ) do sy=1,mapgrid%ny,chunksize cy = minval( (/chunksize, mapgrid%ny-sy+1/) ) edloc = 0. do lxi=1,cx do lyi=1,cy edloc(:,:,lxi,lyi) = density(:,:, (sx + lxi - 1) + (sy + lyi - 2) * mapgrid%nx ) end do end do call nccheck ( nf90_put_var(ncid, varid_density, edloc(:,:,1:cx,1:cy), & (/1,1,sx,sy,ti/), & (/spcgrid%ndir, spcgrid%nfreq, cx, cy,1/)) ) end do end do deallocate(edloc) else do ip=1,msc, chunksize**2 ! abuse sx and cx for size and end index in pointlist sx = minval( (/chunksize**2, msc-ip+1/) ) cx = ip + sx - 1 call nccheck ( nf90_put_var(ncid, varid_density, density(:,:,ip:cx), & (/1,1,ip,ti/), & (/spcgrid%ndir, spcgrid%nfreq, sx, 1/) ) ) end do end if deallocate(density) if (do_scale) deallocate(scale_density) end subroutine agnc_add_spcdata_density subroutine agnc_add_spcdata_3d(ncid, ti_start, energy, theta, spr, spc_as_map) integer, intent(in ) :: ncid, ti_start real(kind=4), dimension(:,:,:), intent(inout) :: energy, theta, spr logical, intent(in ) :: spc_as_map ! local variables integer :: s(3), ti s = shape(energy) do ti = ti_start, ti_start + s(3) - 1 call agnc_add_spcdata_2d(ncid, ti, reshape(energy(:,:, ti), (/s(1), s(2)/)), & reshape( theta(:,:, ti), (/s(1), s(2)/)), & reshape( spr(:,:, ti), (/s(1), s(2)/)), spc_as_map) end do end subroutine agnc_add_spcdata_3d subroutine agnc_add_spcdata_2d(ncid, ti, energy, theta, spr, spc_as_map) integer, intent(in ) :: ncid, ti real(kind=4), dimension(:,:) :: energy, theta, spr logical, intent(in ) :: spc_as_map ! local variables real(kind=4), dimension(:), allocatable :: scale_energy real(kind=4), dimension(:,:,:), allocatable :: eloc, tloc, sloc integer :: varid_scale_energy, & varid_energy, varid_theta, & varid_spr, msc, nt, & ip, ik, ir, & chunksize, sx, cx, sy, cy, & lxi, lyi real(kind=4) :: sf_theta, sf_spr, & ao_theta, ao_spr, & fv_theta, fv_spr type(spcgrid_type) :: spcgrid type(pntgrid_type) :: pntgrid type(mapgrid_type) :: mapgrid type(recordaxe_type) :: recordaxe logical :: do_scale chunksize = 64 do_scale = .false. call agnc_get_griddef(ncid, spcgrid) call agnc_get_recordaxe(ncid, recordaxe) call agnc_collect_spcmeta(ncid, & varid_scale_energy, varid_energy, & varid_theta, ao_theta, sf_theta, fv_theta,& varid_spr, ao_spr, sf_spr, fv_spr) if ( spc_as_map ) then call agnc_get_griddef(ncid, mapgrid) msc = mapgrid%nx * mapgrid%ny else call agnc_get_griddef(ncid, pntgrid) msc = pntgrid%npoints end if if ( msc /= size(energy,2) ) then write(6,*) '#points of this run and the netcdf file is different' STOP 2 end if if ( varid_scale_energy > 0 ) then do_scale = .true. allocate(scale_energy(msc)) ! Note the assumption that AGNC_DUMMY is negative (or zero) scale_energy = maxval(energy, 1) / (2**15-1) do ip=1, msc ! scale or replace energy dummies with _FillValue if (scale_energy(ip) < epsilon(1.) ) scale_energy(ip) = epsilon(1.) do ik=1, spcgrid%nfreq if ( abs(energy(ik,ip) - AGNC_DUMMY) > epsilon(1.)) then energy(ik,ip) = nint(energy(ik,ip) / scale_energy(ip)) else energy(ik,ip) = AGNC_FILL_SHORT end if end do end do if ( spc_as_map ) then ! make a clone of scale_energy. Reshape is acceptable..... call nccheck ( nf90_put_var(ncid, varid_scale_energy, & reshape(scale_energy, (/mapgrid%nx, mapgrid%ny, 1/)), (/1, 1, ti/)) ) else call nccheck ( nf90_put_var(ncid, varid_scale_energy, scale_energy, (/1, ti/)) ) end if else ! replace unscaled energy dummies with _FillValue do ip=1, msc do ik=1, spcgrid%nfreq if ( abs(energy(ik,ip) - AGNC_DUMMY) < 2*epsilon(1.)) & energy(ik,ip) = NF90_FILL_FLOAT end do end do end if ! replace direction and spread dummies with _FillValue do ip=1, msc do ik=1, spcgrid%nfreq if ( abs(theta(ik,ip) - AGNC_DUMMY) > epsilon(1.) ) then theta(ik,ip) = nint((theta(ik,ip) - ao_theta) / sf_theta) else theta(ik,ip) = AGNC_FILL_BYTE end if if ( abs(spr(ik,ip) - AGNC_DUMMY) > epsilon(1.) ) then spr(ik,ip) = nint((spr(ik,ip) - ao_spr) / sf_spr) else spr(ik,ip) = AGNC_FILL_BYTE end if end do end do ! ! chunked write to netcdf ! if ( spc_as_map ) then allocate(eloc(spcgrid%nfreq, chunksize, chunksize)) allocate(tloc(spcgrid%nfreq, chunksize, chunksize)) allocate(sloc(spcgrid%nfreq, chunksize, chunksize)) do sx=1,mapgrid%nx,chunksize cx = minval( (/chunksize, mapgrid%nx-sx+1/) ) do sy=1,mapgrid%ny,chunksize cy = minval( (/chunksize, mapgrid%ny-sy+1/) ) eloc = 0. tloc = 0. sloc = 0. do lxi=1,cx do lyi=1,cy eloc(:,lxi,lyi) = energy(:, (sx + lxi - 1) + (sy + lyi - 2) * mapgrid%nx) tloc(:,lxi,lyi) = theta (:, (sx + lxi - 1) + (sy + lyi - 2) * mapgrid%nx) sloc(:,lxi,lyi) = spr (:, (sx + lxi - 1) + (sy + lyi - 2) * mapgrid%nx) end do end do call nccheck ( nf90_put_var(ncid, varid_energy, eloc(:,1:cx,1:cy), & (/1,sx,sy,ti/), & (/spcgrid%nfreq, cx, cy,1/)) ) call nccheck ( nf90_put_var(ncid, varid_theta , tloc(:,1:cx,1:cy), & (/1,sx,sy,ti/), & (/spcgrid%nfreq, cx, cy,1/)) ) call nccheck ( nf90_put_var(ncid, varid_spr , sloc(:,1:cx,1:cy), & (/1,sx,sy,ti/), & (/spcgrid%nfreq, cx, cy,1/)) ) end do end do deallocate(eloc) deallocate(tloc) deallocate(sloc) else do ip=1,msc, chunksize**2 ! abuse sx and cx for size and end index in pointlist sx = minval( (/chunksize**2, msc-ip+1/) ) cx = ip + sx - 1 call nccheck ( nf90_put_var(ncid, varid_energy, energy(:,ip:cx), & (/1,ip,ti/), & (/spcgrid%nfreq, sx, 1/) )) call nccheck ( nf90_put_var(ncid, varid_theta, theta(:,ip:cx), & (/1,ip,ti/), & (/spcgrid%nfreq, sx, 1/) )) call nccheck ( nf90_put_var(ncid, varid_spr, spr(:,ip:cx), & (/1,ip,ti/), & (/spcgrid%nfreq, sx, 1/) )) end do end if if ( do_scale ) deallocate(scale_energy) end subroutine agnc_add_spcdata_2d subroutine agnc_add_pntdata1d_float(ncid, varname, ti, values, dummyvalue, skip_range_error) ! if the optional dummyvalue is given, the values are searched for this value ! and replaced with the _FillValue ! ti is the index where the data write should be started integer, intent( in) :: ncid, ti character(len=*), intent( in) :: varname real(kind=4), optional, intent( in) :: dummyvalue logical, optional :: skip_range_error real(kind=4), dimension(:) :: values real(kind=4) :: add_offset, scale_factor, & fill_value integer :: varid, nf90_stat, xtype call agnc_get_varid_by_name(ncid, varname, varid) call get_scalies_float(ncid, varid, add_offset, scale_factor, fill_value, xtype) if ( present(dummyvalue) ) then where (abs(values - dummyvalue) < epsilon(1.)) values = fill_value end if ! packed_data_value = nint((unpacked_data_value - add_offset) / scale_factor) if ( xtype == NF90_SHORT .or. xtype == NF90_BYTE .or. xtype == NF90_INT ) then where ( abs(values - fill_value) > epsilon(1.) ) & values = nint( (values - add_offset) / scale_factor) end if nf90_stat = nf90_put_var(ncid, varid, values, (/1, ti/)) if ( nf90_stat /= NF90_NOERR ) then if ( values_in_range_float(ncid, varid, reshape(values, (/size(values)/)), & fill_value, scale_factor, add_offset, varname) ) then call agnc_set_error( nf90_strerror(nf90_stat)) else if (.not. present(skip_range_error)) call agnc_set_error( nf90_strerror(nf90_stat)) end if end if end subroutine agnc_add_pntdata1d_float subroutine agnc_add_pntdata2d_float(ncid, varname, ti, values, dummyvalue, skip_range_error) ! if the optional dummyvalue is given, the values are searched for this value ! and replaced with the _FillValue ! ti is the index where the data write should be started integer, intent( in) :: ncid, ti character(len=*), intent( in) :: varname real(kind=4), optional, intent( in) :: dummyvalue logical, optional :: skip_range_error real(kind=4), dimension(:,:) :: values real(kind=4) :: add_offset, scale_factor, & fill_value integer :: varid, nf90_stat, xtype call agnc_get_varid_by_name(ncid, varname, varid) call get_scalies_float(ncid, varid, add_offset, scale_factor, fill_value, xtype) if ( present(dummyvalue) ) then where (abs(values - dummyvalue) < epsilon(1.)) values = fill_value end if ! packed_data_value = nint((unpacked_data_value - add_offset) / scale_factor) if ( xtype == NF90_SHORT .or. xtype == NF90_BYTE .or. xtype == NF90_INT ) then where ( abs(values - fill_value) > epsilon(1.) ) & values = nint( (values - add_offset) / scale_factor) nf90_stat = nf90_put_var(ncid, varid, values, (/1, ti/)) end if if ( nf90_stat /= NF90_NOERR ) then if ( values_in_range_float(ncid, varid, reshape(values, (/size(values)/)), & fill_value, scale_factor, add_offset, varname) ) then call agnc_set_error( nf90_strerror(nf90_stat)) else if (.not. present(skip_range_error)) call agnc_set_error( nf90_strerror(nf90_stat)) end if end if end subroutine agnc_add_pntdata2d_float subroutine agnc_add_pntdata3d_float(ncid, varname, ti, values, dummyvalue, skip_range_error) ! if the optional dummyvalue is given, the values are searched for this value ! and replaced with the _FillValue ! ti is the index where the data write should be started integer, intent( in) :: ncid, ti character(len=*), intent( in) :: varname real(kind=4), optional, intent( in) :: dummyvalue logical, optional :: skip_range_error real(kind=4), dimension(:,:,:) :: values real(kind=4) :: add_offset, scale_factor, & fill_value integer :: varid, nf90_stat, xtype call agnc_get_varid_by_name(ncid, varname, varid) call get_scalies_float(ncid, varid, add_offset, scale_factor, fill_value, xtype) if ( present(dummyvalue) ) then where (abs(values - dummyvalue) < epsilon(1.)) values = fill_value end if ! packed_data_value = nint((unpacked_data_value - add_offset) / scale_factor) if ( xtype == NF90_SHORT .or. xtype == NF90_BYTE .or. xtype == NF90_INT ) then where ( abs(values - fill_value) > epsilon(1.) ) & values = nint( (values - add_offset) / scale_factor) end if nf90_stat = nf90_put_var(ncid, varid, values, (/1, ti/)) if ( nf90_stat /= NF90_NOERR ) then if ( values_in_range_float(ncid, varid, reshape(values, (/size(values)/)), & fill_value, scale_factor, add_offset, varname) ) then call agnc_set_error( nf90_strerror(nf90_stat)) else if (.not. present(skip_range_error)) call agnc_set_error( nf90_strerror(nf90_stat)) end if end if end subroutine agnc_add_pntdata3d_float subroutine agnc_add_pntdata1d_double(ncid, varname, ti, values, dummyvalue, skip_range_error) ! if the optional dummyvalue is given, the values are searched for this value ! and replaced with the _FillValue ! ti is the index where the data write should be started integer, intent( in) :: ncid, ti character(len=*), intent( in) :: varname real(kind=8), optional, intent( in) :: dummyvalue logical, optional :: skip_range_error real(kind=8), dimension(:) :: values real(kind=8) :: add_offset, scale_factor, & fill_value integer :: varid, nf90_stat, xtype call agnc_get_varid_by_name(ncid, varname, varid) call get_scalies_double(ncid, varid, add_offset, scale_factor, fill_value, xtype) if ( present(dummyvalue) ) then where (abs(values - dummyvalue) < epsilon(1.)) values = fill_value end if ! packed_data_value = nint((unpacked_data_value - add_offset) / scale_factor) if ( xtype == NF90_SHORT .or. xtype == NF90_BYTE .or. xtype == NF90_INT ) then where ( abs(values - fill_value) > epsilon(1.) ) & values = nint( (values - add_offset) / scale_factor) end if nf90_stat = nf90_put_var(ncid, varid, values, (/1, ti/)) if ( nf90_stat /= NF90_NOERR ) then if ( values_in_range_double(ncid, varid, reshape(values, (/size(values)/)), & fill_value, scale_factor, add_offset, varname) ) then call agnc_set_error( nf90_strerror(nf90_stat)) else if (.not. present(skip_range_error)) call agnc_set_error( nf90_strerror(nf90_stat)) end if end if end subroutine agnc_add_pntdata1d_double subroutine agnc_add_pntdata2d_double(ncid, varname, ti, values, dummyvalue, skip_range_error) ! if the optional dummyvalue is given, the values are searched for this value ! and replaced with the _FillValue ! ti is the index where the data write should be started integer, intent( in) :: ncid, ti character(len=*), intent( in) :: varname real(kind=8), optional, intent( in) :: dummyvalue logical, optional :: skip_range_error real(kind=8), dimension(:,:) :: values real(kind=8) :: add_offset, scale_factor, & fill_value integer :: varid, nf90_stat, xtype call agnc_get_varid_by_name(ncid, varname, varid) call get_scalies_double(ncid, varid, add_offset, scale_factor, fill_value, xtype) if ( present(dummyvalue) ) then where (abs(values - dummyvalue) < epsilon(1.)) values = fill_value end if ! packed_data_value = nint((unpacked_data_value - add_offset) / scale_factor) if ( xtype == NF90_SHORT .or. xtype == NF90_BYTE .or. xtype == NF90_INT ) then where ( abs(values - fill_value) > epsilon(1.) ) & values = nint( (values - add_offset) / scale_factor) end if nf90_stat = nf90_put_var(ncid, varid, values, (/1, ti/)) if ( nf90_stat /= NF90_NOERR ) then if ( values_in_range_double(ncid, varid, reshape(values, (/size(values)/)), & fill_value, scale_factor, add_offset, varname) ) then call agnc_set_error( nf90_strerror(nf90_stat)) else if (.not. present(skip_range_error)) call agnc_set_error( nf90_strerror(nf90_stat)) end if end if end subroutine agnc_add_pntdata2d_double subroutine agnc_add_pntdata3d_double(ncid, varname, ti, values, dummyvalue, skip_range_error) ! if the optional dummyvalue is given, the values are searched for this value ! and replaced with the _FillValue ! ti is the index where the data write should be started integer, intent( in) :: ncid, ti character(len=*), intent( in) :: varname real(kind=8), optional, intent( in) :: dummyvalue logical, optional :: skip_range_error real(kind=8), dimension(:,:,:) :: values real(kind=8) :: add_offset, scale_factor, & fill_value integer :: varid, nf90_stat, xtype call agnc_get_varid_by_name(ncid, varname, varid) call get_scalies_double(ncid, varid, add_offset, scale_factor, fill_value, xtype) if ( present(dummyvalue) ) then where (abs(values - dummyvalue) < epsilon(1.)) values = fill_value end if ! packed_data_value = nint((unpacked_data_value - add_offset) / scale_factor) if ( xtype == NF90_SHORT .or. xtype == NF90_BYTE .or. xtype == NF90_INT ) then where ( abs(values - fill_value) > epsilon(1.) ) & values = nint( (values - add_offset) / scale_factor) end if nf90_stat = nf90_put_var(ncid, varid, values, (/1, ti/)) if ( nf90_stat /= NF90_NOERR ) then if ( values_in_range_double(ncid, varid, reshape(values, (/size(values)/)), & fill_value, scale_factor, add_offset, varname) ) then call agnc_set_error( nf90_strerror(nf90_stat)) else if (.not. present(skip_range_error)) call agnc_set_error( nf90_strerror(nf90_stat)) end if end if end subroutine agnc_add_pntdata3d_double ! ! ------------- utils ------------- ! name status ! agnc_get_varid_by_name(ncid, searchfor, varid) done ! close_ncfile(ncid) done ! nccheck done ! axe_is_regular(x, nx, dx) done ! compute_1d_spectra(energy, e1, th1, sp1, undef) done subroutine agnc_get_varid_by_name(ncid, searchfor, varid) ! search for name in: ! 1) name of the variable ! 2) standard_name of the variable ! 3) short_name of the variable integer, intent(in) :: ncid character(len=*), intent(in) :: searchfor integer, intent(out) :: varid integer :: ndims,nvars character(len=NF90_MAX_NAME) :: name call nccheck ( nf90_inquire(ncid, ndims, nvars) ); do varid=1,nvars call nccheck (nf90_inquire_variable(ncid, varid, name) ) if ( trim(name) == trim(searchfor) ) return if ( nf90_get_att(ncid, varid, "standard_name", name) == NF90_NOERR ) then if ( trim(name) == trim(searchfor) ) return end if if (nf90_get_att(ncid, varid, "short_name", name) == NF90_NOERR ) then if ( trim(name) == trim(searchfor) ) return end if if (nf90_get_att(ncid, varid, "long_name", name) == NF90_NOERR ) then if ( trim(name) == trim(searchfor) ) return end if end do if ( varid > nvars ) then varid = -1; end if end subroutine agnc_get_varid_by_name subroutine close_ncfile(ncid) integer , intent(in) :: ncid call nccheck( nf90_close(ncid) ) end subroutine close_ncfile function axe_is_regular_integer(x, nx, dx, fill_check) result (isregular) integer, intent(in) :: x(:),dx,nx logical, intent(in), optional :: fill_check logical :: isregular integer, dimension(:), allocatable :: dx2 isregular = .false. if ( nx > 0 ) then if ( present(fill_check) ) then if ( fill_check .and. all( (x - NF90_FILL_INT) == 0) ) return end if if ( nx > 1 ) then allocate(dx2(nx-1)) dx2 = x(2:nx) - x(1:nx-1) if ( any( (dx2 - dx) /= 0 )) then deallocate(dx2) return end if deallocate(dx2) end if isregular = .true. end if end function axe_is_regular_integer function axe_is_regular_integer64(x, nx, dx, fill_check) result (isregular) integer*8, intent(in) :: x(:) integer, intent(in) :: dx,nx logical, intent(in), optional :: fill_check logical :: isregular integer, dimension(:), allocatable :: dx2 isregular = .false. if ( nx > 0 ) then if ( present(fill_check) ) then if ( fill_check .and. all( (x - NF90_FILL_INT) == 0) ) return end if if ( nx > 1 ) then allocate(dx2(nx-1)) dx2 = x(2:nx) - x(1:nx-1) if ( any( (dx2 - dx) /= 0 )) then deallocate(dx2) return end if deallocate(dx2) end if isregular = .true. end if end function axe_is_regular_integer64 function axe_is_regular_float(x, nx, dx, fill_check) result (isregular) real(kind=4),intent(in) :: x(:),dx integer, intent(in) :: nx logical, intent(in), optional :: fill_check logical :: isregular real(kind=4),dimension(:), allocatable :: dx2 isregular = .false. if ( nx > 0 ) then if ( present(fill_check) ) then if ( fill_check .and. all( (x - NF90_FILL_INT) == 0) ) return end if if ( nx > 1 ) then allocate(dx2(nx-1)) dx2 = x(2:nx) - x(1:nx-1) if ( any( abs(dx2 - dx) > 1e-4 )) then deallocate(dx2) return end if deallocate(dx2) end if isregular = .true. end if end function axe_is_regular_float function axe_is_regular_float_2d(x, nx, dx, fill_check) result (isregular) real(kind=4),intent(in) :: x(:,:),dx integer, intent(in) :: nx logical, intent(in), optional :: fill_check logical :: isregular real(kind=4) :: dx2 ! The reason for 2D meshes is that we can use the SWAN rotated grid feature ! ! I have never used this before so I would really know how to deal with it ! as I'm not sure how x0, x1, and alpha are implemented. isregular = .true. end function axe_is_regular_float_2d function axe_is_regular_double(x, nx, dx, fill_check) result(isregular) real(kind=8),intent(in) :: x(:),dx integer, intent(in) :: nx logical, intent(in), optional :: fill_check logical :: isregular real(kind=8),dimension(:), allocatable :: dx2 isregular = .false. if ( nx > 0 ) then if ( present(fill_check) ) then if ( fill_check .and. all( (x - NF90_FILL_INT) == 0) ) return end if if ( nx > 1 ) then allocate(dx2(nx-1)) dx2 = x(2:nx) - x(1:nx-1) if ( any( abs(dx2 - dx) > 1e-4 )) then deallocate(dx2) return end if deallocate(dx2) end if isregular = .true. end if end function axe_is_regular_double subroutine nccheck(status) integer, intent ( in ) :: status if (status /= nf90_noerr) then write(6,*) trim(nf90_strerror(status)) STOP 1 end if end subroutine nccheck subroutine agnc_set_error(msg) character(len=*), intent (in) :: msg write(6,*) trim(msg) STOP 2 end subroutine agnc_set_error function values_in_range_float(ncid, varid, values, fill_value, scale_factor, add_offset, varname) result (inrange) integer, intent( in) :: ncid, varid real(kind=4), dimension(:), intent( in) :: values real(kind=4), intent( in) :: fill_value, scale_factor, add_offset character(len=*), intent( in) :: varname ! local integer :: xtype, ind1d, fillv_int, value logical :: inrange real(kind=4) :: bmax, bmin inrange = .true. ! Get variable meta information call nccheck ( nf90_inquire_variable(ncid, varid, xtype = xtype)) ! Compute theoretic min and max values given the type and offset / scale select case (xtype) case (NF90_BYTE) bmin = -127. bmax = 127. case (NF90_SHORT) bmin = -32767. bmax = 32767. case (NF90_INT) bmin = -2147483647. bmax = 2147483647. case default return end select if ( (fill_value - bmin - 1) > epsilon(1.) ) then ! This is not the usual NF90_FILL_XXXX value, so compute the integer value of the fill value fillv_int = nint( (fill_value - add_offset) / scale_factor ) else fillv_int = nint(bmin - 1) end if do ind1d=1,size(values) value = nint(values(ind1d)) if ( value /= fillv_int .and. value .lt. bmin .or. value .gt. bmax ) then write(6, '("The value in 1d index ", I12, " lies outside the boundaries provided in nctablemd.f90")') ind1d write(6, *) "fill value int / float: ", fill_value, fill_value * scale_factor + add_offset write(6, *) "Theoretical min / max: ", & bmin * scale_factor + add_offset, & bmax * scale_factor + add_offset write(6, *) "Actual value int / float: ", value, " / ", & values(ind1d) * scale_factor + add_offset inrange = .false. return end if end do end function values_in_range_float function values_in_range_double(ncid, varid, values, fill_value, scale_factor, add_offset, varname) result (inrange) integer, intent( in) :: ncid, varid real(kind=8), dimension(:), intent( in) :: values real(kind=8), intent( in) :: fill_value, scale_factor, add_offset character(len=*), intent( in) :: varname ! local integer :: xtype, ind1d, fillv_int, value logical :: inrange real(kind=8) :: bmax, bmin inrange = .true. ! Get variable meta information call nccheck ( nf90_inquire_variable(ncid, varid, xtype = xtype)) ! Compute theoretic min and max values given the type and offset / scale select case (xtype) case (NF90_BYTE) bmin = -127. bmax = 127. case (NF90_SHORT) bmin = -32767. bmax = 32767. case (NF90_INT) bmin = -2147483647. bmax = 2147483647. case default return end select if ( (fill_value - bmin - 1) > epsilon(1.) ) then ! This is not the usual NF90_FILL_XXXX value, so compute the integer value of the fill value fillv_int = nint( (fill_value - add_offset) / scale_factor ) else fillv_int = nint(bmin - 1) end if do ind1d=1,size(values) value = nint(values(ind1d)) if ( value /= fillv_int .and. value .lt. bmin .or. value .gt. bmax ) then write(6, '("The value in 1d index ", I12, " lies outside the boundaries provided in nctablemd.f90")') ind1d write(6, *) "fill value : ", fill_value write(6, *) "Theoretical min / max: ", & bmin * scale_factor + add_offset, & bmax * scale_factor + add_offset write(6, *) "Actual values int / float: ", value, " / ", & values(ind1d) * scale_factor + add_offset inrange = .false. return end if end do end function values_in_range_double function timeindex32(tarr, t) result (ti) integer, dimension(:), intent( in) :: tarr integer, intent( in) :: t integer :: ti if ( t >= minval(tarr) .and. t <= maxval(tarr) ) then ! somewhere in the existing time axe .and. (t <= max(tarr,1)) ti = minloc( abs(tarr-t), 1) else if ( t < minval(tarr) ) then ! smaller then the existing time axe. This should not be allowed ti = -1 else ! larger. Extend time axe by one ti = 0 end if end if end function timeindex32 function timeindex64(tarr, t) result (ti) integer*8, dimension(:), intent( in) :: tarr integer, intent( in) :: t integer :: ti if ( t >= minval(tarr) .and. t <= maxval(tarr) ) then ! somewhere in the existing time axe .and. (t <= max(tarr,1)) ti = minloc( abs(tarr-t), 1) else if ( t < minval(tarr) ) then ! smaller then the existing time axe. This should not be allowed ti = -1 else ! larger. Extend time axe by one ti = 0 end if end if end function timeindex64 subroutine agnc_collect_spcmeta(ncid, varid_scale_energy, varid_energy, & varid_theta, ao_theta, sf_theta, fv_theta, & varid_spr, ao_spr, sf_spr, fv_spr) integer, intent( in) :: ncid integer, intent(out) :: varid_scale_energy, & varid_energy, varid_theta, & varid_spr real(kind=4), intent(out) :: sf_theta, sf_spr, & ao_theta, ao_spr, & fv_theta, fv_spr call agnc_get_varid_by_name(ncid, "scale_energy_1d", varid_scale_energy) call agnc_get_varid_by_name(ncid, "energy_1d", varid_energy) call agnc_get_varid_by_name(ncid, "theta_1d", varid_theta) call agnc_get_varid_by_name(ncid, "spread_1d", varid_spr) call get_scalies_float(ncid, varid_theta, ao_theta, sf_theta, fv_theta) call get_scalies_float(ncid, varid_spr, ao_spr, sf_spr, fv_spr) end subroutine agnc_collect_spcmeta subroutine agnc_collect_spcmeta_2d(ncid, varid_density, varid_scale_density, fv_density) integer, intent( in) :: ncid integer, intent(out) :: varid_scale_density, & varid_density real(kind=4) :: ao_density, sf_density real(kind=4), intent(out) :: fv_density call agnc_get_varid_by_name(ncid, "density", varid_density) call agnc_get_varid_by_name(ncid, "scale_density", varid_scale_density) call get_scalies_float(ncid, varid_density, ao_density, sf_density, fv_density) end subroutine agnc_collect_spcmeta_2d function datevec_from_epoch( t_in ) result (datevec) integer*8, intent( in) :: t_in integer :: t, year, month ,day, s integer :: datevec(6) ! number of days per month integer, save :: ndpm(12) data ndpm / 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 / t = t_in datevec = (/0, 0, 0, 0, 0, 0/) datevec(6) = mod(t, 60) t=t-datevec(6) datevec(5) = mod(t,3600)/60 t = t - datevec(5)*60 datevec(4) = mod(t,86400)/3600 t = t-datevec(4)*3600 do year = 1970, 2038 if ( (mod(year,400) == 0) .or. ( (mod(year,4) == 0) .and. (mod(year,100) /= 0 )) ) then s = 366*86400 else s = 365*86400 end if if (t < s) then exit else t = t-s end if end do datevec(1) = year do month = 1, 12 if ( ((mod(year,400) == 0) .or. ( (mod(year,4) == 0) & .and. (mod(year,100) /= 0) )) .and. month == 2 ) then s = 29 * 86400 else s = ndpm(month) * 86400 end if if (t < s) then exit else t = t-s end if end do datevec(2) = month datevec(3) = t/86400+1 end function datevec_from_epoch function datevec_from_twoint(date, time) result (datevec) integer, intent( in) :: date, time integer :: datevec(6) datevec(1) = date / 10000 datevec(2) = mod(date,10000) / 100 datevec(3) = mod(date,100) datevec(4) = time / 10000 datevec(5) = mod(time,10000) / 100 datevec(6) = mod(time,100) end function datevec_from_twoint function seconds_since_epoch_twoint ( date, time ) result (seconds_since_epoch) integer, intent( in) :: date, time integer*8 :: seconds_since_epoch seconds_since_epoch = seconds_since_epoch_datevec( datevec_from_twoint(date, time) ) end function seconds_since_epoch_twoint function seconds_since_epoch_datevec ( datevec ) result (seconds_since_epoch) ! convert date time in (int yyyymmdd, hhmmss) format to seconds since 1-jan-1970 integer :: datevec(6) integer*8 :: seconds_since_epoch integer :: year, month, y1, y2, f ! number of days per month integer, save :: ndpm(12) data ndpm / 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 / seconds_since_epoch = 0 y1 = 1970 y2 = 1970 f = 1 if (datevec(2) > 12) then datevec(1) = datevec(1) + datevec(2)/12 datevec(2) = mod(datevec(2),12) end if if (datevec(1) < 1970) then y1 = datevec(1)+1 f = -1 else y2 = datevec(1)-1 end if do year = y1, y2 if ( (mod(year,400) == 0) .or. ( (mod(year,4) == 0) .and. (mod(year,100) /= 0 )) ) then seconds_since_epoch = seconds_since_epoch + 366*86400 * f else seconds_since_epoch = seconds_since_epoch + 365*86400 * f end if end do do month = 1, datevec(2) - 1 if ( ((mod(year,400) == 0) .or. ( (mod(year,4) == 0) & .and. (mod(year,100) /= 0) )) .and. month == 2 ) then seconds_since_epoch = seconds_since_epoch + (ndpm(month) + 1) * 86400 else seconds_since_epoch = seconds_since_epoch + ndpm(month) * 86400 end if end do seconds_since_epoch = seconds_since_epoch + (datevec(3)-1)*86400 + datevec(4)*3600 & + datevec(5) * 60 + datevec(6) end function seconds_since_epoch_datevec subroutine agnc_scan_time(ustr, factor, dvec, emsg) character(len=*), intent( in) :: ustr integer, intent(out) :: factor, dvec(6) character(len=*), intent(out) :: emsg integer :: j, slen, nspace character :: fstr*20, dlim1, dlim2 factor = -1 emsg = ' ' slen = len(trim(ustr)) j = index(ustr, 'since') if ( j < 1 ) then emsg = 'Keyword "since" not found in unit string '//trim(ustr) return end if fstr = trim(ustr(1:j-1)) if ( fstr .eq. 'second' .or. & fstr .eq. 'seconds' .or. & fstr .eq. 'sec' .or. & fstr .eq. 's' ) then factor = 1 elseif ( fstr .eq. 'minute' .or. & fstr .eq. 'minutes' .or. & fstr .eq. 'min' .or. & fstr .eq. 'm' ) then factor = 60 elseif ( fstr .eq. 'day' .or. & fstr .eq. 'days' .or. & fstr .eq. 'd' ) then factor = 86400 elseif ( fstr .eq. 'year' .or. & fstr .eq. 'years' ) then ! int32(365.242198781 * 86400) factor = 31556926 elseif ( fstr .eq. 'week' .or. & fstr .eq. 'weeks' ) then ! int32(365.242198781 * 86400 / 52) factor = 606864 elseif ( fstr .eq. 'month' .or. & fstr .eq. 'months' ) then ! int32(365.242198781 * 86400 / 12) factor = 2629744 else write(emsg,'("Time unit ", A, " not supported")') trim(fstr) return end if ! step to the character after the 'since ', start of yyyy-mm-dd, mm and dd ! do not have to be zero padded... j = j + 6 read(ustr(j:slen),'(I4,A,I2,A,I2)') dvec(1), dlim1, dvec(2), dlim2, dvec(3) ! find the starting charter of the time nspace = index(ustr(j:slen), ' ') if (nspace > 0 .and. nspace < slen) then j = j+nspace read(ustr(j:slen),'(I2,A,I2,A,I2)') dvec(4), dlim1, dvec(5), dlim2, dvec(6) else ! no time provided dvec(4:6) = 0 end if end subroutine agnc_scan_time end module agioncmd