subroutine uEMEP_save_netcdf_file(unit_logfile_in,filename_netcdf,nx,ny,nt_in,val_array_in,x_array,y_array,lon_array,lat_array,name_array,unit_array,title_str,create_file,valid_min,variable_type,scale_factor)
use uEMEP_definitions
use netcdf
implicit none
character(256) filename_netcdf,name_array,unit_array,title_str,temp_name
integer unit_logfile_in
integer nx,ny,nt_in
real val_array(nx,ny,nt_in),val_array_in(nx,ny,nt_in)!,val_array_temp(nx,ny,nt)
real x_array(nx,ny)
real y_array(nx,ny)
real lon_array(nx,ny)
real lat_array(nx,ny) !,lat_array_temp(nx,ny)
!real time_array(nt)
real x_vector(nx)
real y_vector(ny)
logical create_file
real valid_min
character(256) variable_type
real scale_factor
integer ncid
integer y_dimid,x_dimid,time_dimid
integer y_varid,x_varid,lat_varid,lon_varid,val_varid,time_varid,proj_varid
integer dimids3(3),dimids2(2),chunks3(3)
integer n_dims(3)
integer status
integer nf90_type
integer t
integer n_time_total
integer nt
integer(8) time_seconds_output_nc(nt_in) !Need integer 8 for the averaging
integer i_source
character(2) temp_str
integer i
character(256) temp_str2
if (trim(variable_type).eq.'byte') nf90_type=NF90_BYTE
if (trim(variable_type).eq.'short') nf90_type=NF90_SHORT
if (trim(variable_type).eq.'float') nf90_type=NF90_FLOAT
if (trim(variable_type).eq.'double') nf90_type=NF90_DOUBLE
!Assumes x and y are the dimensions
x_vector=x_array(:,1)
y_vector=y_array(1,:)
n_time_total=end_time_nc_index-start_time_nc_index+1
nt=nt_in
val_array=val_array_in
time_seconds_output_nc=time_seconds_output
!Save averages only
if (save_netcdf_average_flag) then
counter_av=counter_av+1
if (counter_av.gt.n_var_av) then
write(unit_logfile_in,*) 'ERROR: Array size for saving averages (n_var_av) not large enough. Stopping'
stop
endif
if (use_single_time_loop_flag) then
val_array_av(:,:,counter_av)=val_array_av(:,:,counter_av)+val_array(:,:,nt) !nt=1 in this case
time_seconds_output_av(counter_av)=time_seconds_output_av(counter_av)+time_seconds_output_nc(nt)
if (t_loop.eq.end_time_loop_index) then
val_array(:,:,nt)=val_array_av(:,:,counter_av)/n_time_total
time_seconds_output_nc(nt)=time_seconds_output_av(counter_av)/n_time_total
endif
!write(unit_logfile_in,'(a,3i)') 'Saving as average single time loop (nt,counter_av):',nt,counter_av,time_seconds_output_nc(nt)
else
val_array_av(:,:,counter_av)=sum(val_array,3)/size(val_array,3)
time_seconds_output_av(counter_av)=sum(time_seconds_output_nc)/size(time_seconds_output_nc,1)
nt=1
val_array(:,:,nt)=val_array_av(:,:,counter_av)
time_seconds_output_nc(nt)=time_seconds_output_av(counter_av)
!write(unit_logfile_in,'(a,3i)') 'Saving as average multiple time loop (nt,counter_av):',nt,counter_av,time_seconds_output_nc(nt)
endif
endif
!Mask the regions if required
if (use_region_select_and_mask_flag) then
do t=1,nt
! NB: Array subgrid_val is no longer allocated or used
!where (use_subgrid_val(:,:,allsource_index).eq.outside_region_index) val_array(:,:,t)=NODATA_value
where (.not.use_subgrid(:,:,allsource_index)) val_array(:,:,t)=NODATA_value
enddo
endif
if (create_file) then
!Create a netcdf file
!call check( nf90_create(filename_netcdf, nf90_clobber, ncid) )
!call check( nf90_create(filename_netcdf, NF90_HDF5, ncid) )
call check( nf90_create(filename_netcdf, IOR(NF90_HDF5, NF90_CLASSIC_MODEL), ncid) ) !New
!Specify global attributes
call check( nf90_put_att(ncid, nf90_global, "Conventions", "CF-1.6" ) )
call check( nf90_put_att(ncid, nf90_global, "title", trim(title_str)) )
call check( nf90_put_att(ncid, nf90_global, "Model", trim(model_version_str) ) )
!Save some model flags for later reference
call check( nf90_put_att(ncid, nf90_global, "EMEP_grid_interpolation_flag", EMEP_grid_interpolation_flag ) )
call check( nf90_put_att(ncid, nf90_global, "no2_chemistry_scheme_flag", no2_chemistry_scheme_flag ) )
call check( nf90_put_att(ncid, nf90_global, "EMEP_grid_interpolation_size", EMEP_grid_interpolation_size ) )
call check( nf90_put_att(ncid, nf90_global, "EMEP_additional_grid_interpolation_size", EMEP_additional_grid_interpolation_size ) )
call check( nf90_put_att(ncid, nf90_global, "no2_background_chemistry_scheme_flag", no2_background_chemistry_scheme_flag ) )
call check( nf90_put_att(ncid, nf90_global, "local_subgrid_method_flag", local_subgrid_method_flag ) )
call check( nf90_put_att(ncid, nf90_global, "EMEP_emission_grid_interpolation_flag", EMEP_emission_grid_interpolation_flag ) )
if (limit_emep_grid_interpolation_region_to_calculation_region) then
call check( nf90_put_att(ncid, nf90_global, "limit_emep_grid_interpolation_region_to_calculation_region", "true" ) )
else
call check( nf90_put_att(ncid, nf90_global, "limit_emep_grid_interpolation_region_to_calculation_region", "false" ) )
endif
call check( nf90_put_att(ncid, nf90_global, "n_local_fraction_grids", n_local_fraction_grids ) )
do i=1,n_local_fraction_grids
write (temp_str,'(i0)') i
temp_str2="local_fraction_grid_size("//trim(temp_str)//")"
call check( nf90_put_att(ncid, nf90_global, trim(temp_str2), local_fraction_grid_size(i)) )
enddo
call check( nf90_put_att(ncid, nf90_global, "local_fraction_grid_for_EMEP_grid_interpolation", local_fraction_grid_for_EMEP_grid_interpolation ) )
call check( nf90_put_att(ncid, nf90_global, "local_fraction_grid_for_EMEP_additional_grid_interpolation", local_fraction_grid_for_EMEP_additional_grid_interpolation ) )
if (.not.use_GNFR_emissions_from_EMEP_flag) then
call check( nf90_put_att(ncid, nf90_global, "use_GNFR_emissions_from_EMEP_flag", "false" ) )
endif
if (use_GNFR_emissions_from_EMEP_flag.and..not.use_GNFR19_emissions_from_EMEP_flag) then
call check( nf90_put_att(ncid, nf90_global, "use_GNFR13_emissions_from_EMEP_flag", "true" ) )
endif
if (use_GNFR19_emissions_from_EMEP_flag) then
call check( nf90_put_att(ncid, nf90_global, "use_GNFR19_emissions_from_EMEP_flag", "true" ) )
endif
!Write out sources
i=0
do i_source=1,n_source_index
if (calculate_source(i_source)) then
i=i+1
write (temp_str,'(i0)') i
temp_str2="uEMEP_source("//trim(temp_str)//")"
call check( nf90_put_att(ncid, nf90_global, trim(temp_str2), trim(source_file_str(i_source)) ) )
endif
enddo
call check( nf90_put_att(ncid, nf90_global, "n_uEMEP_sources", i ))
i=0
do i_source=1,n_source_index
if (calculate_EMEP_source(i_source)) then
i=i+1
write (temp_str,'(i0)') i
temp_str2="EMEP_source("//trim(temp_str)//")"
call check( nf90_put_att(ncid, nf90_global, trim(temp_str2), trim(source_file_str(i_source)) ) )
endif
enddo
call check( nf90_put_att(ncid, nf90_global, "n_EMEP_sources", i ))
!Projection data
if (projection_type.eq.UTM_projection_index) then
call check( nf90_def_var(ncid, "projection_utm", NF90_int, proj_varid) )
call check( nf90_put_att(ncid, proj_varid, "semi_major_axis", 6378137.0 ) )
call check( nf90_put_att(ncid, proj_varid, "inverse_flattening", 298.257222101 ) )
call check( nf90_put_att(ncid, proj_varid, "grid_mapping_name", "transverse_mercator" ) )
call check( nf90_put_att(ncid, proj_varid, "scale_factor_at_central_meridian", 0.9996 ) )
call check( nf90_put_att(ncid, proj_varid, "latitude_of_projection_origin", 0 ) )
call check( nf90_put_att(ncid, proj_varid, "false_easting", 500000. ) )
if (utm_zone.ge.0) then
call check( nf90_put_att(ncid, proj_varid, "false_northing", 0. ) )
else
call check( nf90_put_att(ncid, proj_varid, "false_northing", 10000000. ) )
endif
call check( nf90_put_att(ncid, proj_varid, "longitude_of_central_meridian", utm_lon0 ) )
endif
if (projection_type.eq.LTM_projection_index) then
call check( nf90_def_var(ncid, "projection_tm", NF90_int, proj_varid) )
call check( nf90_put_att(ncid, proj_varid, "semi_major_axis", 6378137.0 ) )
call check( nf90_put_att(ncid, proj_varid, "inverse_flattening", 298.257222101 ) )
call check( nf90_put_att(ncid, proj_varid, "grid_mapping_name", "transverse_mercator" ) )
call check( nf90_put_att(ncid, proj_varid, "scale_factor_at_central_meridian", 0.9996 ) )
call check( nf90_put_att(ncid, proj_varid, "latitude_of_projection_origin", 0 ) )
call check( nf90_put_att(ncid, proj_varid, "false_easting", 500000. ) )
if (utm_zone.ge.0) then
call check( nf90_put_att(ncid, proj_varid, "false_northing", 0. ) )
else
call check( nf90_put_att(ncid, proj_varid, "false_northing", 10000000. ) )
endif
call check( nf90_put_att(ncid, proj_varid, "longitude_of_central_meridian", ltm_lon0 ) )
endif
if (projection_type.eq.RDM_projection_index) then
!Not properly assigned, same as UTM
call check( nf90_def_var(ncid, "projection_RDM", NF90_int, proj_varid) )
call check( nf90_put_att(ncid, proj_varid, "semi_major_axis", 6378137.0 ) )
call check( nf90_put_att(ncid, proj_varid, "inverse_flattening", 298.257222101 ) )
call check( nf90_put_att(ncid, proj_varid, "grid_mapping_name", "transverse_mercator" ) )
call check( nf90_put_att(ncid, proj_varid, "scale_factor_at_central_meridian", 0.9996 ) )
call check( nf90_put_att(ncid, proj_varid, "latitude_of_projection_origin", 0 ) )
call check( nf90_put_att(ncid, proj_varid, "false_easting", 500000. ) )
call check( nf90_put_att(ncid, proj_varid, "false_northing", 0. ) )
call check( nf90_put_att(ncid, proj_varid, "longitude_of_central_meridian", utm_lon0 ) )
endif
if (projection_type.eq.LAEA_projection_index) then
call check( nf90_def_var(ncid, "projection_ETRS89_LAEA", NF90_int, proj_varid) )
!call check( nf90_put_att(ncid, proj_varid, "semi_major_axis", 6378137.0 ) )
!call check( nf90_put_att(ncid, proj_varid, "inverse_flattening", 298.257222101 ) )
!https://github.com/mdsumner/rasterwise/blob/master/README.md
!EPSG:3035
!int ETRS89-LAEA ;
! ETRS89-LAEA:missing_value = -1. ;
! ETRS89-LAEA:grid_mapping_name = "lambert_azimuthal_equal_area" ;
! ETRS89-LAEA:longitude_of_projection_origin = 10. ;
! ETRS89-LAEA:latitude_of_projection_origin = 52. ;
! ETRS89-LAEA:false_easting = 4321000. ;
! ETRS89-LAEA:false_northing = 3210000. ;
! ETRS89-LAEA:inverse_flattening = 298.257222101 ;
! ETRS89-LAEA:semi_major_axis = 6378137. ;
!call check( nf90_put_att(ncid, proj_varid, "grid_mapping_name", "lambert_azimuthal_equal_area" ) )
!call check( nf90_put_att(ncid, proj_varid, "semi_major_axis", 6378137.0 ) )
!call check( nf90_put_att(ncid, proj_varid, "inverse_flattening", 298.257222101 ) )
!call check( nf90_put_att(ncid, proj_varid, "scale_factor_at_central_meridian", 0.9996 ) )
!call check( nf90_put_att(ncid, proj_varid, "latitude_of_projection_origin", 52. ) )
!call check( nf90_put_att(ncid, proj_varid, "false_easting", 4321000. ) )
!call check( nf90_put_att(ncid, proj_varid, "false_northing", 3210000. ) )
!call check( nf90_put_att(ncid, proj_varid, "longitude_of_projection_origin", 10.) )
call check( nf90_put_att(ncid, proj_varid, "grid_mapping_name", "lambert_azimuthal_equal_area" ) )
call check( nf90_put_att(ncid, proj_varid, "semi_major_axis", projection_attributes(5) ) )
call check( nf90_put_att(ncid, proj_varid, "inverse_flattening", projection_attributes(6) ) )
call check( nf90_put_att(ncid, proj_varid, "latitude_of_projection_origin", projection_attributes(2) ) )
call check( nf90_put_att(ncid, proj_varid, "false_easting", projection_attributes(3) ) )
call check( nf90_put_att(ncid, proj_varid, "false_northing", projection_attributes(4) ) )
call check( nf90_put_att(ncid, proj_varid, "longitude_of_projection_origin", projection_attributes(1)) )
endif
!Define the dimensions
call check( nf90_def_dim(ncid,"time",NF90_UNLIMITED, time_dimid) )
call check( nf90_def_dim(ncid, "y", ny, y_dimid) )
call check( nf90_def_dim(ncid, "x", nx, x_dimid) )
!Define the dimension variables
!call check( nf90_def_var(ncid, "time", NF90_DOUBLE, time_dimid, time_varid) )
call check( nf90_def_var(ncid, "time", NF90_INT, time_dimid, time_varid) )
call check( nf90_def_var(ncid, "y", NF90_REAL, y_dimid, y_varid) )
call check( nf90_def_var(ncid, "x", NF90_REAL, x_dimid, x_varid) )
!Define the values
dimids3 = (/ x_dimid, y_dimid, time_dimid /)
dimids2 = (/ x_dimid, y_dimid /)
call check( nf90_def_var(ncid, "lat", NF90_REAL, dimids2, lat_varid) )
call check( nf90_def_var(ncid, "lon", NF90_REAL, dimids2, lon_varid) )
!Specify the units
call check( nf90_put_att(ncid, lat_varid, "units", "degrees_north") )
call check( nf90_put_att(ncid, lon_varid, "units", "degrees_east") )
call check( nf90_put_att(ncid, y_varid, "units", "m") )
call check( nf90_put_att(ncid, x_varid, "units", "m") )
call check( nf90_put_att(ncid, time_varid, "units", trim(unit_dim_nc(time_dim_nc_index))) )
!Specify other dimension attributes
call check( nf90_put_att(ncid, y_varid, "standard_name", "projection_y_coordinate") )
call check( nf90_put_att(ncid, x_varid, "standard_name", "projection_x_coordinate") )
call check( nf90_put_att(ncid, y_varid, "axis", "Y") )
call check( nf90_put_att(ncid, x_varid, "axis", "X") )
!Close the definitions
call check( nf90_enddef(ncid) )
!call check( nf90_put_var(ncid, time_varid, val_dim_nc(1:dim_length_nc(time_dim_nc_index),time_dim_nc_index)) )
!call check( nf90_put_var(ncid, time_varid, time_seconds_output(1:dim_length_nc(time_dim_nc_index))) )
call check( nf90_put_var(ncid, time_varid, time_seconds_output_nc(1:nt)) )
call check( nf90_put_var(ncid, y_varid, y_vector) )
call check( nf90_put_var(ncid, x_varid, x_vector) )
call check( nf90_put_var(ncid, lat_varid, lat_array) )
call check( nf90_put_var(ncid, lon_varid, lon_array) )
call check( nf90_close(ncid) )
endif
!Do not save any average data for single time loops until the last time step where it will save the average
if (save_netcdf_average_flag.and.use_single_time_loop_flag.and.t_loop.ne.end_time_loop_index) then
!call check( nf90_close(ncid) )
return
endif
!Add to the existing file
call check( nf90_open(filename_netcdf, NF90_WRITE, ncid) )
!Get the dimensions id from the existing file
call check( nf90_inq_dimid(ncid,"time",time_dimid) )
call check( nf90_inq_dimid(ncid, "y", y_dimid) )
call check( nf90_inq_dimid(ncid, "x", x_dimid) )
dimids3 = (/ x_dimid, y_dimid, time_dimid /)
chunks3 = (/ nx, ny, 1 /) !New
call check( nf90_inquire_dimension(ncid, dimids3(1), temp_name, n_dims(1)) )
call check( nf90_inquire_dimension(ncid, dimids3(2), temp_name, n_dims(2)) )
call check( nf90_inquire_dimension(ncid, dimids3(3), temp_name, n_dims(3)) )
status=nf90_inq_varid(ncid, trim(name_array), val_varid)
if (status.ne.nf90_NoErr) then
call check( nf90_redef(ncid) )
!if the variable does not exist then create a new one
!write(*,*) 'Creating new: ',trim(name_array)
call check( nf90_def_var(ncid, trim(name_array), nf90_type, dimids3, val_varid) )
! gzip level 3 compression and shuffling
! optional _FillValue for values which never have been written, unpacked value
call check( nf90_def_var_chunking(ncid, val_varid, NF90_CHUNKED, chunks3) ) !New
call check( nf90_def_var_deflate(ncid, val_varid, 1, 1, 3) ) !New
call check( nf90_put_att(ncid, val_varid, "units", trim(unit_array)) )
!Specify other variable attributes
if (nf90_type.eq.NF90_byte) then
call check( nf90_put_att(ncid, val_varid, "missing_value", int(NODATA_value,kind=1) ) ) !New
call check( nf90_put_att(ncid, val_varid, "valid_min", int(valid_min,kind=1)) )
elseif (nf90_type.eq.NF90_short) then
call check( nf90_put_att(ncid, val_varid, "missing_value", int(NODATA_value,kind=2) ) ) !New
call check( nf90_put_att(ncid, val_varid, "valid_min", int(valid_min,kind=2)) )
else
call check( nf90_put_att(ncid, val_varid, "missing_value", NODATA_value ) ) !New
call check( nf90_put_att(ncid, val_varid, "valid_min", valid_min) )
endif
if (projection_type.eq.UTM_projection_index) then
call check( nf90_put_att(ncid, val_varid, "grid_mapping", "projection_utm") )
elseif (projection_type.eq.LTM_projection_index) then
call check( nf90_put_att(ncid, val_varid, "grid_mapping", "projection_tm") )
elseif (projection_type.eq.LAEA_projection_index) then
call check( nf90_put_att(ncid, val_varid, "grid_mapping", "projection_ETRS89_LAEA") )
elseif (projection_type.eq.RDM_projection_index) then
call check( nf90_put_att(ncid, val_varid, "grid_mapping", "projection_RDM") )
endif
call check( nf90_put_att(ncid, val_varid, "coordinates", "lon lat") )
if (scale_factor.ne.1.) call check( nf90_put_att(ncid, val_varid, "scale_factor", scale_factor) )
!Close the definitions
call check( nf90_enddef(ncid) )
endif
if (use_single_time_loop_flag) then
!Add time to the time dimension
call check( nf90_inq_varid(ncid, "time", time_varid) )
!call check( nf90_inquire_dimension(ncid, time_dimid, temp_name, n_dims(3)) )
!n_dims(3)=n_dims(3)+1
if (save_netcdf_average_flag) then
n_dims(3)=nt
else
n_dims(3)=t_loop
endif
call check( nf90_put_var(ncid, time_varid, time_seconds_output_nc(1), start = (/n_dims(3)/) ) )
!write(*,*) n_dims(3),val_dim_nc(1,time_dim_nc_index)
!write(*,*) n_dims
!Add dimension and array to existing
call check( nf90_inq_varid(ncid, trim(name_array), val_varid) )
if (nf90_type.eq.NF90_byte) then
call check( nf90_put_var(ncid, val_varid, int(val_array(:,:,1:nt),kind=1), start=(/1,1,n_dims(3)/), count=(/n_dims(1),n_dims(2),1/)) )
elseif (nf90_type.eq.NF90_short) then
call check( nf90_put_var(ncid, val_varid, int(val_array(:,:,1:nt),kind=2), start=(/1,1,n_dims(3)/), count=(/n_dims(1),n_dims(2),1/)) )
else
call check( nf90_put_var(ncid, val_varid, val_array, start=(/1,1,n_dims(3)/), count=(/n_dims(1),n_dims(2),1/)) )
endif
else
!Write the whole variable to file. Default is float
if (nf90_type.eq.NF90_byte) then
call check( nf90_put_var(ncid, val_varid, int(val_array(:,:,1:nt),kind=1)) )
elseif (nf90_type.eq.NF90_short) then
call check( nf90_put_var(ncid, val_varid, int(val_array(:,:,1:nt),kind=2)) )
else
call check( nf90_put_var(ncid, val_varid, val_array(:,:,1:nt)) )
endif
endif
call check( nf90_close(ncid) )
end subroutine uEMEP_save_netcdf_file