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, open_namelist_file, check_nml_error use mpp_mod, only : mpp_npes, mpp_sync, mpp_root_pe, mpp_send, mpp_recv, mpp_max, & mpp_get_current_pelist, input_nml_file, & COMM_TAG_1, COMM_TAG_2, COMM_TAG_3, COMM_TAG_4 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, ierr 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 #ifdef INTERNAL_FILE_NML read (input_nml_file, station_data_nml, iostat=io_status) ierr = check_nml_error(io_status, 'station_data_nml') #else iunit = open_namelist_file () ierr=1; do while (ierr /= 0) read (iunit, nml=station_data_nml, iostat=io_status, end=10) ierr = check_nml_error(io_status, 'station_data_nml') enddo 10 call close_file (iunit) #endif logunit = stdlog() write(logunit, station_data_nml) 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(),tag=COMM_TAG_1) 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),tag=COMM_TAG_1) 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(),tag=COMM_TAG_2) 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(),tag=COMM_TAG_3) call mpp_send(output_fields(field)%buffer(1,1),plen=size(output_fields(field)%buffer),& to_pe=mpp_root_pe(),tag=COMM_TAG_4) 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),tag=COMM_TAG_2) 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),tag=COMM_TAG_3) call mpp_recv(tmp_buffer(1,1),glen=size(tmp_buffer),from_pe=pelist(pe),tag=COMM_TAG_4) 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