subroutine uEMEP_save_for_EMEP_netcdf_file(filename_netcdf,nx,ny,nt,val_array,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 nx,ny,nt
real val_array(nx,ny,nt)!,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
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,:)
!write(*,*) x_vector
!write(*,*) y_vector
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.4" ) )
call check( nf90_put_att(ncid, nf90_global, "title", trim(title_str)) )
call check( nf90_put_att(ncid, nf90_global, "Model", "uEMEP emissions for EMEP" ) )
!Projection data
if (trim(save_emissions_for_EMEP_projection).eq.'lambert') then
call check( nf90_def_var(ncid, "projection_lambert", NF90_int, proj_varid) )
call check( nf90_put_att(ncid, proj_varid, "standard_parallel", EMEP_projection_attributes(1:2) ) )
call check( nf90_put_att(ncid, proj_varid, "longitude_of_central_meridian", EMEP_projection_attributes(3) ) )
call check( nf90_put_att(ncid, proj_varid, "grid_mapping_name", "lambert_conformal_conic" ) )
call check( nf90_put_att(ncid, proj_varid, "latitude_of_projection_origin", EMEP_projection_attributes(4) ) )
call check( nf90_put_att(ncid, proj_varid, "earth_radius", EMEP_projection_attributes(5) ) )
endif
!Define the dimensions
call check( nf90_def_dim(ncid,"time",NF90_UNLIMITED, time_dimid) )
!if (trim(save_emissions_for_EMEP_projection).eq.'lambert') then
call check( nf90_def_dim(ncid, "y", ny, y_dimid) )
call check( nf90_def_dim(ncid, "x", nx, x_dimid) )
!else
!call check( nf90_def_dim(ncid, "lat", ny, y_dimid) )
!call check( nf90_def_dim(ncid, "lon", nx, x_dimid) )
!endif
!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) )
!if (trim(save_emissions_for_EMEP_projection).eq.'lambert') then
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) )
!else
!call check( nf90_def_var(ncid, "lat", NF90_REAL, y_dimid, y_varid) )
!call check( nf90_def_var(ncid, "lon", NF90_REAL, x_dimid, x_varid) )
!endif
!Define the values
dimids3 = (/ x_dimid, y_dimid, time_dimid /)
dimids2 = (/ x_dimid, y_dimid /)
!if (trim(save_emissions_for_EMEP_projection).eq.'lambert') then
call check( nf90_def_var(ncid, "lat", NF90_REAL, dimids2, lat_varid) )
call check( nf90_def_var(ncid, "lon", NF90_REAL, dimids2, lon_varid) )
!endif
!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") )
if (trim(save_emissions_for_EMEP_projection).eq.'lambert') then
call check( nf90_put_att(ncid, y_varid, "units", "m") )
call check( nf90_put_att(ncid, x_varid, "units", "m") )
else
call check( nf90_put_att(ncid, y_varid, "units", "degrees_north") )
call check( nf90_put_att(ncid, x_varid, "units", "degrees_east") )
endif
call check( nf90_put_att(ncid, time_varid, "units", trim(unit_dim_nc(time_dim_nc_index))) )
!Specify other dimension attributes
!if (trim(save_emissions_for_EMEP_projection).eq.'lambert') then
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") )
!else
!call check( nf90_put_att(ncid, y_varid, "standard_name", "latitude") )
!call check( nf90_put_att(ncid, x_varid, "standard_name", "longitude") )
!endif
!Close the definitions
call check( nf90_enddef(ncid) )
!write(*,*) 'here6',shape(val_dim_nc(1:dim_length_nc(time_dim_nc_index),time_dim_nc_index))
!write(*,*) val_dim_nc(1:dim_length_nc(time_dim_nc_index),time_dim_nc_index)
call check( nf90_put_var(ncid, time_varid, val_dim_nc(save_emissions_start_index:save_emissions_end_index,time_dim_nc_index)) )
!call check( nf90_put_var(ncid, time_varid, time_seconds_output(1:dim_length_nc(time_dim_nc_index))) )
!if (trim(save_emissions_for_EMEP_projection).eq.'lambert') then
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) )
!else
!call check( nf90_put_var(ncid, y_varid, y_vector) )
!call check( nf90_put_var(ncid, x_varid, x_vector) )
!endif
call check( nf90_close(ncid) )
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) )
!if (trim(save_emissions_for_EMEP_projection).eq.'lambert') then
call check( nf90_inq_dimid(ncid, "y", y_dimid) )
call check( nf90_inq_dimid(ncid, "x", x_dimid) )
!else
!call check( nf90_inq_dimid(ncid, "lat", y_dimid) )
!call check( nf90_inq_dimid(ncid, "lon", x_dimid) )
!endif
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)) )
!write(*,*) 'here7'
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
!write(*,*) 'here8'
if (trim(save_emissions_for_EMEP_projection).eq.'lambert') then
call check( nf90_put_att(ncid, val_varid, "grid_mapping", "projection_lambert") )
else
call check( nf90_put_att(ncid, val_varid, "grid_mapping", "latitude_longitude") )
endif
!write(*,*) 'here9'
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) )
!write(*,*) 'here10'
endif
if (use_single_time_loop_flag) then
write(unit_logfile,'(A)') 'ERROR: Saving emissions for EMEP will not work when use_single_time_loop_flag=.true. Set to false'
stop
!Add time to the time dimension
call check( nf90_inq_varid(ncid, "time", time_varid) )
!write(*,*) 'here11a'
!call check( nf90_inquire_dimension(ncid, time_dimid, temp_name, n_dims(3)) )
!n_dims(3)=n_dims(3)+1
n_dims(3)=t_loop
!write(*,*) n_dims(3),val_dim_nc(1,time_dim_nc_index)
call check( nf90_put_var(ncid, time_varid, val_dim_nc(1,time_dim_nc_index), start = (/n_dims(3)/) ) )
!write(*,*) n_dims(3),val_dim_nc(1,time_dim_nc_index)
!write(*,*) n_dims
!write(*,*) 'here11b'
!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,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,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
!write(*,*) 'here11'
else
!Write the variable to file. Default is float
if (nf90_type.eq.NF90_byte) then
call check( nf90_put_var(ncid, val_varid, int(val_array,kind=1)) )
elseif (nf90_type.eq.NF90_short) then
call check( nf90_put_var(ncid, val_varid, int(val_array,kind=2)) )
else
!write(*,*) ncid, val_varid, shape(val_array)
call check( nf90_put_var(ncid, val_varid, val_array) )
endif
endif
call check( nf90_close(ncid) )
end subroutine uEMEP_save_for_EMEP_netcdf_file