!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| module diagnostics !BOP ! !MODULE: diagnostics ! !DESCRIPTION: ! This module contains routines and data for producing diagnostic ! output, including global diagnostics, cfl checks and transport ! diagnostics. ! ! !REVISION HISTORY: ! SVN:$Id: diagnostics.F90 73619 2015-09-28 23:51:23Z mvertens $ ! ! !USES: use POP_KindsMod use POP_IOUnitsMod use POP_SolversMod use domain use constants use prognostic use time_management use io use broadcast use global_reductions use grid use forcing use forcing_fields use timers use exit_mod use tavg, only: define_tavg_field, & accumulate_tavg_field, accumulate_tavg_now, & tavg_method_avg, tavg_method_max, tavg_method_min use movie, only: define_movie_field, movie_requested, update_movie_field ! QL, 150526, new mixed layer depth HMXL_DR use vmix_kpp, only: HMXL, HMXL_DR, KPP_HBLT use registry use io_tools use gather_scatter implicit none private save ! !PUBLIC MEMBER FUNCTIONS: public :: init_diagnostics, & diag_init_sums, & diag_global_preupdate, & diag_global_afterupdate, & diag_print, & diag_transport, & cfl_advect, & cfl_vdiff, & cfl_hdiff, & cfl_check, & check_KE, & diag_velocity ! !PUBLIC DATA MEMBERS: logical (log_kind), public :: & ldiag_global, &! time to compute global diagnostics ldiag_cfl, &! time to compute cfl diagnostics ldiag_transport, &! time to compute transport diagnostics ldiag_velocity, &! compute velocity diagnostics cfl_all_levels, &! writes cfl diags for all vert levels diag_all_levels ! writes some diags for all vert levels !*** arrays for holding various diagnostic results !*** public for now as they are modified directly in baroclinic real (r8), dimension(nx_block,ny_block,max_blocks_clinic), & public :: & DIAG_KE_ADV_2D, &! DIAG_KE_HMIX_2D, &! DIAG_KE_VMIX_2D, &! DIAG_KE_PRESS_2D, &! DIAG_PE_2D real (r8), dimension(nx_block,ny_block,nt,max_blocks_clinic), & public :: & DIAG_TRACER_ADV_2D, &! DIAG_TRACER_HDIFF_2D, &! DIAG_TRACER_VDIFF_2D, &! DIAG_TRACER_SOURCE_2D, &! DIAG_TRACER_SFC_FLX !----------------------------------------------------------------------- ! ! global budget diagnostics ! !----------------------------------------------------------------------- real (r8), dimension(nt),public :: & tracer_mean_initial, &! SUM [volume*tracer] at the beginning ! of a time-averaging interval tracer_mean_final ! SUM [volume*tracer] at the end of a ! time-averaging interval real (r8), public :: & volume_t_initial, & ! T-point volume at the beginning of a ! time-averaging interval volume_t_final ! T-point volume at the end of a ! time-averaging interval !EOP !BOC !----------------------------------------------------------------------- ! ! variables for controlling frequency and output of diagnostics ! !----------------------------------------------------------------------- integer (int_kind) :: & diag_unit, &! i/o unit for output diagnostic file trans_unit, &! i/o unit for output transport file velocity_unit ! i/o unit for output velocity file integer (int_kind) :: & diag_global_flag, &! time flag id for global diags diag_cfl_flag, &! time flag id for cfl diags diag_transp_flag ! time flag id for transp diags character (char_len) :: & diag_outfile, &! current filename for diagnostic output diag_transport_outfile, &! current filename for transport output diag_velocity_outfile ! current filename for velocity output !----------------------------------------------------------------------- ! ! variables specific to ccsm coupled code ! !----------------------------------------------------------------------- logical (log_kind) :: & lccsm = .false. character (char_len) :: & ccsm_diag_date !----------------------------------------------------------------------- ! ! global diagnostics ! !----------------------------------------------------------------------- real (r8) :: & diag_ke, &! mean KE at new time diag_ke_adv, &! KE change due to advection diag_ke_free, &! advective KE change due to free surface diag_ke_hmix, &! KE change due to horizontal diffusion diag_ke_vmix, &! KE change due to vertical diffusion diag_ke_press, &! KE change due to pressure gradient diag_pe, &! change in potential energy diag_press_free, &! pressure work contribution from free surface diag_ke_psfc, &! KE change due to surface pressure diag_ws, &! KE change due to wind stress diag_sealevel ! global mean seal level real (r8), dimension(nt) :: & dtracer_avg, &! avg tracer change for this time step dtracer_abs, &! abs tracer change for this time step diag_tracer_adv, &! tracer change due to advection diag_tracer_hdiff, &! tracer change due to horiz diffusion diag_tracer_vdiff, &! tracer change due to vert diffusion diag_tracer_source, &! tracer change due to source terms avg_tracer, &! global average tracer at new time sfc_tracer_flux ! sfc tracer flux (Fw*tracer) or ! resid sfc flux (w*tracer) real (r8), dimension(km,nt) :: & avg_tracer_k ! mean tracer at every level !----------------------------------------------------------------------- ! ! cfl diagnostics ! !----------------------------------------------------------------------- real (r8), dimension(max_blocks_clinic,km) :: & cfl_advuk_block, &! zonal advective cfl number for level k cfl_advvk_block, &! merid advective cfl number for level k cfl_advwk_block, &! vert advective cfl number for level k cfl_advtk_block, &! total advective cfl number for level k cfl_vdifftk_block, &! tracer vert diff cfl num for level k cfl_vdiffuk_block, &! momentum vert diff cfl num for level k cfl_hdifftk_block, &! tracer horiz diff cfl num for level k cfl_hdiffuk_block ! momentum horiz diff cfl num for level k integer (int_kind), dimension(2,max_blocks_clinic,km) :: & cfladd_advuk_block, &! horiz addresses for max cfl numbers cfladd_advvk_block, &! at level k cfladd_advwk_block, & cfladd_advtk_block, & cfladd_vdifftk_block, & cfladd_vdiffuk_block, & cfladd_hdifftk_block, & cfladd_hdiffuk_block !----------------------------------------------------------------------- ! ! transport diagnostics ! !----------------------------------------------------------------------- type :: transport character(char_len) :: name ! name for transport integer (int_kind) :: type ! type (meridional or zonal) integer (int_kind), dimension(6,max_blocks_clinic) :: & add ! address range for computing transp real (r8) :: mass, &! mass transport heat, &! heat transport salt ! salt transport end type integer (int_kind), parameter :: & transport_type_zonal = 1, &! zonal or meridional transport type transport_type_merid = 2 integer (int_kind) :: & num_transports ! number of transport diagnostics computed type (transport), dimension(:), allocatable :: & transports ! transport info for all requested transports !----------------------------------------------------------------------- ! ! ids for tavg diagnostics computed from diagnostic ! !----------------------------------------------------------------------- integer (int_kind) :: & tavg_HMXL, &! tavg id for average mixed layer depth tavg_HMXL_2, &! tavg id for average mixed layer depth, stream #2 (allows two frequencies) tavg_HMXL_DR, &! tavg id for average mixed layer depth with density criterion, QL, 150526 tavg_XMXL, &! tavg id for maximum mixed layer depth tavg_XMXL_2, &! tavg id for maximum mixed layer depth, stream #2 tavg_TMXL, &! tavg id for minimum mixed layer depth tavg_HBLT, &! tavg id for average boundary layer depth tavg_XBLT, &! tavg id for maximum boundary layer depth tavg_TBLT ! tavg id for minimum boundary layer depth !----------------------------------------------------------------------- ! ! movie ids ! !----------------------------------------------------------------------- integer (int_kind) :: & movie_HMXL, &! movie id for average mixed layer depth movie_HBLT ! movie id for average boundary layer depth !----------------------------------------------------------------------- ! ! solver diagnostics ! !----------------------------------------------------------------------- integer (POP_i4) :: & iterationCount ! current solver iteration count real (POP_r8) :: & rmsResidual ! mean solver residual !----------------------------------------------------------------------- ! ! velocity diagnostics ! !----------------------------------------------------------------------- integer (int_kind), parameter :: & num_vel_loc = 5 integer (int_kind), dimension (num_vel_loc) :: & i_loc, j_loc real (r8), dimension (num_vel_loc) :: & vlat_loc, vlon_loc !EOC !*********************************************************************** contains !*********************************************************************** !BOP ! !IROUTINE: init_diagnostics ! !INTERFACE: subroutine init_diagnostics ! !DESCRIPTION: ! Initializes diagnostic module quantities. ! ! !REVISION HISTORY: ! same as module !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (int_kind) :: & nu, &! unit for contents input file iblock, n, &! dummy loop indices nloc, nreach, &! ditto nml_error, &! namelist i/o error flag ib,ie,jb,je ! beg,end indices for block physical domain integer (int_kind), dimension(6) :: & tmp_add ! transport global addresses ! (ibeg,iend,jbeg,jend,kbeg,kend) character (char_len) :: & diag_global_freq_opt, &! choice for freq of global diagnostics diag_cfl_freq_opt, &! choice for freq of cfl diagnostics diag_transp_freq_opt, &! choice for freq of transport diagnostics diag_transport_file, &! filename for choosing fields for output transport_ctype, &! type of transport (zonal,merid) outfile_tmp, &! temp for appending to outfile name string ! temp for creating outfile name integer (int_kind) :: & diag_global_freq_iopt, &! freq option for computing global diags diag_global_freq, &! freq for computing global diagnostics diag_cfl_freq_iopt, &! freq option for computing cfl diags diag_cfl_freq, &! freq for computing cfl diagnostics diag_transp_freq_iopt, &! freq option for comp transport diags diag_transp_freq ! freq for computing transp diagnostics character (10) :: & cdate ! character date character (47), parameter :: &! output formats for freq options gfreq_fmt = "('Global diagnostics computed every ',i8,a8)", & cfreq_fmt = "('CFL diagnostics computed every ',i8,a8)", & tfreq_fmt = "('Transport diagnostics computed every ',i8,a8)" namelist /diagnostics_nml/diag_global_freq_opt, diag_global_freq, & diag_cfl_freq_opt, diag_cfl_freq, & diag_transp_freq_opt, diag_transp_freq, & diag_transport_file, & diag_all_levels, cfl_all_levels, & diag_outfile, diag_transport_outfile, & diag_velocity_outfile, & ldiag_velocity type (block) :: & this_block ! block information for current block !----------------------------------------------------------------------- ! ! local variables for velocity diagnostics lat/lon search ! !----------------------------------------------------------------------- integer (int_kind), parameter :: search_error = -999 real (r8) :: & diag_vel_search_reach, &! size of lat/lon search length diag_vel_search_inc, &! increment in " " lat_reach, &! local search-length names lon_reach, & min_distance, &!smallest distance from desired lat/lon point distance !actual distance from desired lat/lon point integer (int_kind) :: & i,j, &! dummy indices diag_vel_search_max real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: WORK real (r8), allocatable, dimension (:,:) :: & ULAT_G,ULON_G integer (int_kind), allocatable, dimension(:,:) :: & REGION_MASK_G real (r8), parameter, dimension (num_vel_loc) :: & lon_loc = (/ 220.0_r8, 220.0_r8, 335.0_r8, 335.0_r8, 250.0_r8 /) real (r8), parameter, dimension (num_vel_loc) :: & lat_loc = (/ 0.0_r8, 2.0_r8, 0.0_r8, 2.0_r8, 0.0_r8 /) !----------------------------------------------------------------------- ! ! determine if this is a ccsm coupled run ! !----------------------------------------------------------------------- lccsm = registry_match('lccsm') !----------------------------------------------------------------------- ! ! read diagnostic file output frequency and filenames from namelist ! !----------------------------------------------------------------------- !*** set default values diag_global_freq_opt = 'never' diag_cfl_freq_opt = 'never' diag_transp_freq_opt = 'never' diag_transport_file = 'unknown_transport_file' diag_outfile = 'unknown_diag_outfile' diag_transport_outfile = 'unknown_transport_outfile' diag_all_levels = .false. cfl_all_levels = .false. ldiag_velocity = .false. diag_velocity_outfile = 'unknown_velocity_outfile' if (my_task == master_task) then open (nml_in, file=nml_filename, status='old',iostat=nml_error) if (nml_error /= 0) then nml_error = -1 else nml_error = 1 endif do while (nml_error > 0) read(nml_in, nml=diagnostics_nml,iostat=nml_error) end do if (nml_error == 0) close(nml_in) endif call broadcast_scalar(nml_error, master_task) if (nml_error /= 0) then call exit_POP(sigAbort,'ERROR reading diagnostics_nml') endif !----------------------------------------------------------------------- ! ccsm filenames must meet ccsm output-file naming conventions ! actual ccsm filenames will be constructed in diag_print !----------------------------------------------------------------------- if (lccsm) then call ccsm_date_stamp (ccsm_diag_date, 'ymds') string = diag_outfile diag_outfile = trim(string)/& &/'.'/& &/trim(ccsm_diag_date) string = diag_transport_outfile diag_transport_outfile = trim(string)/& &/'.'/& &/trim(ccsm_diag_date) string = diag_velocity_outfile diag_velocity_outfile = trim(string)/& &/'.'/& &/trim(ccsm_diag_date) else !----------------------------------------------------------------------- ! non-ccsm filenames: ! append runid, initial date to output file names ! concatenation operator must be split across lines to avoid problems ! with preprocessors ! !----------------------------------------------------------------------- if (date_separator == ' ') then cdate(1:4) = cyear cdate(5:6) = cmonth cdate(7:8) = cday cdate(9:10)= ' ' else cdate(1:4) = cyear cdate(5:5) = date_separator cdate(6:7) = cmonth cdate(8:8) = date_separator cdate(9:10) = cday endif outfile_tmp = char_blank outfile_tmp = trim(diag_outfile)/& &/'.'/& &/trim(runid)/& &/'.'/& &/trim(cdate) diag_outfile = char_blank diag_outfile = trim(outfile_tmp) outfile_tmp = char_blank outfile_tmp = trim(diag_transport_outfile)/& &/'.'/& &/trim(runid)/& &/'.'/& &/trim(cdate) diag_transport_outfile = char_blank diag_transport_outfile = trim(outfile_tmp) endif ! lccsm !----------------------------------------------------------------------- ! ! set, broadcast and print options to log file ! !----------------------------------------------------------------------- if (my_task == master_task) then write(stdout,blank_fmt) write(stdout,delim_fmt) write(stdout,blank_fmt) write(stdout,'(a18)') 'Diagnostic options' write(stdout,blank_fmt) write(stdout,delim_fmt) select case (diag_global_freq_opt) case ('never') diag_global_freq_iopt = freq_opt_never write(stdout,'(a31)') 'Global diagnostics not computed' case ('nyear') diag_global_freq_iopt = freq_opt_nyear write(stdout,gfreq_fmt) diag_global_freq, ' years ' case ('nmonth') diag_global_freq_iopt = freq_opt_nmonth write(stdout,gfreq_fmt) diag_global_freq, ' months ' case ('nday') diag_global_freq_iopt = freq_opt_nday write(stdout,gfreq_fmt) diag_global_freq, ' days ' case ('nhour') diag_global_freq_iopt = freq_opt_nhour write(stdout,gfreq_fmt) diag_global_freq, ' hours ' case ('nsecond') diag_global_freq_iopt = freq_opt_nsecond write(stdout,gfreq_fmt) diag_global_freq, ' seconds' case ('nstep') diag_global_freq_iopt = freq_opt_nstep write(stdout,gfreq_fmt) diag_global_freq, ' steps ' case default diag_global_freq_iopt = -1000 end select if (diag_global_freq_iopt /= freq_opt_never) then write(stdout,'(a36,a)') & 'Global diagnostics written to file: ', trim(diag_outfile) if (diag_all_levels) & write(stdout,'(a42)') & 'Diagnostics output for all vertical levels' endif if (lccsm) then write (stdout,'(a39,a)') 'Equatorial velocities written to file: ', & trim(diag_velocity_outfile) endif select case (diag_cfl_freq_opt) case ('never') diag_cfl_freq_iopt = freq_opt_never write(stdout,'(a28)') 'CFL diagnostics not computed' case ('nyear') diag_cfl_freq_iopt = freq_opt_nyear write(stdout,cfreq_fmt) diag_cfl_freq, ' years ' case ('nmonth') diag_cfl_freq_iopt = freq_opt_nmonth write(stdout,cfreq_fmt) diag_cfl_freq, ' months ' case ('nday') diag_cfl_freq_iopt = freq_opt_nday write(stdout,cfreq_fmt) diag_cfl_freq, ' days ' case ('nhour') diag_cfl_freq_iopt = freq_opt_nhour write(stdout,cfreq_fmt) diag_cfl_freq, ' hours ' case ('nsecond') diag_cfl_freq_iopt = freq_opt_nsecond write(stdout,cfreq_fmt) diag_cfl_freq, ' seconds' case ('nstep') diag_cfl_freq_iopt = freq_opt_nstep write(stdout,cfreq_fmt) diag_cfl_freq, ' steps ' case default diag_cfl_freq_iopt = -1000 end select if (diag_cfl_freq_iopt /= freq_opt_never) then write(stdout,'(a33,a)') & 'CFL diagnostics written to file: ',trim(diag_outfile) if (cfl_all_levels) then write(stdout,'(a46)') & 'CFL diagnostics output for all vertical levels' endif endif select case (diag_transp_freq_opt) case ('never') diag_transp_freq_iopt = freq_opt_never write(stdout,'(a34)') 'Transport diagnostics not computed' case ('nyear') diag_transp_freq_iopt = freq_opt_nyear write(stdout,tfreq_fmt) diag_transp_freq, ' years ' case ('nmonth') diag_transp_freq_iopt = freq_opt_nmonth write(stdout,tfreq_fmt) diag_transp_freq, ' months ' case ('nday') diag_transp_freq_iopt = freq_opt_nday write(stdout,tfreq_fmt) diag_transp_freq, ' days ' case ('nhour') diag_transp_freq_iopt = freq_opt_nhour write(stdout,tfreq_fmt) diag_transp_freq, ' hours ' case ('nsecond') diag_transp_freq_iopt = freq_opt_nsecond write(stdout,tfreq_fmt) diag_transp_freq, ' seconds' case ('nstep') diag_transp_freq_iopt = freq_opt_nstep write(stdout,tfreq_fmt) diag_transp_freq, ' steps ' case default diag_transp_freq_iopt = -1000 end select endif call broadcast_scalar(diag_global_freq_iopt, master_task) call broadcast_scalar(diag_global_freq, master_task) call broadcast_scalar(diag_cfl_freq_iopt, master_task) call broadcast_scalar(diag_cfl_freq, master_task) call broadcast_scalar(diag_transp_freq_iopt, master_task) call broadcast_scalar(diag_transp_freq, master_task) call broadcast_scalar(diag_all_levels, master_task) call broadcast_scalar(cfl_all_levels, master_task) call broadcast_scalar(ldiag_velocity, master_task) if (diag_global_freq_iopt == -1000) then call exit_POP(sigAbort, & 'ERROR: unknown global diag frequency option') endif if (diag_cfl_freq_iopt == -1000) then call exit_POP(sigAbort, & 'ERROR: unknown cfl diagnostic frequency option') endif if (diag_transp_freq_iopt == -1000) then call exit_POP(sigAbort, & 'ERROR: unknown transport diag frequency option') endif !----------------------------------------------------------------------- ! ! create time flags for computing diagnostics ! !----------------------------------------------------------------------- call init_time_flag('diag_global', diag_global_flag, default=.false., & owner = 'init_diagnostics', & freq_opt = diag_global_freq_iopt, & freq = diag_global_freq) call init_time_flag('diag_cfl', diag_cfl_flag, default=.false., & owner = 'init_diagnostics', & freq_opt = diag_cfl_freq_iopt, & freq = diag_cfl_freq) call init_time_flag('diag_transp', diag_transp_flag, default=.false., & owner = 'init_diagnostics', & freq_opt = diag_transp_freq_iopt, & freq = diag_transp_freq) !----------------------------------------------------------------------- ! ! read transports file to determine which transports to compute ! !----------------------------------------------------------------------- if (diag_transp_freq_iopt /= freq_opt_never) then call get_unit(nu) if (my_task == master_task) then write(stdout,'(a39,a)') & 'Transport diagnostics written to file: ', & trim(diag_transport_outfile) write(stdout,blank_fmt) open(nu, file=diag_transport_file, status='old') read(nu,*) num_transports write(stdout,'(a14,i4,1x,a27)') 'The following ',num_transports, & 'transports will be computed' endif call broadcast_scalar(num_transports, master_task) allocate( transports(num_transports) ) do n=1,num_transports if (my_task == master_task) then read(nu,'(6(i4,1x),a5,2x,a)') tmp_add, & transport_ctype, transports(n)%name write(stdout,'(a2,a)') ' ',trim(transports(n)%name) select case (transport_ctype(1:5)) case ('zonal') transports(n)%type = transport_type_zonal case ('merid') transports(n)%type = transport_type_merid case default transports(n)%type = -1000 end select endif call broadcast_scalar(transports(n)%type,master_task) call broadcast_array(tmp_add,master_task) if (transports(n)%type == -1000) then call exit_POP(sigAbort,'ERROR: unknown transport type') endif !*** !*** compute local addresses for transport !*** !$OMP PARALLEL DO PRIVATE(iblock,this_block,ib,ie,jb,je) do iblock=1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) ib = this_block%ib ie = this_block%ie jb = this_block%jb je = this_block%je transports(n)%add(1,iblock) = nx_block transports(n)%add(2,iblock) = 1 transports(n)%add(3,iblock) = ny_block transports(n)%add(4,iblock) = 1 transports(n)%add(5,iblock) = tmp_add(5) ! beg k index transports(n)%add(6,iblock) = tmp_add(6) ! end k index if (tmp_add(1) <= this_block%i_glob(ie) .and. & tmp_add(2) >= this_block%i_glob(ib) ) then if (tmp_add(1) >= this_block%i_glob(ib)) then transports(n)%add(1,iblock) = & tmp_add(1) - this_block%i_glob(ib) + ib else transports(n)%add(1,iblock) = ib endif if (tmp_add(2) <= this_block%i_glob(ie)) then transports(n)%add(2,iblock) = & tmp_add(2) - this_block%i_glob(ib) + ib else transports(n)%add(2,iblock) = ie endif endif if (tmp_add(3) <= this_block%j_glob(je) .and. & tmp_add(4) >= this_block%j_glob(jb) ) then if (tmp_add(3) >= this_block%j_glob(jb)) then transports(n)%add(3,iblock) = & tmp_add(3) - this_block%j_glob(jb) + jb else transports(n)%add(3,iblock) = jb endif if (tmp_add(4) <= this_block%j_glob(je)) then transports(n)%add(4,iblock) = & tmp_add(4) - this_block%j_glob(jb) + jb else transports(n)%add(4,iblock) = je endif endif end do ! block loop !$OMP END PARALLEL DO end do if (my_task == master_task) close(nu) call release_unit(nu) endif !----------------------------------------------------------------------- ! ! open output files if required, then write a blank and close. ! !----------------------------------------------------------------------- if (diag_global_freq_iopt /= freq_opt_never .or. & diag_cfl_freq_iopt /= freq_opt_never) then call get_unit(diag_unit) if (my_task == master_task) then open(diag_unit, file=diag_outfile, status='unknown') write(diag_unit,*)' ' close(diag_unit) endif endif if (diag_transp_freq_iopt /= freq_opt_never) then call get_unit(trans_unit) if (my_task == master_task) then open(trans_unit, file=diag_transport_outfile,status='unknown') write(trans_unit,*)' ' close(trans_unit) endif endif if ( lccsm ) then call get_unit(velocity_unit) if (my_task == master_task) then open(velocity_unit, file=diag_velocity_outfile, status='unknown') write(velocity_unit,*)' ' close(velocity_unit) endif endif !----------------------------------------------------------------------- ! ! initialize cfl arrays if required ! !----------------------------------------------------------------------- if (diag_cfl_freq_iopt /= freq_opt_never) then cfl_advuk_block = c0 cfl_advvk_block = c0 cfl_advwk_block = c0 cfl_advtk_block = c0 cfl_vdifftk_block = c0 cfl_vdiffuk_block = c0 cfl_hdifftk_block = c0 cfl_hdiffuk_block = c0 cfladd_advuk_block = 0 cfladd_advvk_block = 0 cfladd_advwk_block = 0 cfladd_advtk_block = 0 cfladd_vdifftk_block = 0 cfladd_vdiffuk_block = 0 cfladd_hdifftk_block = 0 cfladd_hdiffuk_block = 0 endif !----------------------------------------------------------------------- ! ! define tavg diagnostic fields ! !----------------------------------------------------------------------- call define_tavg_field(tavg_HMXL,'HMXL',2, & tavg_method=tavg_method_avg, & long_name='Mixed-Layer Depth', & units='centimeter', grid_loc='2110', & coordinates='TLONG TLAT time') call define_tavg_field(tavg_HMXL_2,'HMXL_2',2, & tavg_method=tavg_method_avg, & long_name='Mixed-Layer Depth', & units='centimeter', grid_loc='2110', & coordinates='TLONG TLAT time') ! QL, 150526, mixed layer depth with density criterion call define_tavg_field(tavg_HMXL_DR,'HMXL_DR',2, & tavg_method=tavg_method_avg, & long_name='Mixed-Layer Depth (density)', & units='centimeter', grid_loc='2110', & coordinates='TLONG TLAT time') call define_tavg_field(tavg_XMXL,'XMXL',2, & tavg_method=tavg_method_max, & long_name='Maximum Mixed-Layer Depth', & units='centimeter', grid_loc='2110', & coordinates='TLONG TLAT time') call define_tavg_field(tavg_XMXL_2,'XMXL_2',2, & tavg_method=tavg_method_max, & long_name='Maximum Mixed-Layer Depth', & units='centimeter', grid_loc='2110', & coordinates='TLONG TLAT time') call define_tavg_field(tavg_TMXL,'TMXL',2, & tavg_method=tavg_method_min, & long_name='Minimum Mixed-Layer Depth', & units='centimeter', grid_loc='2110', & coordinates='TLONG TLAT time') call define_tavg_field(tavg_HBLT,'HBLT',2, & tavg_method=tavg_method_avg, & long_name='Boundary-Layer Depth', & units='centimeter', grid_loc='2110', & coordinates='TLONG TLAT time') call define_tavg_field(tavg_XBLT,'XBLT',2, & tavg_method=tavg_method_max, & long_name='Maximum Boundary-Layer Depth', & units='centimeter', grid_loc='2110', & coordinates='TLONG TLAT time') call define_tavg_field(tavg_TBLT,'TBLT',2, & tavg_method=tavg_method_min, & long_name='Minimum Boundary-Layer Depth', & units='centimeter', grid_loc='2110', & coordinates='TLONG TLAT time') !----------------------------------------------------------------------- ! ! define movie diagnostic fields ! !----------------------------------------------------------------------- call define_movie_field(movie_HMXL,'HMXL',0, & long_name='Mixed-Layer Depth', & units='cm', grid_loc='2110') call define_movie_field(movie_HBLT,'HBLT',0, & long_name='Boundary-Layer Depth', & units='cm', grid_loc='2110') !----------------------------------------------------------------------- ! ! find global i,j indices close to specified lat,lon points ! !----------------------------------------------------------------------- if (ldiag_velocity) then j_loc = search_error i_loc = search_error allocate (ULAT_G(nx_global,ny_global), ULON_G(nx_global,ny_global)) allocate (REGION_MASK_G(nx_global,ny_global)) WORK=ULON*radian call gather_global(ULON_G,WORK,master_task,distrb_clinic) WORK=ULAT*radian call gather_global(ULAT_G,WORK,master_task,distrb_clinic) call gather_global(REGION_MASK_G,REGION_MASK,master_task,distrb_clinic) diag_vel_search_reach = c0 diag_vel_search_inc = 0.1_r8 diag_vel_search_max = 10 if (my_task == master_task) then do nloc = 1, num_vel_loc min_distance = c10 lat_reach = diag_vel_search_reach lon_reach = diag_vel_search_reach reach_loop: do nreach = 1,diag_vel_search_max lat_reach = lat_reach + diag_vel_search_inc lon_reach = lon_reach + diag_vel_search_inc do j=1,ny_global do i=1,nx_global if ( (lat_loc(nloc) - lat_reach <= ULAT_G(i,j) ) .and. & (lat_loc(nloc) + lat_reach >= ULAT_G(i,j) ) .and. & (lon_loc(nloc) - lon_reach <= ULON_G(i,j) ) .and. & (lon_loc(nloc) + lon_reach >= ULON_G(i,j) ) ) then !----------------------------------------------------- ! accept point only if it is an ocean point !----------------------------------------------------- !--------------------------------------------------------- !*** note that REGION_MASK_G is on the t-grid; should test !*** on the u-grid !--------------------------------------------------------- if (REGION_MASK_G(i,j) /= 0) then !--------------------------------------------------------- !*** note that this is crude approximation; will be changed !*** in a more permanent version at a later time distance = sqrt( (lat_loc(nloc)-ULAT_G(i,j))**2 + & (lon_loc(nloc)-ULON_G(i,j))**2 ) !--------------------------------------------------------- if (distance < min_distance) then j_loc(nloc) = j i_loc(nloc) = i min_distance = distance vlat_loc(nloc) = ULAT_G(i,j) vlon_loc(nloc) = ULON_G(i,j) endif endif ! ocean point endif enddo ! i enddo ! j enddo reach_loop enddo ! nloc write(stdout,*) ' ' write(stdout,*) ' Velocity diagnostics will be computed at the following locations:' do nloc = 1, num_vel_loc write(stdout,'(2x,a,i3,a,i3,a, 3x, a,F25.15,a,F25.15,a)' ) & '(', i_loc(nloc), ',', j_loc(nloc),')', & '(',vlat_loc(nloc), ',',vlon_loc(nloc),')' enddo call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout) endif call broadcast_array(j_loc , master_task) call broadcast_array(i_loc , master_task) do nloc = 1, num_vel_loc if ( j_loc(nloc) == search_error .or. i_loc(nloc) == search_error) & call exit_POP (sigAbort, 'velocity diagnostics lat/lon search failed') enddo deallocate (ULAT_G, ULON_G) deallocate (REGION_MASK_G) endif !ldiag_velocity if (my_task == master_task) then write(stdout,blank_fmt) write(stdout,delim_fmt) write(stdout,blank_fmt) write(stdout,'(a18)') 'End Diagnostic options' write(stdout,blank_fmt) write(stdout,delim_fmt) call POP_IOUnitsFlush(POP_stdout); call POP_IOUnitsFlush(stdout) ! temporary endif !EOC end subroutine init_diagnostics !*********************************************************************** !BOP ! !IROUTINE: diag_init_sums ! !INTERFACE: subroutine diag_init_sums ! !DESCRIPTION: ! Sets diagnostic flags and initializes diagnostic sums at beginning ! of each step. ! ! !REVISION HISTORY: ! same as module !EOP !BOC !----------------------------------------------------------------------- ! ! output formats ! !----------------------------------------------------------------------- character (8), parameter :: & time_fmt1 = '(a15,i8)' character (11), parameter :: & time_fmt2 = '(a15,f14.6)' character (20), parameter :: & date_fmt1 = '(a15,a2,1x,a3,1x,a4)' !----------------------------------------------------------------------- ! ! set logical switches for diagnostic timesteps ! do not compute on first pass of Matsuno ! !----------------------------------------------------------------------- ldiag_global = check_time_flag(diag_global_flag) ldiag_cfl = check_time_flag(diag_cfl_flag) ldiag_transport = check_time_flag(diag_transp_flag) if (mix_pass == 1) then ldiag_global = .false. ldiag_cfl = .false. ldiag_transport = .false. endif !----------------------------------------------------------------------- ! ! initialize arrays for global diagnostic sums. ! !----------------------------------------------------------------------- if (ldiag_global) then DIAG_KE_ADV_2D = c0 DIAG_KE_HMIX_2D = c0 DIAG_KE_VMIX_2D = c0 DIAG_KE_PRESS_2D = c0 DIAG_PE_2D = c0 DIAG_TRACER_ADV_2D = c0 DIAG_TRACER_HDIFF_2D = c0 DIAG_TRACER_VDIFF_2D = c0 DIAG_TRACER_SOURCE_2D = c0 DIAG_TRACER_SFC_FLX = c0 endif !----------------------------------------------------------------------- ! ! write some stuff to log file (stdout) ! !----------------------------------------------------------------------- if ((ldiag_global .or. ldiag_cfl) .and. & my_task == master_task) then write(stdout,blank_fmt) write(stdout,time_fmt1) 'Step number : ',nsteps_total write(stdout,date_fmt1) 'Date : ',cday,cmonth3,cyear write(stdout,time_fmt1) 'Hour : ',ihour write(stdout,time_fmt1) 'Minute : ',iminute write(stdout,time_fmt1) 'Second : ',isecond if (tsecond < seconds_in_hour) then write(stdout,time_fmt2) 'Time(seconds): ', tsecond elseif (tsecond < seconds_in_day) then write(stdout,time_fmt2) 'Time(hours) : ', & tsecond/seconds_in_hour elseif (tsecond < 31536000._r8) then write(stdout,time_fmt2) 'Time(days) : ', tday else write(stdout,time_fmt2) 'Time(years) : ', tyear endif endif !----------------------------------------------------------------------- !EOC end subroutine diag_init_sums !*********************************************************************** !BOP ! !IROUTINE: diag_global_preupdate ! !INTERFACE: subroutine diag_global_preupdate(DH,DHU) ! !DESCRIPTION: ! Computes diagnostic sums for those diagnostics computed from ! current time prognostic variables (before updating to new time) ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: real (r8), dimension(nx_block,ny_block,max_blocks_clinic), & intent(in) :: & DH, DHU ! change in surface height (-Fw) at T,U pts !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (int_kind) :: & i, k, & ! horiz, vertical level index n , & ! tracer index iblock , & ! block index dcount , & ! diag counter ib,ie,jb,je real (r8), dimension(max_blocks_clinic,6*nt+9) :: & local_sums ! array for holding block sums ! of each diagnostic real (r8) :: & sum_tmp ! temp for local sum real (r8), dimension(nx_block,ny_block) :: & WORK1, &! local work space TFACT, &! factor for normalizing sums UFACT ! factor for normalizing sums type (block) :: & this_block ! block information for current block !----------------------------------------------------------------------- ! ! compute these diagnostics only if it is time ! !----------------------------------------------------------------------- if (ldiag_global) then !----------------------------------------------------------------------- ! ! compute local sums of each diagnostic ! do all block sums first for better OpenMP performance ! !----------------------------------------------------------------------- local_sums = c0 !$OMP PARALLEL DO PRIVATE(iblock,this_block,ib,ie,jb,je,i, & !$OMP k,n,dcount,WORK1,UFACT,TFACT) do iblock = 1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) ib = this_block%ib ie = this_block%ie jb = this_block%jb je = this_block%je TFACT = TAREA(:,:,iblock)*RCALCT(:,:,iblock) UFACT = UAREA(:,:,iblock)*RCALCU(:,:,iblock) if (ltripole_grid .and. this_block%jblock == nblocks_y) then do i=1,nx_block if (this_block%i_glob(i) > nx_global/2) UFACT(i,je) = c0 end do endif dcount = 0 dcount = dcount + 1 WORK1 = (UVEL(:,:,1,curtime,iblock)*SMF(:,:,1,iblock) & +VVEL(:,:,1,curtime,iblock)*SMF(:,:,2,iblock))*UFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) dcount = dcount + 1 WORK1 = DIAG_KE_ADV_2D(:,:,iblock)*UFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) dcount = dcount + 1 WORK1 = DIAG_KE_HMIX_2D(:,:,iblock)*UFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) dcount = dcount + 1 WORK1 = DIAG_KE_VMIX_2D(:,:,iblock)*UFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) dcount = dcount + 1 WORK1 = DIAG_KE_PRESS_2D(:,:,iblock)*UFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) dcount = dcount + 1 WORK1 = DIAG_PE_2D(:,:,iblock)*TFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) dcount = dcount + 1 WORK1 = -(UBTROP(:,:,curtime,iblock)* & GRADPX(:,:,curtime,iblock) + & VBTROP(:,:,curtime,iblock)* & GRADPY(:,:,curtime,iblock))* & HU(:,:,iblock)*UFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) dcount = dcount + 1 WORK1 = p5*(UVEL(:,:,1,curtime,iblock)**2 + & VVEL(:,:,1,curtime,iblock)**2)* & DHU(:,:,iblock)*UFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) dcount = dcount + 1 WORK1 = PSURF(:,:,curtime,iblock)*DH(:,:,iblock)*TFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) do n=1,nt dcount = dcount + 1 WORK1 = DIAG_TRACER_ADV_2D(:,:,n,iblock)*TFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) dcount = dcount + 1 WORK1 = DIAG_TRACER_HDIFF_2D(:,:,n,iblock)*TFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) dcount = dcount + 1 WORK1 = DIAG_TRACER_VDIFF_2D(:,:,n,iblock)*TFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) dcount = dcount + 1 WORK1 = DIAG_TRACER_SOURCE_2D(:,:,n,iblock)*TFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) dcount = dcount + 1 do k=1,km if (sfc_layer_type == sfc_layer_varthick .and. & k == 1) then WORK1 = merge(TAREA(:,:,iblock)*( & (dz(1) + PSURF(:,:,newtime,iblock)/grav)* & TRACER(:,:,k,n,newtime,iblock) - & (dz(1) + PSURF(:,:,mixtime,iblock)/grav)* & TRACER(:,:,k,n,oldtime,iblock)), & c0, k <= KMT(:,:,iblock)) else if (partial_bottom_cells) then WORK1 = merge(DZT(:,:,k,iblock)*TAREA(:,:,iblock)*& (TRACER(:,:,k,n,newtime,iblock) - & TRACER(:,:,k,n,oldtime,iblock)), & c0, k <= KMT(:,:,iblock)) else WORK1 = merge(dz(k)*TAREA(:,:,iblock)* & (TRACER(:,:,k,n,newtime,iblock) - & TRACER(:,:,k,n,oldtime,iblock)), & c0, k <= KMT(:,:,iblock)) endif endif local_sums(iblock,dcount ) = & local_sums(iblock,dcount ) + & sum(WORK1(ib:ie,jb:je))/c2dtt(k) local_sums(iblock,dcount+1) = & local_sums(iblock,dcount+1) + & sum(abs(WORK1(ib:ie,jb:je)))/c2dtt(k) end do !k loop dcount = dcount + 1 end do !tracer loop end do ! block loop !$OMP END PARALLEL DO !----------------------------------------------------------------------- ! ! finish up scalar diagnostics ! !----------------------------------------------------------------------- dcount = 1 sum_tmp = sum(local_sums(:,dcount)) diag_ws = global_sum(sum_tmp,distrb_clinic)/volume_t dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) diag_ke_adv = global_sum(sum_tmp,distrb_clinic)/volume_t dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) diag_ke_hmix = global_sum(sum_tmp,distrb_clinic)/volume_t dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) diag_ke_vmix = global_sum(sum_tmp,distrb_clinic)/volume_t dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) diag_ke_press = global_sum(sum_tmp,distrb_clinic)/volume_t dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) diag_pe = global_sum(sum_tmp,distrb_clinic)/volume_t dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) diag_press_free = global_sum(sum_tmp,distrb_clinic)/volume_t dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) diag_ke_free = global_sum(sum_tmp,distrb_clinic)/volume_t dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) diag_ke_psfc = global_sum(sum_tmp,distrb_clinic)/volume_t do n=1,nt dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) diag_tracer_adv(n) = & global_sum(sum_tmp,distrb_clinic)/volume_t dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) diag_tracer_hdiff(n) = & global_sum(sum_tmp,distrb_clinic)/volume_t dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) diag_tracer_vdiff(n) = & global_sum(sum_tmp,distrb_clinic)/volume_t dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) diag_tracer_source(n) = & global_sum(sum_tmp,distrb_clinic)/volume_t dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) dtracer_avg(n) = global_sum(sum_tmp,distrb_clinic)/volume_t dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) dtracer_abs(n) = global_sum(sum_tmp,distrb_clinic)/volume_t end do !----------------------------------------------------------------------- endif ! ldiag_global !----------------------------------------------------------------------- if (allocated(HMXL)) then !$OMP PARALLEL DO PRIVATE(iblock) do iblock=1,nblocks_clinic call accumulate_tavg_field(HMXL(:,:,iblock), tavg_HMXL, iblock, 1) end do !$OMP END PARALLEL DO !$OMP PARALLEL DO PRIVATE(iblock) do iblock=1,nblocks_clinic call accumulate_tavg_field(HMXL(:,:,iblock), tavg_HMXL_2, iblock, 1) end do !$OMP END PARALLEL DO ! QL, 150526, mixed layer depth with density criterion !$OMP PARALLEL DO PRIVATE(iblock) do iblock=1,nblocks_clinic call accumulate_tavg_field(HMXL_DR(:,:,iblock), tavg_HMXL_DR, iblock, 1) end do !$OMP END PARALLEL DO !$OMP PARALLEL DO PRIVATE(iblock) do iblock=1,nblocks_clinic call accumulate_tavg_field(HMXL(:,:,iblock), tavg_XMXL, iblock, 1) end do !$OMP END PARALLEL DO !$OMP PARALLEL DO PRIVATE(iblock) do iblock=1,nblocks_clinic call accumulate_tavg_field(HMXL(:,:,iblock), tavg_XMXL_2, iblock, 1) end do !$OMP END PARALLEL DO !$OMP PARALLEL DO PRIVATE(iblock) do iblock=1,nblocks_clinic call accumulate_tavg_field(HMXL(:,:,iblock), tavg_TMXL, iblock, 1) end do !$OMP END PARALLEL DO end if if (allocated(KPP_HBLT)) then !$OMP PARALLEL DO PRIVATE(iblock) do iblock=1,nblocks_clinic call accumulate_tavg_field(KPP_HBLT(:,:,iblock), tavg_HBLT, iblock, 1) end do !$OMP END PARALLEL DO !$OMP PARALLEL DO PRIVATE(iblock) do iblock=1,nblocks_clinic call accumulate_tavg_field(KPP_HBLT(:,:,iblock), tavg_XBLT, iblock, 1) end do !$OMP END PARALLEL DO !$OMP PARALLEL DO PRIVATE(iblock) do iblock=1,nblocks_clinic call accumulate_tavg_field(KPP_HBLT(:,:,iblock), tavg_TBLT, iblock, 1) end do !$OMP END PARALLEL DO end if !----------------------------------------------------------------------- ! ! accumulate movie diagnostics if requested ! !----------------------------------------------------------------------- if (movie_requested(movie_HMXL) ) then !$OMP PARALLEL DO do iblock=1,nblocks_clinic call update_movie_field(HMXL(:,:,iblock), movie_HMXL, iblock, 1) end do !$OMP END PARALLEL DO endif if (movie_requested(movie_HBLT) ) then !$OMP PARALLEL DO do iblock=1,nblocks_clinic call update_movie_field(KPP_HBLT(:,:,iblock), movie_HBLT, iblock, 1) end do !$OMP END PARALLEL DO endif !EOC end subroutine diag_global_preupdate !*********************************************************************** !BOP ! !IROUTINE: diag_global_afterupdate ! !INTERFACE: subroutine diag_global_afterupdate ! !DESCRIPTION: ! Computes diagnostic sums for those diagnostics computed from ! updated prognostic variables (after updating to new time). ! ! !REVISION HISTORY: ! same as module !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (POP_i4) :: & errorCode ! returned error code integer (int_kind) :: & i, k, & ! horiz, vertical level index n , & ! tracer index iblock , & ! block index dcount , & ! diag counter ib,ie,jb,je real (r8), dimension(max_blocks_clinic,2*nt+2) :: & local_sums ! array for holding block sums ! of each diagnostic real (r8), dimension(max_blocks_clinic,km,nt) :: & local_sums_k ! array for holding block sums ! of each diagnostic at each level real (r8) :: & sum_tmp ! temp for local sum real (r8), dimension(nx_block,ny_block) :: & WORK1, &! local work space TFACT, &! factor for normalizing sums UFACT ! factor for normalizing sums type (block) :: & this_block ! block information for current block !----------------------------------------------------------------------- ! ! compute these diagnostics only if it is time ! !----------------------------------------------------------------------- if (ldiag_global) then !----------------------------------------------------------------------- ! ! compute local block sums of each diagnostic ! !----------------------------------------------------------------------- local_sums = c0 if (diag_all_levels) local_sums_k = c0 !$OMP PARALLEL DO & !$OMP PRIVATE(iblock,this_block,ib,ie,jb,je,i,k,n,dcount, & !$OMP WORK1,TFACT,UFACT) do iblock = 1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) ib = this_block%ib ie = this_block%ie jb = this_block%jb je = this_block%je TFACT = TAREA(:,:,iblock)*RCALCT(:,:,iblock) UFACT = UAREA(:,:,iblock)*RCALCU(:,:,iblock) if (ltripole_grid .and. this_block%jblock == nblocks_y) then do i=1,nx_block if (this_block%i_glob(i) > nx_global/2) UFACT(i,je) = c0 end do endif dcount = 0 !*** global mean sea level dcount = dcount + 1 WORK1 = PSURF(:,:,curtime,iblock)/grav*TFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) !*** avg horiz KE at new time dcount = dcount + 1 do k = 1,km if (partial_bottom_cells) then WORK1 = merge(p5*(UVEL(:,:,k,curtime,iblock)**2 + & VVEL(:,:,k,curtime,iblock)**2)* & UFACT*DZU(:,:,k,iblock), & c0, KMU(:,:,iblock) >= k) else WORK1 = merge(p5*(UVEL(:,:,k,curtime,iblock)**2 + & VVEL(:,:,k,curtime,iblock)**2)* & UFACT*dz(k), & c0, KMU(:,:,iblock) >= k) endif local_sums(iblock,dcount) = local_sums(iblock,dcount) + & sum(WORK1(ib:ie,jb:je)) enddo ! k loop !*** average tracers at new time do n=1,nt dcount = dcount + 1 do k = 1,km if (sfc_layer_type == sfc_layer_varthick .and. & k == 1) then WORK1 = merge(TRACER(:,:,k,n,curtime,iblock)* & TAREA(:,:,iblock)* & (c1 + PSURF(:,:,curtime,iblock)/ & (dz(1)*grav))*dz(1) & ,c0,KMT(:,:,iblock) >= k) else if (partial_bottom_cells) then WORK1 = merge(TRACER(:,:,k,n,curtime,iblock)* & TAREA(:,:,iblock)*DZT(:,:,k,iblock),& c0,KMT(:,:,iblock) >= k) else WORK1 = merge(TRACER(:,:,k,n,curtime,iblock)* & TAREA(:,:,iblock)*dz(k), & c0,KMT(:,:,iblock) >= k) endif endif local_sums(iblock,dcount) = local_sums(iblock,dcount) + & sum(WORK1(ib:ie,jb:je)) if (diag_all_levels) then local_sums_k(iblock,k,n) = sum(WORK1(ib:ie,jb:je)) endif end do ! k loop end do ! tracer loop !*** surface tracer fluxes do n=1,nt dcount = dcount + 1 WORK1 = DIAG_TRACER_SFC_FLX(:,:,n,iblock)*TFACT local_sums(iblock,dcount) = sum(WORK1(ib:ie,jb:je)) end do ! tracer loop end do ! block loop !$OMP END PARALLEL DO !----------------------------------------------------------------------- ! ! finish up by computing global sums of local sums ! !----------------------------------------------------------------------- dcount= 1 sum_tmp = sum(local_sums(:,dcount)) diag_sealevel = global_sum(sum_tmp,distrb_clinic)/area_t dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) diag_ke = global_sum(sum_tmp,distrb_clinic)/volume_u do n=1,nt dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) avg_tracer(n) = global_sum(sum_tmp,distrb_clinic)/volume_t if (diag_all_levels) then do k=1,km sum_tmp = sum(local_sums_k(:,k,n)) avg_tracer_k(k,n) = global_sum(sum_tmp,distrb_clinic)/ & area_t_k(k) end do endif end do do n=1,nt dcount = dcount + 1 sum_tmp = sum(local_sums(:,dcount)) sfc_tracer_flux(n) = global_sum(sum_tmp,distrb_clinic)/area_t end do !----------------------------------------------------------------------- ! ! convert some tracer diagnostics to different units ! !----------------------------------------------------------------------- avg_tracer(2) = avg_tracer(2)*salt_to_ppt ! salinity to ppt if (diag_all_levels) then do k=1,km avg_tracer_k(k,2) = avg_tracer_k(k,2)*salt_to_ppt end do endif sfc_tracer_flux(1) = sfc_tracer_flux(1)/hflux_factor ! flux in W/m2 sfc_tracer_flux(2) = sfc_tracer_flux(2)/salinity_factor ! FW flux !----------------------------------------------------------------------- ! ! retrieve solver diagnostics ! !----------------------------------------------------------------------- call POP_SolversGetDiagnostics(iterationCount, rmsResidual, & errorCode) !----------------------------------------------------------------------- endif ! ldiag_global !----------------------------------------------------------------------- !EOC end subroutine diag_global_afterupdate !*********************************************************************** !BOP ! !IROUTINE: diag_print ! !INTERFACE: subroutine diag_print ! !DESCRIPTION: ! Writes scalar diagnostic info to both stdout and diagnostic output ! file. ! ! !REVISION HISTORY: ! same as module !EOP !BOC !----------------------------------------------------------------------- ! ! local variable ! !----------------------------------------------------------------------- integer (int_kind) :: & k,n, &! dummy indices ier character (25) :: diag_name ! name of output diagnostic character (29), parameter :: & diag_fmt = '(1pe22.15,1x,1pe22.15,2x,a25)' character (14), parameter :: & stdo_fmt = '(a32,1pe22.15)' character (21), parameter :: & trcr_fmt = '(a24,i3,a19,1pe22.15)' character (char_len_long) :: & string !----------------------------------------------------------------------- ! ! print diagnostics if it is time ! !----------------------------------------------------------------------- if (ldiag_global) then if (my_task == master_task) then !*** append diagnostic output to end of diag output file open(diag_unit, file=diag_outfile, status='old', & position='append') !----------------------------------------------------------------------- ! ! print global diagnostics to stdout ! !----------------------------------------------------------------------- write(stdout,blank_fmt) write(stdout,'(a17)') 'K.E. diagnostics:' write(stdout,stdo_fmt) & ' mean K.E. :', diag_ke write(stdout,stdo_fmt) & ' K.E. change from horiz. visc. :', diag_ke_hmix write(stdout,stdo_fmt) & ' K.E. change from vert. visc. :', diag_ke_vmix write(stdout,stdo_fmt) & ' K.E. change from wind stress :', diag_ws write(stdout,blank_fmt) write(stdout,stdo_fmt) & ' advective K.E. balance :', diag_ke_adv write(stdout,stdo_fmt) & ' :', diag_ke_free write(stdout,stdo_fmt) & ' P.E. verses pressure work :', diag_pe write(stdout,stdo_fmt) & ' :', diag_ke_press write(stdout,stdo_fmt) & ' surface-pressure contribution :', diag_ke_psfc write(stdout,stdo_fmt) & ' :', diag_press_free write(stdout,blank_fmt) write(stdout,*) 'Tracer diagnostics:' do n=1,nt write(stdout,blank_fmt) write(stdout,trcr_fmt) & ' mean value of tracer ',n, & ' : ', avg_tracer(n) write(stdout,trcr_fmt) & ' abs change in tracer ',n, & ' : ', dtracer_abs(n) write(stdout,trcr_fmt) & ' mean change in tracer ',n, & ' : ', dtracer_avg(n) write(stdout,trcr_fmt) & ' mean change in tracer ',n, & ' from advection : ', diag_tracer_adv(n) write(stdout,trcr_fmt) & ' mean change in tracer ',n, & ' from horiz. diff.: ', diag_tracer_hdiff(n) write(stdout,trcr_fmt) & ' mean change in tracer ',n, & ' from vert. diff.: ', diag_tracer_vdiff(n) write(stdout,trcr_fmt) & ' mean change in tracer ',n, & ' from sources : ', diag_tracer_source(n) write(stdout,trcr_fmt) & ' surface flux of tracer ',n, & ' : ', sfc_tracer_flux(n) if (diag_all_levels) then do k=1,km write(stdout,'(a22,i3,a10,i3,a3,1pe22.15)') & ' mean value of tracer ', n, & ' at level ', k, & ' : ', avg_tracer_k(k,n) end do endif end do write(stdout,blank_fmt) write(stdout,'(a18)') 'Other diagnostics:' write(stdout,'(a24,1pe22.15)') ' global mean sea level: ', & diag_sealevel write(stdout,'(a19,1pe22.15,2x,i4)') ' residual, scans : ', & rmsResidual, iterationCount !----------------------------------------------------------------------- ! ! print global diagnostics to output file ! !----------------------------------------------------------------------- diag_name = 'mean KE' write(diag_unit,diag_fmt) tday, diag_ke, diag_name do n=1,nt write(diag_name,'(a12,i2)') 'mean tracer ',n write(diag_unit,diag_fmt) tday, avg_tracer(n), diag_name if (diag_all_levels) then do k=1,km write(diag_name,'(a12,i2,a8,i2)') & 'mean tracer ',n,' at lvl ',k write(diag_unit,diag_fmt) tday, avg_tracer_k(k,n), & diag_name end do endif end do diag_name = 'mean sea level' write(diag_unit,diag_fmt) tday, diag_sealevel, diag_name diag_name = 'solver residual' write(diag_unit,diag_fmt) tday, rmsResidual, diag_name diag_name = 'solver iterations' write(diag_unit,diag_fmt) tday, real(iterationCount), diag_name do n=1,nt write(diag_name,'(a7,i3,a11)') 'Tracer ',n,' change avg' write(diag_unit,diag_fmt) tday, dtracer_avg(n), diag_name write(diag_name,'(a7,i3,a11)') 'Tracer ',n,' change abs' write(diag_unit,diag_fmt) tday, dtracer_abs(n), diag_name write(diag_name,'(a7,i3,a14)') 'Tracer ',n,' change advect' write(diag_unit,diag_fmt) tday, diag_tracer_adv(n), diag_name write(diag_name,'(a7,i3,a13)') 'Tracer ',n,' change hdiff' write(diag_unit,diag_fmt) tday, diag_tracer_hdiff(n), diag_name write(diag_name,'(a7,i3,a13)') 'Tracer ',n,' change vdiff' write(diag_unit,diag_fmt) tday, diag_tracer_vdiff(n), diag_name write(diag_name,'(a7,i3,a14)') 'Tracer ',n,' change source' write(diag_unit,diag_fmt) tday, diag_tracer_source(n), diag_name write(diag_name,'(a7,i3,a14)') 'Tracer ',n,' surface flux' write(diag_unit,diag_fmt) tday, sfc_tracer_flux(n), diag_name end do diag_name = 'KE change horiz visc' write(diag_unit,diag_fmt) tday, diag_ke_hmix, diag_name diag_name = 'KE change vert visc' write(diag_unit,diag_fmt) tday, diag_ke_vmix, diag_name diag_name = 'KE change wind stress' write(diag_unit,diag_fmt) tday, diag_ws, diag_name diag_name = 'advect KE balance' write(diag_unit,diag_fmt) tday, diag_ke_adv, diag_name diag_name = 'advect KE bal free sfc' write(diag_unit,diag_fmt) tday, diag_ke_free, diag_name diag_name = 'PE change' write(diag_unit,diag_fmt) tday, diag_pe, diag_name diag_name = 'pressure work' write(diag_unit,diag_fmt) tday, diag_ke_press, diag_name diag_name = 'KE change sfc-pressure' write(diag_unit,diag_fmt) tday, diag_ke_psfc, diag_name diag_name = 'press work free sfc' write(diag_unit,diag_fmt) tday, diag_press_free, diag_name !----------------------------------------------------------------------- ! ! close file and print timer info ! !----------------------------------------------------------------------- close(diag_unit) endif ! master_task call timer_print_all() call document ('diag_print', 'file written: '// trim(diag_outfile)) endif ! ldiag_global !----------------------------------------------------------------------- !EOC end subroutine diag_print !*********************************************************************** !BOP ! !IROUTINE: diag_transport ! !INTERFACE: subroutine diag_transport ! !DESCRIPTION: ! Calculates and prints various user-defined transports. ! ! !REVISION HISTORY: ! same as module !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (int_kind) :: & i,j,k,n,iblock, &! dummy loop indices ib,ie,jb,je, &! beg,end indices for block physical domain ier real (r8), dimension(nx_block,ny_block) :: & MASS_Z, &! zonal transport of mass HEAT_Z, &! zonal transport of heat SALT_Z, &! zonal transport of salt MASS_M, &! merid transport of mass HEAT_M, &! merid transport of heat SALT_M, &! merid transport of salt WORK1, WORK2 ! local temp space real (r8) :: & sum_tmp ! temp for computing block sums real (r8), dimension(:,:), allocatable :: & mass_tran, &! local sum of mass transport heat_tran, &! local sum of heat transport salt_tran ! local sum of salt transport character (char_len_long) :: & string type (block) :: & this_block ! block information for current block !----------------------------------------------------------------------- ! ! only compute if it is time ! initialize transport sums ! !----------------------------------------------------------------------- if (ldiag_transport) then if (my_task == master_task) then !*** re-open file and append current values to previous !*** outputs from this run open(trans_unit, file=diag_transport_outfile, & status='old', position='append') endif !----------------------------------------------------------------------- ! ! calculate transports for each block ! !----------------------------------------------------------------------- allocate(mass_tran(nblocks_clinic,num_transports), & heat_tran(nblocks_clinic,num_transports), & salt_tran(nblocks_clinic,num_transports)) !$OMP PARALLEL DO PRIVATE(iblock,this_block,ib,ie,jb,je,k,j,i,n, & !$OMP WORK1,WORK2,MASS_M,SALT_M,HEAT_M, & !$OMP MASS_Z,SALT_Z,HEAT_Z) do iblock = 1,nblocks_clinic this_block = get_block(blocks_clinic(iblock),iblock) ib = this_block%ib ie = this_block%ie jb = this_block%jb je = this_block%je mass_tran(iblock,:) = c0 heat_tran(iblock,:) = c0 salt_tran(iblock,:) = c0 do k = 1,km if (partial_bottom_cells) then do j=jb,je do i=ib,ie MASS_M(i,j) = p5*( & UVEL(i,j ,k,curtime,iblock)* & DYU(i,j ,iblock)* & DZU(i,j ,k,iblock) + & UVEL(i,j-1,k,curtime,iblock)* & DYU(i,j-1,iblock)* & DZU(i,j-1,k,iblock) ) MASS_Z(i,j) = p5*( & VVEL(i ,j,k,curtime,iblock)* & DXU(i ,j,iblock)* & DZU(i ,j,k,iblock) + & VVEL(i-1,j,k,curtime,iblock)* & DXU(i-1,j,iblock)* & DZU(i-1,j,k,iblock) ) end do end do else do j=jb,je do i=ib,ie MASS_M(i,j) = p5*( & UVEL(i,j ,k,curtime,iblock)* & DYU(i,j ,iblock)*dz(k) & + UVEL(i,j-1,k,curtime,iblock)* & DYU(i,j-1,iblock)*dz(k) ) MASS_Z(i,j) = p5*( & VVEL(i ,j,k,curtime,iblock)* & DXU(i ,j,iblock)*dz(k) & + VVEL(i-1,j,k,curtime,iblock)* & DXU(i-1,j,iblock)*dz(k) ) end do end do endif do j=jb,je do i=ib,ie HEAT_M(i,j) = p5*MASS_M(i,j)* & (TRACER(i+1,j,k,1,curtime,iblock) + & TRACER(i ,j,k,1,curtime,iblock)) SALT_M(i,j) = p5*MASS_M(i,j)* & (TRACER(i+1,j,k,2,curtime,iblock) + & TRACER(i ,j,k,2,curtime,iblock)) HEAT_Z(i,j) = p5*MASS_Z(i,j)* & (TRACER(i,j+1,k,1,curtime,iblock) + & TRACER(i,j ,k,1,curtime,iblock)) SALT_Z(i,j) = p5*MASS_Z(i,j)* & (TRACER(i,j+1,k,2,curtime,iblock) + & TRACER(i,j ,k,2,curtime,iblock)) end do end do !----------------------------------------------------------------------- ! ! compute transport sums for each defined transport ! !----------------------------------------------------------------------- do n=1,num_transports if (k >= transports(n)%add(5,iblock) .and. & k <= transports(n)%add(6,iblock) ) then !*** local sum select case (transports(n)%type) case (transport_type_zonal) do j=transports(n)%add(3,iblock), & transports(n)%add(4,iblock) do i=transports(n)%add(1,iblock), & transports(n)%add(2,iblock) mass_tran(iblock,n) = mass_tran(iblock,n) + & MASS_Z(i,j) heat_tran(iblock,n) = heat_tran(iblock,n) + & HEAT_Z(i,j) salt_tran(iblock,n) = salt_tran(iblock,n) + & SALT_Z(i,j) end do end do case (transport_type_merid) do j=transports(n)%add(3,iblock), & transports(n)%add(4,iblock) do i=transports(n)%add(1,iblock), & transports(n)%add(2,iblock) mass_tran(iblock,n) = mass_tran(iblock,n) + & MASS_M(i,j) heat_tran(iblock,n) = heat_tran(iblock,n) + & HEAT_M(i,j) salt_tran(iblock,n) = salt_tran(iblock,n) + & SALT_M(i,j) end do end do end select endif end do ! transport loop end do ! k loop end do ! block loop !$OMP END PARALLEL DO !----------------------------------------------------------------------- ! ! finish global sums for each transport, convert to useful units ! and write to output file ! !----------------------------------------------------------------------- do n=1,num_transports sum_tmp = sum(mass_tran(:,n)) transports(n)%mass = global_sum(sum_tmp,distrb_clinic)* & mass_to_Sv sum_tmp = sum(heat_tran(:,n)) transports(n)%heat = global_sum(sum_tmp,distrb_clinic)* & heat_to_PW sum_tmp = sum(salt_tran(:,n)) transports(n)%salt = global_sum(sum_tmp,distrb_clinic)* & salt_to_Svppt if (my_task == master_task) then write(trans_unit,'(4(1pe13.6,2x), a)') & tday, transports(n)%mass, & transports(n)%heat, & transports(n)%salt, & trim(transports(n)%name) endif end do ! transport loop if (my_task == master_task) then close(trans_unit) endif ! master_task call document ('diag_transport', 'file written: '//trim(diag_transport_outfile)) deallocate(mass_tran, heat_tran, salt_tran) endif ! ldiag_transport !----------------------------------------------------------------------- !EOC end subroutine diag_transport !*********************************************************************** !BOP ! !IROUTINE: cfl_advect ! !INTERFACE: subroutine cfl_advect(k,iblock,UTE,UTW,VTN,VTS,WTK,WTKB,this_block) ! !DESCRIPTION: ! Calculates maximum advective cfl number for a particular block. ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: integer (int_kind), intent(in) :: & k, &! vertical level index iblock ! local block index in baroclinic distribution real (r8), dimension(nx_block,ny_block), intent(in) :: & UTE, UTW, &! advective fluxes through east, west faces VTN, VTS, &! advective fluxes through north,south faces WTK, WTKB ! advective fluxes through top, bottom faces type (block), intent(in) :: & this_block ! block information for current block !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (int_kind) :: & i,j ! dummy loop indices real (r8), dimension(nx_block,ny_block) :: & CFL_E, CFL_W, &! cfl numbers due to advective fluxes through CFL_N, CFL_S, &! each face CFL_T, CFL_B, &! CFL_TOT ! sum for total cfl number !----------------------------------------------------------------------- ! ! only compute cfl numbers if it is time ! !----------------------------------------------------------------------- if (ldiag_cfl) then !----------------------------------------------------------------------- ! ! advective fluxes through east,west,north,south,top,bottom faces ! !----------------------------------------------------------------------- if (partial_bottom_cells) then CFL_E = abs(UTE*TAREA_R(:,:,iblock)/DZU(:,:,k,iblock))*dt(k) CFL_W = abs(UTW*TAREA_R(:,:,iblock)/DZU(:,:,k,iblock))*dt(k) CFL_N = abs(VTN*TAREA_R(:,:,iblock)/DZU(:,:,k,iblock))*dt(k) CFL_S = abs(VTS*TAREA_R(:,:,iblock)/DZU(:,:,k,iblock))*dt(k) CFL_T = abs(WTK /DZT(:,:,k,iblock))*dt(k) CFL_B = abs(WTKB/DZT(:,:,k,iblock))*dt(k) else CFL_E = abs(UTE*TAREA_R(:,:,iblock))*dt(k) CFL_W = abs(UTW*TAREA_R(:,:,iblock))*dt(k) CFL_N = abs(VTN*TAREA_R(:,:,iblock))*dt(k) CFL_S = abs(VTS*TAREA_R(:,:,iblock))*dt(k) CFL_T = abs(WTK /dz(k))*dt(k) CFL_B = abs(WTKB/dz(k))*dt(k) endif CFL_TOT = p5*(CFL_E + CFL_W + CFL_N + CFL_S + CFL_T + CFL_B) !----------------------------------------------------------------------- ! ! find max cfl numbers and locations for this block ! !----------------------------------------------------------------------- cfl_advuk_block(iblock,k) = -c1 cfl_advvk_block(iblock,k) = -c1 cfl_advwk_block(iblock,k) = -c1 cfl_advtk_block(iblock,k) = -c1 cfladd_advuk_block(:,iblock,k) = 0 cfladd_advvk_block(:,iblock,k) = 0 cfladd_advwk_block(:,iblock,k) = 0 cfladd_advtk_block(:,iblock,k) = 0 do j=this_block%jb,this_block%je do i=this_block%ib,this_block%ie !*** zonal cfl number east flux if (CFL_E(i,j) > cfl_advuk_block(iblock,k)) then cfl_advuk_block(iblock,k) = CFL_E(i,j) cfladd_advuk_block(1,iblock,k) = this_block%i_glob(i) cfladd_advuk_block(2,iblock,k) = this_block%j_glob(j) endif !*** zonal cfl number west flux if (CFL_W(i,j) > cfl_advuk_block(iblock,k)) then cfl_advuk_block(iblock,k) = CFL_W(i,j) cfladd_advuk_block(1,iblock,k) = this_block%i_glob(i) cfladd_advuk_block(2,iblock,k) = this_block%j_glob(j) endif !*** meridional cfl number north flux if (CFL_N(i,j) > cfl_advvk_block(iblock,k)) then cfl_advvk_block(iblock,k) = CFL_N(i,j) cfladd_advvk_block(1,iblock,k) = this_block%i_glob(i) cfladd_advvk_block(2,iblock,k) = this_block%j_glob(j) endif !*** meridional cfl number south flux if (CFL_S(i,j) > cfl_advvk_block(iblock,k)) then cfl_advvk_block(iblock,k) = CFL_S(i,j) cfladd_advvk_block(1,iblock,k) = this_block%i_glob(i) cfladd_advvk_block(2,iblock,k) = this_block%j_glob(j) endif !*** vertical cfl number top flux if (CFL_T(i,j) > cfl_advwk_block(iblock,k)) then cfl_advwk_block(iblock,k) = CFL_T(i,j) cfladd_advwk_block(1,iblock,k) = this_block%i_glob(i) cfladd_advwk_block(2,iblock,k) = this_block%j_glob(j) endif !*** vertical cfl number bottom flux if (CFL_B(i,j) > cfl_advwk_block(iblock,k)) then cfl_advwk_block(iblock,k) = CFL_B(i,j) cfladd_advwk_block(1,iblock,k) = this_block%i_glob(i) cfladd_advwk_block(2,iblock,k) = this_block%j_glob(j) endif !*** total advective CFL number if (CFL_TOT(i,j) > cfl_advtk_block(iblock,k)) then cfl_advtk_block(iblock,k) = CFL_TOT(i,j) cfladd_advtk_block(1,iblock,k) = this_block%i_glob(i) cfladd_advtk_block(2,iblock,k) = this_block%j_glob(j) endif end do end do !----------------------------------------------------------------------- ! ! finished computing advective cfl numbers ! !----------------------------------------------------------------------- endif ! ldiag_cfl !----------------------------------------------------------------------- !EOC end subroutine cfl_advect !*********************************************************************** !BOP ! !IROUTINE: cfl_vdiff ! !INTERFACE: subroutine cfl_vdiff(k,VDC,VVC,this_block) ! !DESCRIPTION: ! Calculates maximum vertical diffusion cfl numbers and their ! locations in i,j,k directions. ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: integer (int_kind), intent(in) :: & k ! vertical level index real (r8), dimension(:,:,:,:), intent(in) :: & VDC ! tracer diffusivity real (r8), dimension(:,:,:), intent(in) :: & VVC ! viscosity type (block), intent(in) :: & this_block ! block information for current block !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (int_kind) :: & kmax, klvl, km1, &! vertical level indices i,j,n, &! dummy index bid ! local block id real (r8) :: & local_max ! temp for holding local block max real (r8), dimension(nx_block,ny_block,max_blocks_clinic), save :: & VDCM1, VVCM1 ! save k-1 values for next k level real (r8), dimension(nx_block,ny_block) :: & WORK1, &! local work space DZ1, DZ2, DZ3 ! grid arrays !----------------------------------------------------------------------- ! ! only compute cfl numbers if it is time ! !----------------------------------------------------------------------- if (ldiag_cfl) then !----------------------------------------------------------------------- ! ! determine size of VDC arrays and set up some grid arrays ! !----------------------------------------------------------------------- bid = this_block%local_id kmax = size(VDC, DIM=3) if (kmax > km) then ! KPP uses VDC(0:km+1) klvl = k+1 km1 = k else klvl = k km1 = k-1 endif if (partial_bottom_cells) then if (k < km) then DZ1 = c2/(DZT(:,:,k,bid) + DZT(:,:,k+1,bid)) else DZ1 = c2/DZT(:,:,k,bid) endif if (k == 1) then DZ2 = c2/DZT(:,:,k,bid) else DZ2 = c2/(DZT(:,:,k,bid) + DZT(:,:,k-1,bid)) endif DZ3 = c1/DZT(:,:,k,bid) else DZ1 = dzwr(k ) DZ2 = dzwr(k-1) DZ3 = dzr (k ) endif !----------------------------------------------------------------------- ! ! if VDC a 2d array, compute cfl number for this level and save ! for next level ! !----------------------------------------------------------------------- if (kmax == 1) then if (k == 1) VDCM1(:,:,bid) = c0 WORK1 = abs(c2*(VDC(:,:,1,1)*DZ1 + & VDCM1(:,:,bid)*DZ2)*DZ3) !*** swap to prepare for next k level VDCM1(:,:,bid) = VDC(:,:,1,1) !----------------------------------------------------------------------- ! ! if VDC a 3d array, use k-1 value directly from array ! !----------------------------------------------------------------------- else if (k == 1) then !*** k=1 case (k-1 values = 0) WORK1 = c0 do n=1,size(VDC,DIM=4) WORK1 = max(WORK1, & abs(c2*VDC(:,:,klvl,n)*dzwr(1)*dzr(1))) end do else !*** all other k levels WORK1 = c0 do n=1,size(VDC,DIM=4) WORK1 = max(WORK1,abs(c2*( & VDC(:,:,klvl,n)*DZ1 + & VDC(:,:,km1 ,n)*DZ2)*DZ3)) end do endif ! k=1 endif ! kmax = 1 !----------------------------------------------------------------------- ! ! find max tracer cfl for this block ! !----------------------------------------------------------------------- cfl_vdifftk_block(bid,k) = c0 local_max = -c1 do j=this_block%jb,this_block%je do i=this_block%ib,this_block%ie if (WORK1(i,j) > local_max) then local_max = WORK1(i,j) cfladd_vdifftk_block(1,bid,k) = this_block%i_glob(i) cfladd_vdifftk_block(2,bid,k) = this_block%j_glob(j) endif end do end do cfl_vdifftk_block(bid,k) = local_max*dt(k) !----------------------------------------------------------------------- ! ! repeat for momentum diffusion ! !----------------------------------------------------------------------- kmax = size(VVC, DIM=3) if (partial_bottom_cells) then if (k < km) then DZ1 = c2/(DZU(:,:,k,bid) + DZU(:,:,k+1,bid)) else DZ1 = c2/DZU(:,:,k,bid) endif if (k == 1) then DZ2 = c2/DZU(:,:,k,bid) else DZ2 = c2/(DZU(:,:,k,bid) + DZU(:,:,k-1,bid)) endif DZ3 = c1/DZU(:,:,k,bid) else DZ1 = dzwr(k ) DZ2 = dzwr(k-1) DZ3 = dzr (k ) endif !----------------------------------------------------------------------- ! ! if VVC a 2d array, compute cfl number for this level and save ! for next level ! !----------------------------------------------------------------------- if (kmax == 1) then if (k == 1) then VVCM1(:,:,bid) = c0 endif WORK1 = abs(c2*(VVC(:,:,1) *DZ1 + & VVCM1(:,:,bid)*DZ2)*DZ3) VVCM1(:,:,bid) = VVC(:,:,1) ! swap to prepare for next k level !----------------------------------------------------------------------- ! ! if VVC a 3d array, use k-1 value directly from array ! !----------------------------------------------------------------------- else if (k == 1) then !*** k=1 case (k-1 values = 0) WORK1 = abs(c2*VVC(:,:,1)*dzwr(1)*dzr(1)) else !*** all other k levels WORK1 = abs(c2*(VVC(:,:,k )*DZ1 + & VVC(:,:,k-1)*DZ2)*DZ3) endif ! k=1 endif ! kmax=1 !----------------------------------------------------------------------- ! ! find max momentum cfl for this block ! !----------------------------------------------------------------------- cfl_vdiffuk_block(bid,k) = c0 local_max = -c1 do j=this_block%jb,this_block%je do i=this_block%ib,this_block%ie if (WORK1(i,j) > local_max) then local_max = WORK1(i,j) cfladd_vdiffuk_block(1,bid,k) = this_block%i_glob(i) cfladd_vdiffuk_block(2,bid,k) = this_block%j_glob(j) endif end do end do cfl_vdiffuk_block(bid,k) = local_max*dtu !----------------------------------------------------------------------- ! ! finished computing vertical diffusion cfl numbers ! !----------------------------------------------------------------------- endif ! ldiag_cfl !----------------------------------------------------------------------- !EOC end subroutine cfl_vdiff !*********************************************************************** !BOP ! !IROUTINE: cfl_hdiff ! !INTERFACE: subroutine cfl_hdiff(k,iblock,HDIFFCFL,t_or_u,this_block) ! !DESCRIPTION: ! Calculates maximum horizontal diffusion cfl numbers and their ! locations in i,j,k directions ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: integer (int_kind), intent(in) :: & k, &! vertical level index iblock, &! local block index in baroclinic distribution t_or_u ! called from hdifft(1) or hdiffu(2) real (r8), dimension(nx_block,ny_block), intent(in) :: & HDIFFCFL ! hdiff cfl number at grid points in block type (block), intent(in) :: & this_block ! block information for current block !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (int_kind) :: & i,j ! dummy loop indices real (r8) :: cfltemp ! temporary for holding local maximum integer (int_kind), dimension(2) :: & cfltemp_add ! temporary for holding address of local max !----------------------------------------------------------------------- ! ! compute max cfl number for this block and level ! !----------------------------------------------------------------------- cfltemp = c0 cfltemp_add(:) = 0 do j=this_block%jb,this_block%je do i=this_block%ib,this_block%ie if (HDIFFCFL(i,j) > cfltemp) then cfltemp = HDIFFCFL(i,j) cfltemp_add(1) = this_block%i_glob(i) cfltemp_add(2) = this_block%j_glob(j) endif end do end do select case (t_or_u) case(1) cfl_hdifftk_block(iblock,k) = cfltemp*dt(k) cfladd_hdifftk_block(:,iblock,k) = cfltemp_add(:) case(2) cfl_hdiffuk_block(iblock,k) = cfltemp*dtu cfladd_hdiffuk_block(:,iblock,k) = cfltemp_add(:) end select !----------------------------------------------------------------------- !EOC end subroutine cfl_hdiff !*********************************************************************** !BOP ! !IROUTINE: diag_velocity ! !INTERFACE: subroutine diag_velocity ! !DESCRIPTION: ! Writes velocity diagnostics info to velocity output file. ! ! !REVISION HISTORY: ! same as module !EOP !BOC !----------------------------------------------------------------------- ! ! local variable ! !----------------------------------------------------------------------- character (15), parameter :: & diag_fmt = '(4(1pe13.6,2x))' character (char_len_long) :: & string integer (int_kind) :: & ier, n real (r8), dimension(nx_global,ny_global) :: & ARRAY_G if ( lccsm ) then call gather_global (ARRAY_G, VVEL(:,:,1,curtime,:), master_task, & distrb_clinic) if ( my_task == master_task ) then !*** append velocity output to end of velocity output file open (velocity_unit, file=diag_velocity_outfile, status='old', & position='append') do n=1,num_vel_loc write (velocity_unit,diag_fmt) tday, vlon_loc(n), vlat_loc(n), & ARRAY_G(i_loc(n),j_loc(n)) enddo close (velocity_unit) endif endif !----------------------------------------------------------------------- !EOC end subroutine diag_velocity !*********************************************************************** !BOP ! !IROUTINE: cfl_check ! !INTERFACE: subroutine cfl_check ! !DESCRIPTION: ! Finishes the calculation of all cfl numbers by finding the ! global maxima and locations from the local block maxima. ! The final results are written to stdout (or log file) ! and diagnostic output file. ! ! !REVISION HISTORY: ! same as module !EOP !BOC !----------------------------------------------------------------------- ! ! local variables ! !----------------------------------------------------------------------- integer (int_kind) :: & k,iblock ! loop indices integer (int_kind), dimension(1) :: & imax ! max block location real (r8) :: & cfl_temp, &! temp for holding local max cfl_advu, &! max total cfl number for zonal advection cfl_advv, &! max total cfl number for merid advection cfl_advw, &! max total cfl number for vert advection cfl_advt, &! max total cfl number for total advection cfl_vdifft, &! max cfl num for vert tracer diffusion cfl_vdiffu, &! max cfl num for vert momentum diffusion cfl_hdifft, &! max cfl num for horiz tracer diffusion cfl_hdiffu ! max cfl num for horiz momentum diffusion real (r8), dimension(km) :: & ! max for each level cfl_advuk, &! max total cfl number for zonal advection cfl_advvk, &! max total cfl number for merid advection cfl_advwk, &! max total cfl number for vert advection cfl_advtk, &! max total cfl number for total advection cfl_vdifftk, &! max cfl num for vert tracer diffusion cfl_vdiffuk, &! max cfl num for vert momentum diffusion cfl_hdifftk, &! max cfl num for horiz tracer diffusion cfl_hdiffuk ! max cfl num for horiz momentum diffusion integer (int_kind), dimension(3) :: & cfladd_temp, &! temp for holding local max address cfladd_advu, &! address of max zonal advection cfl number cfladd_advv, &! address of max merid advection cfl number cfladd_advw, &! address of max vert advection cfl number cfladd_advt, &! address of max total advection cfl number cfladd_vdifft, &! address of max vertical tracer diff cfl num cfladd_vdiffu, &! address of max vertical moment diff cfl num cfladd_hdifft, &! address of max horiz tracer diff cfl num cfladd_hdiffu ! address of max horiz momentum diff cfl num integer (int_kind), dimension(2,km) :: & ! addresses for each level cfladd_advuk, &! address of max zonal advection cfl number cfladd_advvk, &! address of max merid advection cfl number cfladd_advwk, &! address of max vert advection cfl number cfladd_advtk, &! address of max total advection cfl number cfladd_vdifftk, &! address of max vertical tracer diff cfl num cfladd_vdiffuk, &! address of max vertical moment diff cfl num cfladd_hdifftk, &! address of max horiz tracer diff cfl num cfladd_hdiffuk ! address of max horiz momentum diff cfl num character (25) :: & name_temp1, name_temp2 ! temporaries for writing names character (*), parameter :: & cfl_fmt1 = "(' max. cfl number for ',a25,': ',1pe12.4,4x,3i5)", & cfl_fmt2 = '(1pe22.15,1x,1pe22.15,2x,a25)' !----------------------------------------------------------------------- ! ! compute cfl diagnostics only if it is time ! !----------------------------------------------------------------------- if (ldiag_cfl) then !----------------------------------------------------------------------- ! ! initialize max cfl numbers ! !----------------------------------------------------------------------- cfl_advu = -c1 cfl_advv = -c1 cfl_advw = -c1 cfl_advt = -c1 cfl_vdifft = -c1 cfl_vdiffu = -c1 cfl_hdifft = -c1 cfl_hdiffu = -c1 cfl_advuk = -c1 cfl_advvk = -c1 cfl_advwk = -c1 cfl_advtk = -c1 cfl_vdifftk = -c1 cfl_vdiffuk = -c1 cfl_hdifftk = -c1 cfl_hdiffuk = -c1 cfladd_advu = 0 cfladd_advv = 0 cfladd_advw = 0 cfladd_advt = 0 cfladd_vdifft = 0 cfladd_vdiffu = 0 cfladd_hdifft = 0 cfladd_hdiffu = 0 cfladd_advuk = 0 cfladd_advvk = 0 cfladd_advwk = 0 cfladd_advtk = 0 cfladd_vdifftk = 0 cfladd_vdiffuk = 0 cfladd_hdifftk = 0 cfladd_hdiffuk = 0 !----------------------------------------------------------------------- ! ! compute level maximum cfl numbers from block results ! !----------------------------------------------------------------------- do k=1,km !*** zonal advection imax = maxloc(cfl_advuk_block(:,k)) cfl_temp = cfl_advuk_block(imax(1),k) cfl_advuk(k) = global_maxval(cfl_temp) cfladd_temp = 0 if (cfl_temp == cfl_advuk(k)) then cfladd_temp(1:2) = cfladd_advuk_block(:,imax(1),k) endif cfladd_advuk(1,k) = global_maxval(cfladd_temp(1)) cfladd_advuk(2,k) = global_maxval(cfladd_temp(2)) !*** meridional advection imax = maxloc(cfl_advvk_block(:,k)) cfl_temp = cfl_advvk_block(imax(1),k) cfl_advvk(k) = global_maxval(cfl_temp) cfladd_temp = 0 if (cfl_temp == cfl_advvk(k)) then cfladd_temp(1:2) = cfladd_advvk_block(:,imax(1),k) endif cfladd_advvk(1,k) = global_maxval(cfladd_temp(1)) cfladd_advvk(2,k) = global_maxval(cfladd_temp(2)) !*** vertical advection imax = maxloc(cfl_advwk_block(:,k)) cfl_temp = cfl_advwk_block(imax(1),k) cfl_advwk(k) = global_maxval(cfl_temp) cfladd_temp = 0 if (cfl_temp == cfl_advwk(k)) then cfladd_temp(1:2) = cfladd_advwk_block(:,imax(1),k) endif cfladd_advwk(1,k) = global_maxval(cfladd_temp(1)) cfladd_advwk(2,k) = global_maxval(cfladd_temp(2)) !*** total advection imax = maxloc(cfl_advtk_block(:,k)) cfl_temp = cfl_advtk_block(imax(1),k) cfl_advtk(k) = global_maxval(cfl_temp) cfladd_temp = 0 if (cfl_temp == cfl_advtk(k)) then cfladd_temp(1:2) = cfladd_advtk_block(:,imax(1),k) endif cfladd_advtk(1,k) = global_maxval(cfladd_temp(1)) cfladd_advtk(2,k) = global_maxval(cfladd_temp(2)) !*** vertical tracer diffusion imax = maxloc(cfl_vdifftk_block(:,k)) cfl_temp = cfl_vdifftk_block(imax(1),k) cfl_vdifftk(k) = global_maxval(cfl_temp) cfladd_temp = 0 if (cfl_temp == cfl_vdifftk(k)) then cfladd_temp(1:2) = cfladd_vdifftk_block(:,imax(1),k) endif cfladd_vdifftk(1,k) = global_maxval(cfladd_temp(1)) cfladd_vdifftk(2,k) = global_maxval(cfladd_temp(2)) !*** vertical momentum diffusion imax = maxloc(cfl_vdiffuk_block(:,k)) cfl_temp = cfl_vdiffuk_block(imax(1),k) cfl_vdiffuk(k) = global_maxval(cfl_temp) cfladd_temp = 0 if (cfl_temp == cfl_vdiffuk(k)) then cfladd_temp(1:2) = cfladd_vdiffuk_block(:,imax(1),k) endif cfladd_vdiffuk(1,k) = global_maxval(cfladd_temp(1)) cfladd_vdiffuk(2,k) = global_maxval(cfladd_temp(2)) !*** horizontal tracer diffusion imax = maxloc(cfl_hdifftk_block(:,k)) cfl_temp = cfl_hdifftk_block(imax(1),k) cfl_hdifftk(k) = global_maxval(cfl_temp) cfladd_temp = 0 if (cfl_temp == cfl_hdifftk(k)) then cfladd_temp(1:2) = cfladd_hdifftk_block(:,imax(1),k) endif cfladd_hdifftk(1,k) = global_maxval(cfladd_temp(1)) cfladd_hdifftk(2,k) = global_maxval(cfladd_temp(2)) !*** horizontal momentum diffusion imax = maxloc(cfl_hdiffuk_block(:,k)) cfl_temp = cfl_hdiffuk_block(imax(1),k) cfl_hdiffuk(k) = global_maxval(cfl_temp) cfladd_temp = 0 if (cfl_temp == cfl_hdiffuk(k)) then cfladd_temp(1:2) = cfladd_hdiffuk_block(:,imax(1),k) endif cfladd_hdiffuk(1,k) = global_maxval(cfladd_temp(1)) cfladd_hdiffuk(2,k) = global_maxval(cfladd_temp(2)) end do !----------------------------------------------------------------------- ! ! compute global maximum cfl numbers from level maxima ! !----------------------------------------------------------------------- do k = 1,km !*** advective cfl numbers if (cfl_advuk(k) > cfl_advu) then cfl_advu = cfl_advuk(k) cfladd_advu(1) = cfladd_advuk(1,k) cfladd_advu(2) = cfladd_advuk(2,k) cfladd_advu(3) = k endif if (cfl_advvk(k) > cfl_advv) then cfl_advv = cfl_advvk(k) cfladd_advv(1) = cfladd_advvk(1,k) cfladd_advv(2) = cfladd_advvk(2,k) cfladd_advv(3) = k endif if (cfl_advwk(k) > cfl_advw) then cfl_advw = cfl_advwk(k) cfladd_advw(1) = cfladd_advwk(1,k) cfladd_advw(2) = cfladd_advwk(2,k) cfladd_advw(3) = k endif if (cfl_advtk(k) > cfl_advt) then cfl_advt = cfl_advtk(k) cfladd_advt(1) = cfladd_advtk(1,k) cfladd_advt(2) = cfladd_advtk(2,k) cfladd_advt(3) = k endif !*** vertical diffusion cfl numbers if (cfl_vdifftk(k) > cfl_vdifft) then cfl_vdifft = cfl_vdifftk(k) cfladd_vdifft(1) = cfladd_vdifftk(1,k) cfladd_vdifft(2) = cfladd_vdifftk(2,k) cfladd_vdifft(3) = k endif if (cfl_vdiffuk(k) > cfl_vdiffu) then cfl_vdiffu = cfl_vdiffuk(k) cfladd_vdiffu(1) = cfladd_vdiffuk(1,k) cfladd_vdiffu(2) = cfladd_vdiffuk(2,k) cfladd_vdiffu(3) = k endif !*** horizontal diffusion cfl numbers if (cfl_hdifftk(k) > cfl_hdifft) then cfl_hdifft = cfl_hdifftk(k) cfladd_hdifft(1) = cfladd_hdifftk(1,k) cfladd_hdifft(2) = cfladd_hdifftk(2,k) cfladd_hdifft(3) = k endif if (cfl_hdiffuk(k) > cfl_hdiffu) then cfl_hdiffu = cfl_hdiffuk(k) cfladd_hdiffu(1) = cfladd_hdiffuk(1,k) cfladd_hdiffu(2) = cfladd_hdiffuk(2,k) cfladd_hdiffu(3) = k endif end do !----------------------------------------------------------------------- ! ! write results to stdout ! !----------------------------------------------------------------------- if (my_task == master_task) then write(stdout,blank_fmt) name_temp1 = 'zonal advection ' write(stdout,cfl_fmt1) name_temp1, cfl_advu, cfladd_advu(:) name_temp1 = 'meridional advection ' write(stdout,cfl_fmt1) name_temp1, cfl_advv, cfladd_advv(:) name_temp1 = 'vertical advection ' write(stdout,cfl_fmt1) name_temp1, cfl_advw, cfladd_advw(:) name_temp1 = 'total advection ' write(stdout,cfl_fmt1) name_temp1, cfl_advt, cfladd_advt(:) if (cfl_vdifft /= c0) then name_temp1 = 'vertical tracer mixing ' write(stdout,cfl_fmt1) name_temp1, cfl_vdifft, & cfladd_vdifft(:) endif if (cfl_vdifft /= c0) then name_temp1 = 'vertical moment. mixing ' write(stdout,cfl_fmt1) name_temp1, cfl_vdiffu, & cfladd_vdiffu(:) endif name_temp1 = 'horizontal tracer mixing ' write(stdout,cfl_fmt1) name_temp1, cfl_hdifft, cfladd_hdifft(:) name_temp1 = 'horiz. momentum mixing ' write(stdout,cfl_fmt1) name_temp1, cfl_hdiffu, cfladd_hdiffu(:) if (cfl_all_levels) then write(stdout,'(a27)') 'CFL numbers at all levels: ' do k=1,km name_temp1 = 'total advection ' write(stdout,cfl_fmt1) name_temp1, cfl_advtk(k), & cfladd_advtk(:,k), k name_temp1 = 'horizontal tracer mixing ' write(stdout,cfl_fmt1) name_temp1, cfl_hdifftk(k), & cfladd_hdifftk(:,k), k end do endif !----------------------------------------------------------------------- ! ! write results to output file ! !----------------------------------------------------------------------- open(diag_unit, file=diag_outfile, status='old', & position='append') name_temp1 = 'cfl_advu' name_temp2 = 'cfl_advv' write(diag_unit,cfl_fmt2) tday, cfl_advu, name_temp1 write(diag_unit,cfl_fmt2) tday, cfl_advv, name_temp2 name_temp1 = 'cfl_advw' name_temp2 = 'cfl_advt' write(diag_unit,cfl_fmt2) tday, cfl_advw, name_temp1 write(diag_unit,cfl_fmt2) tday, cfl_advt, name_temp2 if (cfl_vdifft /= c0) then name_temp1 = 'cfl_vdifft' name_temp2 = 'cfl_vdiffu' write(diag_unit,cfl_fmt2) tday, cfl_vdifft, name_temp1 write(diag_unit,cfl_fmt2) tday, cfl_vdiffu, name_temp2 endif name_temp1 = 'cfl_hdifft' name_temp2 = 'cfl_hdiffu' write(diag_unit,cfl_fmt2) tday, cfl_hdifft, name_temp1 write(diag_unit,cfl_fmt2) tday, cfl_hdiffu, name_temp2 if (cfl_all_levels) then do k=1,km write(name_temp1,"('cfl_advtk ',i3)") k write(name_temp2,"('cfl_advuk ',i3)") k write(diag_unit,cfl_fmt2) tday,cfl_advtk(k), name_temp1 write(diag_unit,cfl_fmt2) tday,cfl_hdifftk(k),name_temp2 end do endif !----------------------------------------------------------------------- ! ! finished writing max cfl numbers ! !----------------------------------------------------------------------- close(diag_unit) endif ! master task !----------------------------------------------------------------------- ! ! all done ! !----------------------------------------------------------------------- endif ! ldiag_cfl !----------------------------------------------------------------------- !EOC end subroutine cfl_check !*********************************************************************** !BOP ! !IROUTINE: check_KE ! !INTERFACE: function check_KE(ke_thresh) ! !DESCRIPTION: ! Checks the kinetic energy diagnostic to determine whether it ! exceeds a threshold value. Also checks for negative values. ! Exceeding these ranges results in a true value returned. ! ! !REVISION HISTORY: ! same as module ! !INPUT PARAMETERS: real (r8), intent(in) :: & ke_thresh ! threshold value for kinetic energy ! !OUTPUT PARAMETERS: logical (log_kind) :: & check_KE ! true if k.e. exceeds threshold or < 0 !EOP !BOC !----------------------------------------------------------------------- check_KE = .false. ! default value if (diag_ke > ke_thresh .or. diag_ke < 0) check_KE = .true. !----------------------------------------------------------------------- !EOC end function check_KE !*********************************************************************** end module diagnostics !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||