!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module ocean_velocity_diag_mod ! ! Ron Pacanowski ! ! ! S.M. Griffies ! ! ! ! Numerical diagnostics for velocity related quantities. ! ! ! ! This module contains some diagnostics for velocity related quantities. ! ! ! ! ! ! Number of time steps between which compute the diagnostics. ! ! ! Perform energy analysis every n timesteps (1==every time step). ! This diagnostic is expensive, so should be used sparingly during ! production runs. ! ! ! ! Maximum number of land cells where will printout nonzero velocity points. ! Default land_cell_num_max=100. ! ! ! Critical value for Courant number, above which the model will be brought down. ! ! ! Large value for Courant number, above which will write some diagnostics. ! ! ! For printing out lots of information about regions of large Courant numbers. ! ! ! ! Set true to do bitwise exact global sum. When it is false, the global ! sum will be non-bitwise_exact, but will significantly increase efficiency. ! The default value is false. ! ! ! For some debugging purposes ! ! ! #include use constants_mod, only: grav, epsln, c2dbars use diag_manager_mod, only: register_diag_field, need_data, send_data use fms_mod, only: open_namelist_file, check_nml_error, close_file, write_version_number use fms_mod, only: FATAL, stdout, stdlog use mpp_domains_mod, only: mpp_global_sum, BITWISE_EXACT_SUM, NON_BITWISE_EXACT_SUM use mpp_mod, only: mpp_error, mpp_max, mpp_sum, mpp_chksum, mpp_pe use mpp_mod, only: mpp_clock_id, mpp_clock_begin, mpp_clock_end, CLOCK_MODULE use time_manager_mod, only: time_type, increment_time use ocean_bih_friction_mod, only: bih_friction, bih_viscosity_check, bih_reynolds_check use ocean_coriolis_mod, only: coriolis_force use ocean_domains_mod, only: get_local_indices use ocean_form_drag_mod, only: form_drag_accel use ocean_lap_friction_mod, only: lap_friction, lap_viscosity_check, lap_reynolds_check use ocean_momentum_source_mod, only: momentum_source use ocean_operators_mod, only: FAX, FAY, REMAP_BT_TO_BU, BDX_ET, BDY_NT use ocean_parameters_mod, only: DEPTH_BASED, PRESSURE_BASED, missing_value use ocean_parameters_mod, only: ENERGETIC, FINITEVOLUME use ocean_parameters_mod, only: rho0, rho0r use ocean_pressure_mod, only: hydrostatic_pressure, geopotential_anomaly use ocean_types_mod, only: ocean_grid_type, ocean_domain_type use ocean_types_mod, only: ocean_time_type, ocean_time_steps_type, ocean_thickness_type use ocean_types_mod, only: ocean_adv_vel_type, ocean_external_mode_type use ocean_types_mod, only: ocean_velocity_type, ocean_density_type use ocean_util_mod, only: write_timestamp, matrix use ocean_velocity_advect_mod, only: horz_advection_of_velocity, vert_advection_of_velocity use ocean_vert_mix_mod, only: vert_friction use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk1_v, wrk2_v, wrk1_v2d implicit none private #include ! useful constants real :: dbars2Pa real :: p5grav real :: p5grav_rho0r real :: grav_rho0r real :: dtime real :: dtimer real :: acor logical :: used ! for kinetic energy integer :: id_ke_tot=-1 integer :: id_clock_energy_analysis integer :: id_clock_press_conversion integer :: id_clock_press_energy integer :: id_clock_friction_energy integer :: id_clock_vert_dissipation ! for potential energy integer :: id_pe_tot =-1 integer :: id_pe_tot_rel=-1 logical :: potential_energy_first=.true. real, dimension(:), allocatable :: potential_0 ! gravitational potential energy (joules) at initial timestep. ! for pressure energy maps integer :: id_u_dot_grad_pint(2) integer :: id_u_dot_grad_pint_xy integer :: id_ubar_dot_grad_pext(2) integer :: id_ubar_dot_grad_pext_xy ! for friction energy maps integer :: id_u_dot_lapfrict_horz=-1 integer :: id_v_dot_lapfrict_horz=-1 integer :: id_u_dot_bihfrict_horz=-1 integer :: id_v_dot_bihfrict_horz=-1 integer :: id_u_dot_frict_vert=-1 integer :: id_v_dot_frict_vert=-1 integer :: id_ubar_dot_frict_bt=-1 integer :: id_vbar_dot_frict_bt=-1 integer :: id_vert_lap_diss=-1 ! pressure conversion diagnostics written in energy analysis subroutine real :: press_int ! integrated power from u_int dot grad(p) real :: press_ext ! integrated power from u_ext dot grad(p) real :: press_conv ! integrated power converted from u dot grap(p) into other terms real :: press_conv_err ! error in converting pressure work into other terms real :: ucell_mass ! for output integer :: unit=6 type(ocean_grid_type), pointer :: Grd =>NULL() type(ocean_domain_type), pointer :: Dom =>NULL() logical :: module_is_initialized = .FALSE. character(len=128) :: version=& '$Id: ocean_velocity_diag.F90,v 16.0.2.2.84.1.4.2 2009/12/01 21:52:13 smg Exp $' character (len=128) :: tagname = & '$Name: mom4p1_pubrel_dec2009_nnz $' public ocean_velocity_diag_init public ocean_velocity_diagnostics public velocity_change public kinetic_energy public potential_energy public compute_vorticity public pressure_conversion public energy_analysis private pressure_energy private friction_energy private vert_dissipation private velocity_land_cell_check private cfl_check1 private cfl_check2 integer :: global_sum_flag integer :: diag_step = -1 integer :: energy_diag_step = -1 integer :: land_cell_num_max = 100 ! for CFL checks real :: max_cfl_value = 100.0 real :: large_cfl_value = 10.0 logical :: verbose_cfl =.false. logical :: do_bitwise_exact_sum = .false. logical :: debug_this_module = .false. namelist /ocean_velocity_diag_nml/ diag_step, energy_diag_step, do_bitwise_exact_sum, debug_this_module, & max_cfl_value, large_cfl_value, verbose_cfl, land_cell_num_max contains !####################################################################### ! ! ! ! Initialize the ocean_velocity_diag module containing subroutines ! diagnosing velocity related properties of the simulation. These are ! not terms in the equations, but rather they are diagnosed from ! terms. ! ! subroutine ocean_velocity_diag_init(Grid, Domain, Time, Time_steps) type(ocean_grid_type), target, intent(in) :: Grid type(ocean_domain_type), target, intent(in) :: Domain type(ocean_time_type), intent(in) :: Time type(ocean_time_steps_type), intent(in) :: Time_steps integer :: ioun, io_status, ierr integer :: stdoutunit,stdlogunit stdoutunit=stdout();stdlogunit=stdlog() if (module_is_initialized) return module_is_initialized = .TRUE. call write_version_number(version, tagname) ioun = open_namelist_file() read(ioun, ocean_velocity_diag_nml, iostat=io_status) write (stdlogunit, ocean_velocity_diag_nml) write (stdoutunit,'(/)') write (stdoutunit, ocean_velocity_diag_nml) ierr = check_nml_error(io_status,'ocean_velocity_diag_nml') call close_file(ioun) if (diag_step == 0) diag_step = 1 if (energy_diag_step == 0) energy_diag_step = 1 if(do_bitwise_exact_sum) then global_sum_flag = BITWISE_EXACT_SUM else global_sum_flag = NON_BITWISE_EXACT_SUM endif if(debug_this_module) then write(stdoutunit,*) & '==>Note: running ocean_velocity_diag_mod with debug_this_module=.true.' endif #ifndef MOM4_STATIC_ARRAYS call get_local_indices(Domain, isd, ied, jsd, jed, isc, iec, jsc, jec) nk = Grid%nk #endif Dom => Domain Grd => Grid p5grav = 0.5*grav grav_rho0r = grav*rho0r p5grav_rho0r = 0.5*grav*rho0r dbars2Pa = 1.0/c2dbars dtime = Time_steps%dtime_u dtimer = 1.0/dtime acor = Time_steps%acor ! kinetic energy id_ke_tot = register_diag_field ('ocean_model', 'ke_tot', Time%model_time, & 'Total ke', '10^15 Joules', missing_value=missing_value, range=(/0.0,1000.0/)) ! potential energy id_pe_tot = register_diag_field ('ocean_model', 'pe_tot', & Time%model_time, 'potential energy', '10^15 Joules',& missing_value=missing_value, range=(/0.0,1e20/)) id_pe_tot_rel = register_diag_field ('ocean_model', 'pe_tot_rel', & Time%model_time, 'potential energy - initial potential energy', '10^15 Joules',& missing_value=missing_value, range=(/0.0,1e20/)) ! pressure energetics id_u_dot_grad_pint(1) = register_diag_field ('ocean_model', 'u_dot_grad_pint_x', & Grd%vel_axes_uv(1:3), Time%model_time, & 'i-current times 3d i-pressure force', 'Watt', & missing_value=missing_value, range=(/-1e15, 1e15/)) id_u_dot_grad_pint(2) = register_diag_field ('ocean_model', 'u_dot_grad_pint_y', & Grd%vel_axes_uv(1:3), Time%model_time, & 'j-current times 3d j-pressure force', 'Watt', & missing_value=missing_value, range=(/-1e15, 1e15/)) id_u_dot_grad_pint_xy = register_diag_field ('ocean_model', 'u_dot_grad_pint_xy', & Grd%vel_axes_uv(1:3), Time%model_time, & 'current dot 3d pressure force', 'Watt', & missing_value=missing_value, range=(/-1e15, 1e15/)) id_ubar_dot_grad_pext(1) = register_diag_field ('ocean_model', 'ubar_dot_grad_pext_x', & Grd%vel_axes_uv(1:2), Time%model_time, & '2d i-current times 2d i-pressure force', 'Watt', & missing_value=missing_value, range=(/-1e15, 1e15/)) id_ubar_dot_grad_pext(2) = register_diag_field ('ocean_model', 'ubar_dot_grad_pext_y', & Grd%vel_axes_uv(1:2), Time%model_time, & '2d j-current times 2d j-pressure force', 'Watt', & missing_value=missing_value, range=(/-1e15, 1e15/)) id_ubar_dot_grad_pext_xy = register_diag_field ('ocean_model', 'ubar_dot_grad_pext_xy', & Grd%vel_axes_uv(1:2), Time%model_time, & '2d current dot 2d pressure force', 'Watt', & missing_value=missing_value, range=(/-1e15, 1e15/)) ! friction energetics id_u_dot_lapfrict_horz = register_diag_field ('ocean_model', 'u_dot_lapfrict_horz', & Grd%vel_axes_uv(1:3), Time%model_time, & 'i-current times horizontal laplacian friction', 'Watt', & missing_value=missing_value, range=(/-1e15, 1e15/)) id_v_dot_lapfrict_horz = register_diag_field ('ocean_model', 'v_dot_lapfrict_horz', & Grd%vel_axes_uv(1:3), Time%model_time, & 'j-current times horizontal laplacian friction', 'Watt', & missing_value=missing_value, range=(/-1e15, 1e15/)) id_u_dot_bihfrict_horz = register_diag_field ('ocean_model', 'u_dot_bihfrict_horz', & Grd%vel_axes_uv(1:3), Time%model_time, & 'i-current times horizontal biharmonic friction', 'Watt', & missing_value=missing_value, range=(/-1e15, 1e15/)) id_v_dot_bihfrict_horz = register_diag_field ('ocean_model', 'v_dot_bihfrict_horz', & Grd%vel_axes_uv(1:3), Time%model_time, & 'j-current times horizontal biharmonic friction', 'Watt', & missing_value=missing_value, range=(/-1e15, 1e15/)) id_u_dot_frict_vert = register_diag_field ('ocean_model', 'u_dot_frict_vert', & Grd%vel_axes_uv(1:3), Time%model_time, & 'i-current times vertical friction', 'Watt', & missing_value=missing_value, range=(/-1e15, 1e15/)) id_v_dot_frict_vert = register_diag_field ('ocean_model', 'v_dot_frict_vert', & Grd%vel_axes_uv(1:3), Time%model_time, & 'j-current times vertical friction', 'Watt', & missing_value=missing_value, range=(/-1e15, 1e15/)) id_ubar_dot_frict_bt = register_diag_field ('ocean_model', 'ubar_dot_frict_bt', & Grd%vel_axes_uv(1:2), Time%model_time, & 'i-ubar current times horizontal friction applied just to ubar', 'Watt', & missing_value=missing_value, range=(/-1e15, 1e15/)) id_vbar_dot_frict_bt = register_diag_field ('ocean_model', 'vbar_dot_frict_bt', & Grd%vel_axes_uv(1:2), Time%model_time, & 'j-vbar current times horizontal friction applied just to vbar', 'Watt', & missing_value=missing_value, range=(/-1e15, 1e15/)) id_vert_lap_diss = register_diag_field ('ocean_model', 'vert_lap_diss', Grd%vel_axes_uv(1:3),& Time%model_time, 'Energy dissipation from vertical Laplacian friction', 'Watt', & missing_value=missing_value, range=(/-1.e20,1.e20/), & standard_name='ocean_kinetic_energy_dissipation_per_unit_area_due_to_vertical_friction') allocate (potential_0(nk)) id_clock_energy_analysis = mpp_clock_id('(Velocity diag: energy analysis) ' ,grain=CLOCK_MODULE) id_clock_press_conversion = mpp_clock_id('(Velocity diag: press convert) ' ,grain=CLOCK_MODULE) id_clock_press_energy = mpp_clock_id('(Velocity diag: press energy) ' ,grain=CLOCK_MODULE) id_clock_friction_energy = mpp_clock_id('(Velocity diag: friction energy)' ,grain=CLOCK_MODULE) id_clock_vert_dissipation = mpp_clock_id('(Velocity diag: vert dissipation' ,grain=CLOCK_MODULE) end subroutine ocean_velocity_diag_init ! NAME="ocean_velocity_diag_init" !####################################################################### ! ! ! ! Call diagnostics related to the velocity. ! ! subroutine ocean_velocity_diagnostics(Time, Thickness, Dens, Velocity) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_density_type), intent(in) :: Dens type(ocean_velocity_type), intent(inout) :: Velocity type(time_type) :: next_time real :: pe_tot = 0.0 real :: pe_tot_rel = 0.0 real :: ke_tot = 0.0 ! diagnostics that do not send_data to diag_manager if (diag_step > 0) then if(mod(Time%itt,diag_step) == 0) then call lap_viscosity_check call lap_reynolds_check(Time, Velocity) call bih_viscosity_check call bih_reynolds_check(Time, Velocity) call velocity_land_cell_check(Time, Velocity) call write_timestamp(Time%model_time) call kinetic_energy(Time, Thickness, Velocity, ke_tot, .false., .true.) call potential_energy(Time, Thickness, Dens, pe_tot, pe_tot_rel, .false., .true.) call cfl_check1(Time, Thickness, Velocity) call cfl_check2(Time, Thickness, Velocity) endif endif ! diagnostics that send_data to diag_manager next_time = increment_time(Time%model_time, int(dtime), 0) if (need_data(id_ke_tot,next_time)) then call kinetic_energy(Time, Thickness, Velocity, ke_tot, .true., .false.) endif if (need_data(id_pe_tot,next_time)) then call potential_energy(Time, Thickness, Dens, pe_tot, pe_tot_rel, .true., .false.) endif end subroutine ocean_velocity_diagnostics ! NAME="ocean_velocity_diagnostics" !####################################################################### ! ! ! ! ! Compute gravitational potential energy (Joules) relative to z=0 ! taken with respect to the value at the initial time step. ! ! ! subroutine potential_energy (Time, Thickness, Dens, pe_tot, pe_tot_rel, diag_flag, write_flag) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_density_type), intent(in) :: Dens real, intent(inout) :: pe_tot real, intent(inout) :: pe_tot_rel logical, intent(in), optional :: diag_flag logical, intent(in), optional :: write_flag logical :: send_diagnostics logical :: write_diagnostics real, dimension(nk) :: potential integer :: k integer :: taup1 integer :: stdoutunit stdoutunit=stdout() if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error in ocean_velocity_diag_mod(potential_energy): module needs initialization ') endif taup1 = Time%taup1 ! assume send_diagnostics=.true., unless diag_flag says it is false. if (PRESENT(diag_flag)) then send_diagnostics = diag_flag else send_diagnostics = .true. endif ! assume write_diagnostics=.true., unless write_flag says it is false. if (PRESENT(write_flag)) then write_diagnostics = write_flag else write_diagnostics = .true. endif potential(:) = 0.0 do k=1,nk potential(k) = mpp_global_sum(Dom%domain2d, & Grd%dat(:,:)*Thickness%dzt(:,:,k)*Dens%rho(:,:,k,taup1) & *Grd%tmask(:,:,k)*Thickness%depth_zt(:,:,k),global_sum_flag) potential(k) = potential(k)*grav enddo if(potential_energy_first) then potential_energy_first=.false. do k=1,nk potential_0(k) = potential(k) enddo endif ! computing pe_tot relative to initial time step value reduces roundoff errors pe_tot = 0.0 pe_tot_rel = 0.0 do k=1,nk pe_tot = pe_tot + potential(k) pe_tot_rel = pe_tot_rel + (potential(k)-potential_0(k)) enddo if(write_diagnostics) then write(stdoutunit,'(/,1x,a,e20.12)') & 'Gravitational potential energy relative to first time (Joules) = ', pe_tot_rel write(stdoutunit,'(/,1x,a,e20.12)') & 'Gravitational potential energy = ', pe_tot endif if(send_diagnostics) then if (id_pe_tot_rel > 0) used = send_data (id_pe_tot_rel, pe_tot_rel*1e-15, Time%model_time) endif if(send_diagnostics) then if (id_pe_tot > 0) used = send_data (id_pe_tot, pe_tot*1e-15, Time%model_time) endif end subroutine potential_energy ! NAME="potential_energy" !####################################################################### ! ! ! ! Compute global integrated horizontal kinetic energy. ! ! subroutine kinetic_energy (Time, Thickness, Velocity, ke_tot, diag_flag, write_flag) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_velocity_type), intent(in) :: Velocity real, intent(inout) :: ke_tot logical, intent(in), optional :: diag_flag logical, intent(in), optional :: write_flag logical :: send_diagnostics logical :: write_diagnostics real :: mass_u_cell integer :: i, j, k, tau integer :: stdoutunit stdoutunit=stdout() if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error in ocean_velocity_diag_mod(kinetic_energy): module needs initialization ') endif ! assume send_diagnostics=.true., unless diag_flag says it is false. if (PRESENT(diag_flag)) then send_diagnostics = diag_flag else send_diagnostics = .true. endif ! assume write_diagnostics=.true., unless write_flag says it is false. if (PRESENT(write_flag)) then write_diagnostics = write_flag else write_diagnostics = .true. endif tau = Time%tau ! total horizontal kinetic energy(joules) at "tau" ke_tot = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec mass_u_cell = Grd%dau(i,j)*Thickness%rho_dzu(i,j,k,tau)*Grd%umask(i,j,k) ke_tot = ke_tot + 0.5*mass_u_cell*(Velocity%u(i,j,k,1,tau)**2 + Velocity%u(i,j,k,2,tau)**2) enddo enddo enddo call mpp_sum (ke_tot) if(write_diagnostics) then write(stdoutunit,'(1x,a,e20.12)') 'Kinetic energy (Joules) = ', ke_tot endif if(send_diagnostics) then if (id_ke_tot > 0) used = send_data (id_ke_tot, ke_tot*1e-15, Time%model_time) endif end subroutine kinetic_energy ! NAME="kinetic_energy" !####################################################################### ! ! ! ! See if there are any points over land with nonzero ocean velocity ! ! subroutine velocity_land_cell_check(Time, Velocity) implicit none type(ocean_time_type), intent(in) :: Time type(ocean_velocity_type), intent(in) :: Velocity integer :: i, j, k, num, taup1 integer :: stdoutunit stdoutunit=stdout() taup1 = Time%taup1 if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error in ocean_velocity_diag_mod(velocity_land_cell_check): module needs initialization ') endif call write_timestamp(Time%model_time) write (stdoutunit,'(//60x,a/)') ' Land cell summary:' write (stdoutunit,'(1x,a/)')'Locations (if any) where land cell velocity is non-zero...' num = 0 kloop: do k=1,nk do j=jsc,jec do i=isc,iec if (Grd%umask(i,j,k) == 0 .and. (Velocity%u(i,j,k,1,taup1) /= 0.0 .or. Velocity%u(i,j,k,2,taup1) /= 0.0)) then num = num + 1 write (unit,9000) i, j, k, Grd%xu(i,j), Grd%yu(i,j), Grd%zt(k), Velocity%u(i,j,k,1,taup1), Velocity%u(i,j,k,2,taup1) endif if(num > land_cell_num_max) then write (stdoutunit,'(1x,a/)')'More than land_cell_num_max land cells with non-zero velocity. Stop the check.' exit kloop endif enddo enddo enddo kloop call mpp_max(num) if (num > 0) then call mpp_error(FATAL,'==>Error: found nonzero ocean velocity over land points. ') endif 9000 format(/' " =>Error: Land cell at (i,j,k) = ','(',i4,',',i4,',',i4,'),',& ' (lon,lat,dpt) = (',f7.2,',',f7.2,',',f7.0,'m) has u=',e10.3, ' v=',e10.3) end subroutine velocity_land_cell_check ! NAME="velocity_land_cell_check" !####################################################################### ! ! ! ! Determine the number of points that have large single-time step ! changes in the absolute value of the velocity. ! ! subroutine velocity_change(Time, Velocity,velocity_change_max,velocity_change_max_num) implicit none type(ocean_time_type), intent(in) :: Time type(ocean_velocity_type), intent(in) :: Velocity real, intent(in) :: velocity_change_max integer, intent(in) :: velocity_change_max_num real :: abschange(2) integer :: i, j, k, n, num integer :: taum1, tau, taup1 integer :: stdoutunit stdoutunit=stdout() if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error in ocean_velocity_diag_mod(velocity_change): module needs initialization ') endif tau = Time%tau taum1 = Time%taum1 taup1 = Time%taup1 call write_timestamp(Time%model_time) write (stdoutunit,'(//60x,a/)') & ' Velocity change summary (test for leap-frog noise):' write (stdoutunit,'(1x,a,e12.6/)') & 'Locations (if any) where abs(0.5*(u(taup1)+u(taum1))-u(tau)) (m/s) > ',velocity_change_max num = 0 do k=1,nk do j=jsc,jec do i=isc,iec if(Grd%umask(i,j,k) > 0) then do n=1,2 abschange(n) = abs(0.5*(Velocity%u(i,j,k,n,taup1)+Velocity%u(i,j,k,n,taum1))-Velocity%u(i,j,k,n,tau)) enddo if (abschange(1) > velocity_change_max .or. abschange(2) > velocity_change_max) then num = num + 1 write (unit,9000) i, j, k, Grd%xu(i,j), Grd%yu(i,j), Grd%zt(k), abschange(1), abschange(2) endif endif enddo enddo enddo call mpp_max(num) if (num > velocity_change_max_num) then call mpp_error(FATAL, & 'Found too many points where velocity changed too much over single time step.') endif 9000 format(/' " =>Error: Ocean at (i,j,k) = ','(',i4,',',i4,',',i4,'),', & ' (lon,lat,dpt) = (',f7.2,',',f7.2,',',f7.0,'m) has excessive change(u)=',e10.3, & ' or excessive change(v)=',e10.3) end subroutine velocity_change ! NAME="velocity_change" !####################################################################### ! ! ! ! Compute z-component to vorticity on B-grid. ! ! subroutine compute_vorticity(Time, Velocity, vorticity_z) type(ocean_time_type), intent(in) :: Time type(ocean_velocity_type), intent(in) :: Velocity real, dimension(isd:,jsd:,:), intent(inout) :: vorticity_z integer :: i, j, k integer :: tau if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error in ocean_velocity_diag_mod(compute_vorticity): module needs initialization ') endif tau = Time%tau do k=1,nk do j=jsc,jec do i=isc,iec vorticity_z(i,j,k) = 0.5*((Velocity%u(i,j,k,2,tau)-Velocity%u(i-1,j,k,2,tau) )/Grd%dxtn(i,j) & +(Velocity%u(i,j-1,k,2,tau)-Velocity%u(i-1,j-1,k,2,tau))/Grd%dxtn(i,j-1) & -(Velocity%u(i,j,k,1,tau)-Velocity%u(i,j-1,k,1,tau) )/Grd%dxte(i,j) & -(Velocity%u(i-1,j,k,1,tau)-Velocity%u(i-1,j-1,k,1,tau))/Grd%dxte(i,j-1) & ) enddo enddo enddo end subroutine compute_vorticity ! NAME="compute_vorticity" !####################################################################### ! ! ! ! Perform pressure conversion error analysis. This analysis should be ! computed prior to update_ucell_thickness since we need to use dzu(tau) ! here rather than dzu(taup1). ! ! subroutine pressure_conversion (Time, Thickness, rho, Dens, Ext_mode, & Adv_vel, Velocity, vert_coordinate_class) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_density_type), intent(in) :: Dens type(ocean_adv_vel_type), intent(in) :: Adv_vel type(ocean_external_mode_type), intent(in) :: Ext_mode type(ocean_velocity_type), intent(inout) :: Velocity real, dimension(isd:,jsd:,:), intent(in) :: rho integer, intent(in) :: vert_coordinate_class integer :: tau, taum1, taup1 integer :: i, j, k, n, kb real, dimension(isd:ied,jsd:jed,2) :: ubar real :: column_mass, boxarea real :: uint, uext, fx real :: term, term1, term2 if(id_u_dot_grad_pint(1) > 0 .or. id_u_dot_grad_pint(2) > 0 .or. & id_ubar_dot_grad_pext(1) > 0 .or. id_ubar_dot_grad_pext(2) > 0 .or. & id_u_dot_grad_pint_xy > 0 .or. id_ubar_dot_grad_pext_xy > 0) then call pressure_energy(Time, Thickness, Ext_mode, Velocity, vert_coordinate_class) endif if (.not. (energy_diag_step > 0 .and. mod(Time%itt, energy_diag_step) == 0)) return call mpp_clock_begin(id_clock_press_conversion) if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_velocity_diag_mod (pressure_conversion): module must be initialized') endif tau = Time%tau taum1 = Time%taum1 taup1 = Time%taup1 ucell_mass = mpp_global_sum(Dom%domain2d, Grd%dau(:,:)*Thickness%mass_u(:,:,tau), global_sum_flag) ! truncate out the lower order bits if using non_bitwise_reproducible sums if (global_sum_flag .eq. NON_BITWISE_EXACT_SUM) ucell_mass = real(ucell_mass, kind=FLOAT_KIND) press_conv = 0.0 press_int = 0.0 press_ext = 0.0 press_conv_err = 0.0 ! place density anomaly in wrk2 since it is needed frequently do k=1,nk do j=jsd,jed do i=isd,ied wrk2(i,j,k) = Grd%tmask(i,j,k)*(rho(i,j,k) - rho0) enddo enddo enddo ! place mass tendency in wrk3 as it is needed frequently do k=1,nk do j=jsd,jed do i=isd,ied wrk3(i,j,k) = & Grd%tmask(i,j,k)*(Thickness%rho_dzt_tendency(i,j,k)-Thickness%mass_source(i,j,k)) enddo enddo enddo ! place U-cell area into wrk4 do k=1,nk do j=jsd,jed do i=isd,ied wrk4(i,j,k) = Grd%umask(i,j,k)*Grd%dau(i,j) enddo enddo enddo ! vertically averaged horizontal velocity do j=jsd,jed do i=isd,ied ubar(i,j,1) = Grd%umask(i,j,1)*Ext_mode%udrho(i,j,1,tau)/(Thickness%mass_u(i,j,tau)+epsln) ubar(i,j,2) = Grd%umask(i,j,1)*Ext_mode%udrho(i,j,2,tau)/(Thickness%mass_u(i,j,tau)+epsln) enddo enddo ! work done by baroclinic pressure ! press_force [=] Pa [=] kg/(m*s^2) ! press_conv [=] kg*m^2/s^3 = Watt do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) uext = ubar(i,j,n) uint = Velocity%u(i,j,k,n,tau) - uext term = Velocity%press_force(i,j,k,n)*boxarea press_int = press_int + uint*term press_ext = press_ext + uext*term enddo enddo enddo enddo ! work done by depth independent pressure on depth independent flow ! grad_ps [=] Pa/m = kg/(m^2*s^2) ! grad_mb [=] m/s^2 ! press_conv [=] kg*m^2/s^3 = Watt ! as this term also appears in the conversion of u dot grad(p), ! we include it with press_conv. if(vert_coordinate_class==DEPTH_BASED) then do j=jsc,jec do i=isc,iec column_mass = wrk4(i,j,1)*Thickness%mass_u(i,j,tau) term = rho0r*(Ext_mode%grad_ps(i,j,1)*ubar(i,j,1) & + Ext_mode%grad_ps(i,j,2)*ubar(i,j,2)) press_ext = press_ext - term*column_mass press_conv = press_conv - term*column_mass enddo enddo elseif(vert_coordinate_class==PRESSURE_BASED) then do j=jsc,jec do i=isc,iec column_mass = wrk4(i,j,1)*Thickness%mass_u(i,j,tau) term = Ext_mode%grad_anompb(i,j,1)*ubar(i,j,1) & + Ext_mode%grad_anompb(i,j,2)*ubar(i,j,2) press_ext = press_ext - term*column_mass press_conv = press_conv - term*column_mass enddo enddo endif ! compute the conversion terms if(vert_coordinate_class==DEPTH_BASED) then wrk1(:,:,:) = hydrostatic_pressure(Thickness, wrk2(:,:,:)) if(Thickness%method==ENERGETIC) then do j=jsc,jec do i=isc,iec kb = Grd%kmt(i,j) if (kb /= 0) then k = 1 term = -rho0r*Grd%tmask(i,j,k)*Grd%dat(i,j)*wrk1(i,j,k) & *(wrk3(i,j,k)+Adv_vel%wrho_bt(i,j,k-1)) press_conv = press_conv + term fx = Grd%dat(i,j)*p5grav do k=2,kb term1 = -fx*Adv_vel%wrho_bt(i,j,k-1)*Thickness%dzwt(i,j,k-1) & *(wrk2(i,j,k-1) + wrk2(i,j,k)) term2 = -Grd%dat(i,j)*wrk1(i,j,k)*wrk3(i,j,k) press_conv = press_conv + rho0r*Grd%tmask(i,j,k)*(term1 + term2) enddo endif enddo enddo else do j=jsc,jec do i=isc,iec kb = Grd%kmt(i,j) if (kb /= 0) then k = 1 term = -rho0r*Grd%tmask(i,j,k)*Grd%dat(i,j)*wrk1(i,j,k) & *(wrk3(i,j,k)+Adv_vel%wrho_bt(i,j,k-1)) press_conv = press_conv + term fx = Grd%dat(i,j)*grav do k=2,kb term1 = -fx*Adv_vel%wrho_bt(i,j,k-1) & *(Thickness%dztlo(i,j,k-1)*wrk2(i,j,k-1) + Thickness%dztup(i,j,k)*wrk2(i,j,k)) term2 = -Grd%dat(i,j)*wrk1(i,j,k)*wrk3(i,j,k) press_conv = press_conv + rho0r*Grd%tmask(i,j,k)*(term1 + term2) enddo endif enddo enddo endif ! gradient of geopotential (area sum vanishes on k-level if geodepth_zt is independent of (i,j)) do k=1,nk wrk1_v(:,:,k,1) = BDX_ET((Adv_vel%uhrho_et(:,:,k))*FAX(wrk2(:,:,k)))*Grd%dat(:,:) wrk1_v(:,:,k,2) = BDY_NT((Adv_vel%vhrho_nt(:,:,k))*FAY(wrk2(:,:,k)))*Grd%dat(:,:) fx = grav*rho0r do j=jsc,jec do i=isc,iec term = -fx*(wrk1_v(i,j,k,1) + wrk1_v(i,j,k,2))*Thickness%geodepth_zt(i,j,k) press_conv = press_conv + term enddo enddo enddo elseif(vert_coordinate_class==PRESSURE_BASED) then wrk1(:,:,:) = geopotential_anomaly(Thickness, wrk2(:,:,:)) if(Thickness%method==ENERGETIC) then do j=jsc,jec do i=isc,iec kb = Grd%kmt(i,j) if (kb /= 0) then k = 1 term = -Grd%tmask(i,j,k)*Grd%dat(i,j)*wrk1(i,j,k) & *(wrk3(i,j,k)+Adv_vel%wrho_bt(i,j,k-1)) press_conv = press_conv + term fx = Grd%dat(i,j)*p5grav_rho0r do k=2,kb term1 = -fx*Adv_vel%wrho_bt(i,j,k-1)*Thickness%dzwt(i,j,k-1) & *(wrk2(i,j,k-1) + wrk2(i,j,k)) term2 = -Grd%dat(i,j)*wrk1(i,j,k)*wrk3(i,j,k) press_conv = press_conv + Grd%tmask(i,j,k)*(term1 + term2) enddo endif enddo enddo else do j=jsc,jec do i=isc,iec kb = Grd%kmt(i,j) if (kb /= 0) then k = 1 term = -Grd%tmask(i,j,k)*Grd%dat(i,j)*wrk1(i,j,k) & *(wrk3(i,j,k)+Adv_vel%wrho_bt(i,j,k-1)) press_conv = press_conv + term fx = Grd%dat(i,j)*grav_rho0r do k=2,kb term1 = -fx*Adv_vel%wrho_bt(i,j,k-1) & *(Thickness%dztlo(i,j,k-1)*wrk2(i,j,k-1) + & Thickness%dztup(i,j,k) *wrk2(i,j,k)) term2 = -Grd%dat(i,j)*wrk1(i,j,k)*wrk3(i,j,k) press_conv = press_conv + Grd%tmask(i,j,k)*(term1 + term2) enddo endif enddo enddo endif do k=1,nk do j=jsd,jed do i=isd,iec wrk2_v(i,j,k,1) = (Dens%pressure_at_depth(i+1,j,k) & -Dens%pressure_at_depth(i,j,k))*dbars2Pa enddo enddo do j=jsd,jec do i=isd,ied wrk2_v(i,j,k,2) = (Dens%pressure_at_depth(i,j+1,k) & -Dens%pressure_at_depth(i,j,k))*dbars2Pa enddo enddo enddo do k=1,nk wrk1_v(:,:,k,1) = FAY(FAX(wrk2(:,:,k))*wrk2_v(:,:,k,1))*Grd%dyu(:,:) wrk1_v(:,:,k,2) = FAX(FAY(wrk2(:,:,k))*wrk2_v(:,:,k,2))*Grd%dxu(:,:) do j=jsc,jec do i=isc,iec term = rho0r*Thickness%dzu(i,j,k) & *(Velocity%u(i,j,k,1,tau)*wrk1_v(i,j,k,1) + Velocity%u(i,j,k,2,tau)*wrk1_v(i,j,k,2)) press_conv = press_conv + term enddo enddo enddo endif call mpp_sum(press_int) call mpp_sum(press_ext) call mpp_sum(press_conv) press_conv_err = press_int + press_ext - press_conv call mpp_clock_end(id_clock_press_conversion) end subroutine pressure_conversion ! NAME="pressure_conversion" !####################################################################### ! ! ! ! Diagnose u dot grad(p) for diagnostic purposes. These maps ! when summed over all grid points will result in an energy ! that is equal to pint+pext as computed in pressure_conversion. ! The advantage is that here we can look at maps of these ! contributions. ! ! subroutine pressure_energy (Time, Thickness, Ext_mode, Velocity, vert_coordinate_class) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_external_mode_type), intent(in) :: Ext_mode type(ocean_velocity_type), intent(inout) :: Velocity integer, intent(in) :: vert_coordinate_class integer :: tau integer :: i, j, k, n real, dimension(isd:ied,jsd:jed,2) :: ubar real :: column_mass, boxarea real :: term call mpp_clock_begin(id_clock_press_energy) tau = Time%tau wrk1_v(:,:,:,:) = 0.0 wrk1_v2d(:,:,:) = 0.0 ! place U-cell area into wrk4 do k=1,nk do j=jsd,jed do i=isd,ied wrk4(i,j,k) = Grd%umask(i,j,k)*Grd%dau(i,j) enddo enddo enddo ! vertically averaged horizontal velocity do j=jsd,jed do i=isd,ied ubar(i,j,1) = Grd%umask(i,j,1)*Ext_mode%udrho(i,j,1,tau)/(Thickness%mass_u(i,j,tau)+epsln) ubar(i,j,2) = Grd%umask(i,j,1)*Ext_mode%udrho(i,j,2,tau)/(Thickness%mass_u(i,j,tau)+epsln) enddo enddo ! work done by horizontal baroclinic gradients ! press_force [=] Pa [=] kg/(m*s^2) ! wrk1_v [=] wrk2_v [=] kg*m^2/s^3 = Watt do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) term = Velocity%press_force(i,j,k,n)*boxarea wrk1_v(i,j,k,n) = Velocity%u(i,j,k,n,tau)*term enddo enddo enddo enddo ! work done by depth independent pressure on depth independent flow ! grad_ps [=] Pa/m = kg/(m^2*s^2) ! grad_mb [=] m/s^2 ! wrk2_v [=] kg*m^2/s^3 = Watt if(vert_coordinate_class==DEPTH_BASED) then do j=jsc,jec do i=isc,iec column_mass = wrk4(i,j,1)*Thickness%mass_u(i,j,tau) term = rho0r*Ext_mode%grad_ps(i,j,1)*ubar(i,j,1) wrk1_v2d(i,j,1) = -term*column_mass term = rho0r*Ext_mode%grad_ps(i,j,2)*ubar(i,j,2) wrk1_v2d(i,j,2) = -term*column_mass enddo enddo elseif(vert_coordinate_class==PRESSURE_BASED) then do j=jsc,jec do i=isc,iec column_mass = wrk4(i,j,1)*Thickness%mass_u(i,j,tau) term = Ext_mode%grad_anompb(i,j,1)*ubar(i,j,1) wrk1_v2d(i,j,1) = -term*column_mass term = Ext_mode%grad_anompb(i,j,2)*ubar(i,j,2) wrk1_v2d(i,j,2) = -term*column_mass enddo enddo endif if (id_u_dot_grad_pint(1) > 0) used = send_data (id_u_dot_grad_pint(1), wrk1_v(:,:,:,1), & Time%model_time, rmask=Grd%umask(:,:,:), & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) if (id_u_dot_grad_pint(2) > 0) used = send_data (id_u_dot_grad_pint(2), wrk1_v(:,:,:,2), & Time%model_time, rmask=Grd%umask(:,:,:), & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) if (id_u_dot_grad_pint_xy > 0) used = send_data (id_u_dot_grad_pint_xy, & wrk1_v(:,:,:,1)+wrk1_v(:,:,:,2), & Time%model_time, rmask=Grd%umask(:,:,:), & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) if (id_ubar_dot_grad_pext(1) > 0) used = send_data (id_ubar_dot_grad_pext(1), wrk1_v2d(:,:,1), & Time%model_time, rmask=Grd%umask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) if (id_ubar_dot_grad_pext(2) > 0) used = send_data (id_ubar_dot_grad_pext(2), wrk1_v2d(:,:,2), & Time%model_time, rmask=Grd%umask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) if (id_ubar_dot_grad_pext_xy > 0) used = send_data (id_ubar_dot_grad_pext_xy, & wrk1_v2d(:,:,1)+wrk1_v2d(:,:,2), & Time%model_time, rmask=Grd%umask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) call mpp_clock_end(id_clock_press_energy) end subroutine pressure_energy ! NAME="pressure_energy" !####################################################################### ! ! ! ! Diagnose u dot Friction for diagnostic purposes. ! ! NOTE: ! ! A) DO NOT split into baroclinic and barotropic pieces. Just compute ! u dot F using full velocity field u. Otherwise, the calculation emulates ! that done in energy_analysis subroutine. ! ! B) DO NOT remove the effects from bottom drag and from surface stress. ! The reason is that bottom drag and surface stress are incorporated to ! the vertical friction operator, even when doing vertical friction ! implicitly in time. So it is tough to remove these effects in an ! explicitl diagnostic manner. So the u dot vertical friction piece ! includes BOTH surface and bottom stress. ! ! ! subroutine friction_energy (Time, Thickness, Adv_vel, Ext_mode, Velocity, visc_cbu) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_adv_vel_type), intent(in) :: Adv_vel type(ocean_external_mode_type), intent(in) :: Ext_mode type(ocean_velocity_type), intent(inout) :: Velocity real, dimension(isd:,jsd:,:), intent(in) :: visc_cbu integer :: tau integer :: i, j, k, n real, dimension(isd:ied,jsd:jed,2) :: ubar real :: boxarea, term call mpp_clock_begin(id_clock_friction_energy) tau = Time%tau wrk2_v(:,:,:,:) = 0.0 ! for holding u dot friction ! (do NOT use wrk1_v, as that is used in friction modules) wrk1_v2d(:,:,:) = 0.0 ! for holding ubar dot horizontal stress acting JUST on ubar ! place U-cell area into wrk4 do k=1,nk do j=jsd,jed do i=isd,ied wrk4(i,j,k) = Grd%umask(i,j,k)*Grd%dau(i,j) enddo enddo enddo ! vertically averaged horizontal velocity do j=jsd,jed do i=isd,ied ubar(i,j,1) = Grd%umask(i,j,1)*Ext_mode%udrho(i,j,1,tau)/(Thickness%mass_u(i,j,tau)+epsln) ubar(i,j,2) = Grd%umask(i,j,1)*Ext_mode%udrho(i,j,2,tau)/(Thickness%mass_u(i,j,tau)+epsln) enddo enddo ! work done by laplacian friction on horizontal velocity field. ! lap_friction returns thickness weighted and density weighted ! friction with dimensions (kg/m^3)*(m^2/s^2) = N/m^2 = Pa ! wrk1_v [=] kg*m^2/s^3 = Watt call lap_friction(Time, Thickness, Adv_vel, Velocity, energy_analysis_step=.true.) do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) term = Velocity%wrkv(i,j,k,n)*boxarea wrk2_v(i,j,k,n) = Velocity%u(i,j,k,n,tau)*term enddo enddo enddo enddo if (id_u_dot_lapfrict_horz > 0) used = send_data (id_u_dot_lapfrict_horz, wrk2_v(:,:,:,1), & Time%model_time, rmask=Grd%umask(:,:,:), & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) if (id_v_dot_lapfrict_horz > 0) used = send_data (id_v_dot_lapfrict_horz, wrk2_v(:,:,:,2), & Time%model_time, rmask=Grd%umask(:,:,:), & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) ! work done by biharmonic friction ! lap_friction returns thickness weighted and density weighted ! friction with dimensions (kg/m^3)*(m^2/s^2) = N/m^2 = Pa ! eng(5) [=] kg*m^2/s^3 = Watt call bih_friction(Time, Thickness, Velocity, energy_analysis_step=.true.) do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) term = Velocity%wrkv(i,j,k,n)*boxarea wrk2_v(i,j,k,n) = Velocity%u(i,j,k,n,tau)*term enddo enddo enddo enddo if (id_u_dot_bihfrict_horz > 0) used = send_data (id_u_dot_bihfrict_horz, wrk2_v(:,:,:,1), & Time%model_time, rmask=Grd%umask(:,:,:), & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) if (id_v_dot_bihfrict_horz > 0) used = send_data (id_v_dot_bihfrict_horz, wrk2_v(:,:,:,2), & Time%model_time, rmask=Grd%umask(:,:,:), & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) ! contribution from lap_friction acting just on barotropic. ! friction with dimensions (kg/m^3)*(m^2/s^2) = N/m^2 = Pa ! wrk1_v2d [=] kg*m^2/s^3 = Watt do n=1,2 do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,1) term = Velocity%lap_friction_bt(i,j,n)*boxarea wrk1_v2d(i,j,n) = ubar(i,j,n)*term enddo enddo enddo ! contribution from bih_friction acting just on the barotropic. ! friction with dimensions (kg/m^3)*(m^2/s^2) = N/m^2 = Pa ! eng(5) [=] kg*m^2/s^3 = Watt do n=1,2 do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,1) term = Velocity%lap_friction_bt(i,j,n)*boxarea wrk1_v2d(i,j,n) = wrk1_v2d(i,j,n) + ubar(i,j,n)*term enddo enddo enddo if (id_ubar_dot_frict_bt > 0) used = send_data (id_ubar_dot_frict_bt, wrk1_v2d(:,:,1), & Time%model_time, rmask=Grd%umask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) if (id_vbar_dot_frict_bt > 0) used = send_data (id_vbar_dot_frict_bt, wrk1_v2d(:,:,2), & Time%model_time, rmask=Grd%umask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) ! work done by vertical friction, including wind stress and bottom stress. ! implicit acceleration has been saved from previous call to ocean_implicit_accel. ! the implicit piece include contribution from visc_cbu_form_drag, as well as visc_cbu. ! vert_frict returns thickness weighted and density weighted vertical friction (kg/m^3)*(m^2/s^2) ! wrk2_v [=] kg*m^2/s^3 = Watt call vert_friction(Time, Thickness, Velocity, visc_cbu, energy_analysis_step=.true.) do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) term = (Velocity%wrkv(i,j,k,n)+ Velocity%vfrict_impl(i,j,k,n))*boxarea wrk2_v(i,j,k,n) = Velocity%u(i,j,k,n,tau)*term enddo enddo enddo enddo if (id_u_dot_frict_vert > 0) used = send_data (id_u_dot_frict_vert, wrk2_v(:,:,:,1), & Time%model_time, rmask=Grd%umask(:,:,:), & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) if (id_v_dot_frict_vert > 0) used = send_data (id_v_dot_frict_vert, wrk2_v(:,:,:,2), & Time%model_time, rmask=Grd%umask(:,:,:), & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) call mpp_clock_end(id_clock_friction_energy) end subroutine friction_energy ! NAME="friction_energy" !####################################################################### ! ! ! ! Diagnose dissipation from vertical friction due just to viscosity. ! ! Assumptions: ! ! 1/ Ignore bottom drag here...just concerned with viscosity. ! ! 2/ Assume vertical friction is handled implicitly in time. ! ! ! subroutine vert_dissipation (Time, Thickness, Velocity, visc_cbu, visc_cbu_form_drag) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_velocity_type), intent(in) :: Velocity real, dimension(isd:,jsd:,:), intent(in) :: visc_cbu real, dimension(isd:,jsd:,:,:), intent(in) :: visc_cbu_form_drag integer :: itime integer :: i, j, k, n real :: du_dz type(time_type) :: next_time call mpp_clock_begin(id_clock_vert_dissipation) itime = Time%taup1 next_time = increment_time(Time%model_time, int(dtime), 0) if (need_data(id_vert_lap_diss,next_time)) then wrk1(:,:,:) = 0.0 do n=1,2 do k=1,nk-1 do j=jsc,jec do i=isc,iec du_dz = Grd%umask(i,j,k+1) & *(Velocity%u(i,j,k,n,itime)-Velocity%u(i,j,k+1,n,itime))/Thickness%dzwu(i,j,k) wrk1(i,j,k) = wrk1(i,j,k) & -rho0*Thickness%dzu(i,j,k)*Grd%dau(i,j)*du_dz*du_dz & *(visc_cbu(i,j,k)+visc_cbu_form_drag(i,j,k,n)) enddo enddo enddo enddo if(id_vert_lap_diss>0) then used = send_data (id_vert_lap_diss, wrk1(:,:,:), & Time%model_time, rmask=Grd%umask(:,:,:), & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) endif endif call mpp_clock_end(id_clock_vert_dissipation) end subroutine vert_dissipation ! NAME="vert_dissipation" !####################################################################### ! ! ! ! Perform energy analysis by taking scalar product of horizontal ! velocity with the velocity equations and integrating over the ocean volume. ! ! Pressure conversions have already been computed in pressure_conversion ! subroutine. It is necessary to perform that analysis earlier than ! the call to update_ucell_thickness inside ocean_model_mod, whereas ! the remaining elements in the energy analysis can be called at the ! end of the update for velocity. ! ! subroutine energy_analysis (Time, Thickness, Ext_mode, Adv_vel, Dens, & pme, river, upme, uriver, visc_cbu, visc_cbu_form_drag, & Velocity) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_adv_vel_type), intent(in) :: Adv_vel type(ocean_density_type), intent(in) :: Dens type(ocean_external_mode_type), intent(in) :: Ext_mode type(ocean_velocity_type), intent(inout) :: Velocity real, dimension(isd:,jsd:), intent(in) :: pme real, dimension(isd:,jsd:), intent(in) :: river real, dimension(isd:,jsd:,:), intent(in) :: upme real, dimension(isd:,jsd:,:), intent(in) :: uriver real, dimension(isd:,jsd:,:), intent(in) :: visc_cbu real, dimension(isd:,jsd:,:,:), intent(in) :: visc_cbu_form_drag integer :: i, j, k, n integer :: tau, taum1, taup1, tau_m0 real :: boxarea real :: uint, uext, term, vel real :: ke_tot, pe_tot, pe_tot_rel, ke_per_mass real, dimension(isd:ied,jsd:jed,2) :: ubar real, dimension(isd:ied,jsd:jed) :: pme_u real, dimension(isd:ied,jsd:jed) :: river_u real, dimension(13) :: engint ! internal mode energy integral components (Watt) real, dimension(13) :: engext ! external mode energy integral components (Watt) real :: compression ! power from compression of grid cell and/or mass sources real :: enleak ! difference between integrated advection of momentum ! and power from compression and water forcing integer :: stdoutunit stdoutunit=stdout() ! get 3d maps for energetics from friction if(id_u_dot_lapfrict_horz > 0 .or. id_v_dot_lapfrict_horz > 0 .or. & id_u_dot_bihfrict_horz > 0 .or. id_v_dot_bihfrict_horz > 0 .or. & id_u_dot_frict_vert > 0 .or. id_v_dot_frict_vert > 0 .or. & id_ubar_dot_frict_bt > 0 .or. id_vbar_dot_frict_bt > 0) then call friction_energy(Time, Thickness, Adv_vel, Ext_mode, Velocity, visc_cbu) endif ! get energy dissipation from vertical friction call vert_dissipation(Time, Thickness, Velocity, visc_cbu, visc_cbu_form_drag) if (.not. (energy_diag_step > 0 .and. mod(Time%itt, energy_diag_step) == 0)) return call mpp_clock_begin(id_clock_energy_analysis) if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_velocity_diag_mod (energy_analysis): module must be initialized') endif tau = Time%tau taum1 = Time%taum1 taup1 = Time%taup1 tau_m0 = Time%tau_m0 compression = 0.0 engint(:) = 0.0 engext(:) = 0.0 ! place mass tendency in wrk3 as it is needed frequently do k=1,nk do j=jsd,jed do i=isd,ied wrk3(i,j,k) = & Grd%tmask(i,j,k)*(Thickness%rho_dzt_tendency(i,j,k)-Thickness%mass_source(i,j,k)) enddo enddo enddo ! place U-cell area into wrk4 do k=1,nk do j=jsd,jed do i=isd,ied wrk4(i,j,k) = Grd%umask(i,j,k)*Grd%dau(i,j) enddo enddo enddo ! vertically averaged horizontal velocity do j=jsd,jed do i=isd,ied ubar(i,j,1) = Grd%umask(i,j,1)*Ext_mode%udrho(i,j,1,tau)/(Thickness%mass_u(i,j,tau)+epsln) ubar(i,j,2) = Grd%umask(i,j,1)*Ext_mode%udrho(i,j,2,tau)/(Thickness%mass_u(i,j,tau)+epsln) enddo enddo if (debug_this_module) then write(stdoutunit,*) ' ' write(stdoutunit,*) '===chksums in energy analysis===' write(stdoutunit,*) 'dtime(seconds) = ',dtime, 'and dtimer = ',dtimer wrk1(isc:iec,jsc:jec,:) = Thickness%rho_dzu(isc:iec,jsc:jec,:,taum1)*Grd%umask(isc:iec,jsc:jec,:) write(stdoutunit,*) 'Thickness%rho_dzu(taum1) chksum = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:)) wrk1(isc:iec,jsc:jec,:) = Thickness%rho_dzu(isc:iec,jsc:jec,:,tau)*Grd%umask(isc:iec,jsc:jec,:) write(stdoutunit,*) 'Thickness%rho_dzu(tau) chksum = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:)) wrk1(isc:iec,jsc:jec,:) = Thickness%rho_dzu(isc:iec,jsc:jec,:,taup1)*Grd%umask(isc:iec,jsc:jec,:) write(stdoutunit,*) 'Thickness%rho_dzu(taup1) chksum = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:)) wrk1(isc:iec,jsc:jec,1) = Thickness%mass_u(isc:iec,jsc:jec,tau)*Grd%umask(isc:iec,jsc:jec,1) write(stdoutunit,*) 'Thickness%mass_u(tau) chksum = ', mpp_chksum(wrk1(isc:iec,jsc:jec,1)) wrk1(isc:iec,jsc:jec,:) = Velocity%u(isc:iec,jsc:jec,:,1,taum1)*Grd%umask(isc:iec,jsc:jec,:) write(stdoutunit,*) 'Velocity%u(n=1,taum1) chksum = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:)) wrk1(isc:iec,jsc:jec,:) = Velocity%u(isc:iec,jsc:jec,:,1,tau)*Grd%umask(isc:iec,jsc:jec,:) write(stdoutunit,*) 'Velocity%u(n=1,tau) chksum = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:)) wrk1(isc:iec,jsc:jec,:) = Velocity%u(isc:iec,jsc:jec,:,1,taup1)*Grd%umask(isc:iec,jsc:jec,:) write(stdoutunit,*) 'Velocity%u(n=1,taup1) chksum = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:)) wrk1(isc:iec,jsc:jec,:) = Velocity%u(isc:iec,jsc:jec,:,2,taum1)*Grd%umask(isc:iec,jsc:jec,:) write(stdoutunit,*) 'Velocity%u(n=2,taum1) chksum = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:)) wrk1(isc:iec,jsc:jec,:) = Velocity%u(isc:iec,jsc:jec,:,2,tau)*Grd%umask(isc:iec,jsc:jec,:) write(stdoutunit,*) 'Velocity%u(n=2,tau) chksum = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:)) wrk1(isc:iec,jsc:jec,:) = Velocity%u(isc:iec,jsc:jec,:,2,taup1)*Grd%umask(isc:iec,jsc:jec,:) write(stdoutunit,*) 'Velocity%u(n=2,taup1) chksum = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:)) wrk1(isc:iec,jsc:jec,1) = ubar(isc:iec,jsc:jec,1)*Grd%umask(isc:iec,jsc:jec,1) write(stdoutunit,*) 'ubar(n=1) chksum = ', mpp_chksum(wrk1(isc:iec,jsc:jec,1)) wrk1(isc:iec,jsc:jec,1) = ubar(isc:iec,jsc:jec,2)*Grd%umask(isc:iec,jsc:jec,1) write(stdoutunit,*) 'ubar(n=2) chksum = ', mpp_chksum(wrk1(isc:iec,jsc:jec,1)) write(stdoutunit,*)'=================================' endif ! scalar product of u with the time rate of change ! eng(1) [=] kg*m^2/s^3 = Watt do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) uext = ubar(i,j,n) uint = Velocity%u(i,j,k,n,tau) - uext term = (Thickness%rho_dzu(i,j,k,taup1)*Velocity%u(i,j,k,n,taup1) & -Thickness%rho_dzu(i,j,k,taum1)*Velocity%u(i,j,k,n,taum1)) & *dtimer*boxarea engint(1) = engint(1) + uint*term engext(1) = engext(1) + uext*term enddo enddo enddo enddo ! work done by horizontal advection ! horz_advection_of_velocity [=] thickness and density weighted advection ! which has dimensions (kg/m^3)*(m^2/s^2) = Pa ! eng(2) [=] kg*m^2/s^3 = Watt call horz_advection_of_velocity(Time, Thickness, Adv_vel, Velocity, & energy_analysis_step=.true.) do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) uext = ubar(i,j,n) uint = Velocity%u(i,j,k,n,tau) - uext term = Velocity%wrkv(i,j,k,n)*boxarea engint(2) = engint(2) - uint*term engext(2) = engext(2) - uext*term enddo enddo enddo enddo ! work done by vertical advection ! vertical_advection_of_velocity [=] thickness and density weighted advection ! which as dimensions (kg/m^3)*(m^2/s^2) = Pa ! eng(3) [=] kg*m^2/s^3 = Watt call vert_advection_of_velocity(Time, Adv_vel, Velocity, & pme, river, upme, uriver, & energy_analysis_step=.true.) do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) uext = ubar(i,j,n) uint = Velocity%u(i,j,k,n,tau) - uext term = Velocity%wrkv(i,j,k,n)*boxarea engint(3) = engint(3) - uint*term engext(3) = engext(3) - uext*term enddo enddo enddo enddo ! work done by laplacian friction ! lap_friction returns thickness weighted and density weighted ! friction with dimensions (kg/m^3)*(m^2/s^2) = N/m^2 = Pa ! eng(4) [=] kg*m^2/s^3 = Watt call lap_friction(Time, Thickness, Adv_vel, Velocity, energy_analysis_step=.true.) do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) uext = ubar(i,j,n) uint = Velocity%u(i,j,k,n,tau) - uext term = Velocity%wrkv(i,j,k,n)*boxarea engint(4) = engint(4) + uint*term engext(4) = engext(4) + uext*term enddo enddo enddo enddo ! add contribution from lap_friction acting just on the barotropic. ! friction with dimensions (kg/m^3)*(m^2/s^2) = N/m^2 = Pa ! eng(4) [=] kg*m^2/s^3 = Watt do n=1,2 do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,1) uext = ubar(i,j,n) term = Velocity%lap_friction_bt(i,j,n)*boxarea engext(4) = engext(4) + uext*term enddo enddo enddo ! work done by biharmonic friction ! lap_friction returns thickness weighted and density weighted ! friction with dimensions (kg/m^3)*(m^2/s^2) = N/m^2 = Pa ! eng(5) [=] kg*m^2/s^3 = Watt call bih_friction(Time, Thickness, Velocity, energy_analysis_step=.true.) do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) uext = ubar(i,j,n) uint = Velocity%u(i,j,k,n,tau) - uext term = Velocity%wrkv(i,j,k,n)*boxarea engint(5) = engint(5) + uint*term engext(5) = engext(5) + uext*term enddo enddo enddo enddo ! add contribution from bih_friction acting just on the barotropic. ! friction with dimensions (kg/m^3)*(m^2/s^2) = N/m^2 = Pa ! eng(5) [=] kg*m^2/s^3 = Watt do n=1,2 do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,1) uext = ubar(i,j,n) term = Velocity%bih_friction_bt(i,j,n)*boxarea engext(5) = engext(5) + uext*term enddo enddo enddo ! work done by vertical friction (contributions from wind and bottom drag removed below) ! implicit acceleration has been saved from previous call to ocean_implicit_accel. ! Greatbatch form drag implemented through vertical viscosity is part of implict accel. ! vert_frict returns thickness weighted and density weighted vertical friction (kg/m^3)*(m^2/s^2) ! eng(6) [=] kg*m^2/s^3 = Watt call vert_friction(Time, Thickness, Velocity, visc_cbu, energy_analysis_step=.true.) do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) uext = ubar(i,j,n) uint = Velocity%u(i,j,k,n,tau) - uext term = (Velocity%wrkv(i,j,k,n)+ Velocity%vfrict_impl(i,j,k,n))*boxarea engint(6) = engint(6) + uint*term engext(6) = engext(6) + uext*term enddo enddo enddo enddo ! work done by wind stress and bottom drag ! smf [=] N/m^2 = Pa ! bmf [=] N/m^2 = Pa ! eng(7) [=] kg*m^2/s^3 = Watt do n=1,2 ! wind stress k = 1 do j=jsc,jec do i=isc,iec uext = ubar(i,j,n) uint = Velocity%u(i,j,k,n,tau) - uext term = wrk4(i,j,k)*Velocity%smf(i,j,n) engint(7) = engint(7) + uint*term engext(7) = engext(7) + uext*term enddo enddo ! bottom stress do j=jsc,jec do i=isc,iec k = Grd%kmu(i,j) if (k /= 0) then uext = ubar(i,j,n) uint = Velocity%u(i,j,k,n,tau) - uext term = wrk4(i,j,k)*Velocity%bmf(i,j,n) engint(8) = engint(8) - uint*term engext(8) = engext(8) - uext*term endif enddo enddo enddo ! subtract wind stress and bottom stress contributions ! from vertical friction to facilitate diagnostics engint(6) = engint(6) - engint(7) - engint(8) engext(6) = engext(6) - engext(7) - engext(8) ! work due to Coriolis force (zero on B_grid only when acor=0.0) ! coriolis returns thickness weighted and density weighted Coriolis force (kg/m^3)*(m^2/s^2) ! eng(9) [=] kg*m^2/s^3 = Watt call coriolis_force(Time, Thickness, Velocity, energy_analysis_step=.true.) do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) uext = ubar(i,j,n) uint = Velocity%u(i,j,k,n,tau) - uext term = Velocity%wrkv(i,j,k,n)*boxarea engint(9) = engint(9) + uint*term engext(9) = engext(9) + uext*term enddo enddo enddo enddo ! power (Watt) from surface water fluxes entering (+power) or leaving (-power) ocean pme_u(:,:) = REMAP_BT_TO_BU(pme) river_u(:,:) = REMAP_BT_TO_BU(river) k = 1 do n=1,2 do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) vel = Velocity%u(i,j,k,n,tau) uext = ubar(i,j,n) uint = vel - uext term = pme_u(i,j)*upme(i,j,n) + river_u(i,j)*uriver(i,j,n) term = term*boxarea engint(10) = engint(10) + uint*term engext(10) = engext(10) + uext*term enddo enddo enddo ! contribution (Watt) to kinetic energy conversion from surface water fluxes k = 1 do n=1,2 do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) vel = Velocity%u(i,j,k,n,tau) uext = ubar(i,j,n) uint = vel - uext term = pme_u(i,j)*upme(i,j,n)+river_u(i,j)*uriver(i,j,n) term = 0.5*term*boxarea engint(11) = engint(11) + uint*term engext(11) = engext(11) + uext*term enddo enddo enddo ! contribution (Watt) to kinetic energy from momentum sources call momentum_source(Time, Thickness, Velocity, energy_analysis_step=.true.) do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) uext = ubar(i,j,n) uint = Velocity%u(i,j,k,n,tau) - uext term = Velocity%wrkv(i,j,k,n)*boxarea engint(12) = engint(12) + uint*term engext(12) = engext(12) + uext*term enddo enddo enddo enddo ! contribution (Watt) to kinetic energy from parameterized form drag via Aiki scheme call form_drag_accel(Time, Thickness, Velocity, energy_analysis_step=.true.) do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) uext = ubar(i,j,n) uint = Velocity%u(i,j,k,n,tau) - uext term = Velocity%wrkv(i,j,k,n)*boxarea engint(13) = engint(13) + uint*term engext(13) = engext(13) + uext*term enddo enddo enddo enddo ! rho_dzu_tendency and mass source are converted from kinetic energy advection. ! these terms arise from the ability of the grid cell volume/mass to change due ! to changes in its thickness or through the advent of volume/mass sources. ! note that some of the surface water work is encompassed here, hence the need ! to compute engint(11) and engext(11). do k=1,nk wrk1(:,:,k) = Grd%umask(:,:,k)*REMAP_BT_TO_BU(wrk3(:,:,k)) do j=jsc,jec do i=isc,iec boxarea = wrk4(i,j,k) ke_per_mass = 0.5*(Velocity%u(i,j,k,1,tau)*Velocity%u(i,j,k,1,tau) & +Velocity%u(i,j,k,2,tau)*Velocity%u(i,j,k,2,tau)) term = boxarea*ke_per_mass*wrk1(i,j,k) compression = compression + term enddo enddo enddo call mpp_sum(compression) call mpp_sum(engint, size(engint(:))) call mpp_sum(engext, size(engext(:))) enleak = engint(2) + engint(3) - engint(11) + engext(2) + engext(3) - engext(11)- compression call kinetic_energy(Time, Thickness, Velocity, ke_tot, .false., .false.) call potential_energy(Time, Thickness, Dens, pe_tot, pe_tot_rel, .false., .false.) write(stdoutunit,*) ' ' write (stdoutunit,'(//60x,a)') & ' Globally integrated energy analysis over model time step ' call write_timestamp(Time%model_time) if(acor > 0) then write (stdoutunit,'(1x,/a)') & ' ==> Note: acor > 0 means Coriolis force will not conserve energy.' endif if(Grd%tripolar) then write (stdoutunit,'(/1x,/a)') & ' ==> NOTE: Energy conversion errors are nontrivial when using the tripolar grid. ' endif write(stdoutunit,'(/,1x,a,e20.12)') & 'Potential energy relative to first time (Joules) = ', pe_tot_rel write(stdoutunit,'(1x,a,e20.12)') & 'Potential energy (Joules) = ', pe_tot write(stdoutunit,'(1x,a,e20.12)') & 'Kinetic energy (Joules) = ', ke_tot write (stdoutunit,9100) & 'Globally integrated U dot momentum eqns (Watt)', ucell_mass, Grd%ucella(1) write (stdoutunit,9101) ' time rate of change ', engint(1)+engext(1), engint(1), engext(1) write (stdoutunit,9101) ' horizontal advection ', engint(2)+engext(2), engint(2), engext(2) write (stdoutunit,9101) ' vertical advection ', engint(3)+engext(3), engint(3), engext(3) write (stdoutunit,9101) ' pressure force ', press_int+press_ext, press_int, press_ext write (stdoutunit,9101) ' laplacian friction ', engint(4)+engext(4), engint(4), engext(4) write (stdoutunit,9101) ' biharmonic friction ', engint(5)+engext(5), engint(5), engext(5) write (stdoutunit,9101) ' vertical friction ', engint(6)+engext(6), engint(6), engext(6) write (stdoutunit,9101) ' work by wind ', engint(7)+engext(7), engint(7), engext(7) write (stdoutunit,9101) ' work by bottom drag ', engint(8)+engext(8), engint(8), engext(8) write (stdoutunit,9101) ' coriolis forces ', engint(9)+engext(9), engint(9), engext(9) write (stdoutunit,9101) ' work by water flux ', engint(10)+engext(10), engint(10), engext(10) write (stdoutunit,9101) ' work by momentum src ', engint(12)+engext(12), engint(12), engext(12) write (stdoutunit,9101) ' work by Aiki drag ', engint(13)+engext(13), engint(13), engext(13) write (stdoutunit,9110) press_conv, press_int+press_ext, press_conv_err, compression, enleak 9100 format(/1x,a,1x,/,' ocean mass =',e16.9,' kg',/, ' ocean surface area =',e16.9,& ' m^2',//,32x,'total(Watt) internal(Watt) external(Watt)') 9101 format(a21,3(3x,es24.17)) 9102 format(a21,3x,es24.17) 9110 format(/ ' power contributed by pressure conversion and mass tendencies = ',es24.17,& /,' power from horizontal pressure force -u dot force(p) = ',es24.17,& /,' power mismatch between -u dot grad(p) and its conversions = ',es24.17,& /,' power from grid cell mass changes due to compresibility and/or sources = ',es24.17,& /,' power mismatch between -u dot grad (v u) and its conversions = ',es24.17/) call mpp_clock_end(id_clock_energy_analysis) end subroutine energy_analysis ! NAME="energy_analysis" !####################################################################### ! ! ! ! Perform the first of two CFL checks on horizontal velocity. ! ! Vectorized version from Russell.Fiedler@csiro.au computes cfl ! values at a single latitude. The location of the maximum at this ! latitude is calculated via the maxloc() intrinsic. The maximum ! value for this processor is then updated if necessary. ! ! subroutine cfl_check1(Time, Thickness, Velocity) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_velocity_type), intent(in) :: Velocity real :: cflup,cflvp ! percent of cfl criteria reached by velocity component real :: cflum,cflvm ! velocity component which comes closest to its cfl criteria integer :: icflu,jcflu,kcflu ! "i","j","k" coordinate of "cflum" integer :: icflv,jcflv,kcflv ! "i","j","k" coordinate of "cflvm" real :: dtmax real :: cflup0, cflvp0 real :: fudge integer :: i, j, k, tau ! work array containing cfl values real , dimension(isc:iec) :: tempcfl ! array required for maxloc() intrinsic function integer, dimension(1) :: itemp integer :: stdoutunit stdoutunit=stdout() if (.not.module_is_initialized) then call mpp_error(FATAL,& '==>Error in ocean_velocity_diag_mod (cfl_check1): module needs initialization ') endif tau = Time%tau write (stdoutunit,'(/60x,a/)') ' Horizontal CFL summary I: ' write (stdoutunit,'(1x,a/)')'Location of largest horizontal CFL.' ! to distinguish processors fudge = 1 + 1.e-12*mpp_pe() icflu=isc; jcflu=jsc; kcflu=1; cflup=epsln; cflum=0.0 icflv=isc; jcflv=jsc; kcflv=1; cflvp=epsln; cflvm=0.0 dtmax = dtime do k=1,nk do j=jsc,jec do i=isc,iec tempcfl(i)=abs(Velocity%u(i,j,k,1,tau)/Grd%dxu(i,j)) enddo itemp=maxloc(tempcfl)+isc-1 if(tempcfl(itemp(1))>cflup) then cflup=tempcfl(itemp(1)) icflu=itemp(1) jcflu=j kcflu=k endif do i=isc,iec tempcfl(i)=abs(Velocity%u(i,j,k,2,tau)/Grd%dyu(i,j)) enddo itemp=maxloc(tempcfl)+isc-1 if(tempcfl(itemp(1))>cflvp) then cflvp=tempcfl(itemp(1)) icflv=itemp(1) jcflv=j kcflv=k endif enddo enddo ! multiply by 100.0 to convert to percentages ! multiply by dtmax to convert to dimensionless CFL number ! multipby cflwtp by rho0r for dimensional reasons as well cflup = 100.0*cflup*dtmax cflvp = 100.0*cflvp*dtmax cflum = Velocity%u(icflu,jcflu,kcflu,1,tau) cflvm = Velocity%u(icflv,jcflv,kcflv,2,tau) cflup = cflup*fudge cflvp = cflvp*fudge cflup0 = cflup cflvp0 = cflvp call mpp_max(cflup) call mpp_max(cflvp) if (cflup == cflup0) then write (unit,'(a,g10.3,a,f7.2,a,g10.3,a,i4,a,i4,a,i3,a,a,f12.3,a,f12.3,a,f10.3,a)')& ' u (',cflum,' m/s) is ',cflup/fudge,' % of CFL (',abs(100.0*cflum/(cflup/fudge)), & ' m/s) at (i,j,k) = (',icflu+Dom%ioff,',',jcflu+Dom%joff,',',kcflu,'),', & ' (lon,lat,thk) = (',Grd%xu(icflu,jcflu),',',Grd%yu(icflu,jcflu),',', & Thickness%depth_zu(icflu,jcflu,kcflu),' m)' write (unit,'(a,e12.3,a,e12.3,a,f10.3)')& ' where grid is (m) (dxu,dyu,dzu) = (', & Grd%dxu(icflu,jcflu),',',Grd%dyu(icflu,jcflu),',',Thickness%dzu(icflu,jcflu,kcflu) endif if (cflvp == cflvp0) then write (unit,'(a,g10.3,a,f7.2,a,g10.3,a,i4,a,i4,a,i3,a,a,f12.3,a,f12.3,a,f10.3,a)') & ' v (',cflvm,' m/s) is ',cflvp/fudge,' % of CFL (',abs(100.0*cflvm/(cflvp/fudge)), & ' m/s) at (i,j,k) = (',icflv+Dom%ioff,',',jcflv+Dom%joff,',',kcflv,'),', & ' (lon,lat,thk) = (',Grd%xu(icflv,jcflv),',',Grd%yu(icflv,jcflv),',', & Thickness%depth_zu(icflv,jcflv,kcflv),' m)' write (unit,'(a,e12.3,a,e12.3,a,f10.3)') & ' where grid is (m) (dxu,dyu,dzu) = (', & Grd%dxu(icflv,jcflv),',',Grd%dyu(icflv,jcflv),',',& Thickness%dzu(icflv,jcflv,kcflv) endif end subroutine cfl_check1 ! NAME="cfl_check1" !####################################################################### ! ! ! ! Perform the second of two CFL checks on horizontal velocity. ! ! Bring the model down if too many large Courant numbers detected. ! subroutine cfl_check2(Time, Thickness, Velocity) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_velocity_type), intent(in) :: Velocity real :: dtmax, cflu, cflv real :: umax, pcflu, vmax, pcflv, scl real, dimension(isd:ied,jsd:jed) :: u, v integer :: i, j, k integer :: is, ie, js, je integer :: tau integer :: stdoutunit stdoutunit=stdout() if (.not.module_is_initialized) then call mpp_error(FATAL,& '==>Error ocean_velocity_diag_mod (cfl_check2): module needs initialization ') endif tau = Time%tau write (stdoutunit,'(/60x,a/)') ' Horizontal CFL summary II: ' write (stdoutunit,'(1x,a,f5.2/)') & 'Locations (if any) where horizontal Courant number exceeds ', large_cfl_value dtmax = dtime do k=1,nk u(:,:) = Velocity%u(:,:,k,1,tau) v(:,:) = Velocity%u(:,:,k,2,tau) do j=jsc,jec do i=isc,iec cflu = abs((dtmax/Grd%dxu(i,j))*u(i,j)) cflv = abs((dtmax/Grd%dyu(i,j))*v(i,j)) if (cflu >= large_cfl_value .or. cflv >= large_cfl_value) then write (unit,'(/,a,i4,a1,i3,a,i3,a,f6.3)') & ' Note: CFL horizontal velocity limit exceeded at (i,j,k) = (',& i+Dom%ioff, ',',j+Dom%joff,',',k,') by factor =',large_cfl_value umax = Grd%dxu(i,j)/dtmax pcflu = abs(100.0*u(i,j)/umax) vmax = Grd%dyu(i,j)/dtmax pcflv = abs(100.0*v(i,j)/vmax) write (unit,'(a,f8.2,a,g15.8,a)') & ' u reached ', pcflu,' % of the local CFL limit (',umax,' m/s)' write (unit,'(a,f8.2,a,g15.8,a)') & ' v reached ', pcflv,' % of the local CFL limit (',vmax,' m/s)' write (unit,'(a,e12.3,a,e12.3,a,f10.3)') & ' where the grid cell has dimensions (metres) (dxu,dyu,dzu) = (',& Grd%dxu(i,j),',',Grd%dyu(i,j),',',Thickness%dzu(i,j,k) endif if (cflu >= max_cfl_value .or. cflv >= max_cfl_value) then write (unit,'(/,a,i4,a1,i3,a,i3,a,f6.3)') & ' Note: CFL horizontal velocity limit exceeded at (i,j,k) = (',& i+Dom%ioff, ',',j+Dom%joff,',',k,') by factor =',max_cfl_value write(*,'(/a)') & '==>Error in ocean_velocity_diag: max_cfl_value exceeded for horizontal velocity.' write(*,'(a)') & ' The model is bringing itself down in order for the user to determine ' write(*,'(a)') & ' whether the simulation is going to produce Inf or NaN, or whether it ' write(*,'(a)') & ' will successfully run through the present phase with large Courant numbers.' write(*,'(a)') & ' It is often the case that early spin-ups produce very large Courant' write(*,'(a)') & ' numbers, with the model later adjusting to stable behaviour with modest' write(*,'(a)') & ' Courant numbers. To test if the model will reach stable behaviour,' write(*,'(a)') & ' try significally increasing "max_cfl_value" in ocean_velocity_diag_nml.' if(verbose_cfl) then write(*,'(/a,f6.3/)') & ' Information about region where horizontal Courant number is larger than',& max_cfl_value is = max(isc,i-3) ie = min(iec,i+3) js = max(jsc,j-3) je = min(jec,j+3) scl = 0.0 write (*,9100)'u (m/s)', Time%itt, k, Thickness%depth_zu(is,js,k), & Grd%xu(is,j), Grd%xu(ie,j), Grd%yu(i,js), Grd%yu(i,je), scl call matrix (u(is:ie,js:je), is, ie, js, je, scl) write (*,9100) 'v (m/s)', Time%itt, k, Thickness%depth_zu(is,js,k), & Grd%xu(is,j), Grd%xu(ie,j), Grd%yu(i,js), Grd%yu(i,je), scl call matrix (v(is:ie,js:je), is, ie, js, je, scl) endif call mpp_error(FATAL,& '==>Error in ocean_velocity_diag_mod: max_cfl_value exceeded.') endif enddo enddo enddo 9100 format(1x,a12,1x,'ts=',i10,1x,',k=',i3,', z(m)=',f8.1,', lon(deg):',f6.2,' --> ',f6.2,& ', lat(deg):',f6.2,' --> ',f6.2,', scaling=',1pg10.3) end subroutine cfl_check2 ! NAME="cfl_check2" end module ocean_velocity_diag_mod