!*********************************************************************** !* Copyright 2012,2013 * !* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, * !* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,* !* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, * !* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, * !* Adam Dingwell * !* * !* This file is part of FLEXPART WRF * ! * ! FLEXPART is free software: you can redistribute it and/or modify * ! it under the terms of the GNU General Public License as published by* ! the Free Software Foundation, either version 3 of the License, or * ! (at your option) any later version. * ! * ! FLEXPART is distributed in the hope that it will be useful, * ! but WITHOUT ANY WARRANTY; without even the implied warranty of * ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * ! GNU General Public License for more details. * ! * ! You should have received a copy of the GNU General Public License * ! along with FLEXPART. If not, see . * !********************************************************************** subroutine write_ncheader(itime,nesting_level) !***************************************************************************** ! * ! This routine perdefines a netcdf ouput file with information on flexpart * ! settings, releases and topography. * ! * ! Author: A. Dingwell * ! * ! 27 May 2013 * ! * ! Modifications * ! June 5 2013: J. Brioude: generate a header*nc * !***************************************************************************** use point_mod use outg_mod use com_mod implicit none include 'netcdf.inc' integer :: itime,stat ! seconds since simulation start integer :: nesting_level ! 0 for main grid (mother) 1 for nest (child) ! this is written to be easy to expand is additional ! are desired in the future real(kind=dp) :: jul ! Julian date integer :: jjjjmmdd,ihmmss ! date & time as integer character :: adate*8,atime*6 ! date and time strings, used for filename ! Grid related variables real :: xp1,yp1,xp2,yp2 ! temporary coordinates real :: xsw,xne,ysw,yne,tmpx,tmpy,tmplon,tmplat,xl2,yl2 integer :: ncgrid_nx,ncgrid_ny ! nx,ny of current grid integer :: ncgrid_dx,ncgrid_dy ! dx,dy of current grid in m or latlon real :: ncgrid_swlon,ncgrid_swlat ! SW corner of current grid in latlon real :: ncgrid_nelon,ncgrid_nelat ! NE corner of current grid in latlon real :: ncgrid_xm0,ncgrid_ym0 ! lower-left grid coord in metres real :: ncgrid_lon0,ncgrid_lat0 ! lower-left grid coord in latlon ! Grid related 2D-variables (reassigning these here is a bit inefficient but ! it lets us keep a consistent structure of the code, besides it's only once ! per output real,allocatable,dimension (:,:) :: ncgrid_oro,ncgrid_area ! of current grid ! Iterators integer i,j,ix,jy ! NETCDF file related variables integer nclvlid,nclonid,nclatid,ncrecid,ncspcid,ncageid !outgrid dimension ids integer ncrelid,ncrseid ! release points dimension ids integer ncrnvid,ncrmvid,ncspvid ! release points: number,mass,species ids integer ncrtvid,ncrxvid,ncryvid,ncrzvid ! release points: t,x,y,z min/max limits integer nctovid,ncarvid ! Topography and grid area variable-ids integer ncstr1id,ncstr2id,ncstr3id ! decrtiption string length dimid integer nclvlvid,nclonvid,nclatvid,ncspcvid,ncagevid ! outgrid dimension variables integer nclonvid2,nclatvid2 integer ncdimsid3(6),ncdimsid2(5) ! arrays of dimension ids for outgrid 3D & 2D ! NETCDF filename & attribute related variables character descr*11,units*5,ncname*29,coord*11,coordxy*10 integer coordxylen character unit2d*10 ! unit for deposition fields integer unit2dlen ! length of character string ! NETCDF misc variables integer ncid ! local container for netcdf file-id (either ncout or ncoutn) integer ncret ! Return-value of calls to nf_* utils integer :: deflate_level=1 ! compression level integer :: shuffle=1 ! shuffle ! integer :: chunks(2) ! shuffle ! Attribute notation: descr = 'description' units = 'units' coord = 'coordinates' coordxy = 'XLONG XLAT' coordxylen = 10 ! Determine current calendar date, needed for the file name !********************************************************** jul=bdate+real(itime,kind=dp)/86400._dp call caldate(jul,jjjjmmdd,ihmmss) write(adate,'(i8.8)') jjjjmmdd write(atime,'(i6.6)') ihmmss !************************ ! Open header output file !************************ ! write(ncname,'(A8,I2.2,A1,I8.8,A1,I6.6,A3)') & ! 'flxout_d',nesting_level+1,'_',jjjjmmdd,'_',ihmmss,'.nc' ! filename write(ncname,'(A8,I2.2,A3)') & 'header_d',nesting_level+1,'.nc' ! filename ! call nf_set_log_level(3) if (option_verbose.ge.1) write(*,*) & 'write_ncheader: creating file: ',path(1)(1:length(1))//ncname ! call nf_set_chunk_cache(32000000) ! ncret = nf_create(path(1)(1:length(1))//ncname, nf_clobber,ncid) ncret = nf_create(path(1)(1:length(1))//ncname, NF_NETCDF4,ncid) call check_ncerror(ncret) ! Determine which nest/outfile we just created so we can set up the grid !*********************************************************************** if (nesting_level.eq.0) then ! current grid is main grid ncout = ncid ! copy current file handle to ncout ncgrid_nx = numxgrid ncgrid_ny = numygrid ncgrid_nelon = outgrid_nelon ncgrid_nelat = outgrid_nelat ncgrid_swlon = outgrid_swlon ncgrid_swlat = outgrid_swlat allocate(ncgrid_oro(ncgrid_nx,ncgrid_ny),stat=stat) allocate(ncgrid_area(ncgrid_nx,ncgrid_ny),stat=stat) ncgrid_oro = oroout(0:ncgrid_nx-1,0:ncgrid_ny-1) ncgrid_area = area(0:ncgrid_nx-1,0:ncgrid_ny-1) if (outgrid_option.eq.1) then ! input was in latlon ncgrid_dx = dxoutl ncgrid_dy = dyoutl ncgrid_lon0 = outlon0 ncgrid_lat0 = outlat0 else ! input was in metres ncgrid_dx = dxout ncgrid_dy = dyout ncgrid_xm0 = out_xm0 ncgrid_ym0 = out_ym0 endif elseif (nesting_level.eq.1) then ! current grid is nested ncoutn = ncid ! copy current file handle to ncoutn ncgrid_nx = numxgridn ncgrid_ny = numygridn ncgrid_nelon = outgridn_nelon ncgrid_nelat = outgridn_nelat ncgrid_swlon = outgridn_swlon ncgrid_swlat = outgridn_swlat allocate(ncgrid_oro(ncgrid_nx,ncgrid_ny),stat=stat) allocate(ncgrid_area(ncgrid_nx,ncgrid_ny),stat=stat) ncgrid_oro = orooutn(0:ncgrid_nx-1,0:ncgrid_ny-1) ncgrid_area = arean(0:ncgrid_nx-1,0:ncgrid_ny-1) if (outgrid_option.eq.1) then ! input was in latlon ncgrid_dx = dxoutln ncgrid_dy = dyoutln ncgrid_lon0 = outlon0n ncgrid_lat0 = outlat0n else ! input was in metres ncgrid_dx = dxoutn ncgrid_dy = dyoutn ncgrid_xm0 = out_xm0n ncgrid_ym0 = out_ym0n endif endif if (option_verbose.ge.10) & write(*,*) 'write_ncheader: ncout,ncoutn=',ncout,ncoutn ! Write the header information !***************************** !ncret = nf_put_att_text(ncout,nf_global,'TITLE',20,version) !call check_ncerror(ncret) if (ldirect.eq.1) then ! Forward simulation if (option_verbose.ge.10) write(*,10) 'forward simulation attributes' ncret = nf_put_att_int(ncid,nf_global,'SIMULATION_START_DATE',nf_int,1,ibdate) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'SIMULATION_START_TIME',nf_int,1,ibtime) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'SIMULATION_END_DATE',nf_int,1,iedate) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'SIMULATION_END_TIME',nf_int,1,ietime) call check_ncerror(ncret) else ! Backward simulation if (option_verbose.ge.10) write(*,10) 'backward simulation attributes' ncret = nf_put_att_int(ncid,nf_global,'SIMULATION_START_DATE',nf_int,1,iedate) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'SIMULATION_START_TIME',nf_int,1,ietime) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'SIMULATION_END_DATE',nf_int,1,ibdate) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'SIMULATION_END_TIME',nf_int,1,ibtime) call check_ncerror(ncret) endif if (option_verbose.ge.10) write(*,10) 'map projection attributes' if (outgrid_option .eq. 1) then ncret = & nf_put_att_text(ncid,nf_global,'OUTPUT_PROJECTION',20,'Regular Latit/Longit') call check_ncerror(ncret) else if (map_proj_id.eq.1) then ncret = & nf_put_att_text(ncid,nf_global,'OUTPUT_PROJECTION',17,'Lambert conformal') call check_ncerror(ncret) elseif (map_proj_id.eq.2) then ncret = & nf_put_att_text(ncid,nf_global,'OUTPUT_PROJECTION',13,'stereographic') call check_ncerror(ncret) elseif (map_proj_id.eq.3) then ncret = & nf_put_att_text(ncid,nf_global,'OUTPUT_PROJECTION',8,'mercator') call check_ncerror(ncret) elseif (map_proj_id.eq.4) then ncret = & nf_put_att_text(ncid,nf_global,'OUTPUT_PROJECTION',6,'global') call check_ncerror(ncret) endif endif ! Write info common model settings !********************************* if (option_verbose.ge.10) write(*,10) 'common model attributes' if (option_verbose.ge.10) write(*,10) 'OUTPUT_INTERVAL' ncret = nf_put_att_int(ncid,nf_global,'OUTPUT_INTERVAL',nf_int,1,loutstep) call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'AVERAGING_TIME' ncret = nf_put_att_int(ncid,nf_global,'AVERAGING_TIME',nf_int,1,loutaver) call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'AVERAGE_SAMPLING' ncret = nf_put_att_int(ncid,nf_global,'AVERAGE_SAMPLING',nf_int,1,loutsample) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'NSPEC',nf_int,1,nspec) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'NUMRECEPTOR',nf_int,1,numreceptor) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'NAGECLASS',nf_int,1,nageclass) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'NUMRELEASES',nf_int,1,numpoint) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'DISPERSION_METHOD',nf_int,1,method) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'SUBGRID_TOPOGRAPHY',nf_int,1,lsubgrid) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'CONVECTION_PARAM',nf_int,1,lconvection) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'SUBGRID_TOPOGRAPHY',nf_int,1,lsubgrid) call check_ncerror(ncret) ! Write information on output grid setup !*************************************** if (option_verbose.ge.10) write(*,10) 'WEST-EAST_GRID_DIMENSION' ncret = nf_put_att_int(ncid,nf_global,'WEST-EAST_GRID_DIMENSION', & nf_int,1,ncgrid_nx) call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'SOUTH-NORTH_GRID_DIMENSION' ncret = nf_put_att_int(ncid,nf_global,'SOUTH-NORTH_GRID_DIMENSION', & nf_int,1,ncgrid_ny) call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'BOTTOM-TOP_GRID_DIMENSION' ncret = nf_put_att_int(ncid,nf_global,'BOTTOM-TOP_GRID_DIMENSION', & nf_int,1,numzgrid) if (option_verbose.ge.10) write(*,10) 'DX and DY' ncret = nf_put_att_int(ncid,nf_global,'DX',nf_int,1,ncgrid_dx) call check_ncerror(ncret) ncret = nf_put_att_int(ncid,nf_global,'DY',nf_int,1,ncgrid_dy) call check_ncerror(ncret) ! Set up netcdf dimensions !************************* if (option_verbose.ge.10) write(*,10) 'main grid dimensions' ncret = nf_def_dim(ncid,'Time',nf_unlimited,ncrecid) call check_ncerror(ncret) ncret = nf_def_dim(ncid,'DateStrLen',15,ncstr3id) !TODO: WRF format call check_ncerror(ncret) ncret = nf_def_dim(ncid,'west_east',ncgrid_nx,nclonid) call check_ncerror(ncret) ncret = nf_def_dim(ncid,'south_north',ncgrid_ny,nclatid) call check_ncerror(ncret) ncret = nf_def_dim(ncid,'bottom_top',numzgrid,nclvlid) call check_ncerror(ncret) ncret = nf_def_dim(ncid,'species',nspec,ncspcid) call check_ncerror(ncret) ncret = nf_def_dim(ncid,'SpeciesStrLen',10,ncstr1id) call check_ncerror(ncret) ncret = nf_def_dim(ncid,'ageclass',nageclass,ncageid) call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'release point dimensions' ncret = nf_def_dim(ncid,'releases',numpoint,ncrelid) call check_ncerror(ncret) ncret = nf_def_dim(ncid,'ReleaseStrLen',45,ncstr2id) call check_ncerror(ncret) ncret = nf_def_dim(ncid,'ReleaseStartEnd',2,ncrseid) call check_ncerror(ncret) ! Select which dimensions to use for main output grids ncdimsid3(1) = nclonid ! X ncdimsid3(2) = nclatid ! Y ncdimsid3(3) = nclvlid ! Z if (ldirect.eq.1) ncdimsid3(4) = ncspcid ! species if (ldirect.eq.-1) ncdimsid3(4) = ncrelid ! points ncdimsid3(5) = ncageid ! ageclass ncdimsid3(6) = ncrecid ! t ncdimsid2(1) = nclonid ! X ncdimsid2(2) = nclatid ! Y if (ldirect.eq.1) ncdimsid2(3) = ncspcid ! species if (ldirect.eq.-1) ncdimsid2(3) = ncrelid ! points ncdimsid2(4) = ncageid ! ageclass ncdimsid2(5) = ncrecid ! t ! Set up dimension variables !*************************** ! XLONG if (option_verbose.ge.10) write(*,10) 'XLONG dimension variable' ncret = nf_def_var(ncid,'XLONG',nf_real,2,ncdimsid2(1:2),nclonvid) ! Turn on deflate compression, fletcher32 checksum. ncret = NF_DEF_VAR_deflate(ncid,nclonvid, shuffle, 1, deflate_level) ! if (ncret .ne. nf_noerr) call handle_err(retval) ! ncret = NF_DEF_VAR_FLETCHER32(ncid, nclonvid, NF_FLETCHER32) ! if (ncret .ne. nf_noerr) call handle_err(retval) ! ncret = nf_def_var_deflate(ncid,'XLONG',nf_real,2,ncdimsid2(1:2),nclonvid,deflate_level=deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,nclonvid,descr,42,'Longitude of center grid, west is negative') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,nclonvid,units,11,'degree_east') call check_ncerror(ncret) ncret = nf_def_var(ncid,'XLONG_CORNER',nf_real,2,ncdimsid2(1:2),nclonvid2) ncret = NF_DEF_VAR_deflate(ncid,nclonvid2, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,nclonvid2,descr,57,'Longitude of lower left corner of grids, west is negative') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,nclonvid2,units,11,'degree_east') call check_ncerror(ncret) ! XLAT if (option_verbose.ge.10) write(*,10) 'XLAT dimension variable' ncret = nf_def_var(ncid,'XLAT',nf_real,2,ncdimsid2(1:2),nclatvid) ! ncret = nf_def_var_deflate(ncid,'XLAT',nf_real,2,ncdimsid2(1:2),nclatvid,deflate_level=deflate_level) ncret = NF_DEF_VAR_deflate(ncid,nclatvid, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,nclatvid,descr,42,'Latitude of center grid, south is negative') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,nclatvid,units,12,'degree_north') call check_ncerror(ncret) ncret = nf_def_var(ncid,'XLAT_CORNER',nf_real,2,ncdimsid2(1:2),nclatvid2) ncret = NF_DEF_VAR_deflate(ncid,nclatvid2, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,nclatvid2,descr,57,'Latitude of lower left corner of grids, south is negative') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,nclatvid2,units,12,'degree_north') call check_ncerror(ncret) ! ZTOP if (option_verbose.ge.10) write(*,10) 'ZTOP dimension variable' ncret = nf_def_var(ncid,'ZTOP',nf_real,1,ncdimsid3(3),nclvlvid) ncret = NF_DEF_VAR_deflate(ncid,nclvlvid, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,nclvlvid,descr,32, & 'UPPER BOUNDARY OF MODEL LAYER') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,nclvlvid,units,1,'m') call check_ncerror(ncret) ! SPECIES if (option_verbose.ge.10) write(*,10) 'SPECIES dimension variable' ncret = nf_def_var(ncid,'SPECIES',nf_char,2,(/ncstr1id,ncspcid/),ncspcvid) ncret = NF_DEF_VAR_deflate(ncid,ncspcvid, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncspcvid,descr,15,'NAME OF SPECIES') call check_ncerror(ncret) ! AGECLASSES if (option_verbose.ge.10) write(*,10) 'AGECLASSES dimension variable' ncret = nf_def_var(ncid,'AGECLASS',nf_int,1,ncageid,ncagevid) ncret = NF_DEF_VAR_deflate(ncid,ncagevid, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncagevid,descr,27,'MAX AGE OF SPECIES IN CLASS') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncagevid,units,1,'s') call check_ncerror(ncret) ! TIMES if (option_verbose.ge.10) write(*,10) 'TIMES dimension variable' ncret = nf_def_var(ncid,'Times',nf_char,2,(/ncstr3id,ncrecid/),ncrecvid) ncret = NF_DEF_VAR_deflate(ncid,ncrecvid, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncrecvid,descr,42, & 'TIME OF OUTPUT (END OF AVERAGING INTERVAL)') ! Release related variables if (option_verbose.ge.10) write(*,10) 'ReleaseName variable' ncret = nf_def_var(ncid,'ReleaseName',nf_char,2,(/ncstr2id,ncrelid/),ncrnvid) ncret = NF_DEF_VAR_deflate(ncid,ncrnvid, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncrnvid,descr,25,'RELEASE IDENTIFIER/COMMENT') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncrnvid,units,1,'-') call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'ReleaseTstart_end variable' ncret = nf_def_var(ncid,'ReleaseTstart_end', & nf_int,2,(/ncrseid,ncrelid/),ncrtvid) ncret = NF_DEF_VAR_deflate(ncid,ncrtvid, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncrtvid,descr,32, & 'BEGINNING/ENDING TIME OF RELEASE (SECONDS SINCE RUN START)') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncrtvid,units,1,'s') call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'ReleaseXstart_end variable' ncret = nf_def_var(ncid,'ReleaseXstart_end', & nf_float,2,(/ncrseid,ncrelid/),ncrxvid) ncret = NF_DEF_VAR_deflate(ncid,ncrxvid, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncrxvid,descr,32, & 'WEST/EAST BOUNDARIES OF SOURCE') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncrxvid,units,12,'degree_north') call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'ReleaseYstart_end variable' ncret = nf_def_var(ncid,'ReleaseYstart_end', & nf_float,2,(/ncrseid,ncrelid/),ncryvid) ncret = NF_DEF_VAR_deflate(ncid,ncryvid, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncryvid,descr,32, & 'SOUTH/NORTH BOUNDARIES OF SOURCE') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncryvid,units,12,'degree_north') call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'ReleaseZstart_end variable' ncret = nf_def_var(ncid,'ReleaseZstart_end', & nf_float,2,(/ncrseid,ncrelid/),ncrzvid) ncret = NF_DEF_VAR_deflate(ncid,ncrzvid, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncrzvid,descr,31, & 'BOTTOM/TOP BOUNDARIES OF SOURCE') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncrzvid,units,1,'m') call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'ReleaseNP variable' ncret = nf_def_var(ncid,'ReleaseNP',nf_int,1,ncrelid,ncspvid) ncret = NF_DEF_VAR_deflate(ncid,ncspvid, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncspvid,descr,34, & 'TOTAL NUMBER OF PARTICLES RELEASED') ncret = nf_put_att_text(ncid,ncspvid,units,1,'-') call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'ReleaseXMass variable' ncret = nf_def_var(ncid,'ReleaseXMass',nf_real,2,(/ncspcid,ncrelid/),ncrmvid) ncret = NF_DEF_VAR_deflate(ncid,ncrmvid, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncrmvid,descr,18,'TOTAL MASS RELEASED') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncrmvid,units,2,'kg') call check_ncerror(ncret) ! Since we need to exit define mode before we can insert ! variable data, we will include the last file attributes and ! define the last variables here. ! DIRECTION INDEPENDENT OUTPUT VARIABLES if (option_verbose.ge.10) write(*,10) 'TOPOGRAPHY variable' ncret = nf_def_var(ncid,'TOPOGRAPHY',NF_real,2,ncdimsid2(1:2),nctovid) ! ncret = nf_def_var_deflate(ncid,'TOPOGRAPHY',NF_real,2,ncdimsid2(1:2),nctovid,deflate_level=deflate_level) ncret = NF_DEF_VAR_deflate(ncid,nctovid, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,nctovid,descr,33, & 'TERRAIN ELEVATION ABOVE SEA LEVEL') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,nctovid,units,1,'m') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,nctovid,coord,coordxylen,coordxy) call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'GRIDAREA variable' ncret = nf_def_var(ncid,'GRIDAREA',NF_real,2,ncdimsid2(1:2),ncarvid) ! ncret = nf_def_var_deflate(ncid,'GRIDAREA',NF_real,2,ncdimsid2(1:2),ncarvid,deflate_level=deflate_level) ncret = NF_DEF_VAR_deflate(ncid,ncarvid, shuffle, 1, deflate_level) call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncarvid,descr,30, & 'SURFACE AREA OF EACH GRID CELL') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncarvid,units,2,'m2') call check_ncerror(ncret) ncret = nf_put_att_text(ncid,ncarvid,coord,coordxylen,coordxy) call check_ncerror(ncret) ! MAIN OUTPUT VARIABLES ! if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then ! CONCENTRATION ! if (option_verbose.ge.10) write(*,10) 'CONC variable' ! ncret = nf_def_var(ncid,'CONC',NF_REAL,6,ncdimsid3,nccovid) !! chunks(1) = ncgrid_nx !! chunks(2) = ncgrid_ny !!! ncret = NF_DEF_VAR_CHUNKING(ncid, nccovid, NF_CHUNKED, chunks) !! if (ncret .ne. nf_noerr) call check_ncerror(ncret) ! ncret = NF_DEF_VAR_deflate(ncid,nccovid, shuffle, 1, deflate_level) ! ! call check_ncerror(ncret) ! ncret = nf_put_att_text(ncid,nccovid,descr,33, & ! 'CONCENTRATION OF AIRBORNE SPECIES') ! call check_ncerror(ncret) ! ncret = nf_put_att_text(ncid,nccovid,coord,coordxylen,coordxy) ! call check_ncerror(ncret) ! endif ! if ((iout.eq.2).or.(iout.eq.3)) then ! MIXING RATIO ! if (option_verbose.ge.10) write(*,10) 'MIXINGRATIO variable' ! ncret = nf_def_var(ncid,'MIXINGRATIO',NF_REAL,6,ncdimsid3,ncravid) ! ncret = NF_DEF_VAR_deflate(ncid,ncravid, shuffle, 1, deflate_level) ! call check_ncerror(ncret) ! ncret = nf_put_att_text(ncid,ncravid,descr,37, & ! 'MASS MIXING RATIO OF AIRBORNE SPECIES') ! call check_ncerror(ncret) ! ncret = nf_put_att_text(ncid,ncravid,coord,coordxylen,coordxy) ! call check_ncerror(ncret) ! endif ! ! if (ldirect.eq.1) then ! Forward run ! unit2d = 'pg m-2' ! unit2dlen = 6 ! ! if (option_verbose.ge.10) write(*,10) 'DRYDEP variable' ! write(*,*) ncdimsid2 ! ncret = nf_def_var(ncid,'DRYDEP',NF_REAL,5,ncdimsid2,ncddvid) ! ncret = NF_DEF_VAR_deflate(ncid,ncddvid, shuffle, 1, deflate_level) ! call check_ncerror(ncret) ! ncret = nf_put_att_text(ncid,ncddvid,descr,32, & ! 'ACCUMULATED TOTAL DRY DEPOSITION') ! call check_ncerror(ncret) ! ncret = nf_put_att_text(ncid,ncddvid,units,unit2dlen,unit2d) ! call check_ncerror(ncret) ! ncret = nf_put_att_text(ncid,ncddvid,coord,coordxylen,coordxy) ! call check_ncerror(ncret) ! ! if (option_verbose.ge.10) write(*,10) 'WETDEP variable' ! ncret = nf_def_var(ncid,'WETDEP',NF_REAL,5,ncdimsid2,ncwdvid) ! ncret = NF_DEF_VAR_deflate(ncid,ncwdvid, shuffle, 1, deflate_level) ! call check_ncerror(ncret) ! ncret = nf_put_att_text(ncid,ncwdvid,descr,32, & ! 'ACCUMULATED TOTAL WET DEPOSITION') ! call check_ncerror(ncret) ! ncret = nf_put_att_text(ncid,ncwdvid,units,unit2dlen,unit2d) ! call check_ncerror(ncret) ! ncret = nf_put_att_text(ncid,ncwdvid,coord,coordxylen,coordxy) ! call check_ncerror(ncret) ! ! ! Add unit attr to mixing ratio and concentration fields ! ncret = nf_put_att_text(ncid,nccovid,units,6,'ng m-3') !CONC ! call check_ncerror(ncret) ! ncret = nf_put_att_text(ncid,ncravid,units,3,'ppt') !MIX ! call check_ncerror(ncret) ! else ! Backward run ! if (ind_rel.eq.1) then ! release in mass ! write(*,*) 'A' ! !Concentration field should be in 's' (?) ! ncret = nf_put_att_text(ncid,nccovid,units,1,'s') !CONC ! !call check_ncerror(ncret) ! !Mixing ratio field should be in 's kg m-3' (?) ! !ncret = nf_put_att_text(ncid,ncravid,units,8,'s kg m-3') !RATIO ! !call check_ncerror(ncret) ! else ! release in mass mix ! !Concentration should be in 's m3 kg-1' (?) ! write(*,*) 'B' ! ncret = nf_put_att_text(ncid,nccovid,units,9,'s m3 kg-1') !CONC ! !call check_ncerror(ncret) ! !Mixing ratio should be in 's' (?) ! ncret = nf_put_att_text(ncid,ncravid,units,1,'s') !RATIO ! !call check_ncerror(ncret) ! endif ! endif ! Backward/Forward run ! EXIT DEFINE MODE, ENTER DATA MODE ncret = nf_enddef(ncid) call check_ncerror(ncret) ! DIMENSION VARIABLES if (option_verbose.ge.10) write(*,10) 'ZTOP data' ncret = nf_put_var_real(ncid,nclvlvid,outheight) call check_ncerror(ncret) ! X,Y - Lon,Lat if (option_verbose.ge.10) write(*,10) 'XLAT,XLONG data' if (outgrid_option.eq.0) then ! irregular do jy=1,ncgrid_ny do ix=1,ncgrid_nx tmpx=ncgrid_xm0+(float(ix)-0.5)*ncgrid_dx tmpy=ncgrid_ym0+(float(jy)-0.5)*ncgrid_dy call xymeter_to_ll_wrf(tmpx,tmpy,tmplon,tmplat) ncret = nf_put_vara_real(ncid,nclonvid,(/ix,jy/),(/1,1/),tmplon) call check_ncerror(ncret) ncret = nf_put_vara_real(ncid,nclatvid,(/ix,jy/),(/1,1/),tmplat) call check_ncerror(ncret) tmpx=ncgrid_xm0+(float(ix)-1.)*ncgrid_dx tmpy=ncgrid_ym0+(float(jy)-1.)*ncgrid_dy call xymeter_to_ll_wrf(tmpx,tmpy,tmplon,tmplat) ncret = nf_put_vara_real(ncid,nclonvid2,(/ix,jy/),(/1,1/),tmplon) call check_ncerror(ncret) ncret = nf_put_vara_real(ncid,nclatvid2,(/ix,jy/),(/1,1/),tmplat) call check_ncerror(ncret) enddo enddo else do jy=1,ncgrid_ny do ix=1,ncgrid_nx call ll_to_xymeter_wrf(ncgrid_swlon,ncgrid_swlat,xsw,ysw) call ll_to_xymeter_wrf(ncgrid_nelon,ncgrid_nelat,xne,yne) tmpx=xsw+(xne-xsw)*float(ix-1)/float(ncgrid_nx-1) tmpy=ysw+(yne-ysw)*float(jy-1)/float(ncgrid_ny-1) call xymeter_to_ll_wrf(tmpx,tmpy,tmplon,tmplat) xl2=ncgrid_lon0+(float(ix)-0.5)*dxoutl !long yl2=ncgrid_lat0+(float(jy)-0.5)*dyoutl !lat ncret = nf_put_vara_real(ncid,nclonvid,(/ix,jy/),(/1,1/),xl2) call check_ncerror(ncret) ncret = nf_put_vara_real(ncid,nclatvid,(/ix,jy/),(/1,1/),yl2) call check_ncerror(ncret) xl2=ncgrid_lon0+(float(ix)-1.)*dxoutl !long yl2=ncgrid_lat0+(float(jy)-1.)*dyoutl !lat ncret = nf_put_vara_real(ncid,nclonvid2,(/ix,jy/),(/1,1/),xl2) call check_ncerror(ncret) ncret = nf_put_vara_real(ncid,nclatvid2,(/ix,jy/),(/1,1/),yl2) call check_ncerror(ncret) enddo enddo endif ! outgrid_option ! Write information on release points: total number, then for each point: ! start, end, coordinates, # of particles, name, mass !************************************************************************ do i=1,numpoint xp1=xpoint1(i)*dx+xlon0 ! This is probably wrong, but it seems to be yp1=ypoint1(i)*dy+ylat0 ! the same in writeheader*.f90, so I'll leave xp2=xpoint2(i)*dx+xlon0 ! it for now... //AD yp2=ypoint2(i)*dy+ylat0 ! if (option_verbose.ge.10) write(*,10) 'ReleaseTstart_end data' ncret = nf_put_vara_int(ncid,ncrtvid, & ! ReleaseTstart_end (/1,i/),(/2,1/),(/ireleasestart(i),ireleaseend(i)/)) call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'ReleaseXstart_end data' ncret = nf_put_vara_real(ncid,ncrxvid, & ! ReleaseXstart_end (/1,i/),(/2,1/),(/xp1,xp2/)) call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'ReleaseYstart_end data' ncret = nf_put_vara_real(ncid,ncryvid, & !ReleaseYstart_end (/1,i/),(/2,1/),(/yp1,yp2/)) call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'ReleaseZstart_end data' ncret = nf_put_vara_real(ncid,ncrzvid, & !ReleaseZstart_end (/1,i/),(/2,1/),(/zpoint1(i),zpoint2(i)/)) call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'ReleaseXMass data' ncret = nf_put_vara_real(ncid,ncrmvid, & !ReleaseXMass (/1,i/),(/nspec,1/),xmass(i,1:nspec)) call check_ncerror(ncret) if (option_verbose.ge.10) write(*,10) 'ReleaseNP data' ncret = nf_put_vara_int(ncid,ncspvid, & !ReleaseNP i,1,npart(i)) call check_ncerror(ncret) !Release Name/Comment j=1 ! Find the length of each release point comment/name do while( j.lt.45.and.compoint(i)(j+1:j+1).ne." ") j=j+1 enddo if (option_verbose.ge.10) write(*,10) 'ReleaseName data' ncret = nf_put_vara_text(ncid,ncrnvid,(/1,i/),(/j,1/),compoint(i)(1:j)) call check_ncerror(ncret) enddo ! Write age class information !**************************** if (option_verbose.ge.10) write(*,10) 'AGECLASSES data' ncret = nf_put_var_int(ncid,ncagevid,lage(1:nageclass)) call check_ncerror(ncret) ! Write topography to output file !******************************** if (option_verbose.ge.10) write(*,10) 'TOPOGRAPHY data' ! do ix=0,ncgrid_nx-1 do ix=1,ncgrid_nx ncret = nf_put_vara_real(ncid,nctovid, & ! (/ix+1,1/),(/1,ncgrid_ny/),ncgrid_oro(ix,0:ncgrid_ny-1)) (/ix,1/),(/1,ncgrid_ny/),ncgrid_oro(ix,1:ncgrid_ny)) call check_ncerror(ncret) enddo ! Write grid cell surface area !***************************** if (option_verbose.ge.10) write(*,10) 'GRIDAREA data' ! do ix=0,ncgrid_nx-1 do ix=1,ncgrid_nx ncret = nf_put_vara_real(ncid,ncarvid, & ! (/ix+1,1/),(/1,ncgrid_ny/),ncgrid_area(ix,0:ncgrid_ny-1)) (/ix,1/),(/1,ncgrid_ny/),ncgrid_area(ix,1:ncgrid_ny)) call check_ncerror(ncret) enddo ! SAVE CREATED NETCDF TO FILE !**************************** if (option_verbose.ge.1) write(*,*) 'write_ncheader: writing to disk' ncret = nf_sync(ncid) call check_ncerror(ncret) return 10 format('write_ncheader: Setting up ',A) ncret=nf_close(ncid) deallocate(ncgrid_oro,ncgrid_area) end subroutine write_ncheader