!################################################################# ! ! contains output for specific stations ! !### macro's ##################################################### ! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" ! !################################################################# module user_output_station use GO, only : gol, goErr, goPr use binas, only : ae, twopi, grav use dims, only : im, jm, lm, dx, dy, xref, yref, xbeg, ybeg, xend, yend use dims, only : nregions, region_name, itaur, xcyc use chem_param, only : ntrace, ntracet, names, fscale use file_hdf, only : THdfFile, Init, Done, TSds, SetDim use ParTools, only : myid, root, tracer_active, tracer_loc, ntracet_ar implicit none private public :: read_stationlist, stat, stations public :: init_station_output public :: evaluate_stationconc public :: reset_stationconc_accumulator public :: write_Stationfiles, StationlistFFilename , output_stationconc public :: free_stationfields, write_stationconc public :: n_stat, io_nrecords #ifdef with_optics public :: sds_station_aod,sds_station_alpha #endif #ifdef with_pm public :: sds_station_pm #endif ! and now a very advanced construction thanks to internet #if defined(with_pm) || defined(with_optics) public :: n_stat2 public :: io_hdf ! Optics wants to write wavelength information here. ! Optics can do it more easily than user_output. ! Optics uses user_output_station. This means that user_output_station ! cannot use optics. ! PM also uses this. #endif ! MODIFIED by Peter Beramaschi. ! PARALLEL version Maarten Krol Jul 2003 ! jun 2005: switch to HDF output file ! jul 2007: paralel writing to increase performance.. ! 10 Jul 2012 - Ph. Le Sager - merged with chem base version to account ! for with_optics, with_pm character(len=100) :: StationlistFFilename ! list of 'regular' stations character(len=100) :: StationlistFFilename2 ! list of 'validation' stations ! character(len=100) :: StationOutFFilename character(len=100) :: InputDir = './' logical :: output_stations_meteo = .true. logical :: write_Stationfiles = .true. character(len=300) :: StationSpeclist integer :: StationlistFileFormat ! 1 : old format ! 2: new format integer,save :: n_stat ! number of 'regular' stations (used for 4DVAR) integer :: n_stat2 ! number of 'validation' stations (not used for 4DVAR) integer, dimension(ntrace) :: station_speclist ! list of species for output (may differ per PE) integer :: station_nspec ! number of species for output (may differ per PE) integer, parameter :: n_metdat = 7 ! max number of meteo data (for output in station output file) integer, parameter :: n_interpol = 3 ! number of station interpolation methods integer, parameter :: n_neighbors = 6 ! number of neighboring grid cells used for estimate of represent. error character(len=4), dimension(n_interpol), parameter :: c_interpol = (/'_grd','_slp','_int'/) character(len=12),dimension(n_metdat ), parameter :: C_metdat = & (/'Temperature', & 'Pressure ',& 'uwind ',& 'vwind ',& 'windspeed ',& 'winddir ',& 'blh '/) character(len=12),dimension(n_metdat ), parameter :: U_metdat = & (/'K ', & 'hPa ', & 'm/s ', & 'm/s ', & 'm/s ', & 'degrees wrtN', & 'm '/) type(THdfFile) :: io_hdf ! contains HDF output file integer :: io_record = 0 integer,save :: io_nrecords ! 'open' datasets of the io_hdf output file: type(TSds), dimension(:,:), allocatable :: Sds type(TSds), dimension(n_metdat) :: Sds_meteo type(TSds) :: Sds_idate type stations character(len=13) :: ident ! station identifier character(len=50) :: name ! station name real :: lat ! latitude of station real :: lon ! longitude of station real :: height ! height of station real :: height_sample ! sampling height real :: height_surface ! height of model surface real :: height_sample_above_surf ! height_sample-height_surface integer :: index ! index of station integer :: region integer :: is integer :: js integer :: ls character(len=2) :: type ! sampling type (FM = flask; CM = continuous) integer :: land ! land or ocean site integer :: nonsurf ! 1: nonsurface sites (e.g. towers); 0: 'regular' sites real,dimension(n_metdat) :: metdat ! meteodata (interpolated) at stations real,dimension(:,:),pointer :: c_mean ! c_mean(species, interpolation) real,dimension(:,:),pointer :: c_std ! c_std(species, interpolation) integer :: ncount ! counter ! In the chem version, the following comment was found: ! "No one is estimating representativeness errors. You can compare the ! different interpolations to acquire something similar." And the code ! with *_nb was commented. ! That code is kept live here in the user_output version until further decision. ! grid neighbors (for estimate of representativeness error) real,dimension(:,:),pointer :: c_mean_nb ! mean concentrations of neighbors (species, neighb. index) real,dimension(:,:),pointer :: c_std_nb ! standard deviation (species, neighb. index) integer,dimension(:,:),pointer :: ncount_nb ! counter (species, neighb. index) real,dimension(:),pointer :: dc_mod ! estimate of representativeness error (species) !(based on 3D gradient to grid cell neighbors) end type stations type(stations), dimension(:), allocatable :: stat #ifdef with_optics Type(TSds) :: sds_station_aod Type(TSds) :: sds_station_alpha #endif #ifdef with_pm Type(TSds) :: sds_station_pm #endif ! --- const ------------------------- character(len=*), parameter :: mname = 'user_output_station' contains !=========================================================================================================== !=========================================================================================================== subroutine free_stationfields(status) implicit none !__IO____________________________________________________________________ integer, intent(inout) :: status !__LOCAL_VARIABLES_______________________________________________________ integer :: i_stat character(len=*), parameter :: rname = mname//'/free_stationfields' !__START_SUBROUTINE______________________________________________________ if(write_Stationfiles) then call close_stationfiles( status ) IF_NOTOK_RETURN(status=1) endif write(gol,*) '************************************', station_nspec ; call goPr !cmk : the following lines do not work on huygens???? if (station_nspec > 0 ) then do i_stat = 1, n_stat+n_stat2 deallocate(stat(i_stat)%c_mean) deallocate(stat(i_stat)%c_std ) enddo do i_stat = 1, n_stat deallocate(stat(i_stat)%c_mean_nb) deallocate(stat(i_stat)%c_std_nb) deallocate(stat(i_stat)%ncount_nb) deallocate(stat(i_stat)%dc_mod) enddo deallocate( stat ) endif !ok status = 0 end subroutine free_stationfields !=========================================================================================================== !=========================================================================================================== subroutine read_stationlist( status ) use GO, only : TrcFile, Init, Done, ReadRc use global_data, only : rcfile use chem_param, only : ntracet, names implicit none !__IO____________________________________________________________________ integer, intent(inout) :: status !__LOCAL_VARIABLES_______________________________________________________ integer :: slen1, slen2 !string length character(len=200) :: dummy integer :: fstatus integer :: i_stat integer :: region integer,dimension(6) :: idate_sim integer :: itau_sim, isim, ispec integer :: ipos, i, n ! x,y resolution (in degrees) for current region real :: dxr, dyr character(len=8) :: name_spec !CMK changed from 6-->8 integer :: sunit1=110 integer :: sunit2=111 type(TrcFile) :: rcF character(len=*), parameter :: rname = mname//'/read_stationlist' !__START_SUBROUTINE______________________________________________________ call Init( rcF, rcfile, status) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'inputdir.station', InputDir,status) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'input.station.filename', StationlistFFilename,status) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'input.station.filename2', StationlistFFilename2,status) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'input.station.fileformat', StationlistFileFormat,status) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'station.speclist', StationSpeclist,status) IF_NOTOK_RETURN(status=1) call Done( rcF ,status) IF_NOTOK_RETURN(status=1) n_stat=0 n_stat2=0 slen1=len_trim(StationlistFFilename) if ( slen1 > 0 ) then !========================================= ! open StationFile (with list of stations) !========================================= open(UNIT=sunit1, FORM='FORMATTED', & FILE=trim(InputDir)//trim(StationlistFFilename), & STATUS='OLD', ERR=1000) read(sunit1, '(a)') dummy ! count number of stations do read(sunit1, '(a)', end=100) dummy !if((dummy(1:1) /= ' ') .and. (dummy(1:7) /= 'species')) then n_stat = n_stat+1 enddo 100 continue end if print*,n_stat slen2=len_trim(StationlistFFilename2) if ( slen2 > 0 ) then !========================================= ! open StationFile (with list of stations) !========================================= open(UNIT=sunit2, FORM='FORMATTED', & FILE=trim(InputDir)//trim(StationlistFFilename2), & STATUS='OLD', ERR=1001) read(sunit2, '(a)') dummy ! count number of stations do read(sunit2, '(a)', END=200) dummy if((dummy(1:1) /= ' ') .and. (dummy(1:7) /= 'species')) then n_stat2 = n_stat2+1 else goto 200 endif enddo 200 continue end if ! Determine Which Pe handles the output of the required species: station_nspec = 0 ipos = 1 do i = index(StationSpeclist(ipos:),' ') if (i == 0) then name_spec = trim(StationSpeclist(ipos:)) else name_spec = trim(StationSpeclist(ipos:ipos+i-1)) ipos = ipos + i endif do n=1,ntracet if (trim(name_spec).eq.'all'.and.tracer_active(n) ) then station_nspec = station_nspec + 1 station_speclist(station_nspec) = tracer_loc(n) else if ((trim(name_spec).eq.trim(names(n)) .and. tracer_active(n)) ) then station_nspec = station_nspec + 1 station_speclist(station_nspec) = tracer_loc(n) endif enddo if(trim(name_spec).eq.'') exit enddo write(gol,*) trim(mname),':station_nspec; station_speclist', station_nspec, station_speclist(1:station_nspec) call goPr if (station_nspec > 0) then allocate(stat(n_stat+n_stat2)) do i_stat = 1, n_stat+n_stat2 allocate(stat(i_stat)%c_mean(station_nspec,n_interpol)) allocate(stat(i_stat)%c_std(station_nspec,n_interpol)) enddo do i_stat = 1, n_stat allocate(stat(i_stat)%c_mean_nb(station_nspec,n_neighbors)) allocate(stat(i_stat)%c_std_nb(station_nspec,n_neighbors)) allocate(stat(i_stat)%ncount_nb(station_nspec,n_neighbors)) allocate(stat(i_stat)%dc_mod(station_nspec)) enddo dyr = dy/yref(1) dxr = dx/xref(1) if ( slen1 > 0 ) then !========================= ! read 'regular' stations !========================= rewind(sunit1) read(sunit1, '(a)') dummy do i_stat=1, n_stat select case(StationlistFileFormat) case(1) read(sunit1, '(a7,1x,f8.2,f8.2,f8.1,1x,a50)') & stat(i_stat)%ident, & stat(i_stat)%lat, & stat(i_stat)%lon, & stat(i_stat)%height, & stat(i_stat)%name stat(i_stat)%nonsurf = 0 case(2) read(sunit1, '(a11,f8.2,4x,f8.2,4x,f8.1,1x,a2,1x,i2,1x,a50)') & stat(i_stat)%ident, & stat(i_stat)%lat, & stat(i_stat)%lon, & stat(i_stat)%height, & stat(i_stat)%type, & stat(i_stat)%nonsurf, & stat(i_stat)%name case(3) read(sunit1, '(a13,f6.2,4x,f8.2,4x,f8.0,4x,i1,5x,i1,a50)') & stat(i_stat)%ident, & stat(i_stat)%lat, & stat(i_stat)%lon, & stat(i_stat)%height, & stat(i_stat)%land, & stat(i_stat)%nonsurf, & stat(i_stat)%name stat(i_stat)%ident = adjustl(stat(i_stat)%ident) endselect stat(i_stat)%index = i_stat ! assume global region as default for average mixing ratios stat(i_stat)%region = 1 stat(i_stat)%is = & int((stat(i_stat)%lon-float(xbeg(1)))/dxr + 0.99999) stat(i_stat)%js = & int((stat(i_stat)%lat-float(ybeg(1)))/dyr + 0.99999) end do ! close stationfile close(sunit1) endif if ( slen2 > 0 ) then !========================= ! read 'validation' stations !========================= rewind(sunit2) read(sunit2, '(a)') dummy do i_stat=n_stat+1, n_stat+n_stat2 select case(StationlistFileFormat) case(1) read(sunit2, '(a7,1x,f8.2,f8.2,f8.1,1x,a50)') & stat(i_stat)%ident, & stat(i_stat)%lat, & stat(i_stat)%lon, & stat(i_stat)%height, & stat(i_stat)%name stat(i_stat)%nonsurf = 0 case(2) read(sunit2, '(a11,f8.2,f8.2,f8.1,1x,a2,1x,i2,1x,a50)') & stat(i_stat)%ident, & stat(i_stat)%lat, & stat(i_stat)%lon, & stat(i_stat)%height, & stat(i_stat)%type, & stat(i_stat)%nonsurf, & stat(i_stat)%name case(3) read(sunit1, '(a13,f6.2,4x,f8.2,4x,f8.0,4x,i1,5x,i1,a50)') & stat(i_stat)%ident, & stat(i_stat)%lat, & stat(i_stat)%lon, & stat(i_stat)%height, & stat(i_stat)%land, & stat(i_stat)%nonsurf, & stat(i_stat)%name stat(i_stat)%ident = adjustl(stat(i_stat)%ident) endselect stat(i_stat)%index = i_stat ! assume global region as default for average mixing ratios stat(i_stat)%region = 1 stat(i_stat)%is = & int((stat(i_stat)%lon-float(xbeg(1)))/dxr + 0.99999) stat(i_stat)%js = & int((stat(i_stat)%lat-float(ybeg(1)))/dyr + 0.99999) end do ! close stationfile close(sunit2) end if !================================================= ! determine region with highest refinement factor !================================================= do i_stat=1, n_stat+n_stat2 do region=1, nregions dyr = dy/yref(region) dxr = dx/xref(region) if ( (stat(i_stat)%lon .ge. xbeg(region) .and. & stat(i_stat)%lon .le. xend(region)) .and. & (stat(i_stat)%lat .ge. ybeg(region) .and. & stat(i_stat)%lat .le. yend(region) ) ) then !===================== ! station is in region !===================== ! determine region with hightes refinement factor if (xref(region) > xref(stat(i_stat)%region)) then stat(i_stat)%region = region stat(i_stat)%is = & int((stat(i_stat)%lon-float(xbeg(region)))/dxr + 0.99999) stat(i_stat)%js = & int((stat(i_stat)%lat-float(ybeg(region)))/dyr + 0.99999) end if end if end do end do if (myid == 0) then write(gol,*) '==========================' ; call goPr write(gol,*) ' regular stations ' ; call goPr write(gol,*) '==========================' ; call goPr do i_stat=1, n_stat write(gol, '(a13,1x,f8.2,f8.2,f8.1,1x,i4,1x,i2,1x,a50)') & stat(i_stat)%ident, & stat(i_stat)%lat, & stat(i_stat)%lon, & stat(i_stat)%height, & stat(i_stat)%region, & stat(i_stat)%nonsurf, & stat(i_stat)%name ; call goPr end do write(gol,*) '==========================' ; call goPr write(gol,*) ' validation stations ' ; call goPr write(gol,*) '==========================' ; call goPr do i_stat=n_stat+1, n_stat+n_stat2 write(gol, '(a11,1x,f8.2,f8.2,f8.1,1x,i4,1x,i2,1x,a50)') & stat(i_stat)%ident, & stat(i_stat)%lat, & stat(i_stat)%lon, & stat(i_stat)%height, & stat(i_stat)%region, & stat(i_stat)%nonsurf, & stat(i_stat)%name ; call goPr end do endif endif ! station nspec > 0 call reset_stationconc_accumulator ! ok status = 0 return ! error handling 1000 write(gol,*) 'read_stationlist: could not open ' // trim(InputDir)//trim(StationlistFFilename) ; call goErr 1001 write(gol,*) 'read_stationlist: could not open ' // trim(InputDir)//trim(StationlistFFilename2) ; call goErr status = 1 end subroutine read_stationlist !=========================================================================================================== !=========================================================================================================== subroutine init_station_output(status) !===================================== ! HDF output for stations ! (instantaneous model results) !===================================== use Dims, only : itaui, itaue, ndyn_max, idatei, idatee use file_hdf, only : Init, WriteAttribute, WriteData, TSds, Done use GO, only : TrcFile, Init, Done, ReadRc use global_data, only : rcfile, outdir implicit none ! --- in/out ----------------------------- integer, intent(out) :: status !__LOCAL_VARIABLES_______________________________________________________ integer :: i_stat integer :: io_status type(TrcFile) :: rcF character(len=12) :: xname integer :: i, ii, j, offsetn type(TSds) :: sds_temp character(len=200) :: station_filename character(len=3) :: smyid character(len=10) :: sidatei, sidatee ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/init_station_output' !__START_SUBROUTINE______________________________________________________ offsetn = sum(ntracet_ar(0:myid-1)) if (write_Stationfiles .and. station_nspec > 0)then call Init( rcF, rcfile ,status) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'output.station.filename', station_filename,status) IF_NOTOK_RETURN(status=1) call Done( rcF ,status) IF_NOTOK_RETURN(status=1) #ifdef MPI i = index(station_filename,'.hdf') write(smyid,'(a1,i2.2)') '_',myid write(sidatei,'(i4.4,3i2.2)') idatei(1:4) write(sidatee,'(i4.4,3i2.2)') idatee(1:4) if(i /= 0) then station_filename = station_filename(1:i-1)//smyid//'_'//sidatei//'_'//sidatee//'.hdf' else station_filename = trim(station_filename)//smyid//'_'//sidatei//'_'//sidatee//'.hdf' endif #else i = index(station_filename,'.hdf') write(sidatei,'(i4.4,3i2.2)') idatei(1:4) write(sidatee,'(i4.4,3i2.2)') idatee(1:4) if(i /= 0) then station_filename = station_filename(1:i-1)//'_'//sidatei//'_'//sidatee//'.hdf' else station_filename = trim(station_filename)//'_'//sidatei//'_'//sidatee//'.hdf' endif #endif call Init( io_hdf, trim(outdir)//'/'//trim(station_filename),'create', status ) IF_NOTOK_RETURN(status=1) io_record = 0 io_nrecords =abs(itaue-itaui)/ndyn_max call WriteAttribute( io_hdf, 'io_nrecords', io_nrecords,status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'n_stat', n_stat,status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'n_stat2', n_stat2,status) ! validation stations IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'n_spec', station_nspec,status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'station_speclist', & station_speclist(1:station_nspec), status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'station_specnames', & names(offsetn+station_speclist(1:station_nspec)), status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'station_lat', stat(:)%lat, status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'station_lon', stat(:)%lon, status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'station_height', stat(:)%height, status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'station_index', stat(:)%index, status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'station_nonsurf', stat(:)%nonsurf, status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'station_region', stat(:)%region, status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'station_is', stat(:)%is, status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'station_js', stat(:)%js, status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'station_names', stat(:)%name, status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'station_ident', stat(:)%ident, status) IF_NOTOK_RETURN(status=1) allocate(Sds(station_nspec, n_interpol)) do i=1, station_nspec do ii= 1,n_interpol xname = trim(names(offsetn + station_speclist(i)))//c_interpol(ii) call init(Sds(i,ii), io_hdf, xname, (/n_stat+n_stat2, io_nrecords/), 'real(4)', status) IF_NOTOK_RETURN(status=1) call SetDim( sds(i,ii), 0, 'n_stations', 'station number', (/ (j, j=1,n_stat+n_stat2) /) , status) IF_NOTOK_RETURN(status=1) call SetDim( sds(i,ii), 1, 'n_records ', 'time slot #', (/ (j, j=1,io_nrecords) /) , status) IF_NOTOK_RETURN(status=1) call WriteAttribute(Sds(i,ii), 'Unit', 'Volume mixing ratio (ppb)', status) IF_NOTOK_RETURN(status=1) enddo enddo if (myid /= 0) output_stations_meteo = .false. if (output_stations_meteo) then do i=1,n_metdat call init(sds_meteo(i), io_hdf, c_metdat(i), (/n_stat+n_stat2, io_nrecords/), 'real(4)', status) IF_NOTOK_RETURN(status=1) call SetDim( sds_meteo(i), 0, 'n_stations', 'station number', (/ (j, j=1,n_stat+n_stat2) /) , status) IF_NOTOK_RETURN(status=1) call SetDim( sds_meteo(i), 1, 'n_records ', 'time slot #', (/ (j, j=1,io_nrecords) /) , status) IF_NOTOK_RETURN(status=1) call WriteAttribute(Sds_meteo(i), 'Unit', u_metdat(i), status) IF_NOTOK_RETURN(status=1) enddo endif call init(sds_idate, io_hdf, 'idate', (/6, io_nrecords/), 'integer(2)', status) IF_NOTOK_RETURN(status=1) call SetDim( sds_idate, 1, 'n_records ', 'time slot #', (/ (j, j=1,io_nrecords) /), status ) IF_NOTOK_RETURN(status=1) call WriteAttribute(Sds_idate, 'Unit', 'year month day hour minutes seconds', status) IF_NOTOK_RETURN(status=1) end if ! write stationfiles = T ! ok status = 0 end subroutine init_station_output !=========================================================================================================== !=========================================================================================================== subroutine close_stationfiles( status ) use file_hdf, only: Done implicit none !__IO____________________________________________________________________ integer, intent(inout) :: status !__LOCAL_VARIABLES_______________________________________________________ integer :: i, ii character(len=*), parameter :: rname = mname//'/close_station_files' !__START_SUBROUTINE______________________________________________________ if (write_Stationfiles .and. station_nspec > 0)then do i=1, station_nspec do ii= 1,n_interpol call Done(sds(i,ii), status) IF_NOTOK_RETURN(status=1) enddo enddo if (myid /= 0) output_stations_meteo = .false. if (output_stations_meteo) then do i=1,n_metdat call Done(sds_meteo(i),status) IF_NOTOK_RETURN(status=1) enddo endif call Done(sds_idate, status) IF_NOTOK_RETURN(status=1) deallocate(Sds) call Done(io_hdf, status) IF_NOTOK_RETURN(status=1) write(gol,*) 'Station file closed ' ; call goPr endif ! ok status = 0 end subroutine close_stationfiles !=========================================================================================================== !=========================================================================================================== subroutine output_stationconc(region, status) use global_data, only : mass_dat #ifndef without_convection use global_data, only : conv_dat #endif use Meteo, only : m_dat, gph_dat, temper_dat, pu_dat, pv_dat, phlb_dat use datetime, only : tau2date implicit none !__IO____________________________________________________________________ integer,intent(in) :: region integer,intent(inout):: status !__LOCAL_VARIABLES_______________________________________________________ real,dimension(:,:,:,:),pointer :: rm ! tracer mass real,dimension(:,:,:,:),pointer :: rxm,rym,rzm ! slopes of tracer mass real,dimension(:,:,:) ,pointer :: m ! air mass real,dimension(:,:,:) ,pointer :: gph ! geopotential height real, dimension(:,:,:,:), allocatable, target :: rmg ! MPI arrays to gather fields. real, dimension(:,:,:,:), allocatable, target :: rxmg, rymg, rzmg real, dimension(:,:,: ), allocatable, target :: mg integer :: i_interpol ! interpolation procedures: real :: sampleheight ! 1: no interpolation ! 2: interpolation based on slopes ! 3: 3D linear interpolations integer :: i,j,l ! grid cell indices integer :: is,js,ls ! i,j,l index of grid cell in which station lis located integer :: isn,jsn,lsn ! i,j,l index of grid cell which is taken as neighbour for linear interpolation integer :: lst, lstn ! MPI implementation integer :: js_north, js_south ! j index of grid cell for neighbours for windspeed_v interpolation integer :: is_west ! i index of grid cell for neighbour for windspeed_u interpolation integer :: n ! index of tracer integer :: i_stat ! index of station real :: ris,rjs,rls ! deviation of station from center of the grid cell (-0.5 ... +0.5) real :: wcx,wcy,wcz ! x,y,z-weight of grid cell (1.0 ... 0.5) real :: dxr, dyr ! x,y resolution (in degrees) for current region real,dimension(0:lm(region)) :: height ! height of vertical grid boundaries integer, dimension(6) :: idate_f ! current idate for region integer :: sunit ! unit number for station output file real :: rm_interpol ! tracer mass, interpolated to station location integer :: ispec ! index over activated tracers integer :: offsetn ! for paralel processing !============= ! meteo output !============= real,dimension(n_stat,n_metdat) :: metdat ! meteodata (interpolated) at stations real,dimension(:,:,:), pointer :: T ! temperature real,dimension(:,:,:), pointer :: phlb ! pressure grid boundaries real,dimension(:,:,:), pointer :: pu ! mass flux x-direction [kg/s] real,dimension(:,:,:), pointer :: pv ! mass flux y-direction [kg/s] real,dimension(:,:), pointer :: blh ! boundary layer height [m] real :: T_interpol ! temperature, interpolated to station location real :: p_interpol ! pressure, interpolated to station location real :: windspeed_u_interpol ! wind speed [m/s] x-direction , interpolated to station location real :: windspeed_v_interpol ! wind speed [m/s] y-direction , interpolated to station location real :: windspeed_interpol ! wind speed [m/s] , interpolated to station location real :: winddir_interpol ! wind direction , interpolated to station location real :: blh_interpol ! boundary layer height , interpolated to station location real :: lat_j ! latitude of grid cell center js [degrees] real :: lat_jn ! latitude of grid cell center jsn [degrees] real :: ddx_j ! x-extension of grid cell js [m] real :: ddx_jn ! x-extension of grid cell jsn [m] real :: ddy ! y-extension of grid cell js [m] real :: vmr ! volume mi integer :: communicator,root_id,lmr,imr,jmr character(len=*), parameter :: rname = mname//'/output_station_conc' !__START_SUBROUTINE______________________________________________________ if (station_nspec == 0) then status = 0 return endif imr = im(region) ; jmr = jm(region) ; lmr = lm(region) m => m_dat(region)%data rm => mass_dat(region)%rm_t rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t ! assign pointers (METEO) gph => gph_dat(region)%data t => temper_dat(region)%data pu => pu_dat(region)%data pv => pv_dat(region)%data #ifndef without_convection blh => conv_dat(region)%blh #endif phlb => phlb_dat(region)%data ! perform now on all PEs ! x and y resolution [degrees] for current region dyr = dy/yref(region) dxr = dx/xref(region) ! idate for current region call tau2date(itaur(region),idate_f) !==================== ! loop over stations !==================== offsetn = sum(ntracet_ar(0:myid-1)) do i_stat=1, n_stat+n_stat2 if (stat(i_stat)%region == region)then !avoid edges! ris = (stat(i_stat)%lon-float(xbeg(region)))/dxr + 0.99999 rjs = (stat(i_stat)%lat-float(ybeg(region)))/dyr + 0.99999 is = int(ris) ! i-index of grid cell in which station is located js = int(rjs) ! j-index of grid cell in which station is located ! x-deviation from the center of the grid cell (-0.5...+0.5) ris = ris-is-0.5 ! y-deviation from the center of the grid cell (-0.5...+0.5) rjs = rjs-js-0.5 !================================= !the neighbour for x interpolation !================================= if ( ris .gt. 0 ) then isn = is+1 else isn = is-1 endif !================================= !the neighbour for y interpolation !================================= if ( rjs .gt. 0 ) then jsn = js+1 else jsn = js-1 endif ! x- / y-weighting of grid cell in which station is located wcx = (1.0-abs(ris)) ! 1.0 ... 0.5 wcy = (1.0-abs(rjs)) ! 1.0 ... 0.5 !================================================================= ! if index of neighbour is exceeding range of region set ! neighbour = current cell (i.e. no interpolation) ! in case of cyclic x-boundaries take corresponding cyclic i index !================================================================= if ( jsn == 0) jsn=1 if ( jsn == jm(region)+1 ) jsn=jm(region) ! isn-->jsn (wouter Peters) if ( xcyc(region) == 0 ) then ! non-cyclic boundaries if ( isn == 0) isn=1 if ( isn == im(region)+1 ) isn=im(region) else ! cyclic x-boundaries if ( isn == 0 ) isn=im(region) if ( isn == im(region)+1 ) isn=1 end if !============================================ ! interpolate the gph to xy position of station !============================================ !ls = 1 !layer do l=0,lm(region) !CMK note: gph now from 1-->lm+1 (dec 2002) height(l) = wcx * wcy* gph(is,js,l+1) + & (1.0-wcx)* wcy* gph(isn,js,l+1) + & wcx *(1.0-wcy)* gph(is,jsn,l+1) + & (1.0-wcx)*(1.0-wcy)* gph(isn,jsn,l+1) end do !=========================== ! determine layer of station !=========================== if(stat(i_stat)%nonsurf == 1) then sampleheight=height(0)+stat(i_stat)%height else sampleheight=stat(i_stat)%height endif ! do not allow sampleheight below model surface sampleheight = max(sampleheight, height(0)) do l=0,lm(region) if(height(l) > sampleheight) exit end do select case(l) case(0) !PB this should no longer occur !! ! force station to model surface ls = 1 rls = -0.5 case default ls = l ! station's layer ! the off-set from the center of the layer (-0.5--->+0.5) ! (interpolation is in (m)) rls = (sampleheight - height(l-1)) / & (height(l)-height(l-1)) - 0.5 end select stat(i_stat)%ls = ls ! store for output to file stat(i_stat)%height_sample = sampleheight ! store for output to file stat(i_stat)%height_surface = height(0) ! store for output to file stat(i_stat)%height_sample_above_surf = sampleheight-height(0) ! store for output to file !================================= !the neighbour for z interpolation !================================= if ( rls > 0 ) then lsn = ls+1 else lsn = ls-1 end if ! z-weighting of grid cell in which station is located wcz = (1.0-abs(rls)) !.0 ... 0.5 !========================================================= ! if vertical neighbor is 0 (which does not exist) ! take vertical layer with l=2 for EXTRApolation to ground !========================================================= if(lsn == 0) then lsn=2 wcz=1.0-rls ! 1.0 ... 1.5 endif if(lsn == lmr+1) then !========================================================= ! if vertical neighbor is lmr+1 (which does not exist) ! -> no interpolation !========================================================= lsn=lmr ! no interpolation wcz=1.0 endif do ispec=1,station_nspec n = station_speclist(ispec) lst = ls lstn = lsn !======================================== ! value of grid box without interpolation !======================================== vmr = rm(is,js,lst,n) / m(is,js,lst) * fscale(offsetn+n) stat(i_stat)%c_mean(ispec,1) = & stat(i_stat)%c_mean(ispec,1) + vmr stat(i_stat)%c_std(ispec,1) = & stat(i_stat)%c_std(ispec,1) + vmr*vmr if(i_stat <= n_stat) then ! for 'regular' stations only !================================================= ! values of neigboring grid boxes ! (used as estimate for representativeness error !================================================= ! lower neighbor if(ls > 1)then vmr = rm(is,js,ls-1,n) / m(is,js,ls-1) * fscale(offsetn+n) stat(i_stat)%c_mean_nb(ispec,1) = stat(i_stat)%c_mean_nb(ispec,1) + vmr stat(i_stat)%c_std_nb(ispec,1) = stat(i_stat)%c_std_nb(ispec,1) + vmr*vmr stat(i_stat)%ncount_nb(ispec,1) = stat(i_stat)%ncount_nb(ispec,1) + 1 endif ! upper neighbor if(ls < lm(region)) then vmr = rm(is,js,ls+1,n) / m(is,js,ls+1) * fscale(offsetn+n) stat(i_stat)%c_mean_nb(ispec,2) = stat(i_stat)%c_mean_nb(ispec,2) + vmr stat(i_stat)%c_std_nb(ispec,2) = stat(i_stat)%c_std_nb(ispec,2) + vmr*vmr stat(i_stat)%ncount_nb(ispec,2) = stat(i_stat)%ncount_nb(ispec,2) + 1 endif ! western neighbor if(is > 1) then vmr = rm(is-1,js,ls,n) / m(is-1,js,ls) * fscale(offsetn+n) stat(i_stat)%c_mean_nb(ispec,3) = stat(i_stat)%c_mean_nb(ispec,3) + vmr stat(i_stat)%c_std_nb(ispec,3) = stat(i_stat)%c_std_nb(ispec,3) + vmr*vmr stat(i_stat)%ncount_nb(ispec,3) = stat(i_stat)%ncount_nb(ispec,3) + 1 endif ! eastern neighbor if(is < im(region)) then vmr = rm(is+1,js,ls,n) / m(is+1,js,ls) * fscale(offsetn+n) stat(i_stat)%c_mean_nb(ispec,4) = stat(i_stat)%c_mean_nb(ispec,4) + vmr stat(i_stat)%c_std_nb(ispec,4) = stat(i_stat)%c_std_nb(ispec,4) + vmr*vmr stat(i_stat)%ncount_nb(ispec,4) = stat(i_stat)%ncount_nb(ispec,4) + 1 endif ! southern neighbor if(js > 1) then vmr = rm(is,js-1,ls,n) / m(is,js-1,ls) * fscale(offsetn+n) stat(i_stat)%c_mean_nb(ispec,5) = stat(i_stat)%c_mean_nb(ispec,5) + vmr stat(i_stat)%c_std_nb(ispec,5) = stat(i_stat)%c_std_nb(ispec,5) + vmr*vmr stat(i_stat)%ncount_nb(ispec,5) = stat(i_stat)%ncount_nb(ispec,5) + 1 endif ! northern neighbor if(js < jm(region)) then vmr = rm(is,js+1,ls,n) / m(is,js+1,ls) * fscale(offsetn+n) stat(i_stat)%c_mean_nb(ispec,6) = stat(i_stat)%c_mean_nb(ispec,6) + vmr stat(i_stat)%c_std_nb(ispec,6) = stat(i_stat)%c_std_nb(ispec,6) + vmr*vmr stat(i_stat)%ncount_nb(ispec,6) = stat(i_stat)%ncount_nb(ispec,6) + 1 endif end if !======================================== ! interpolation using slopes !======================================== !rm-value is obtained from rm + slopes. !slope = rxm = (rm*dX/dx *deltaX/2) rm_interpol = rm(is,js,lst,n) + 2.0*(ris*rxm(is,js,lst,n) + & rjs*rym(is,js,lst,n) + & rls*rzm(is,js,lst,n) ) rm_interpol = max(0.0,rm_interpol) vmr = rm_interpol/m(is,js,lst) *fscale(offsetn+n) stat(i_stat)%c_mean(ispec,2) = & stat(i_stat)%c_mean(ispec,2) + vmr stat(i_stat)%c_std(ispec,2) = & stat(i_stat)%c_std(ispec,2) + vmr*vmr !======================================== ! 3D linear interpolation !======================================== ! this PE handles also the neighbour layer if (lstn <= lmr) then ! lstn = 0 is already excluded (apply vertical EXTRApolation with upper neighbor; see above) rm_interpol = wcx *wcy *wcz *rm(is,js,lst,n) /m(is,js,lst) + & (1.0-wcx)*wcy *wcz *rm(isn,js,lst,n) /m(isn,js,lst) + & wcx *(1.0-wcy)*wcz *rm(is,jsn,lst,n) /m(is,jsn,lst) + & (1.0-wcx)*(1.0-wcy)*wcz *rm(isn,jsn,lst,n) /m(isn,jsn,lst) + & wcx *wcy *(1.0-wcz) *rm(is,js,lstn,n) /m(is,js,lstn) + & (1.0-wcx)*wcy *(1.0-wcz) *rm(isn,js,lstn,n) /m(isn,js,lstn) + & wcx *(1.0-wcy)*(1.0-wcz) *rm(is,jsn,lstn,n) /m(is,jsn,lstn) + & (1.0-wcx)*(1.0-wcy)*(1.0-wcz) *rm(isn,jsn,lstn,n)/m(isn,jsn,lstn) else !other PE calculates neighbouring contribution !PB do NOT apply factor wcz here (no interpolation in z direction in this case) !!!!! rm_interpol = wcx *wcy *rm(is,js,lst,n) /m(is,js,lst) + & (1.0-wcx)*wcy *rm(isn,js,lst,n) /m(isn,js,lst) + & wcx *(1.0-wcy) *rm(is,jsn,lst,n) /m(is,jsn,lst) + & (1.0-wcx)*(1.0-wcy) *rm(isn,jsn,lst,n) /m(isn,jsn,lst) endif vmr = rm_interpol*fscale(offsetn+n) stat(i_stat)%c_mean(ispec,3) = & stat(i_stat)%c_mean(ispec,3) + vmr stat(i_stat)%c_std(ispec,3) = & stat(i_stat)%c_std(ispec,3) + vmr*vmr end do !ispec !====================== !meteo data at stations !====================== if (output_stations_meteo) then !==================================== ! 3D linear interpolation Temperature !==================================== T_interpol = & wcx *wcy *wcz *T(is,js,ls) + & (1.0-wcx)*wcy *wcz *T(isn,js,ls) + & wcx *(1.0-wcy)*wcz *T(is,jsn,ls) + & (1.0-wcx)*(1.0-wcy)*wcz *T(isn,jsn,ls) + & wcx *wcy *(1.0-wcz) *T(is,js,lsn) + & (1.0-wcx)*wcy *(1.0-wcz) *T(isn,js,lsn) + & wcx *(1.0-wcy)*(1.0-wcz) *T(is,jsn,lsn) + & (1.0-wcx)*(1.0-wcy)*(1.0-wcz) *T(isn,jsn,lsn) stat(i_stat)%metdat(1) = & stat(i_stat)%metdat(1) + T_interpol !==================================== ! 3D linear interpolation pressure !==================================== p_interpol =(((0.5-rls) * phlb(is,js,ls) + (0.5+rls) * phlb(is,js,ls+1)) * wcx * wcy + & ((0.5-rls) * phlb(isn,js,ls) + (0.5+rls) * phlb(isn,js,ls+1)) * (1.0-wcx) * wcy + & ((0.5-rls) * phlb(is,jsn,ls) + (0.5+rls) * phlb(is,jsn,ls+1)) * wcx * (1.0-wcy) + & ((0.5-rls) * phlb(isn,jsn,ls)+ (0.5+rls) * phlb(isn,jsn,ls+1))* (1.0-wcx) * (1.0-wcy) ) *0.01 ![Pa] -> [hPa] stat(i_stat)%metdat(2) = & stat(i_stat)%metdat(2) + p_interpol !==================================== ! 3D windspeed_u_interpol (x-direction) !==================================== ! pu is eastward flux through east grid box boundary ! windspeed_u_interpol wind component from east ! latitude of grid cell center js [degrees] lat_j =ybeg(region)+(js -0.5)*dyr ! latitude of grid cell center jsn [degrees] lat_jn=ybeg(region)+(jsn-0.5)*dyr ! x-extension of grid cell js [m] ddx_j =ae * twopi * dxr / 360.0 * cos(lat_j*twopi/360.0) ! x-extension of grid cell jsn [m] ddx_jn=ae * twopi * dxr / 360.0 * cos(lat_jn*twopi/360.0) is_west=is-1 if ( xcyc(region) == 0 ) then ! non-cyclic boundaries if ( is_west == 0 ) is_west=1 else ! cyclic x-boundaries if ( is_west == 0 ) is_west=im(region) end if !west border !east border windspeed_u_interpol=((0.5-ris) * pu(is_west,js,ls)/m(is_west,js,ls) + & (0.5+ris) * pu(is,js,ls)/ m(is,js,ls) ) * & ddx_j * wcy * wcz + & ((0.5-ris) * pu(is_west,jsn,ls)/m(is_west,jsn,ls) + & (0.5+ris) * pu(is,jsn,ls)/m(is,jsn,ls) ) * & ddx_jn * (1.0-wcy)* wcz + & ((0.5-ris) * pu(is_west,js,lsn)/m(is_west,js,lsn) + & (0.5+ris) * pu(is,js,lsn)/ m(is,js,lsn)) * & ddx_j * wcy * (1.0-wcz) + & ((0.5-ris) * pu(is_west,jsn,lsn)/m(is_west,js,lsn) + & (0.5+ris) * pu(is,jsn,lsn)/ m(is,jsn,lsn) ) * & ddx_jn * (1.0-wcy)* (1.0-wcz) !==================================== ! 3D windspeed_v_interpol (y-direction) !==================================== ! pv northward flux through south grid box boundary ! windspeed_v_interpol wind component from north ddy =ae * twopi * dyr / 360.0 ! y-extension of grid cell [m] js_south = js-1 if (js_south==0) js_south=1 js_north = js+1 if (js_north==(jm(region)+1)) js_north=jm(region) !south border !north border windspeed_v_interpol=(( & (0.5-rjs)*pv(is,js,ls)/m(is,js_south,ls) + & (0.5+rjs)*pv(is,js_north,ls)/ m(is,js,ls) ) * & wcx * wcz + & ((0.5-rjs)*pv(isn,js,ls)/m(isn,js_south,ls) + & (0.5+rjs)*pv(isn,js_north,ls)/ m(isn,js,ls) ) * & (1.0 - wcx)* wcz + & ((0.5-rjs)*pv(is,js,lsn)/m(is,js_south,lsn) + & (0.5+rjs)*pv(is,js_north,lsn)/ m(is,js,lsn) ) * & wcx * (1.0-wcz) + & ((0.5-rjs)*pv(isn,js,lsn)/m(isn,js_south,lsn) + & (0.5+rjs)*pv(isn,js_north,lsn)/ m(isn,js,lsn) ) * & (1.0 - wcx)* (1.0-wcz) ) * & ddy windspeed_interpol = & sqrt(windspeed_u_interpol**2 + windspeed_v_interpol**2) stat(i_stat)%metdat(3) = & stat(i_stat)%metdat(3) + windspeed_u_interpol stat(i_stat)%metdat(4) = & stat(i_stat)%metdat(4) + windspeed_v_interpol stat(i_stat)%metdat(5) = & stat(i_stat)%metdat(5) + windspeed_interpol #ifndef without_convection !============================================ ! interpolate the blh to xy position of station !============================================ blh_interpol = & wcx * wcy* blh(is,js) + & (1.0-wcx)* wcy* blh(isn,js) + & wcx *(1.0-wcy)* blh(is,jsn) + & (1.0-wcx)*(1.0-wcy)* blh(isn,jsn) stat(i_stat)%metdat(7) = & stat(i_stat)%metdat(7) + blh_interpol #endif endif ! meteo out... stat(i_stat)%ncount = stat(i_stat)%ncount + 1 end if end do ! end i_stat loop nullify(m) nullify(rm) nullify(rxm) nullify(rym) nullify(rzm) nullify(gph) nullify(t) nullify(pu) nullify(pv) nullify(phlb) nullify(blh) ! ok status = 0 end subroutine output_stationconc !=========================================================================================================== !=========================================================================================================== subroutine reset_stationconc_accumulator !================================================================== ! initialize accumulation arrays for average station concentration !================================================================== implicit none integer i_stat if (station_nspec > 0) then do i_stat=1, n_stat+n_stat2 stat(i_stat)%c_mean(:,:) = 0.0 stat(i_stat)%c_std(:,:) = 0.0 if (output_stations_meteo) stat(i_stat)%metdat(:) = 0.0 stat(i_stat)%ncount = 0 end do endif end subroutine reset_stationconc_accumulator !=========================================================================================================== !=========================================================================================================== subroutine evaluate_stationconc(status) !================================================================== ! evaluate mean concentration and std !================================================================== use dims, only : idate use file_hdf, only : Init, WriteData, TSds, Done, select use file_hdf, only : CheckInfo implicit none !__IO____________________________________________________________________ integer, intent(inout) :: status !__LOCAL_VARIABLES_______________________________________________________ integer i_stat, i, ii, io_status integer n_data, isds, nsds, itp real temp character(len=12) :: xname real, dimension(:,:), allocatable :: cout ! to get output integer, dimension(:,:), allocatable :: iidate real :: windspeed_u_interpol ! wind speed [m/s] x-direction , interpolated to station location real :: windspeed_v_interpol ! wind speed [m/s] y-direction , interpolated to station location real :: windspeed_interpol ! wind speed [m/s] , interpolated to station location real :: winddir_interpol ! wind direction , interpolated to station location character(len=*), parameter :: rname = mname//'/evaluate_stationconc' !__START_SUBROUTINE______________________________________________________ !if(myid /= root) return if (station_nspec > 0) then do i_stat=1, n_stat+n_stat2 n_data = stat(i_stat)%ncount if (n_data > 0) then stat(i_stat)%c_mean(:,:) = stat(i_stat)%c_mean(:,:) / n_data if(output_stations_meteo) stat(i_stat)%metdat(:) = stat(i_stat)%metdat(:) / n_data if (n_data > 1) then do i=1, station_nspec do ii= 1, n_interpol temp = ((stat(i_stat)%c_std(i,ii) - n_data * ((stat(i_stat)%c_mean(i,ii))**2) ) / (n_data-1) ) if(temp > 0) then stat(i_stat)%c_std(i,ii) = sqrt(temp) else stat(i_stat)%c_std(i,ii) = 0.0 endif enddo enddo else stat(i_stat)%c_std(:,:)=-1.0 end if else stat(i_stat)%c_mean(:,:)= -1.0 stat(i_stat)%c_std(:,:) = -1.0 if(output_stations_meteo) stat(i_stat)%metdat(:) = -999. end if if(output_stations_meteo) then windspeed_u_interpol = stat(i_stat)%metdat(3) windspeed_v_interpol = stat(i_stat)%metdat(4) windspeed_interpol = stat(i_stat)%metdat(5) !==================================== ! wind direction !==================================== ! definition of winddirection: ! east : 90 ! south : 180 ! west : 270 ! north : 360 if ((windspeed_u_interpol < 0) .and. (windspeed_v_interpol < 0))then ! wind from north east (0...90) winddir_interpol = atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi else if ((windspeed_u_interpol < 0) .and. (windspeed_v_interpol > 0)) then ! wind from south east (90...180) winddir_interpol = 180.0+atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi else if ((windspeed_u_interpol < 0) .and. (windspeed_v_interpol == 0)) then ! wind from east (90) winddir_interpol = 90.0 else if ((windspeed_u_interpol > 0) .and. (windspeed_v_interpol < 0)) then ! wind from north west (270... 360) winddir_interpol = 360.0+atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi else if ((windspeed_u_interpol > 0) .and. (windspeed_v_interpol > 0)) then ! wind from south west (180...270) winddir_interpol = 180.0+atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi else if ((windspeed_u_interpol > 0) .and. (windspeed_v_interpol == 0)) then ! wind from west (270) winddir_interpol = 270.0 else if ((windspeed_u_interpol == 0) .and. (windspeed_v_interpol < 0)) then ! wind from north (360) winddir_interpol = 360.0 else if ((windspeed_u_interpol == 0) .and. (windspeed_v_interpol > 0)) then ! wind from south (180) winddir_interpol = 180.0 else if ((windspeed_u_interpol == 0) .and. (windspeed_v_interpol == 0)) then ! no wind winddir_interpol = -1.0 else write(gol,*)'output_stationconc: error wind direction'; call goErr TRACEBACK status=1; return end if stat(i_stat)%metdat(6) = winddir_interpol endif end do endif !ok status = 0 end subroutine evaluate_stationconc !=========================================================================================================== !=========================================================================================================== subroutine write_stationconc(status) !================================================================== ! write mean concentration !================================================================== use dims, only : idate use file_hdf, only : Init, WriteData, TSds, Done use file_hdf, only : WriteAttribute implicit none !__IO____________________________________________________________________ integer, intent(inout) :: status !__LOCAL_VARIABLES_______________________________________________________ integer i, ii, io_status, i_stat integer n_data, isds, nsds, itp real, dimension(:,:), allocatable :: cout ! to get output integer, dimension(:,:), allocatable :: iidate character(len=*), parameter :: rname = mname//'/write_stationconc' !__START_SUBROUTINE______________________________________________________ ! output if (station_nspec > 0 ) then if (io_record == 0) then call WriteAttribute( io_hdf, 'station_ls', stat(:)%ls, status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'height_sample', stat(:)%height_sample, status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'height_surface', stat(:)%height_surface, status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'height_sample_above_surf', stat(:)%height_sample_above_surf, status) IF_NOTOK_RETURN(status=1) endif allocate(cout(n_stat+n_stat2,1)) do i=1, station_nspec do ii= 1,n_interpol do i_stat = 1, n_stat+n_stat2 cout(i_stat,1) = stat(i_stat)%c_mean(i,ii) enddo call writedata(sds(i,ii), cout, status, start=(/0,io_record/)) IF_NOTOK_RETURN(status=1) enddo enddo if (output_stations_meteo) then do i=1, n_metdat do i_stat = 1, n_stat+n_stat2 cout(i_stat,1) = stat(i_stat)%metdat(i) enddo call writedata(sds_meteo(i), cout, status, start=(/0,io_record/)) IF_NOTOK_RETURN(status=1) enddo endif deallocate(cout) allocate(iidate(6,1)) iidate(:,1) = idate call writedata(sds_idate, iidate, status, start=(/0,io_record/) ) IF_NOTOK_RETURN(status=1) deallocate(iidate) io_record = io_record + 1 endif ! ok status = 0 end subroutine write_stationconc !=========================================================================================================== !=========================================================================================================== end module user_output_station