!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS 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 FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module station_data_mod ! ! Giang Nong ! ! ! This module is used for outputing model results in a list ! of stations (not gridded arrrays). The user needs to supply ! a list of stations with lat, lon values of each station. ! Data at a single point (I,J) that is closest to each station will ! be written to a file. No interpolation is made when a station ! is between two or more grid points.
! In the output file, a 3D field will have a format of array(n1,n2) and ! a 2D field is array(n1) where n1 is number of stations and n2 is number ! of vertical levels or depths. ! ! ! Here are some basic steps of how to use station_data_mod
!1/Call data_station_init
! user needs to supply 2 tables: list_stations and station_data_table as follows:
! example of list of stations (# sign means comment)
! # station_id lat lon
! station_1 20.4 100.8
! example of station_data_table (# sign means comment)
! # General descriptor
! Am2p14 station data
!# start time (should be the same as model's initial time)
! 19800101
!# file inforamtion
!# filename, output_frequency, frequency_unit, time_axis_unit
! "ocean_day" 1 "days" "hours"
!# field information
!# module field_name filename time_method pack
! Ice_mod temperature ocean_day . true. 2
! Ice_mod pressure ocean_day .false. 2
! 2/ ! Call register_station_field to register each field that needs to be written to a file, the call ! register_station_field returns a field_id that will be used later in send_station_data
! 3/ ! Call send_station_data will send data at each station in the list ! to a file
! 4/ Finally, call station_data_end after the last time step.
! use axis_utils_mod, only: nearest_index use mpp_io_mod, only : mpp_open, MPP_RDONLY, MPP_ASCII, mpp_close,MPP_OVERWR,MPP_NETCDF, & mpp_write_meta, MPP_SINGLE, mpp_write, fieldtype,mpp_flush use fms_mod, only : error_mesg, FATAL, WARNING, stdlog, write_version_number,& mpp_pe, lowercase, stdout, close_file use mpp_mod, only : mpp_npes, mpp_sync, mpp_root_pe, mpp_send, mpp_recv, mpp_max, & mpp_get_current_pelist use mpp_domains_mod,only: domain2d, mpp_get_compute_domain use diag_axis_mod, only : diag_axis_init use diag_output_mod,only: write_axis_meta_data, write_field_meta_data,diag_fieldtype,done_meta_data use diag_manager_mod,only : get_date_dif, DIAG_SECONDS, DIAG_MINUTES, DIAG_HOURS, & DIAG_DAYS, DIAG_MONTHS, DIAG_YEARS use diag_util_mod,only : diag_time_inc use time_manager_mod, only: operator(>),operator(>=),time_type,get_calendar_type,NO_CALENDAR,set_time, & set_date, increment_date, increment_time implicit none private integer, parameter :: max_fields_per_file = 150 integer, parameter :: max_files = 31 integer :: num_files = 0 integer :: num_stations = 0 integer :: max_stations = 20 integer :: max_output_fields = 100 integer :: num_output_fields = 0 real :: EMPTY = 0.0 real :: MISSING = 1.E20 logical :: module_is_initialized = .false. logical :: need_write_axis = .true. integer :: base_year, base_month, base_day, base_hour, base_minute, base_second type (time_type) :: base_time character (len=10) :: time_unit_list(6) = (/'seconds ', 'minutes ', & 'hours ', 'days ', 'months ', 'years '/) integer, parameter :: EVERY_TIME = 0 integer, parameter :: END_OF_RUN = -1 character(len=128) :: version = '' character(len=128) :: tagname = '' character(len=256) :: global_descriptor character (len = 7) :: avg_name = 'average' integer :: total_pe integer :: lat_axis, lon_axis integer, allocatable:: pelist(:) character(len=32) :: pelist_name type file_type character(len=128) :: name integer :: output_freq integer :: output_units integer :: time_units integer :: fields(max_fields_per_file) integer :: num_fields integer :: file_unit integer :: time_axis_id, time_bounds_id type (time_type) :: last_flush type(fieldtype) :: f_avg_start, f_avg_end, f_avg_nitems, f_bounds end type file_type type station_type character(len=128) :: name real :: lat, lon integer :: id integer :: global_i, global_j ! index of global grid integer :: local_i, local_j ! index on the current PE logical :: need_compute ! true if the station present in this PE end type station_type type group_field_type integer :: output_file integer :: num_station ! number of stations on this PE integer, pointer :: station_id(:) =>null() ! id of station on this PE character(len=128) :: output_name, module_name,long_name, units logical :: time_average,time_max,time_min, time_ops, register integer :: pack, axes(2), num_axes character(len=8) :: time_method !could be: true, false, max, min, mean, ... real, pointer :: buffer(:, :)=>null() integer :: counter,nlevel type(time_type) :: last_output, next_output type(fieldtype) :: f_type end type group_field_type type global_field_type real, pointer :: buffer(:,:)=>null() integer :: counter end type global_field_type type(global_field_type),save :: global_field type (file_type),save :: files(max_files) type(group_field_type),allocatable,save :: output_fields(:) type (station_type),allocatable :: stations(:) type(diag_fieldtype),save :: diag_field public register_station_field, send_station_data, station_data_init, station_data_end interface register_station_field module procedure register_station_field2d module procedure register_station_field3d end interface interface send_station_data module procedure send_station_data_2d module procedure send_station_data_3d end interface contains ! ! ! ! read in lat. lon of each station
! create station_id based on lat, lon
! read station_data_table, initialize output_fields and output files
! ! subroutine station_data_init() character(len=128) :: station_name real :: lat, lon integer :: iunit,nfiles,nfields,time_units,output_freq_units,j,station_id,io_status,logunit logical :: init_verbose character(len=128) :: record type file_part_type character(len=128) :: name integer :: output_freq character(len=10) :: output_freq_units integer :: format ! must always be 1 for netcdf files character(len=10) :: time_unit end type file_part_type type field_part_type character(len=128) :: module_name,field_name,file_name character(len=8) :: time_method integer :: pack end type field_part_type type(file_part_type) :: record_files type(field_part_type) :: record_fields namelist /station_data_nml/ max_output_fields, max_stations,init_verbose if (module_is_initialized) return init_verbose = .false. total_pe = mpp_npes() allocate(pelist(total_pe)) call mpp_get_current_pelist(pelist, pelist_name) ! read namelist call mpp_open(iunit, 'input.nml',form=MPP_ASCII,action=MPP_RDONLY) read(iunit,station_data_nml,iostat=io_status) logunit = stdlog() write(logunit, station_data_nml) if (io_status > 0) then call error_mesg('station_data_init', 'Error reading station_data_nml',FATAL) endif call mpp_close (iunit) allocate(output_fields(max_output_fields), stations(max_stations)) ! read list of stations if(init_verbose) then logunit = stdout() write(logunit, *) ' ' write(logunit, *) '****** Summary of STATION information from list_stations ********' write(logunit, *) ' ' write(logunit, *) 'station name ', ' latitude ', ' longitude ' write(logunit, *) ' ' endif call mpp_open(iunit, 'list_stations',form=MPP_ASCII,action=MPP_RDONLY) do while (num_stations 0) then stations(station_id)%name = station_name else call error_mesg('station_data_init','station DUPLICATED in file list_stations', FATAL) endif logunit = stdout() if( init_verbose.and. mpp_pe() == mpp_root_pe()) & write(logunit,1)stations(station_id)%name,stations(station_id)%lat,stations(station_id)%lon 1 format(1x,A18, 1x,F8.2,4x,F8.2) 75 continue enddo call error_mesg('station_data_init','max_stations exceeded, increase it via namelist', FATAL) 76 continue call mpp_close (iunit) logunit = stdout() if(init_verbose) write(logunit, *)'*****************************************************************' ! read station_data table call mpp_open(iunit, 'station_data_table',form=MPP_ASCII,action=MPP_RDONLY) ! Read in the global file labeling string read(iunit, *, end = 99, err=99) global_descriptor ! Read in the base date read(iunit, *, end = 99, err = 99) base_year, base_month, base_day, & base_hour, base_minute, base_second if (get_calendar_type() /= NO_CALENDAR) then base_time = set_date(base_year, base_month, base_day, base_hour, & base_minute, base_second) else ! No calendar - ignore year and month base_time = set_time(base_hour*3600+base_minute*60+base_second, base_day) base_year = 0 base_month = 0 end if nfiles=0 do while (nfiles <= max_files) read(iunit,'(a)',end=86,err=85) record if (record(1:1) == '#') cycle read(record,*,err=85,end=85)record_files%name,record_files%output_freq, & record_files%output_freq_units,record_files%format,record_files%time_unit if(record_files%format /= 1) cycle !avoid reading field part time_units = 0 output_freq_units = 0 do j = 1, size(time_unit_list(:)) if(record_files%time_unit == time_unit_list(j)) time_units = j if(record_files%output_freq_units == time_unit_list(j)) output_freq_units = j end do if(time_units == 0) & call error_mesg('station_data_init',' check time unit in station_data_table',FATAL) if(output_freq_units == 0) & call error_mesg('station_data_init',', check output_freq in station_data_table',FATAL) call init_file(record_files%name,record_files%output_freq, output_freq_units,time_units) 85 continue enddo call error_mesg('station_data_init','max_files exceeded, increase max_files', FATAL) 86 continue rewind(iunit) nfields=0 do while (nfields <= max_output_fields) read(iunit,'(a)',end=94,err=93) record if (record(1:1) == '#') cycle read(record,*,end=93,err=93) record_fields if (record_fields%pack .gt. 8 .or.record_fields%pack .lt. 1) cycle !avoid reading file part nfields=nfields+1 call init_output_field(record_fields%module_name,record_fields%field_name, & record_fields%file_name,record_fields%time_method,record_fields%pack) 93 continue enddo call error_mesg('station_data_init','max_output_fields exceeded, increase it via nml ', FATAL) 94 continue call close_file(iunit) call check_duplicate_output_fields call write_version_number (version, tagname) module_is_initialized = .true. return 99 continue call error_mesg('station_data_init','error reading station_datatable',FATAL) end subroutine station_data_init !---------------------------------------------------------------------- subroutine check_duplicate_output_fields() ! pair(output_name and output_file) should be unique in data_station_table, ERROR1 ! pair(module_name and output_name) should be unique in data_station_table, ERROR2 integer :: i, j, tmp_file character(len=128) :: tmp_name, tmp_module if(num_output_fields <= 1) return do i = 1, num_output_fields-1 tmp_name = trim(output_fields(i)%output_name) tmp_file = output_fields(i)%output_file tmp_module = trim(output_fields(i)%module_name) do j = i+1, num_output_fields if((tmp_name == trim(output_fields(j)%output_name)).and. & (tmp_file == output_fields(j)%output_file)) & call error_mesg (' ERROR1 in station_data_table:', & &' module/field '//tmp_module//'/'//tmp_name//' duplicated', FATAL) if((tmp_name == trim(output_fields(j)%output_name)).and. & (tmp_module == trim(output_fields(j)%module_name))) & call error_mesg (' ERROR2 in station_data_table:', & &' module/field '//tmp_module//'/'//tmp_name//' duplicated', FATAL) enddo enddo end subroutine check_duplicate_output_fields !---------------------------------------------------------------------- function get_station_id(lat,lon) integer :: get_station_id, i real, intent(in):: lat,lon ! each station should have distinct lat and lon get_station_id = -1 do i = 1, num_stations if(stations(i)%lat == lat .and. stations(i)%lon == lon) return enddo num_stations = num_stations + 1 stations(num_stations)%id = num_stations stations(num_stations)%lat = lat stations(num_stations)%lon = lon stations(num_stations)%need_compute = .false. stations(num_stations)%global_i = -1; stations(num_stations)%global_j = -1 stations(num_stations)%local_i = -1 ; stations(num_stations)%local_j = -1 get_station_id = num_stations end function get_station_id !---------------------------------------------------------------------- subroutine init_file(filename, output_freq, output_units, time_units) character(len=*), intent(in) :: filename integer, intent(in) :: output_freq, output_units, time_units character(len=128) :: time_units_str real, dimension(1) :: tdata num_files = num_files + 1 if(num_files >= max_files) & call error_mesg('station_data, init_file', ' max_files exceeded, incease max_files', FATAL) files(num_files)%name = trim(filename) files(num_files)%output_freq = output_freq files(num_files)%output_units = output_units files(num_files)%time_units = time_units files(num_files)%num_fields = 0 files(num_files)%last_flush = base_time files(num_files)%file_unit = -1 !---- register axis_id and time boundaries id write(time_units_str, 11) trim(time_unit_list(files(num_files)%time_units)), base_year, & base_month, base_day, base_hour, base_minute, base_second 11 format(a, ' since ', i4.4, '-', i2.2, '-', i2.2, ' ', i2.2, ':', i2.2, ':', i2.2) files(num_files)%time_axis_id = diag_axis_init ('Time', tdata, time_units_str, 'T', & 'Time' , set_name=trim(filename)) files(num_files)%time_bounds_id = diag_axis_init('nv',(/1.,2./),'none','N','vertex number',& set_name=trim(filename)) end subroutine init_file !-------------------------------------------------------------------------- subroutine init_output_field(module_name,field_name,file_name,time_method,pack) character(len=*), intent(in) :: module_name, field_name, file_name character(len=*), intent(in) :: time_method integer, intent(in) :: pack integer :: out_num, file_num,num_fields, method_selected, l1 character(len=8) :: t_method ! Get a number for this output field num_output_fields = num_output_fields + 1 if(num_output_fields > max_output_fields) & call error_mesg('station_data', 'max_output_fields exceeded, increase it via nml', FATAL) out_num = num_output_fields file_num = find_file(file_name) if(file_num < 0) & call error_mesg('station_data,init_output_field', 'file '//trim(file_name) & //' is NOT found in station_data_table', FATAL) ! Insert this field into list of fields of this file files(file_num)%num_fields = files(file_num)%num_fields + 1 if(files(file_num)%num_fields > max_fields_per_file) & call error_mesg('station_data, init_output_field', 'max_fields_per_file exceeded ', FATAL) num_fields = files(file_num)%num_fields files(file_num)%fields(num_fields) = out_num output_fields(out_num)%output_name = trim(field_name) output_fields(out_num)%module_name = trim(module_name) output_fields(out_num)%counter = 0 output_fields(out_num)%output_file = file_num output_fields(out_num)%pack = pack output_fields(out_num)%time_average = .false. output_fields(out_num)%time_min = .false. output_fields(out_num)%time_max = .false. output_fields(out_num)%time_ops = .false. output_fields(out_num)%register = .false. t_method = lowercase(time_method) select case (trim(t_method)) case('.true.') output_fields(out_num)%time_average = .true. output_fields(out_num)%time_method = 'mean' case('mean') output_fields(out_num)%time_average = .true. output_fields(out_num)%time_method = 'mean' case('average') output_fields(out_num)%time_average = .true. output_fields(out_num)%time_method = 'mean' case('avg') output_fields(out_num)%time_average = .true. output_fields(out_num)%time_method = 'mean' case('.false.') output_fields(out_num)%time_average = .false. output_fields(out_num)%time_method = 'point' case ('max') call error_mesg('station_data, init_output_field','time_method MAX is not supported',& FATAL) output_fields(out_num)%time_max = .true. output_fields(out_num)%time_method = 'max' l1 = len_trim(output_fields(out_num)%output_name) if(output_fields(out_num)%output_name(l1-2:l1) /= 'max') & output_fields(out_num)%output_name = trim(field_name)//'_max' case ('min') call error_mesg('station_data, init_output_field','time_method MIN is not supported',& FATAL) output_fields(out_num)%time_min = .true. output_fields(out_num)%time_method = 'min' l1 = len_trim(output_fields(out_num)%output_name) if(output_fields(out_num)%output_name(l1-2:l1) /= 'min') & output_fields(out_num)%output_name = trim(field_name)//'_min' case default call error_mesg('station_data, init_output_field', 'error in time_method of field '& //trim(field_name), FATAL) end select if (files(file_num)%output_freq == EVERY_TIME) & output_fields(out_num)%time_average = .false. output_fields(out_num)%time_ops = output_fields(out_num)%time_min.or.output_fields(out_num)%time_max & .or.output_fields(out_num)%time_average output_fields(out_num)%time_method = trim(time_method) end subroutine init_output_field !-------------------------------------------------------------------------- function find_file(name) integer :: find_file character(len=*), intent(in) :: name integer :: i find_file = -1 do i = 1, num_files if(trim(files(i)%name) == trim(name)) then find_file = i return end if end do end function find_file ! ! ! ! This function is similar to register_diag_field of diag_manager_mod. All arguments ! are inputs that user needs to supply, some are optional. The names of input args are ! self-describing.
levels is absent for 2D fields.
! Note that pair (module_name, fieldname) must be unique in the ! station_data_table or a fatal error will occur.
! A field id is returned from this call that will be used later in send_station_data.
! ! !-------------------------------------------------------------------------- function register_station_field2d (module_name,fieldname,glo_lat,glo_lon,init_time, & domain,longname,units) integer :: register_station_field2d character(len=*), intent(in) :: module_name, fieldname real,dimension(:), intent(in) :: glo_lat,glo_lon type(domain2d), intent(in) :: domain type(time_type), intent(in) :: init_time character(len=*), optional, intent(in) :: longname, units real :: levels(1:1) levels = 0. register_station_field2d = register_station_field3d (module_name,fieldname,glo_lat,glo_lon,& levels,init_time,domain,longname,units) end function register_station_field2d !-------------------------------------------------------------------------- function register_station_field3d (module_name,fieldname,glo_lat,glo_lon,levels,init_time, & domain,longname,units) ! write field meta data on ROOT PE only ! allocate buffer integer :: register_station_field3d character(len=*), intent(in) :: module_name, fieldname real,dimension(:), intent(in) :: glo_lat,glo_lon,levels !in X,Y,Z direction respectively type(domain2d), intent(in) :: domain type(time_type), intent(in) :: init_time character(len=*), optional, intent(in) :: longname,units integer :: i,ii, nlat, nlon,nlevel, isc, iec, jsc, jec character(len=128) :: error_msg integer :: local_num_stations ! number of stations on this PE integer :: out_num ! position of this field in array output_fields integer :: file_num, freq, output_units, outunit real, allocatable :: station_values(:), level_values(:) character(len=128) :: longname2,units2 if(PRESENT(longname)) then longname2 = longname else longname2 = fieldname endif if(PRESENT(units)) then units2 = units else units2 = "none" endif nlat = size(glo_lat); nlon = size(glo_lon); nlevel=size(levels) allocate(station_values(num_stations), level_values(nlevel)) do i = 1, nlevel level_values(i) = real(i) enddo ! determine global index of this field in all stations outunit = stdout() do i = 1,num_stations station_values(i) = real(i) if(stations(i)%latglo_lat(nlat)) then write(error_msg,'(F9.3)') stations(i)%lat write(outunit,*) 'Station with latitude '//trim(error_msg)//' outside global latitude values' call error_mesg ('register_station_field', 'latitude out of range', FATAL) endif if(stations(i)%longlo_lon(nlon)) then write(error_msg,'(F9.3)') stations(i)%lon write(outunit,*) 'Station with longitude '//trim(error_msg)//' outside global longitude values' call error_mesg ('register_station_field', 'longitude out of range', FATAL) endif stations(i)%global_i = nearest_index(stations(i)%lon, glo_lon) stations(i)%global_j = nearest_index(stations(i)%lat, glo_lat) if(stations(i)%global_i<0 .or. stations(i)%global_j<0) & call error_mesg ('register_station_field', 'Error in global index of station',FATAL) enddo ! determine local index of this field in all stations , local index starts from 1 call mpp_get_compute_domain(domain, isc,iec,jsc,jec) local_num_stations = 0 do i = 1,num_stations if(isc<=stations(i)%global_i .and. iec>= stations(i)%global_i .and. & jsc<=stations(i)%global_j .and. jec>= stations(i)%global_j) then stations(i)%need_compute = .true. stations(i)%local_i = stations(i)%global_i - isc + 1 stations(i)%local_j = stations(i)%global_j - jsc + 1 local_num_stations = local_num_stations +1 endif enddo ! get the position of this field in the array output_fields out_num = find_output_field(module_name, fieldname) if(out_num < 0 .and. mpp_pe() == mpp_root_pe()) then call error_mesg ('register_station_field', & 'module/field_name '//trim(module_name)//'/'//& trim(fieldname)//' NOT found in station_data table', WARNING) register_station_field3d = out_num return endif if(local_num_stations>0) then allocate(output_fields(out_num)%station_id(local_num_stations)) allocate(output_fields(out_num)%buffer(local_num_stations,nlevel)) output_fields(out_num)%buffer = EMPTY ! fill out list of available stations in this PE ii=0 do i = 1,num_stations if(stations(i)%need_compute) then ii = ii+ 1 if(ii>local_num_stations) call error_mesg ('register_station_field', & 'error in determining local_num_station', FATAL) output_fields(out_num)%station_id(ii)=stations(i)%id endif enddo endif output_fields(out_num)%num_station = local_num_stations if( mpp_pe() == mpp_root_pe()) then allocate(global_field%buffer(num_stations,nlevel)) global_field%buffer = MISSING endif output_fields(out_num)%register = .true. output_fields(out_num)%output_name = fieldname file_num = output_fields(out_num)%output_file output_fields(out_num)%last_output = init_time freq = files(file_num)%output_freq output_units = files(file_num)%output_units output_fields(out_num)%next_output = diag_time_inc(init_time, freq, output_units) register_station_field3d = out_num output_fields(out_num)%long_name = longname2 output_fields(out_num)%units = units2 output_fields(out_num)%nlevel = nlevel ! deal with axes output_fields(out_num)%axes(1) = diag_axis_init('Stations',station_values,'station number', 'X') if(nlevel == 1) then output_fields(out_num)%num_axes = 1 else output_fields(out_num)%num_axes = 2 output_fields(out_num)%axes(2) = diag_axis_init('Levels',level_values,'level number', 'Y' ) endif if(need_write_axis) then lat_axis = diag_axis_init('Latitude', stations(1:num_stations)%lat,'station latitudes', 'n') lon_axis = diag_axis_init('Longitude',stations(1:num_stations)%lon, 'station longitudes', 'n') endif need_write_axis = .false. ! call mpp_sync() end function register_station_field3d !------------------------------------------------------------------------- function find_output_field(module_name, field_name) integer find_output_field character(len=*), intent(in) :: module_name, field_name integer :: i find_output_field = -1 do i = 1, num_output_fields if(trim(output_fields(i)%module_name) == trim(module_name) .and. & lowercase(trim(output_fields(i)%output_name)) == & lowercase(trim(field_name))) then find_output_field = i return endif end do end function find_output_field !------------------------------------------------------------------------- subroutine opening_file(file) ! open file, write axis meta_data for all files (only on ROOT PE, ! do nothing on other PEs) integer, intent(in) :: file character(len=128) :: time_units integer :: j,field_num,num_axes,axes(5),k logical :: time_ops integer :: time_axis_id(1),time_bounds_id(1) write(time_units, 11) trim(time_unit_list(files(file)%time_units)), base_year, & base_month, base_day, base_hour, base_minute, base_second 11 format(a, ' since ', i4.4, '-', i2.2, '-', i2.2, ' ', i2.2, ':', i2.2, ':', i2.2) call mpp_open(files(file)%file_unit, files(file)%name, action=MPP_OVERWR, & form=MPP_NETCDF, threading=MPP_SINGLE, fileset=MPP_SINGLE) call mpp_write_meta (files(file)%file_unit, 'title', cval=trim(global_descriptor)) time_ops = .false. do j = 1, files(file)%num_fields field_num = files(file)%fields(j) if(output_fields(field_num)%time_ops) then time_ops = .true. exit endif enddo !write axis meta data do j = 1, files(file)%num_fields field_num = files(file)%fields(j) num_axes = output_fields(field_num)%num_axes axes(1:num_axes) = output_fields(field_num)%axes(1:num_axes) do k = 1,num_axes if(axes(k)<0) & call error_mesg ('station_data opening_file','output_name '// & trim(output_fields(field_num)%output_name)// & ' has axis_id = -1', FATAL) enddo axes(num_axes + 1) = lat_axis axes(num_axes + 2) = lon_axis axes(num_axes + 3) = files(file)%time_axis_id ! need in write_axis: name, unit,long_name call write_axis_meta_data(files(file)%file_unit, axes(1:num_axes + 3), time_ops) if(time_ops) then axes(num_axes + 4) = files(file)%time_bounds_id call write_axis_meta_data(files(file)%file_unit, axes(1:num_axes + 4)) endif end do !write field meta data do j = 1, files(file)%num_fields field_num = files(file)%fields(j) num_axes = output_fields(field_num)%num_axes axes(1:num_axes) = output_fields(field_num)%axes(1:num_axes) num_axes = num_axes + 1 axes(num_axes) = files(file)%time_axis_id diag_field = write_field_meta_data(files(file)%file_unit,output_fields(field_num)%output_name, & axes(1:num_axes),output_fields(field_num)%units, & output_fields(field_num)%long_name,time_method=output_fields(field_num)%time_method,& pack=output_fields(field_num)%pack) output_fields(field_num)%f_type = diag_field%Field end do if(time_ops) then time_axis_id(1) = files(file)%time_axis_id time_bounds_id(1) = files(file)%time_bounds_id diag_field=write_field_meta_data(files(file)%file_unit, avg_name // '_T1',time_axis_id, & time_units,"Start time for average period", pack=1) files(file)%f_avg_start = diag_field%Field diag_field=write_field_meta_data(files(file)%file_unit,avg_name // '_T2' ,time_axis_id, & time_units,"End time for average period", pack=1) files(file)%f_avg_end = diag_field%Field diag_field=write_field_meta_data(files(file)%file_unit,avg_name // '_DT' ,time_axis_id, & time_units,"Length of average period", pack=1) files(file)%f_avg_nitems = diag_field%Field diag_field=write_field_meta_data(files(file)%file_unit, 'Time_bounds', (/time_bounds_id,time_axis_id/), & trim(time_unit_list(files(file)%time_units)), & 'Time axis boundaries', pack=1) files(file)%f_bounds = diag_field%Field endif call done_meta_data(files(file)%file_unit) end subroutine opening_file ! ! ! ! data should have the size of compute domain(isc:iec,jsc:jec)
! time is model's time
! field_id is returned from register_station_field
! only data at stations will be be sent to root_pe which, in turn, sends to output file ! ! !------------------------------------------------------------------------- subroutine send_station_data_2d(field_id, data, time) integer, intent(in) :: field_id real, intent(in) :: data(:,:) type(time_type), intent(in) :: time real :: data3d(size(data,1),size(data,2),1) data3d(:,:,1) = data call send_station_data_3d(field_id, data3d, time) end subroutine send_station_data_2d !------------------------------------------------------------------------- subroutine send_station_data_3d(field_id, data, time) integer, intent(in) :: field_id real, intent(in) :: data(:,:,:) type(time_type), intent(in) :: time integer :: freq,units,file_num,local_num_stations,i,ii, max_counter integer :: index_x, index_y, station_id integer, allocatable :: station_ids(:) ! root_pe only, to receive local station_ids real, allocatable :: tmp_buffer(:,:) ! root_pe only, to receive local buffer from each PE if (.not.module_is_initialized) & call error_mesg ('send_station_data_3d',' station_data NOT initialized', FATAL) if(field_id < 0) return file_num = output_fields(field_id)%output_file if( mpp_pe() == mpp_root_pe() .and. files(file_num)%file_unit < 0) then call opening_file(file_num) endif freq = files(file_num)%output_freq units = files(file_num)%output_units ! compare time with next_output if (time > output_fields(field_id)%next_output .and. freq /= END_OF_RUN) then ! time to write out ! ALL PEs, including root PE, must send data to root PE call mpp_send(output_fields(field_id)%num_station,plen=1,to_pe=mpp_root_pe()) if(output_fields(field_id)%num_station > 0) then call mpp_send(output_fields(field_id)%station_id(1),plen=size(output_fields(field_id)%station_id),& to_pe=mpp_root_pe()) call mpp_send(output_fields(field_id)%buffer(1,1),plen=size(output_fields(field_id)%buffer),& to_pe=mpp_root_pe()) endif ! get max_counter if the field is averaged if(output_fields(field_id)%time_average) then max_counter = output_fields(field_id)%counter call mpp_max(max_counter, pelist) endif ! receive local data from all PEs if(mpp_pe() == mpp_root_pe()) then do i = 1,size(pelist) call mpp_recv(local_num_stations,glen=1,from_pe=pelist(i)) if(local_num_stations> 0) then allocate(station_ids(local_num_stations)) allocate(tmp_buffer(local_num_stations,output_fields(field_id)%nlevel)) call mpp_recv(station_ids(1), glen=size(station_ids), from_pe=pelist(i)) call mpp_recv(tmp_buffer(1,1),glen=size(tmp_buffer), from_pe=pelist(i)) do ii = 1,local_num_stations global_field%buffer(station_ids(ii),:) = tmp_buffer(ii,:) enddo deallocate(station_ids, tmp_buffer) endif enddo ! send global_buffer content to file if(output_fields(field_id)%time_average) then if(max_counter == 0 ) & call error_mesg ('send_station_data','counter=0 for averaged field '// & output_fields(field_id)%output_name, FATAL) global_field%buffer = global_field%buffer/real(max_counter) endif ! check if global_field contains any missing values if(any(global_field%buffer == MISSING)) & call error_mesg ('send_station_data','Global_field contains MISSING, field '// & output_fields(field_id)%output_name, FATAL) call station_data_out(file_num,field_id,global_field%buffer,output_fields(field_id)%next_output) global_field%buffer = MISSING endif call mpp_sync() ! clear buffer, increment next_output time and reset counter on ALL PEs if(output_fields(field_id)%num_station>0) output_fields(field_id)%buffer = EMPTY output_fields(field_id)%last_output = output_fields(field_id)%next_output output_fields(field_id)%next_output = diag_time_inc(output_fields(field_id)%next_output,& freq, units) output_fields(field_id)%counter = 0; max_counter = 0 endif ! accumulate buffer only do i = 1 , output_fields(field_id)%num_station station_id = output_fields(field_id)%station_id(i) index_x = stations(station_id)%local_i; index_y = stations(station_id)%local_j if(index_x>size(data,1) .or. index_y>size(data,2)) & call error_mesg ('send_station_data','local index out of range for field '// & output_fields(field_id)%output_name, FATAL) if(output_fields(field_id)%time_average) then output_fields(field_id)%buffer(i,:) = output_fields(field_id)%buffer(i,:) + & data(index_x,index_y,:) ! accumulate buffer else ! not average output_fields(field_id)%buffer(i,:) = data(index_x,index_y,:) endif enddo if(output_fields(field_id)%time_average) & output_fields(field_id)%counter = output_fields(field_id)%counter + 1 end subroutine send_station_data_3d !------------------------------------------------------------------------ subroutine station_data_out(file, field, data, time,final_call_in) integer, intent(in) :: file, field real, intent(inout) :: data(:, :) type(time_type), intent(in) :: time logical, optional, intent(in):: final_call_in logical :: final_call integer :: i, num real :: dif, time_data(2, 1, 1), dt_time(1, 1, 1), start_dif, end_dif final_call = .false. if(present(final_call_in)) final_call = final_call_in dif = get_date_dif(time, base_time, files(file)%time_units) call mpp_write(files(file)%file_unit,output_fields(field)%f_type, data, dif) start_dif = get_date_dif(output_fields(field)%last_output, base_time,files(file)%time_units) end_dif = dif do i = 1, files(file)%num_fields num = files(file)%fields(i) if(output_fields(num)%time_ops) then if(num == field) then time_data(1, 1, 1) = start_dif call mpp_write(files(file)%file_unit, files(file)%f_avg_start, & time_data(1:1,:,:), dif) time_data(2, 1, 1) = end_dif call mpp_write(files(file)%file_unit, files(file)%f_avg_end, & time_data(2:2,:,:), dif) dt_time(1, 1, 1) = end_dif - start_dif call mpp_write(files(file)%file_unit, files(file)%f_avg_nitems, & dt_time(1:1,:,:), dif) ! Include boundary variable for CF compliance call mpp_write(files(file)%file_unit, files(file)%f_bounds, & time_data(1:2,:,:), dif) exit endif end if end do if(final_call) then if(time >= files(file)%last_flush) then call mpp_flush(files(file)%file_unit) files(file)%last_flush = time endif else if(time > files(file)%last_flush) then call mpp_flush(files(file)%file_unit) files(file)%last_flush = time endif endif end subroutine station_data_out ! ! ! ! Must be called after the last time step to write the buffer content ! ! !----------------------------------------------------------------------------- subroutine station_data_end(time) type(time_type), intent(in) :: time !model's time integer :: freq, max_counter, local_num_stations integer :: file, nfield, field, pe, col integer, allocatable :: station_ids(:) ! root_pe only, to receive local station_ids real, allocatable :: tmp_buffer(:,:) ! root_pe only, to receive local buffer from each PE do file = 1, num_files freq = files(file)%output_freq do nfield = 1, files(file)%num_fields field = files(file)%fields(nfield) if(.not. output_fields(field)%register) cycle if(time >= output_fields(field)%next_output .or. freq == END_OF_RUN) then ! ALL PEs, including root PE, must send data to root PE call mpp_send(output_fields(field)%num_station,plen=1,to_pe=mpp_root_pe()) if(output_fields(field)%num_station > 0) then call mpp_send(output_fields(field)%station_id(1),plen=size(output_fields(field)%station_id),& to_pe=mpp_root_pe()) call mpp_send(output_fields(field)%buffer(1,1),plen=size(output_fields(field)%buffer),& to_pe=mpp_root_pe()) endif ! get max_counter if the field is averaged if(output_fields(field)%time_average) then max_counter = output_fields(field)%counter call mpp_max(max_counter, pelist) endif ! only root PE receives local data from all PEs if(mpp_pe() == mpp_root_pe()) then do pe = 1,size(pelist) call mpp_recv(local_num_stations,glen=1,from_pe=pelist(pe)) if(local_num_stations> 0) then allocate(station_ids(local_num_stations)) allocate(tmp_buffer(local_num_stations,output_fields(field)%nlevel)) call mpp_recv(station_ids(1), glen=size(station_ids), from_pe=pelist(pe)) call mpp_recv(tmp_buffer(1,1),glen=size(tmp_buffer),from_pe=pelist(pe)) do col = 1,local_num_stations global_field%buffer(station_ids(col),:) = tmp_buffer(col,:) enddo deallocate(station_ids, tmp_buffer) endif enddo ! send global_buffer content to file if(output_fields(field)%time_average)then if(max_counter == 0 )& call error_mesg ('send_station_end','counter=0 for averaged field '// & output_fields(field)%output_name, FATAL) global_field%buffer = global_field%buffer/real(max_counter) endif ! check if global_field contains any missing values if(any(global_field%buffer == MISSING)) & call error_mesg ('send_station_end','Global_field contains MISSING, field '// & output_fields(field)%output_name, FATAL) call station_data_out(file,field,global_field%buffer,output_fields(field)%next_output,.true.) global_field%buffer = MISSING endif call mpp_sync() endif !deallocate field buffer if(output_fields(field)%num_station>0) & deallocate(output_fields(field)%buffer, output_fields(field)%station_id) enddo ! nfield enddo ! file if(mpp_pe() == mpp_root_pe()) deallocate(global_field%buffer) end subroutine station_data_end !----------------------------------------------------------------------------------- end module station_data_mod