module diag_integral_mod #include ! ! fil ! ! ! ! ! diag_integral_mod computes and outputs global and / or ! hemispheric physics integrals. ! ! ! ! shared modules: use time_manager_mod, only: time_type, get_time, set_time, & time_manager_init, & operator(+), operator(-), & operator(==), operator(>=), & operator(/=) use mpp_mod, only: input_nml_file use fms_mod, only: open_file, file_exist, error_mesg, & open_namelist_file, check_nml_error, & fms_init, & mpp_pe, mpp_root_pe,& WARNING,& FATAL, write_version_number, & stdlog, close_file use constants_mod, only: radius, constants_init use mpp_mod, only: mpp_sum, mpp_init !-------------------------------------------------------------------- implicit none private !---------------------------------------------------------------------- ! diag_integral_mod computes and outputs global and / or ! hemispheric physics integrals. !---------------------------------------------------------------------- !--------------------------------------------------------------------- !----------- version number for this module ------------------- character(len=128) :: version = '$Id: diag_integral.F90,v 20.0 2013/12/13 23:17:05 fms Exp $' character(len=128) :: tagname = '$Name: tikal $' !--------------------------------------------------------------------- !------ interfaces ------ public & diag_integral_init, diag_integral_field_init, & sum_diag_integral_field, diag_integral_output, & diag_integral_end interface sum_diag_integral_field module procedure sum_field_2d, & sum_field_2d_hemi, & sum_field_3d, & sum_field_wght_3d end interface private & ! from diag_integral_init: set_axis_time, & ! from diag_integral_field_init and sum_diag_integral_field: get_field_index, & ! from diag_integral_output and diag_integral_end: write_field_averages, & ! from write_field_averages: format_text_init, format_data_init, & get_axis_time, & ! from diag_integral_output: diag_integral_alarm, & ! from sum_diag_integral_field: vert_diag_integral !--------------------------------------------------------------------- !------ namelist ------- integer, parameter :: & mxch = 64 ! maximum number of characters in ! the optional output file name real :: & output_interval = -1.0 ! time interval at which integrals ! are to be output character(len=8) :: & time_units = 'hours' ! time units associated with ! output_interval character(len=mxch) :: & file_name = ' ' ! optional integrals output file name logical :: & print_header = .true. ! print a header for the integrals ! file ? integer :: & fields_per_print_line = 4 ! number of fields to write per line ! of output namelist /diag_integral_nml/ output_interval, time_units, & file_name, print_header, & fields_per_print_line !--------------------------------------------------------------------- !------- public data ------ !--------------------------------------------------------------------- !------- private data ------ !--------------------------------------------------------------------- ! variables associated with the determination of when integrals ! are to be written. ! Next_alarm_time next time at which integrals are to be ! written ! Alarm_interval time interval between writing integrals ! Zero_time time_type variable set to (0,0); used as ! flag to indicate integrals are not being ! output ! Time_init_save initial time associated with experiment; ! used as a base for defining time !--------------------------------------------------------------------- type (time_type) :: Next_alarm_time, Alarm_interval, Zero_time type (time_type) :: Time_init_save !--------------------------------------------------------------------- ! variables used in determining weights associated with each ! contribution to the integrand. ! area area of each grid box ! idim x dimension of grid on local processor ! jdim y dimension of grid on local processor ! field_size number of columns on global domain ! sum_area surface area of globe !--------------------------------------------------------------------- real, allocatable, dimension(:,:) :: area integer :: idim, jdim, field_size real :: sum_area !--------------------------------------------------------------------- ! variables used to define the integral fields: ! max_len_name maximum length of name associated with integral ! max_num_field maximum number of integrals allowed ! num_field number of integrals that have been activated ! field_name(i) name associated with integral i ! field_format(i) output format for integral i ! field_sum(i) integrand for integral i ! field_count(i) number of values in integrand i !--------------------------------------------------------------------- integer, parameter :: max_len_name = 12 integer, parameter :: max_num_field = 32 integer :: num_field = 0 character(len=max_len_name) :: field_name (max_num_field) character(len=16) :: field_format (max_num_field) real :: field_sum (max_num_field) integer :: field_count (max_num_field) !--------------------------------------------------------------------- ! variables defining output formats. ! format_text format statement for header ! format_data format statement for data output ! do_format_data a data format needs to be generated ? ! nd number of characters in data format statement ! nt number of characters in text format statement !--------------------------------------------------------------------- character(len=160) :: format_text, format_data logical :: do_format_data = .true. integer :: nd, nt !-------------------------------------------------------------------- ! miscellaneous variables. !--------------------------------------------------------------------- integer :: diag_unit = 0 ! unit number for output file logical :: module_is_initialized = .false. ! module is initialized ? !----------------------------------------------------------------------- !----------------------------------------------------------------------- contains !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! PUBLIC SUBROUTINES ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !#################################################################### ! ! ! diag_integral_init is the constructor for diag_integral_mod. ! ! ! diag_integral_init is the constructor for diag_integral_mod. ! ! ! ! Initial time to start the integral ! ! ! current time ! ! ! array of model latitudes at cell boundaries [radians] ! ! ! array of model longitudes at cell boundaries [radians] ! ! ! subroutine diag_integral_init (Time_init, Time, blon, blat, area_in) !-------------------------------------------------------------------- ! diag_integral_init is the constructor for diag_integral_mod. !-------------------------------------------------------------------- type (time_type), intent(in), optional :: Time_init, Time real,dimension(:,:), intent(in), optional :: blon, blat, area_in !-------------------------------------------------------------------- ! intent(in),optional variables: ! ! Time_init ! Time ! blon ! blat ! !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! local variables: real :: rsize integer :: unit, io, ierr, nc, logunit integer :: field_size_local real :: sum_area_local !--------------------------------------------------------------------- ! local variables: ! ! r2 ! rsize ! unit ! io ! ierr ! seconds ! nc ! i,j ! !-------------------------------------------------------------------- !--------------------------------------------------------------------- ! if routine has already been executed, exit. !--------------------------------------------------------------------- if (module_is_initialized) return !--------------------------------------------------------------------- ! verify that modules used by this module that are not called later ! have already been initialized. !--------------------------------------------------------------------- call fms_init call mpp_init call constants_init call time_manager_init !---------------------------------------------------------------------- ! if this is the initialization call, proceed. if this was simply ! a verification of previous initialization, return. !-------------------------------------------------------------------- if (present(Time_init) .and. present(Time) .and. & present(blon) .and. present(blat) ) then !----------------------------------------------------------------------- ! read namelist. !----------------------------------------------------------------------- if ( file_exist('input.nml')) then #ifdef INTERNAL_FILE_NML read (input_nml_file, nml=diag_integral_nml, iostat=io) ierr = check_nml_error(io,'diag_integral_nml') #else unit = open_namelist_file ( ) ierr=1; do while (ierr /= 0) read (unit, nml=diag_integral_nml, iostat=io, end=10) ierr = check_nml_error(io,'diag_integral_nml') end do 10 call close_file (unit) #endif endif !--------------------------------------------------------------------- ! write version number and namelist to logfile. !--------------------------------------------------------------------- call write_version_number(version, tagname) logunit = stdlog() if (mpp_pe() == mpp_root_pe() ) & write (logunit, nml=diag_integral_nml) !-------------------------------------------------------------------- ! save the initial time to time-stamp the integrals which will be ! calculated. !--------------------------------------------------------------------- Time_init_save = Time_init !--------------------------------------------------------------------- ! define the model grid sizes and the total number of columns on ! the processor. sum over all processors and store the global ! number of columns in field_size. !--------------------------------------------------------------------- idim = size(blon,1) - 1 jdim = size(blon,2) - 1 field_size_local = idim*jdim rsize = real(field_size_local) call mpp_sum (rsize) field_size = nint(rsize) !--------------------------------------------------------------------- ! define an array to hold the surface area of each grid column ! so that the integrals may be weighted properly. sum over the ! processor, and then over all processors, storing the total ! global surface area in sum_area. !--------------------------------------------------------------------- allocate (area(idim,jdim)) area = area_in sum_area_local = sum(area) sum_area = sum_area_local call mpp_sum (sum_area) !-------------------------------------------------------------------- ! if integral output is to go to a file, open the file on unit ! diag_unit. !-------------------------------------------------------------------- if (file_name(1:1) /= ' ' ) then nc = len_trim(file_name) diag_unit = open_file (file_name(1:nc), action='write') endif !--------------------------------------------------------------------- ! define the variables needed to control the time interval of ! output. Zero time is a flag indicating that the alarm is not set, ! i.e., integrals are not desired. otherwise set the next time to ! output integrals to be at the value of nml variable ! output_interval from now. !--------------------------------------------------------------------- Zero_time = set_time (0,0) if (output_interval >= -0.01) then Alarm_interval = set_axis_time (output_interval, time_units) Next_alarm_time = Time + Alarm_interval else Alarm_interval = Zero_time endif Next_alarm_time = Time + Alarm_interval !-------------------------------------------------------------------- ! deallocate the local array and mark the module as initialized. !-------------------------------------------------------------------- module_is_initialized = .true. endif ! (present optional arguments) !----------------------------------------------------------------------- end subroutine diag_integral_init !###################################################################### ! ! ! diag_integral_field_init registers and intializes an integral field ! ! ! diag_integral_field_init registers and intializes an integral field ! ! ! ! Name of the field to be integrated ! ! ! Output format of the field to be integrated ! ! ! subroutine diag_integral_field_init (name, format) !--------------------------------------------------------------------- ! !--------------------------------------------------------------------- character(len=*), intent(in) :: name, format !--------------------------------------------------------------------- ! intent(in) variables: ! ! name ! format ! !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! local variables: integer :: field ! index assigned to the current integral !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! note: no initialization is required for this interface. all needed ! variables are initialized in the source. !--------------------------------------------------------------------- !-------------------------------------------------------------------- ! make sure the integral name is not too long. !-------------------------------------------------------------------- if (len(name) > max_len_name ) then call error_mesg ('diag_integral_mod', & ' integral name too long', FATAL) endif !--------------------------------------------------------------------- ! check to be sure the integral name has not already been ! initialized. !--------------------------------------------------------------------- field = get_field_index (name) if (field /= 0) then call error_mesg ('diag_integral_mod', & 'integral name already exists', FATAL) endif !------------------------------------------------------------------- ! prepare to register the integral. make sure that there are not ! more integrals registered than space was provided for; if so, exit. !---------------------------------------------------------------------- num_field = num_field + 1 if (num_field > max_num_field) then call error_mesg ('diag_integral_mod', & 'too many fields initialized', FATAL) endif !-------------------------------------------------------------------- ! register the name and output format desired for the given integral. ! initialize its value and the number of grid points that have been ! counted to zero. !-------------------------------------------------------------------- field_name (num_field) = name field_format (num_field) = format field_sum (num_field) = 0.0 field_count (num_field) = 0 !---------------------------------------------------------------------- end subroutine diag_integral_field_init !##################################################################### ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! INTERFACE SUM_DIAG_INTEGRAL_FIELD ! ! call sum_diag_integral_field (name, data, is, js) ! or ! call sum_diag_integral_field (name, data, wt, is, js) ! or ! call sum_diag_integral_field (name, data, is, ie, js, je) ! ! in the first option data may be either ! real, intent(in) :: data(:,:) [ sum_field_2d ] ! or ! real, intent(in) :: data(:,:,:) [ sum_field_3d ] ! !------------------------------------------------------------------- ! intent(in) arguments: ! ! character(len=*), intent(in) :: name ! real, intent(in) :: wt(:,:,:) ! integer, optional, intent(in) :: is, ie, js, je ! !-------------------------------------------------------------------- ! intent(in) arguments: ! ! name name associated with integral ! data field of integrands to be summed over ! wt vertical weighting factor to be applied to integrands ! when summing ! is,ie,js,je starting/ending i,j indices over which summation is ! to occur ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !##################################################################### ! ! ! Perform a 2 dimensional summation of named field ! ! ! Perform a 2 dimensional summation of named field ! ! ! ! Name of the field to be integrated ! ! ! field of integrands to be summed over ! ! ! starting i,j indices over which summation is ! to occur ! ! ! subroutine sum_field_2d (name, data, is, js) character(len=*), intent(in) :: name real, intent(in) :: data(:,:) integer, optional, intent(in) :: is, js !--------------------------------------------------------------------- ! local variables: integer :: field ! index of desired integral integer :: i1, j1, i2, j2 ! location indices of current data in ! processor-global coordinates !---------------------------------------------------------------------- ! be sure module has been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized ) then call error_mesg ('diag_integral_mod', & 'module has not been initialized', FATAL ) endif !--------------------------------------------------------------------- ! obtain the index of the current integral. make certain it is valid. !--------------------------------------------------------------------- field = get_field_index (name) if (field == 0) then call error_mesg ('diag_integral_mod', & 'field does not exist', FATAL) endif !--------------------------------------------------------------------- ! define the processor-global indices of the current data. use the ! value 1 for the initial grid points, if is and js are not input. !--------------------------------------------------------------------- i1 = 1; if (present(is)) i1 = is j1 = 1; if (present(js)) j1 = js i2 = i1 + size(data,1) - 1 j2 = j1 + size(data,2) - 1 !--------------------------------------------------------------------- ! increment the count of points toward this integral and add the ! values at this set of grid points to the accumulation array. !--------------------------------------------------------------------- !$OMP CRITICAL field_count (field) = field_count(field) + & size(data,1)*size(data,2) field_sum (field) = field_sum (field) + & sum (data * area(i1:i2,j1:j2)) !$OMP END CRITICAL !-------------------------------------------------------------------- end subroutine sum_field_2d !####################################################################### ! ! ! Perform a 3 dimensional summation of named field ! ! ! Perform a 3 dimensional summation of named field ! ! ! ! Name of the field to be integrated ! ! ! field of integrands to be summed over ! ! ! starting i,j indices over which summation is ! to occur ! ! ! subroutine sum_field_3d (name, data, is, js) character(len=*), intent(in) :: name real, intent(in) :: data(:,:,:) integer, optional, intent(in) :: is, js !--------------------------------------------------------------------- ! local variables: real, dimension (size(data,1), & size(data,2)) :: data2 integer :: field integer :: i1, j1, i2, j2 !--------------------------------------------------------------------- ! local variables: ! ! data2 ! field ! index of desired integral ! i1, j1, i2, j2 ! location indices of current data in ! processor-global coordinates ! !-------------------------------------------------------------------- !---------------------------------------------------------------------- ! be sure module has been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized ) then call error_mesg ('diag_integral_mod', & 'module has not been initialized', FATAL ) endif !--------------------------------------------------------------------- ! obtain the index of the current integral. make certain it is valid. !--------------------------------------------------------------------- field = get_field_index (name) if (field == 0) then call error_mesg ('diag_integral_mod', & 'field does not exist', FATAL) endif !--------------------------------------------------------------------- ! define the processor-global indices of the current data. use the ! value 1 for the initial grid points, if is and js are not input. !--------------------------------------------------------------------- i1 = 1; if (present(is)) i1 = is j1 = 1; if (present(js)) j1 = js i2 = i1 + size(data,1) - 1 j2 = j1 + size(data,2) - 1 !--------------------------------------------------------------------- ! increment the count of points toward this integral. sum first ! in the vertical and then add the values at this set of grid points ! to the accumulation array. !--------------------------------------------------------------------- !$OMP CRITICAL field_count (field) = field_count (field) + & size(data,1)*size(data,2) data2 = sum(data,3) field_sum (field) = field_sum (field) + & sum (data2 * area(i1:i2,j1:j2)) !$OMP END CRITICAL !--------------------------------------------------------------------- end subroutine sum_field_3d !####################################################################### ! ! ! Perform a 3 dimensional weighted summation of named field ! ! ! Perform a 3 dimensional weighted summation of named field ! ! ! ! Name of the field to be integrated ! ! ! field of integrands to be summed over ! ! ! the weight function to be evaluated at summation ! ! ! starting i,j indices over which summation is ! to occur ! ! ! subroutine sum_field_wght_3d (name, data, wt, is, js) character(len=*), intent(in) :: name real, intent(in) :: data(:,:,:), wt(:,:,:) integer, optional, intent(in) :: is, js !--------------------------------------------------------------------- ! local variables: real, dimension (size(data,1),size(data,2)) :: data2 integer :: field, i1, j1, i2, j2 !--------------------------------------------------------------------- ! local variables: ! ! data2 ! field ! index of desired integral ! i1, j1, i2, j2 ! location indices of current data in ! processor-global coordinates ! !-------------------------------------------------------------------- !---------------------------------------------------------------------- ! be sure module has been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized ) then call error_mesg ('diag_integral_mod', & 'module has not been initialized', FATAL ) endif !--------------------------------------------------------------------- ! obtain the index of the current integral. make certain it is valid. !--------------------------------------------------------------------- field = get_field_index (name) if (field == 0) then call error_mesg ('diag_integral_mod', & 'field does not exist', FATAL) endif !--------------------------------------------------------------------- ! define the processor-global indices of the current data. use the ! value 1 for the initial grid points, if is and js are not input. !--------------------------------------------------------------------- i1 = 1; if (present(is)) i1 = is j1 = 1; if (present(js)) j1 = js i2 = i1 + size(data,1) - 1 j2 = j1 + size(data,2) - 1 !--------------------------------------------------------------------- ! increment the count of points toward this integral. sum first ! in the vertical (including a vertical weighting factor) and then ! add the values at this set of grid points to the accumulation ! array. !--------------------------------------------------------------------- !$OMP CRITICAL field_count (field) = field_count (field) + & size(data,1)*size(data,2) data2 = vert_diag_integral (data, wt) field_sum(field) = field_sum (field) + & sum (data2 * area(i1:i2,j1:j2)) !$OMP END CRITICAL !---------------------------------------------------------------------- end subroutine sum_field_wght_3d !####################################################################### ! ! ! Perform a 2 dimensional hemispherical summation of named field ! ! ! Perform a 2 dimensional hemispherical summation of named field ! ! ! ! Name of the field to be integrated ! ! ! field of integrands to be summed over ! ! ! starting/ending i,j indices over which summation is ! to occur ! ! ! subroutine sum_field_2d_hemi (name, data, is, ie, js, je) character(len=*), intent(in) :: name real, intent(in) :: data(:,:) integer, intent(in) :: is, js, ie, je !--------------------------------------------------------------------- ! local variables: integer :: field, i1, j1, i2, j2 !--------------------------------------------------------------------- ! local variables: ! ! field ! index of desired integral ! i1, j1, i2, j2 ! location indices of current data in ! processor-global coordinates ! !-------------------------------------------------------------------- !---------------------------------------------------------------------- ! be sure module has been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized ) then call error_mesg ('diag_integral_mod', & 'module has not been initialized', FATAL ) endif !--------------------------------------------------------------------- ! obtain the index of the current integral. make certain it is valid. !--------------------------------------------------------------------- field = get_field_index (name) if (field == 0) then call error_mesg ('diag_integral_mod', & 'field does not exist', FATAL) endif !---------------------------------------------------------------------- ! define the processor-global indices of the current data. this form ! is needed to handle case of 2d domain decomposition with physics ! window smaller than processor domain size. !---------------------------------------------------------------------- i1 = mod ( (is-1), size(data,1) ) + 1 i2 = i1 + size(data,1) - 1 !-------------------------------------------------------------------- ! for a hemispheric sum, sum one jrow at a time in case a processor ! has data from both hemispheres. !-------------------------------------------------------------------- j1 = mod ( (js-1) ,size(data,2) ) + 1 j2 = j1 !---------------------------------------------------------------------- ! increment the count of points toward this integral. include hemi- ! spheric factor of 2 in field_count. add the data values at this ! set of grid points to the accumulation array. !---------------------------------------------------------------------- !$OMP CRITICAL field_count (field) = field_count (field) + 2* (i2-i1+1)*(j2-j1+1) field_sum (field) = field_sum (field) + & sum (data(i1:i2,j1:j2)*area(is:ie,js:je)) !$OMP END CRITICAL !--------------------------------------------------------------------- end subroutine sum_field_2d_hemi !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! END INTERFACE SUM_DIAG_INTEGRAL_FIELD ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !################################################################## ! ! ! diag_integral_output determines if this is a timestep on which ! integrals are to be written. if not, it returns; if so, it calls ! write_field_averages. ! ! ! diag_integral_output determines if this is a timestep on which ! integrals are to be written. if not, it returns; if so, it calls ! write_field_averages. ! ! ! ! integral time stamp at the current time ! ! ! subroutine diag_integral_output (Time) !--------------------------------------------------------------------- ! diag_integral_output determines if this is a timestep on which ! integrals are to be written. if not, it returns; if so, it calls ! write_field_averages. !--------------------------------------------------------------------- type (time_type), intent(in) :: Time !----------------------------------------------------------------------- ! intent(in) variables: ! ! Time integral time stamp at the current time ! [ time_type ] ! !--------------------------------------------------------------------- !---------------------------------------------------------------------- ! be sure module has been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized ) then call error_mesg ('diag_integral_mod', & 'module has not been initialized', FATAL ) endif !--------------------------------------------------------------------- ! see if integral output is desired at this time. !--------------------------------------------------------------------- if ( diag_integral_alarm(Time) ) then !--------------------------------------------------------------------- ! write the integrals by calling write_field_averages. upon return ! reset the alarm to the next diagnostics time. !--------------------------------------------------------------------- call write_field_averages (Time) Next_alarm_time = Next_alarm_time + Alarm_interval endif !----------------------------------------------------------------------- end subroutine diag_integral_output !####################################################################### ! ! ! diag_integral_end is the destructor for diag_integral_mod. ! ! ! diag_integral_end is the destructor for diag_integral_mod. ! ! ! ! integral time stamp at the current time ! ! ! subroutine diag_integral_end (Time) !-------------------------------------------------------------------- ! diag_integral_end is the destructor for diag_integral_mod. !-------------------------------------------------------------------- type (time_type), intent(in) :: Time !---------------------------------------------------------------------- ! be sure module has been initialized. !--------------------------------------------------------------------- if (.not. module_is_initialized ) then call error_mesg ('diag_integral_mod', & 'module has not been initialized', FATAL ) endif !--------------------------------------------------------------------- ! if the alarm interval was set to Zero_time (meaning no integral ! output during the model run) call write_field_averages to output ! the integrals valid over the entire period of integration. !--------------------------------------------------------------------- if (Alarm_interval == Zero_time ) then ! if (Alarm_interval /= Zero_time ) then ! else call write_field_averages (Time) endif !--------------------------------------------------------------------- ! deallocate module variables. !--------------------------------------------------------------------- deallocate (area) !--------------------------------------------------------------------- ! mark the module as uninitialized. !--------------------------------------------------------------------- module_is_initialized = .false. !-------------------------------------------------------------------- end subroutine diag_integral_end !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! PRIVATE SUBROUTINES ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !####################################################################### ! ! ! Function to convert input time to a time_type ! ! ! Function to convert input time to a time_type ! ! ! ! integral time stamp at the current time ! ! ! input units, not used ! ! ! function set_axis_time (atime, units) result (Time) !-------------------------------------------------------------------- ! !-------------------------------------------------------------------- real, intent(in) :: atime character(len=*), intent(in) :: units type(time_type) :: Time !--------------------------------------------------------------------- ! intent(in) variables: ! ! atime ! units ! ! result: ! ! Time ! !-------------------------------------------------------------------- !--------------------------------------------------------------------- ! local variables: integer :: sec ! seconds corresponding to the input ! variable atime integer :: day = 0 ! day component of time_type variable !-------------------------------------------------------------------- ! convert the input time to seconds, regardless of input units. !-------------------------------------------------------------------- if (units(1:3) == 'sec') then sec = int(atime + 0.5) else if (units(1:3) == 'min') then sec = int(atime*60. + 0.5) else if (units(1:3) == 'hou') then sec = int(atime*3600. + 0.5) else if (units(1:3) == 'day') then sec = int(atime*86400. + 0.5) else call error_mesg('diag_integral_mod', & 'Invalid units sent to set_axis_time', FATAL) endif !-------------------------------------------------------------------- ! convert the time in seconds to a time_type variable. !-------------------------------------------------------------------- Time = set_time (sec, day) end function set_axis_time !###################################################################### ! ! ! get_field_index returns returns the index associated with an ! integral name. ! ! ! get_field_index returns returns the index associated with an ! integral name. ! ! ! ! Name associated with an integral ! ! ! function get_field_index (name) result (index) !--------------------------------------------------------------------- ! get_field_index returns returns the index associated with an ! integral name. !--------------------------------------------------------------------- character(len=*), intent(in) :: name integer :: index !-------------------------------------------------------------------- ! intent(in) variables: ! ! name ! ! result: ! ! index ! !-------------------------------------------------------------------- !--------------------------------------------------------------------- ! local variables: integer :: nc integer :: i !--------------------------------------------------------------------- ! !-------------------------------------------------------------------- nc = len_trim (name) if (nc > max_len_name) then call error_mesg ('diag_integral_mod', & 'name too long', FATAL) endif !-------------------------------------------------------------------- ! search each field name for the current string. when found exit ! with the index. if not found index will be 0 upon return, which ! initiates error condition. !-------------------------------------------------------------------- index = 0 do i = 1, num_field if (name(1:nc) == & field_name(i) (1:len_trim(field_name(i))) ) then index = i exit endif end do !--------------------------------------------------------------------- end function get_field_index !##################################################################### ! ! ! Subroutine to sum multiple fields, average them and then write the result ! to an output file. ! ! ! Subroutine to sum multiple fields, average them and then write the result ! to an output file. ! ! ! ! integral time stamp at the current time ! ! ! subroutine write_field_averages (Time) !--------------------------------------------------------------------- ! !--------------------------------------------------------------------- type (time_type), intent(in) :: Time !-------------------------------------------------------------------- ! intent(in) variables: ! ! Time ! !-------------------------------------------------------------------- !-------------------------------------------------------------------- ! local variables: real :: field_avg(max_num_field) logical :: field_skip(max_num_field) real :: xtime, rcount integer :: nn, ninc, nst, nend, fields_to_print integer :: i, kount integer(LONG_KIND) :: icount !-------------------------------------------------------------------- ! local variables: ! ! field_avg ! field_skip ! xtime ! rcount ! nn ! ninc ! nst ! nend ! fields_to_print ! i ! kount ! !-------------------------------------------------------------------- !-------------------------------------------------------------------- ! each header and data format may be different and must be generated ! as needed. !---------------------------------------------------------------------- fields_to_print = 0 do i = 1, num_field !-------------------------------------------------------------------- ! increment the fields_to_print counter. sum the integrand and the ! number of data points contributing to it over all processors. !-------------------------------------------------------------------- fields_to_print = fields_to_print + 1 rcount = real(field_count(i)) call mpp_sum (rcount) call mpp_sum (field_sum(i)) icount = rcount field_skip(i) = .false. !-------------------------------------------------------------------- ! verify that all the data expected for an integral has been ! obtained. !-------------------------------------------------------------------- if (icount == 0 ) then call error_mesg & ('diag_integral_mod', & 'field_count equals zero for field_name ' // & field_name(i)(1:len_trim(field_name(i))), WARNING ) field_skip(i) = .true. end if kount = icount/field_size if ((field_size*1.0)*kount /= rcount) then print*,"name,pe,kount,field_size,icount,rcount=",trim(field_name(i)),mpp_pe(),kount,field_size,icount,rcount call error_mesg & ('diag_integral_mod', & 'field_count not a multiple of field_size', WARNING ) field_skip(i) = .true. endif !---------------------------------------------------------------------- ! define the global integral for field i. reinitialize the point ! and data accumulators. !---------------------------------------------------------------------- if (field_skip(i)) then field_avg(fields_to_print) = 9e33 ! missing value else field_avg(fields_to_print) = field_sum(i)/ & (sum_area*float(kount)) end if field_sum (i) = 0.0 field_count(i) = 0 end do !-------------------------------------------------------------------- ! only the root pe will write out data. !-------------------------------------------------------------------- if ( mpp_pe() /= mpp_root_pe() ) return !--------------------------------------------------------------------- ! define the time associated with the integrals just calculated. !--------------------------------------------------------------------- xtime = get_axis_time (Time-Time_init_save, time_units) !--------------------------------------------------------------------- ! generate the new header and data formats. !--------------------------------------------------------------------- nst = 1 nend = fields_per_print_line ninc = (num_field-1)/fields_per_print_line + 1 do nn=1, ninc nst = 1 + (nn-1)*fields_per_print_line nend = MIN (nn*fields_per_print_line, num_field) if (print_header) call format_text_init (nst, nend) call format_data_init (nst, nend) if (diag_unit /= 0) then write (diag_unit,format_data(1:nd)) & xtime, (field_avg(i),i=nst,nend) else write (*, format_data(1:nd)) & xtime, (field_avg(i),i=nst,nend) endif end do !----------------------------------------------------------------------- end subroutine write_field_averages !####################################################################### ! ! ! format_text_init generates the header records to be output in the ! integrals file. ! ! ! format_text_init generates the header records to be output in the ! integrals file. ! ! ! ! starting/ending integral index which will be included ! in this format statement ! ! ! subroutine format_text_init (nst_in, nend_in) !---------------------------------------------------------------------- ! format_text_init generates the header records to be output in the ! integrals file. !---------------------------------------------------------------------- integer, intent(in), optional :: nst_in, nend_in !--------------------------------------------------------------------- ! intent(in),optional variables: ! ! nst_in starting integral index which will be included ! in this format statement ! nend_in ending integral index which will be included ! in this format statement ! !-------------------------------------------------------------------- !-------------------------------------------------------------------- ! local variables: integer :: i, nc, nst, nend !-------------------------------------------------------------------- ! local variables: ! ! i ! nc ! nst ! nend ! !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! only the root pe need execute this routine, since only it will ! be outputting integrals. !--------------------------------------------------------------------- if (mpp_pe() /= mpp_root_pe()) return !---------------------------------------------------------------------- ! define the starting and ending integral indices that will be ! included in this format statement. !---------------------------------------------------------------------- if (present (nst_in) ) then nst = nst_in nend = nend_in else nst = 1 nend = num_field endif !-------------------------------------------------------------------- ! define the first 11 characters in the format statement. !-------------------------------------------------------------------- nt = 11 format_text(1:nt) = "('# time" !-------------------------------------------------------------------- ! generate the rest of the format statement, which will cover ! integral indices nst to nend. if satndard printout is desired, ! cycle through the loop. !-------------------------------------------------------------------- do i=nst,nend nc = len_trim(field_name(i)) format_text(nt+1:nt+nc+5) = ' ' // field_name(i)(1:nc) nt = nt+nc+5 end do !--------------------------------------------------------------------- ! include the end of the format statement. !--------------------------------------------------------------------- format_text(nt+1:nt+2) = "')" nt = nt+2 !-------------------------------------------------------------------- ! write the format statement to either an output file or to stdout. !-------------------------------------------------------------------- if (diag_unit /= 0) then write (diag_unit, format_text(1:nt)) else write (*, format_text(1:nt)) endif !--------------------------------------------------------------------- end subroutine format_text_init !####################################################################### ! ! ! format_text_init generates the format to be output in the ! integrals file. ! ! ! format_text_init generates the format to be output in the ! integrals file. ! ! ! ! starting/ending integral index which will be included ! in this format statement ! ! ! subroutine format_data_init (nst_in, nend_in) !--------------------------------------------------------------------- ! format_data_init generates the format that will write out the ! integral data. !--------------------------------------------------------------------- integer, intent(in), optional :: nst_in, nend_in !-------------------------------------------------------------------- ! intent(in),optional variables: ! ! nst_in starting integral index which will be included ! in this format statement ! nend_in ending integral index which will be included ! in this format statement ! !-------------------------------------------------------------------- !-------------------------------------------------------------------- ! local variables: integer :: i, nc, nst, nend !-------------------------------------------------------------------- ! local variables: ! ! i ! nc ! nst ! nend ! !--------------------------------------------------------------------- !-------------------------------------------------------------------- ! define the start of the format, which covers the time stamp of the ! integrals. this section is 9 characters long. !-------------------------------------------------------------------- nd = 9 format_data(1:nd) = '(1x,f10.2' !-------------------------------------------------------------------- ! define the indices of the integrals that are to be written by this ! format statement. !-------------------------------------------------------------------- if ( present (nst_in) ) then nst = nst_in nend = nend_in else nst = 1 nend = num_field endif !------------------------------------------------------------------- ! complete the data format. use the format defined for the ! particular integral in setting up the format statement. !------------------------------------------------------------------- do i=nst,nend nc = len_trim(field_format(i)) format_data(nd+1:nd+nc+5) = ',1x,' // field_format(i)(1:nc) nd = nd+nc+5 end do !------------------------------------------------------------------- ! close the format statement. !------------------------------------------------------------------- format_data(nd+1:nd+1) = ')' nd = nd + 1 !------------------------------------------------------------------- end subroutine format_data_init !####################################################################### ! ! ! Function to convert the time_type input variable into units of ! units and returns it in atime. ! ! ! Function to convert the time_type input variable into units of ! units and returns it in atime. ! ! ! ! integral time stamp ! ! ! input units of time_type ! ! ! function get_axis_time (Time, units) result (atime) !--------------------------------------------------------------------- ! !--------------------------------------------------------------------- type(time_type), intent(in) :: Time character(len=*), intent(in) :: units real :: atime !---------------------------------------------------------------------- ! intent(in) variables: ! ! Time ! units ! ! result: ! ! atime ! !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! local variables: integer :: sec, day ! components of time_type variable !------------------------------------------------------------------- ! get_axis_time converts the time_type input variable into units of ! units and returns it in atime. !------------------------------------------------------------------- call get_time (Time, sec, day) if (units(1:3) == 'sec') then atime = float(sec) + 86400.*float(day) else if (units(1:3) == 'min') then atime = float(sec)/60. + 1440.*float(day) else if (units(1:3) == 'hou') then atime = float(sec)/3600. + 24.*float(day) else if (units(1:3) == 'day') then atime = float(sec)/86400. + float(day) endif !-------------------------------------------------------------------- end function get_axis_time !##################################################################### ! ! ! Function to check if it is time to write integrals. ! if not writing integrals, return. ! ! ! Function to check if it is time to write integrals. ! if not writing integrals, return. ! ! ! ! current time ! ! ! function diag_integral_alarm (Time) result (answer) !-------------------------------------------------------------------- ! !-------------------------------------------------------------------- type (time_type), intent(in) :: Time logical :: answer !--------------------------------------------------------------------- ! intent(in) variables: ! ! Time ! ! result: ! ! answer ! !--------------------------------------------------------------------- !-------------------------------------------------------------------- ! check if it is time to write integrals. if not writing integrals, ! return. !-------------------------------------------------------------------- answer = .false. if (Alarm_interval == Zero_time) return if (Time >= Next_alarm_time) answer = .true. !-------------------------------------------------------------------- end function diag_integral_alarm !####################################################################### ! ! ! Function to perform a weighted integral in the vertical ! direction of a 3d data field ! ! ! Function to perform a weighted integral in the vertical ! direction of a 3d data field ! ! ! ! integral field data arrays ! ! ! integral field weighting functions ! ! ! function vert_diag_integral (data, wt) result (data2) !---------------------------------------------------------------------- ! !---------------------------------------------------------------------- real, dimension (:,:,:), intent(in) :: data, wt real, dimension (size(data,1),size(data,2)) :: data2 !--------------------------------------------------------------------- ! intent(in) variables; ! ! data ! wt ! ! result: ! data2 ! !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! local variables: real, dimension(size(data,1),size(data,2)) :: wt2 !--------------------------------------------------------------------- ! local variables: ! ! wt2 ! !--------------------------------------------------------------------- !-------------------------------------------------------------------- wt2 = sum(wt,3) if (count(wt2 == 0.) > 0) then call error_mesg ('diag_integral_mod', & 'vert sum of weights equals zero', FATAL) endif data2 = sum(data*wt,3) / wt2 !--------------------------------------------------------------------- end function vert_diag_integral !####################################################################### end module diag_integral_mod