!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! 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_mod ! ! R.C. Pacanowski ! ! ! A. Rosati ! ! ! ! S.M. Griffies ! ! ! ! M.J. Harrison ! ! ! ! Time step velocity field. ! ! ! ! This module steps the velocity field forward in time. ! ! ! ! ! ! Durran, Numerical Methods for Wave Equations in Geophysical ! Fluid Dynamics (1999). ! ! ! ! R.C. Pacanowski and S.M. Griffies, The MOM3 Manual (1999). ! NOAA/Geophysical Fluid Dynamics Laboratory ! ! ! ! S.M. Griffies, R.C. Pacanowski, R.M. Schmidt, and V. Balaji ! Tracer Conservation with an Explicit Free Surface Method for ! Z-coordinate Ocean Models ! Monthly Weather Review (2001) vol 129 pages 1081--1098 ! ! ! ! S.M. Griffies, M.J. Harrison, R.C. Pacanowski, and ! A. Rosati, A Technical Guide to MOM4 (2004). ! NOAA/Geophysical Fluid Dynamics Laboratory ! ! ! ! S.M. Griffies, Fundamentals of Ocean Climate Models (2004). ! Princeton University Press. ! ! ! ! S.M. Griffies (2005) ! Elements of mom4p1 ! ! ! ! ! ! ! ! For debugging. ! ! ! Set true to write a restart. False setting only for rare ! cases where wish to benchmark model without measuring the cost ! of writing restarts and associated chksums. ! Default is write_a_restart=.true. ! ! ! ! When updating the velocity, this method first computes uprime ! as the updated velocity minus the barotropic pressure gradient. ! This approach is motivated from the rigid lid approach, in which ! the surface pressure was never used to update the barotropic ! fields. ! With the explicit free surface, we have the choice to ! update the full velocity field, with the barotropic contributions ! to the pressure field resulting from a time average in the ! external mode algorithm. This approach is for testing only, ! and it has been found to be unstable for many cases. ! ! update_velocity_via_uprime=.true. uses the older aproach, ! in which the udrho,vdrho fields are taken from the external ! mode module. ! ! update_velocity_via_uprime=.false. only takes the time averaged ! pressure from the external mode, and thus updates the full ! velocity and so recomputes the udrho,vdrho fields. ! ! Default update_velocity_via_uprime=.true. ! The case of update_velocity_via_uprime=.false. is for testing only! ! ! ! ! For running with time independent constant velocity. ! For use with idealized cases. Default=.false. ! ! ! For running with use_constant_velocity. ! Set the i-velocity component to this value. ! Default constant_u=0.0 ! ! ! For running with use_constant_velocity. ! Set the j-velocity component to this value. ! Default constant_v=0.0 ! ! ! ! For debugging. Will freeze the baroclinic velocity fields. ! ! ! For debugging. Will not use explicit-a part of the tendency. ! ! ! For debugging. Will not use explicit-b part of the tendency. ! ! ! For debugging. Will not use implicit part of the tendency. ! ! ! ! Truncate the baroclinic velocity to a maximum value. Useful for cases where ! the initial spin-up initiates spuriously large model velocities that would ! otherwise cause the model to blow-up. Also can be used as a very simple ! "polar filter" in cases where have spherical coordinates and wish to avoid ! using the traditional polar filters. ! ! ! Speed above which will truncate the baroclinic velocity ! ! ! Latitude poleward of which we truncate the velocity. Useful in cases ! when wish to truncate the velocity only in polar regions. Default is 0.0 ! ! ! For verbose printout ! ! ! ! Maximum internal gravity wave speed--used for diagnosing conservative ! estimate of stable time steps. ! ! ! ! Dimensionless parameter for 2nd order Adams-Bashforth implementation of ! velocity advection. Values between 0.5 and 1.0 are recommended. ! Value of 0.5 leads to second order accurate, but it is formally ! weakly unstable (Durran, Section 2.3.4). ! ! ! For a third order treatment of the velocity advection. ! This is stable and so needs no temporal dissipation ! (Section 2.3.6 of Durran). This is the model default. ! ! use constants_mod, only: epsln, grav use diag_manager_mod, only: register_diag_field, send_data use fms_mod, only: write_version_number, check_nml_error, close_file, open_namelist_file use fms_mod, only: file_exist, FATAL, WARNING, NOTE use fms_io_mod, only: field_size, restart_file_type, register_restart_field use fms_io_mod, only: save_restart, restore_state, reset_field_pointer use mpp_domains_mod, only: mpp_update_domains, mpp_global_sum, BGRID_NE, BITWISE_EXACT_SUM use mpp_mod, only: mpp_chksum, mpp_error, mpp_pe, mpp_min, mpp_max, mpp_sum use mpp_mod, only: mpp_broadcast, stdout, stdlog use ocean_bih_friction_mod, only: bih_friction use ocean_coriolis_mod, only: coriolis_force, coriolis_force_implicit 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 use ocean_momentum_source_mod, only: momentum_source use ocean_obc_mod, only: ocean_obc_update_boundary use ocean_operators_mod, only: DIV_UD use ocean_parameters_mod, only: TWO_LEVEL, THREE_LEVEL use ocean_parameters_mod, only: missing_value, rho0r use ocean_parameters_mod, only: DEPTH_BASED use ocean_pressure_mod, only: pressure_force use ocean_types_mod, only: ocean_domain_type, ocean_grid_type, ocean_thickness_type use ocean_types_mod, only: ocean_time_type, ocean_time_steps_type use ocean_types_mod, only: ocean_density_type, ocean_thickness_type, ocean_velocity_type use ocean_types_mod, only: ocean_adv_vel_type, ocean_external_mode_type, ocean_options_type use ocean_util_mod, only: write_timestamp use ocean_velocity_advect_mod, only: horz_advection_of_velocity, vert_advection_of_velocity use ocean_velocity_diag_mod, only: compute_vorticity, kinetic_energy, potential_energy use ocean_vert_mix_mod, only: vert_friction, vert_friction_implicit use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk1_v implicit none private ! for diagnostics integer :: id_u(2) =-1 integer :: id_u_on_depth(2) =-1 integer :: id_usurf(2) =-1 integer :: id_ubott(2) =-1 integer :: id_vorticity_z =-1 integer :: id_speed =-1 integer :: id_accel(2) =-1 integer :: id_converge_rho_ud_t =-1 logical :: used !for restart file integer :: id_restart_u = 0 integer :: id_restart_v = 0 integer :: id_restart_advu = 0 integer :: id_restart_advv = 0 type(restart_file_type), save :: Vel_restart type(restart_file_type), save :: Adv_restart integer :: unit=6 character(len=128) :: & version='$Id: ocean_velocity.F90,v 16.0.2.2.18.1.54.1 2009/10/10 00:42:04 nnz Exp $' character (len=128) :: tagname = & '$Name: mom4p1_pubrel_dec2009_nnz $' logical :: have_obc = .false. #include real :: dtime integer :: time_index integer :: tendency type(ocean_domain_type), pointer :: Dom =>NULL() type(ocean_grid_type), pointer :: Grd =>NULL() public ocean_velocity_init public ocean_velocity_end public ocean_explicit_accel_a public ocean_explicit_accel_b public ocean_implicit_accel public update_ocean_velocity public ocean_velocity_restart private ocean_velocity_chksum private check_gravity_wave_cfl private velocity_truncate private remap_s_to_depth ! Adams-Bashforth treatment of velocity advection logical :: adams_bashforth_third = .true. real :: adams_bashforth_epsilon = 0.6 real :: abtau_m0 = 23.0/12.0 real :: abtau_m1 = -16.0/12.0 real :: abtau_m2 = 5.0/12.0 logical :: update_velocity_via_uprime = .true. logical :: module_is_initialized = .false. logical :: zero_tendency = .false. logical :: zero_tendency_explicit_a = .false. logical :: zero_tendency_explicit_b = .false. logical :: zero_tendency_implicit = .false. logical :: truncate_velocity = .false. logical :: truncate_verbose = .false. real :: truncate_velocity_lat = 0.0 real :: truncate_velocity_value = 1.0 ! for idealized constant velocity situations logical :: use_constant_velocity = .false. real :: constant_u = 0.0 real :: constant_v = 0.0 ! for CFL check on gravity waves real :: max_cgint = 2.0 ! speed (m/s) used to compute CFL check for internal gravity waves real :: max_dt_for_cgint ! maximum time step (s) allowed by CFL to resolve internal gravity waves ! for debugging logical :: debug_this_module = .false. ! for writing restarts logical :: write_a_restart = .true. namelist /ocean_velocity_nml/ debug_this_module, write_a_restart, max_cgint, & zero_tendency, zero_tendency_explicit_a, zero_tendency_explicit_b, zero_tendency_implicit, & truncate_velocity, truncate_verbose, truncate_velocity_lat, truncate_velocity_value, & adams_bashforth_third, adams_bashforth_epsilon, & use_constant_velocity, constant_u, constant_v, & update_velocity_via_uprime contains !####################################################################### ! ! ! ! Initialize terms for the velocity equation. ! ! subroutine ocean_velocity_init (Grid, Domain, Time, Time_steps, Ocean_options, & Velocity, obc, debug) 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 type(ocean_options_type), intent(inout) :: Ocean_options type(ocean_velocity_type), intent(inout) :: Velocity logical, intent(in) :: obc logical, optional, intent(in) :: debug integer :: ioun, io_status, ierr integer :: m, n integer :: taum1, tau, taup1 integer :: tau_m0, tau_m1 integer, dimension(4) :: siz character(len=128) :: filename integer :: stdoutunit,stdlogunit stdoutunit=stdout();stdlogunit=stdlog() if ( module_is_initialized ) then call mpp_error(FATAL, & '==>Error from ocean_velocity_mod (ocean_velocity_init): module already initialized') endif module_is_initialized = .TRUE. call write_version_number( version, tagname ) have_obc = obc tendency = Time_steps%tendency dtime = Time_steps%dtime_u ! provide for namelist over-ride of defaults ioun = open_namelist_file() read (ioun, ocean_velocity_nml,iostat=io_status) write (stdoutunit,'(/)') write (stdoutunit, ocean_velocity_nml) write (stdlogunit, ocean_velocity_nml) ierr = check_nml_error(io_status, 'ocean_velocity_nml') call close_file(ioun) if (PRESENT(debug) .and. .not. debug_this_module) then debug_this_module = debug endif if(debug_this_module) then write(stdoutunit,'(a)') & '==>Note: running with debug_this_module=.true. so will print lots of checksums.' endif if(.not. write_a_restart) then write(stdoutunit,'(a)') '==>Note: running ocean_velocity with write_a_restart=.false.' write(stdoutunit,'(a)') ' Will NOT write restart file, and so cannot restart the run.' endif Dom => Domain Grd => Grid #ifndef MOM4_STATIC_ARRAYS call get_local_indices(Domain, isd, ied, jsd, jed, isc, iec, jsc, jec) nk=Grid%nk allocate (Velocity%u(isd:ied,jsd:jed,nk,2,3)) allocate (Velocity%smf(isd:ied,jsd:jed,2)) allocate (Velocity%bmf(isd:ied,jsd:jed,2)) allocate (Velocity%gamma(isd:ied,jsd:jed)) allocate (Velocity%rossby_radius(isd:ied,jsd:jed)) allocate (Velocity%press_force(isd:ied,jsd:jed,nk,2)) allocate (Velocity%accel(isd:ied,jsd:jed,nk,2)) allocate (Velocity%source(isd:ied,jsd:jed,nk,2)) allocate (Velocity%wrkv(isd:ied,jsd:jed,nk,2)) allocate (Velocity%advection(isd:ied,jsd:jed,nk,2,3)) allocate (Velocity%vfrict_impl(isd:ied,jsd:jed,nk,2)) allocate (Velocity%lap_friction_bt(isd:ied,jsd:jed,2)) allocate (Velocity%bih_friction_bt(isd:ied,jsd:jed,2)) #endif Velocity%smf = 0.0 Velocity%bmf = 0.0 Velocity%gamma = 0.0 Velocity%rossby_radius = 1.e5 Velocity%accel = 0.0 Velocity%source = 0.0 Velocity%wrkv = 0.0 Velocity%press_force = 0.0 Velocity%advection = 0.0 Velocity%vfrict_impl = 0.0 Velocity%lap_friction_bt = 0.0 Velocity%bih_friction_bt = 0.0 ! register fields for diagnostic output id_u(1) = register_diag_field ('ocean_model', 'u', Grd%vel_axes_uv(1:3), Time%model_time, & 'i-current', 'm/sec', missing_value=missing_value, range=(/-10.0,10.0/), & standard_name='sea_water_x_velocity') id_u(2) = register_diag_field ('ocean_model', 'v', Grd%vel_axes_uv(1:3), Time%model_time, & 'j-current', 'm/sec', missing_value=missing_value, range=(/-10.0,10.0/), & standard_name='sea_water_y_velocity') id_u_on_depth(1) = register_diag_field ('ocean_model', 'u_on_depth', Grd%vel_axes_uv_depth(1:3), Time%model_time, & 'i-current mapped to depth surface', 'm/sec', missing_value=missing_value, range=(/-10.0,10.0/)) id_u_on_depth(2) = register_diag_field ('ocean_model', 'v_on_depth', Grd%vel_axes_uv_depth(1:3), Time%model_time, & 'j-current mapped to depth surface', 'm/sec', missing_value=missing_value, range=(/-10.0,10.0/)) id_usurf(1) = register_diag_field ('ocean_model', 'usurf', Grd%vel_axes_uv(1:2), Time%model_time, & 'i-surface current', 'm/sec', missing_value=missing_value, range=(/-10.0,10.0/)) id_usurf(2) = register_diag_field ('ocean_model', 'vsurf', Grd%vel_axes_uv(1:2), Time%model_time, & 'j-surface current', 'm/sec', missing_value=missing_value, range=(/-10.0,10.0/)) id_ubott(1) = register_diag_field ('ocean_model', 'ubott', Grd%vel_axes_uv(1:2), Time%model_time, & 'i-bottom current', 'm/sec', missing_value=missing_value, range=(/-10.0,10.0/)) id_ubott(2) = register_diag_field ('ocean_model', 'vbott', Grd%vel_axes_uv(1:2), Time%model_time, & 'j-bottom current', 'm/sec', missing_value=missing_value, range=(/-10.0,10.0/)) id_vorticity_z = register_diag_field ('ocean_model', 'vorticity_z', Grd%tracer_axes(1:3), Time%model_time, & 'vertical vorticity component: v_x-u_y', '1/sec', missing_value=missing_value, range=(/-1e6,1e6/)) id_speed = register_diag_field ('ocean_model', 'speed', Grd%vel_axes_uv(1:3), Time%model_time, & 'speed of horizontal current', 'm/sec', missing_value=missing_value, range=(/-10.0,10.0/)) id_accel(1) = register_diag_field ('ocean_model', 'accel_u', Grd%vel_axes_uv(1:3), Time%model_time, & 'baroclinic forcing of rho*u', '(kg/m^3)*(m^2/s^2)', missing_value=missing_value, range=(/-1e9,1e9/)) id_accel(2) = register_diag_field ('ocean_model', 'accel_v', Grd%vel_axes_uv(1:3), Time%model_time, & 'baroclinic forcing of rho*v', '(kg/m^3)*(m^2/s^2)', missing_value=missing_value, range=(/-1e9,1e9/)) id_converge_rho_ud_t = register_diag_field ('ocean_model', 'converge_rho_ud_t', Grd%tracer_axes(1:2), & Time%model_time,'convergence rho*ud on T cells computed in velocity_mod', & '(kg/m^3)*(m/s)', missing_value=missing_value, range=(/-10.0,10.0/)) if (truncate_velocity) then write (stdoutunit,'(/a)')& '==>Warning: truncate_velocity=.true. Model will truncate baroclinc velocity' write (stdoutunit,'(a,f5.3)')& ' so that max horz velocity magnitude is (m/s) ',truncate_velocity_value write (stdoutunit,'(a,f8.3)') & ' Truncation occurs for regions poleward of the latitude',truncate_velocity_lat write (stdoutunit,'(a)')& ' This option may be useful during the spin-up phase of an experiment.' write (stdoutunit,'(a/)')& ' It is also of use for polar filtering. ' endif if(tendency==THREE_LEVEL) then abtau_m0= 1.0 abtau_m1= 0.0 abtau_m2= 0.0 endif if(tendency==TWO_LEVEL) then write(stdoutunit,'(/a)')'==>Note from ocean_velocity_mod: use of twolevel time_tendency' write(stdoutunit,'(a)') ' necessitates an Adams-Bashforth treatment of velocity advection.' if(adams_bashforth_third) then write(stdoutunit,'(/a)')' Using 3rd order Adams-Bashforth for velocity advection.' write(stdoutunit,'(a/)')' This is the mom4 default. ' abtau_m0= 23.0/12.0 abtau_m1=-16.0/12.0 abtau_m2= 5.0/12.0 else write(stdoutunit,'(/a,f8.3/)')' Using 2nd order Adams-Bashforth with dimensionless AB parameter = ', & adams_bashforth_epsilon abtau_m0= 1.0 + adams_bashforth_epsilon abtau_m1= 0.0 - adams_bashforth_epsilon abtau_m2= 0.0 if(adams_bashforth_epsilon < 0.5 .or. adams_bashforth_epsilon > 1.0) then call mpp_error(FATAL,& '==>Error from ocean_velocity_mod: adams_bashforth parameter must be between 0.5 & 1.0') endif endif endif if(zero_tendency) then call mpp_error(NOTE, & '==>zero_tendency=.true. so will not time step the baroclinic velocity field ') Ocean_options%baroclinic_tendency = 'Did NOT time step the baroclinic velocity.' else Ocean_options%baroclinic_tendency = 'Time stepped the baroclinic velocity.' endif if(zero_tendency_explicit_a) then call mpp_error(NOTE, & '==>zero_tendency_explicit_a=.true. so will not include explicit_a velocity tendency ') endif if(zero_tendency_explicit_b) then call mpp_error(NOTE, & '==>zero_tendency_explicit_b=.true. so will not include explicit_b velocity tendency ') endif if(zero_tendency_implicit) then call mpp_error(NOTE, & '==>zero_tendency_implicit=.true. so will not include implicit part of the velocity tendency ') endif if(update_velocity_via_uprime) then call mpp_error(NOTE, & '==>update_velocity_via_uprime=.true., so keep udrho from external mode solver.') else call mpp_error(WARNING, & '==>update_velocity_via_uprime=.false., so recompute udrho after updating full velocity. For testing only!') endif if(use_constant_velocity) then call mpp_error(NOTE, & '==>use_constant_velocity=.true. so will hold the velocity to a constant.') zero_tendency=.true. Ocean_options%baroclinic_tendency = 'Did NOT time step velocity; kept it = constant.' write(stdoutunit,'(a,f6.2,a,f6.2,a)') & '==>Note from ocean_velocity_init: velocity fixed at (u,v)(m/s) = (',constant_u,',',constant_v,')' endif call check_gravity_wave_cfl() taum1 = Time%taum1 tau = Time%tau taup1 = Time%taup1 tau_m1 = Time%tau_m1 tau_m0 = Time%tau_m0 ! initialize to small but nonzero velocity do m=1,3 do n=1,2 Velocity%u(:,:,:,n,m) = epsln*Grd%umask(:,:,:) enddo enddo if(use_constant_velocity) then do m=1,3 Velocity%u(:,:,:,1,m) = constant_u*Grd%umask(:,:,:) Velocity%u(:,:,:,2,m) = constant_v*Grd%umask(:,:,:) enddo endif filename = 'ocean_velocity.res.nc' if(tendency==THREE_LEVEL) then id_restart_u = register_restart_field(Vel_restart, filename,'u',Velocity%u(:,:,:,1,tau), & Velocity%u(:,:,:,1,taup1), Dom%domain2d) id_restart_v = register_restart_field(Vel_restart, filename,'v',Velocity%u(:,:,:,2,tau), & Velocity%u(:,:,:,2,taup1), Dom%domain2d ) elseif (tendency==TWO_LEVEL) then id_restart_u = register_restart_field(Vel_restart, filename,'u',Velocity%u(:,:,:,1,taup1), & Dom%domain2d ) id_restart_v = register_restart_field(Vel_restart, filename,'v',Velocity%u(:,:,:,2,taup1), & Dom%domain2d ) filename = 'ocean_velocity_advection.res.nc' id_restart_advu = register_restart_field(Adv_restart, filename,'advectionu',Velocity%advection(:,:,:,1,tau_m1),& Velocity%advection(:,:,:,1,tau_m0), Dom%domain2d ) id_restart_advv = register_restart_field(Adv_restart, filename,'advectionv',Velocity%advection(:,:,:,2,tau_m1),& Velocity%advection(:,:,:,2,tau_m0), Dom%domain2d ) end if filename = 'INPUT/ocean_velocity.res.nc' if (file_exist(trim(filename))) then !--- read restart file data by calling restore_state call restore_state(Vel_restart) if(tendency==THREE_LEVEL) then write (stdoutunit,'(a)') ' Reading THREE_LEVEL restart for velocity from INPUT/ocean_velocity.res.nc' write (stdoutunit,'(a)') ' Expecting two time records for each restart field.' call mpp_update_domains(Velocity%u(:,:,:,1,tau), Velocity%u(:,:,:,2,tau),& Dom%domain2d,gridtype=BGRID_NE) call mpp_update_domains(Velocity%u(:,:,:,1,taup1), Velocity%u(:,:,:,2,taup1),& Dom%domain2d,gridtype=BGRID_NE) if(have_obc) then call ocean_obc_update_boundary(Velocity%u(:,:,:,1,taup1), 'M','n') call ocean_obc_update_boundary(Velocity%u(:,:,:,2,taup1), 'M','t') call ocean_obc_update_boundary(Velocity%u(:,:,:,1,taup1), 'Z','t') call ocean_obc_update_boundary(Velocity%u(:,:,:,2,taup1), 'Z','n') call ocean_obc_update_boundary(Velocity%u(:,:,:,1,tau), 'M','n') call ocean_obc_update_boundary(Velocity%u(:,:,:,2,tau), 'M','t') call ocean_obc_update_boundary(Velocity%u(:,:,:,1,tau), 'Z','t') call ocean_obc_update_boundary(Velocity%u(:,:,:,2,tau), 'Z','n') endif elseif (tendency==TWO_LEVEL) then write (stdoutunit,'(/a)') ' Reading TWO_LEVEL restart for velocity from INPUT/ocean_velocity.res.nc' write (stdoutunit,'(a)') ' Expecting only one time record for each restart field.' call field_size(filename,'u', siz) if (siz(4) > 1) then write(stdoutunit,'(/a)') & '==>ERROR: Attempt to read ocean_velocity.res.nc from a 3-level time scheme (2 time records)' write(stdoutunit,'(a)') & ' when running mom4 with 2-level timestepping (only need 1 time record in restart).' write(stdoutunit,'(a)') & ' Reduce restart file to a single time record in order to avoid confusion.' call mpp_error(FATAL, & 'Reading 3-time lev ocean_velocity.res.nc (w/ 2 time records) while using 2-lev (needs only 1 record)') endif call mpp_update_domains(Velocity%u(:,:,:,1,taup1), Velocity%u(:,:,:,2,taup1),& Dom%domain2d,gridtype=BGRID_NE) write (stdoutunit,'(a)') ' Finished reading restart for velocity field' if(.not. Time%init) then write (stdoutunit,'(/a)') ' Reading TWO_LEVEL restart from INPUT/ocean_velocity_advection.res.nc' write (stdoutunit,'(a)') ' Expecting two time records for the advection velocity fields.' call restore_state(Adv_restart) write (stdoutunit,'(a)') ' Finished reading restart for velocity advection operator' call mpp_update_domains(Velocity%advection(:,:,:,1,tau_m1),Dom%domain2d) call mpp_update_domains(Velocity%advection(:,:,:,1,tau_m0),Dom%domain2d) call mpp_update_domains(Velocity%advection(:,:,:,2,tau_m1),Dom%domain2d) call mpp_update_domains(Velocity%advection(:,:,:,2,tau_m0),Dom%domain2d) endif if(have_obc) then call ocean_obc_update_boundary(Velocity%u(:,:,:,1,taup1), 'M','n') call ocean_obc_update_boundary(Velocity%u(:,:,:,2,taup1), 'M','t') call ocean_obc_update_boundary(Velocity%u(:,:,:,1,taup1), 'Z','t') call ocean_obc_update_boundary(Velocity%u(:,:,:,2,taup1), 'Z','n') endif endif else if (Time%init) then call mpp_error(NOTE,& '==> From ocean_velocity_mod: Initializing velocity to zero since Time%init=.true. & &and did not find INPUT/ocean_velocity.res.nc.') endif if (.NOT. Time%init) then call mpp_error(FATAL,& '==> Error from ocean_velocity_mod: Expecting INPUT/ocean_velocity.res.nc to exist.& &This file was not found and Time%init=.false.') endif endif write(stdoutunit,*)' ' write(stdoutunit,*) '===Initial velocity checksums ==>' call write_timestamp(Time%model_time) if(tendency==THREE_LEVEL) then write(stdoutunit,*) 'From ocean_velocity_mod: initial velocity chksum (tau)' call ocean_velocity_chksum(Velocity, tau, write_advection=.false.) write(stdoutunit,*) 'From ocean_velocity_mod: initial velocity chksum (taup1)' call ocean_velocity_chksum(Velocity, taup1, write_advection=.false.) else write(stdoutunit,*) 'From ocean_velocity_mod: initial velocity chksum (taup1)' call ocean_velocity_chksum(Velocity, taup1, write_advection=.true.) endif end subroutine ocean_velocity_init ! NAME="ocean_velocity_init" !####################################################################### ! ! ! ! Check CFL for internal gravity waves. ! ! subroutine check_gravity_wave_cfl() logical :: cfl_error=.false. real :: gridsp, dtcg, max_dt_for_cgint0 real :: cfl_grid_factor, dtuv integer :: i, j, icg, jcg integer :: stdlogunit stdlogunit=stdlog() ! estimate the maximum timestep allowable for resolving internal gravity waves if(tendency==TWO_LEVEL) then cfl_grid_factor=1.0 dtuv = dtime elseif(tendency==THREE_LEVEL) then cfl_grid_factor=0.5 dtuv = 0.5*dtime endif icg = isc; jcg = jsc; gridsp = 1.0e20; max_dt_for_cgint = 1.0e6 do j=jsc,jec do i=isc,iec if (Grd%kmu(i,j) > 0) then gridsp = 1.0/(Grd%dxu(i,j)*Grd%dxu(i,j)) + 1.0/(Grd%dyu(i,j)*Grd%dyu(i,j)) gridsp = sqrt(1.0/gridsp) dtcg = cfl_grid_factor*gridsp/(epsln+max_cgint) if (dtcg < max_dt_for_cgint) then max_dt_for_cgint = dtcg; icg = i; jcg = j endif endif enddo enddo max_dt_for_cgint = nint(max_dt_for_cgint) max_dt_for_cgint = max_dt_for_cgint + 0.001*mpp_pe() ! to separate redundancies max_dt_for_cgint0 = max_dt_for_cgint call mpp_min (max_dt_for_cgint) ! show the most unstable location for baroclinic gravity waves if (max_dt_for_cgint == max_dt_for_cgint0 .and. nint(dtuv) > 0) then if (dtuv <= max_dt_for_cgint) then write (unit,'(/a,i4,a,i4,a,f6.2,a,f6.2,a)') & ' Baroclinic time step stability most nearly violated at U-cell (i,j) = (',& icg+Dom%ioff,',',jcg+Dom%joff,'), (lon,lat) = (',Grd%xu(icg,jcg),',',Grd%yu(icg,jcg),').' write(unit,'(a,i6)') ' The number of kmu-levels at this point is ',Grd%kmu(icg,jcg) write(unit,'(a,e12.6)') ' The dxu grid distance (m) at this point is ',Grd%dxu(icg,jcg) write(unit,'(a,e12.6)') ' The dyu grid distance (m) at this point is ',Grd%dyu(icg,jcg) cfl_error=.false. else write (unit,'(/a,i4,a,i4,a,f6.2,a,f6.2,a)') & '==>Error: Baroclinic time step stability violated at U-cell (i,j) = (',& icg+Dom%ioff,',',jcg+Dom%joff,'), (lon,lat) = (',Grd%xu(icg,jcg),',',Grd%yu(icg,jcg),').' write(unit,'(a,i6)') ' The number of kmu-levels at this point is ',Grd%kmu(icg,jcg) write(unit,'(a,e12.6)') ' The dxu grid distance (m) at this point is ',Grd%dxu(icg,jcg) write(unit,'(a,e12.6)') ' The dyu grid distance (m) at this point is ',Grd%dyu(icg,jcg) cfl_error=.true. endif write (unit,'(a,f5.2,a/a,f6.0,a,f6.0,a)')& ' Due to a specified maximum baroclinic gravity wave speed of ',max_cgint,' m/s.',& ' "dtuv" must be less than ',max_dt_for_cgint,' sec. "dtuv" = ', dtuv,' sec.' if(cfl_error) then call mpp_error(FATAL, & '==>Error: time step instability detected for baroclinic gravity waves in ocean_model_mod') endif write (stdlogunit,'(/a/)') & ' Note: A more appropriate maximum baroclinic gravity wave speed can be specified via namelist.' endif max_dt_for_cgint = nint(max_dt_for_cgint) end subroutine check_gravity_wave_cfl ! NAME="check_gravity_wave_cfl" !####################################################################### ! ! ! ! ! Time explicit contributions to thickness weighted and density ! weighted acceleration. ! ! Omit here the Coriolis force and vertical friction. ! They are omitted in order to facilitate the construction of the ! vertically integrated forcing of the barotropic dynamics. ! They will be added later in ocean_explicit_accel_b. ! ! ! subroutine ocean_explicit_accel_a(Velocity, Time, Adv_vel, Thickness, Dens, & rho, pme, river, upme, uriver) type(ocean_velocity_type), intent(inout) :: Velocity type(ocean_time_type), intent(in) :: Time type(ocean_adv_vel_type), intent(in) :: Adv_vel type(ocean_thickness_type), intent(in) :: Thickness type(ocean_density_type), intent(in) :: Dens real, dimension(isd:,jsd:,:), intent(in) :: rho 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 integer :: i, j, k, n integer :: taum1, tau, taup1 integer :: tau_m0, tau_m1, tau_m2 integer :: stdoutunit stdoutunit=stdout() taum1 = Time%taum1 tau = Time%tau taup1 = Time%taup1 tau_m2 = Time%tau_m2 tau_m1 = Time%tau_m1 tau_m0 = Time%tau_m0 ! initialize to zero do n=1,2 do k=1,nk do j=jsd,jed do i=isd,ied Velocity%accel(i,j,k,n) = 0.0 Velocity%source(i,j,k,n) = 0.0 Velocity%press_force(i,j,k,n) = 0.0 Velocity%advection(i,j,k,n,tau_m0) = 0.0 enddo enddo enddo enddo if (.not. zero_tendency .and. .not. zero_tendency_explicit_a) then ! fill Velocity%advection(tau_m0) call horz_advection_of_velocity(Time, Thickness, Adv_vel, Velocity, energy_analysis_step=.false.) call vert_advection_of_velocity(Time, Adv_vel, Velocity, & pme, river, upme, uriver, energy_analysis_step=.false.) do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec Velocity%accel(i,j,k,n) = -abtau_m0*Velocity%advection(i,j,k,n,tau_m0) & -abtau_m1*Velocity%advection(i,j,k,n,tau_m1) & -abtau_m2*Velocity%advection(i,j,k,n,tau_m2) enddo enddo enddo enddo ! include contribution from horizontal pressure gradient call pressure_force(Time, Thickness, Dens, Velocity, rho(:,:,:)) do n=1,2 do k=1,nk do j=jsc,jec do i=isc,iec Velocity%accel(i,j,k,n) = Velocity%accel(i,j,k,n) + Velocity%press_force(i,j,k,n) enddo enddo enddo enddo ! compute horizontal friction and add to Velocity%accel call lap_friction(Time, Thickness, Adv_vel, Velocity, energy_analysis_step=.false.) call bih_friction(Time, Thickness, Velocity, energy_analysis_step=.false.) ! compute momentum source/sinks and add to Velocity%accel call momentum_source(Time, Thickness, Velocity, energy_analysis_step=.false.) ! compute parameterized form drag and add to Velocity%accel call form_drag_accel(Time, Thickness, Velocity, energy_analysis_step=.false.) endif if (debug_this_module) then write(stdoutunit,*) ' ' write(stdoutunit,*) 'From ocean_velocity_mod: acceleration chksums' call write_timestamp(Time%model_time) write(stdoutunit,*) 'explicit accel_a(1) = ',mpp_chksum(Velocity%accel(isc:iec,jsc:jec,:,1)) write(stdoutunit,*) 'explicit accel_a(2) = ',mpp_chksum(Velocity%accel(isc:iec,jsc:jec,:,2)) endif end subroutine ocean_explicit_accel_a ! NAME="ocean_explicit_accel_a" !####################################################################### ! ! ! ! ! Add the contributions from the Coriolis force, computed explicitly ! in time, and those from explicit vertical friction. Add these to ! the thickness weighted and density weighted acceleration. ! ! Note: no visc_cbu_form_drag is included here, since it must ! be handled via implicit vertical friction to maintain stability ! for general cases. ! ! ! subroutine ocean_explicit_accel_b(visc_cbu, Time, Thickness, Velocity) real, dimension(isd:,jsd:,:), intent(in) :: visc_cbu type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_velocity_type), intent(inout) :: Velocity integer :: stdoutunit stdoutunit=stdout() if (.not. zero_tendency .and. .not. zero_tendency_explicit_b) then call vert_friction(Time, Thickness, Velocity, visc_cbu, energy_analysis_step=.false.) call coriolis_force(Time, Thickness, Velocity, energy_analysis_step=.false.) if (debug_this_module) then write(stdoutunit,*) ' ' write(stdoutunit,*) 'From ocean_velocity_mod: explicit acceleration chksums' call write_timestamp(Time%model_time) write(stdoutunit,*) 'accel_b(1) = ',mpp_chksum(Velocity%accel(isc:iec,jsc:jec,:,1)) write(stdoutunit,*) 'accel_b(2) = ',mpp_chksum(Velocity%accel(isc:iec,jsc:jec,:,2)) endif endif end subroutine ocean_explicit_accel_b ! NAME="ocean_explicit_accel_b" !####################################################################### ! ! ! ! ! Add the time implicit contributions from the Coriolis force ! and vertical friction. Add these to the thickness weighted ! and density weighted acceleration. ! ! Note the contribution from visc_cbu_form_drag is used for ! implicit vertical friction. ! ! ! subroutine ocean_implicit_accel(visc_cbu, visc_cbu_form_drag, Time, Thickness, Velocity) real, dimension(isd:,jsd:,:), intent(in) :: visc_cbu real, dimension(isd:,jsd:,:,:), intent(in) :: visc_cbu_form_drag type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_velocity_type), intent(inout) :: Velocity integer :: stdoutunit stdoutunit=stdout() if (.not. zero_tendency .and. .not. zero_tendency_implicit) then call vert_friction_implicit(visc_cbu, visc_cbu_form_drag, Time, Thickness, Velocity) call coriolis_force_implicit(Time, Velocity) if (debug_this_module) then write(stdoutunit,*) ' ' write(stdoutunit,*) 'From ocean_velocity_mod: acceleration chksums after implicit update' call write_timestamp(Time%model_time) write(stdoutunit,*) 'accel(1) = ',mpp_chksum(Velocity%accel(isc:iec,jsc:jec,:,1)) write(stdoutunit,*) 'accel(2) = ',mpp_chksum(Velocity%accel(isc:iec,jsc:jec,:,2)) endif endif ! since implicit_accel is the last piece of accel computed, it is ! time to now output the full acceleration to diagnostics manager. if (id_accel(1) > 0) used = send_data(id_accel(1), Velocity%accel(:,:,:,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_accel(2) > 0) used = send_data(id_accel(2), Velocity%accel(:,:,:,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) end subroutine ocean_implicit_accel ! NAME="ocean_implicit_accel" !####################################################################### ! ! ! ! Update velocity components. ! There are two general methods available. ! ! 1/ update baroclinic velocity; then add (udrho,vdrho) from external ! mode algorithm to get the full velocity field. This method is ! analogous to the older approach with the rigid lid. ! ! 2/ update the full velocity, so there is no need to introduce the ! intermediate step with the baroclinic velocity. To remain stable, ! we must use the time filtered barotropic pressure gradient. ! We then diagnose the vertically integrated horizontal momentum ! (udrho,vdrho) and its convergence, since these fields are needed ! elsewhere. ! ! ! subroutine update_ocean_velocity(Time, Thickness, barotropic_split, & vert_coordinate_class, Ext_mode, Velocity) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness integer, intent(in) :: barotropic_split integer, intent(in) :: vert_coordinate_class type(ocean_external_mode_type), intent(inout) :: Ext_mode type(ocean_velocity_type), intent(inout) :: Velocity real, dimension(isd:ied,jsd:jed) :: tmp real, dimension(isd:ied,jsd:jed) :: tmpu real, dimension(isd:ied,jsd:jed) :: tmpv integer :: i, j, k, kbot, n integer :: taum1, tau, taup1 integer :: stdoutunit stdoutunit=stdout() taum1 = Time%taum1 tau = Time%tau taup1 = Time%taup1 tmp = 0.0 if (barotropic_split > 1 .and. update_velocity_via_uprime) then do n=1,2 ! compute uprime do k=1,nk do j=jsc,jec do i=isc,iec Velocity%u(i,j,k,n,taup1) = Thickness%rho_dzur(i,j,k) & *(Thickness%rho_dzu(i,j,k,taum1)*Velocity%u(i,j,k,n,taum1) & +dtime*Velocity%accel(i,j,k,n)*Grd%umask(i,j,k) ) enddo enddo enddo ! sum uprime*rho_dzu tmp=0.0 do k=1,nk do j=jsc,jec do i=isc,iec tmp(i,j) = tmp(i,j) & +Velocity%u(i,j,k,n,taup1)*Thickness%rho_dzu(i,j,k,taup1)*Grd%umask(i,j,k) enddo enddo enddo ! (vertical mean) - (value from barotropic integration) do j=jsc,jec do i=isc,iec tmp(i,j) = Grd%umask(i,j,1)*(tmp(i,j)-Ext_mode%udrho(i,j,n,taup1)) & /(Thickness%mass_u(i,j,taup1)+epsln) enddo enddo ! replace vertical mean with value from barotropic integration to update velocity do k=1,nk do j=jsc,jec do i=isc,iec Velocity%u(i,j,k,n,taup1) = (Velocity%u(i,j,k,n,taup1)-tmp(i,j))*Grd%umask(i,j,k) enddo enddo enddo enddo ! finish the n=1,2 do-loop ! update full velocity directly with the full pressure gradient. ! barotropic contribution to pressure gradient is time averaged ! if barotropic_split>1. Unless dtuv is very small, this approach ! can be rather unstable. User beware! else do n=1,2 if(vert_coordinate_class==DEPTH_BASED) then ! include gradient of surface pressure do k=1,nk do j=jsc,jec do i=isc,iec Velocity%u(i,j,k,n,taup1) = & (Thickness%rho_dzu(i,j,k,taum1)*Velocity%u(i,j,k,n,taum1) & +dtime*(Velocity%accel(i,j,k,n) & -rho0r*Thickness%rho_dzu(i,j,k,tau)*Ext_mode%grad_ps(i,j,n)*Grd%umask(i,j,k))) & *Thickness%rho_dzur(i,j,k) enddo enddo enddo else ! include gradient of bottom pressure do k=1,nk do j=jsc,jec do i=isc,iec Velocity%u(i,j,k,n,taup1) = & (Thickness%rho_dzu(i,j,k,taum1)*Velocity%u(i,j,k,n,taum1) & +dtime*(Velocity%accel(i,j,k,n) & -rho0r*Thickness%rho_dzu(i,j,k,tau)*Ext_mode%grad_anompb(i,j,n)*Grd%umask(i,j,k))) & *Thickness%rho_dzur(i,j,k) enddo enddo enddo endif ! diagnose vertically integrated horizontal momentum Ext_mode%udrho(:,:,n,taup1) = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec Ext_mode%udrho(i,j,n,taup1) = Ext_mode%udrho(i,j,n,taup1) & +Velocity%u(i,j,k,n,taup1)*Thickness%rho_dzu(i,j,k,taup1) enddo enddo enddo enddo ! end for n=1,2 ! fill halos for vertically integrated horizontal momentum call mpp_update_domains (Ext_mode%udrho(:,:,1,taup1),Ext_mode%udrho(:,:,2,taup1), & Dom%domain2d, gridtype=BGRID_NE) ! construct convergence of udrho at (tau+1) Ext_mode%conv_rho_ud_t(:,:,taup1) = -DIV_UD(Ext_mode%udrho(:,:,:,taup1)) call mpp_update_domains (Ext_mode%conv_rho_ud_t(:,:,taup1), Dom%domain2d) endif ! endif for splitting or uprime method ! override above update of velocity in the case of zero tendency if(zero_tendency) then Velocity%u(isc:iec,jsc:jec,:,:,taup1) = Velocity%u(isc:iec,jsc:jec,:,:,tau) endif ! for debugging do n=1,2 if(truncate_velocity) then call velocity_truncate(n, taup1, Velocity) endif if (debug_this_module .and. n==2) then write(stdoutunit,*) ' ' write(stdoutunit,*) 'From update_ocean_velocity at the u-prime step' call write_timestamp(Time%model_time) write(stdoutunit,*) 'accel(1) = ',mpp_chksum(Velocity%accel(isc:iec,jsc:jec,:,1)) write(stdoutunit,*) 'accel(2) = ',mpp_chksum(Velocity%accel(isc:iec,jsc:jec,:,2)) write(stdoutunit,*) 'rho_dzu(taup1) = ',mpp_chksum(Thickness%rho_dzu(isc:iec,jsc:jec,:,taup1)) write(stdoutunit,*) 'rho_dzu(taum1) = ',mpp_chksum(Thickness%rho_dzu(isc:iec,jsc:jec,:,taum1)) write(stdoutunit,*) 'rho_dzur = ',mpp_chksum(Thickness%rho_dzur(isc:iec,jsc:jec,:)) write(stdoutunit,*) 'udrho(1) = ',mpp_chksum(Ext_mode%udrho(isc:iec,jsc:jec,1,taup1)) write(stdoutunit,*) 'udrho(2) = ',mpp_chksum(Ext_mode%udrho(isc:iec,jsc:jec,2,taup1)) write(stdoutunit,*) 'mass_u = ',mpp_chksum(Thickness%mass_u(isc:iec,jsc:jec,taup1)) call ocean_velocity_chksum(Velocity, taum1, write_advection=.false.) call ocean_velocity_chksum(Velocity, taup1, write_advection=.false.) endif enddo !--- update velocity at the global halo points to make the gradient 0 across boundary if(have_obc) then call ocean_obc_update_boundary(Velocity%u(:,:,:,1,taup1), 'M','n') call ocean_obc_update_boundary(Velocity%u(:,:,:,2,taup1), 'M','t') call ocean_obc_update_boundary(Velocity%u(:,:,:,1,taup1), 'Z','t') call ocean_obc_update_boundary(Velocity%u(:,:,:,2,taup1), 'Z','n') endif if (debug_this_module) then write(stdoutunit,*) ' ' write(stdoutunit,*) 'From update_ocean_velocity: velocity(taup1) chksum' call write_timestamp(Time%model_time) call ocean_velocity_chksum(Velocity, taup1, write_advection=.false.) endif ! send diagnostics to diagnostics manager if (id_u(1) > 0) used = send_data (id_u(1), Velocity%u(:,:,:,1,tau), & 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(2) > 0) used = send_data (id_u(2), Velocity%u(:,:,:,2,tau), & 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_usurf(1) > 0) used = send_data (id_usurf(1), Velocity%u(:,:,1,1,tau), & Time%model_time, rmask=Grd%umask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) if (id_usurf(2) > 0) used = send_data (id_usurf(2), Velocity%u(:,:,1,2,tau), & Time%model_time, rmask=Grd%umask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) if (id_converge_rho_ud_t > 0) used = send_data (id_converge_rho_ud_t, Ext_mode%conv_rho_ud_t(:,:,tau), & Time%model_time, rmask=Grd%tmask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) if(id_speed > 0) then do k=1,nk do j=jsd,jed do i=isd,ied wrk1(i,j,k) = sqrt(Velocity%u(i,j,k,1,tau)**2 + Velocity%u(i,j,k,2,tau)**2) enddo enddo enddo used = send_data (id_speed, 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 if(id_ubott(1) > 0 .or. id_ubott(2) > 0) then tmpu = 0.0 ; tmpv = 0.0 do j=jsd,jed do i=isd,ied kbot = Grd%kmu(i,j) if(kbot > 1) then tmpu(i,j) = Velocity%u(i,j,kbot,1,tau) tmpv(i,j) = Velocity%u(i,j,kbot,2,tau) endif enddo enddo if (id_ubott(1) > 0) used = send_data (id_ubott(1), tmpu(:,:), & Time%model_time, rmask=Grd%umask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) if (id_ubott(2) > 0) used = send_data (id_ubott(2), tmpv(:,:), & Time%model_time, rmask=Grd%umask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) endif if(id_vorticity_z > 0) then wrk1=0.0 call compute_vorticity(Time, Velocity, wrk1) used = send_data (id_vorticity_z, wrk1(:,:,:), & Time%model_time, rmask=Grd%tmask(:,:,:), & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) endif if(id_u_on_depth(1) > 0) then call remap_s_to_depth(Thickness, Time, Velocity%u(:,:,:,1,tau), 1) endif if(id_u_on_depth(2) > 0) then call remap_s_to_depth(Thickness, Time, Velocity%u(:,:,:,2,tau), 2) endif end subroutine update_ocean_velocity ! NAME="update_ocean_velocity" !####################################################################### ! ! ! Write out restart files registered through register_restart_file ! subroutine ocean_velocity_restart(Time, Velocity, time_stamp) type(ocean_time_type), intent(in) :: Time type(ocean_velocity_type), intent(in) :: Velocity character(len=*), intent(in), optional :: time_stamp integer :: tau, taup1, tau_m1, tau_m0 tau = Time%tau taup1 = Time%taup1 tau_m1 = Time%tau_m1 tau_m0 = Time%tau_m0 ! resetting the pointer to data if(tendency == THREE_LEVEL) then call reset_field_pointer(Vel_restart, id_restart_u, Velocity%u(:,:,:,1,tau), Velocity%u(:,:,:,1,taup1) ) call reset_field_pointer(Vel_restart, id_restart_v, Velocity%u(:,:,:,2,tau), Velocity%u(:,:,:,2,taup1) ) elseif (tendency==TWO_LEVEL) then call reset_field_pointer(Vel_restart, id_restart_u, Velocity%u(:,:,:,1,taup1) ) call reset_field_pointer(Vel_restart, id_restart_v, Velocity%u(:,:,:,2,taup1) ) end if call save_restart(Vel_restart, time_stamp) if(tendency==TWO_LEVEL) then call reset_field_pointer(Adv_restart, id_restart_advu, Velocity%advection(:,:,:,1,tau_m1), & Velocity%advection(:,:,:,1,tau_m0) ) call reset_field_pointer(Adv_restart, id_restart_advv, Velocity%advection(:,:,:,2,tau_m1), & Velocity%advection(:,:,:,2,tau_m0) ) call save_restart(Adv_restart, time_stamp) end if return end subroutine ocean_velocity_restart ! NAME="ocean_velocity_restart" !####################################################################### ! ! ! ! Write the velocity field to a restart ! subroutine ocean_velocity_end(Time, Velocity) type(ocean_time_type), intent(in) :: Time type(ocean_velocity_type), intent(in) :: Velocity integer :: tau, taup1 integer :: tau_m0, tau_m1 integer :: stdoutunit stdoutunit=stdout() tau = Time%tau taup1 = Time%taup1 tau_m1 = Time%tau_m1 tau_m0 = Time%tau_m0 call ocean_velocity_restart(Time, Velocity) if(tendency==THREE_LEVEL) then write(stdoutunit,*)' ' write(stdoutunit,*) 'From ocean_velocity_mod: ending velocity chksum (tau)' call write_timestamp(Time%model_time) call ocean_velocity_chksum(Velocity, tau, write_advection=.false.) write(stdoutunit,*)' ' call write_timestamp(Time%model_time) write(stdoutunit,*) 'From ocean_velocity_mod: ending velocity chksum (taup1)' call ocean_velocity_chksum(Velocity, taup1, write_advection=.false.) elseif(tendency==TWO_LEVEL) then write(stdoutunit,*)' ' write(stdoutunit,*) 'From ocean_velocity_mod: ending velocity chksum (taup1)' call write_timestamp(Time%model_time) call ocean_velocity_chksum(Velocity, taup1, write_advection=.true.) endif return end subroutine ocean_velocity_end ! NAME="ocean_velocity_end" !####################################################################### ! ! ! ! Truncate velocity so that either component ! has magnitude no larger than nml specified value. ! ! subroutine velocity_truncate(n, taup1, Velocity) integer, intent(in) :: n integer, intent(in) :: taup1 type(ocean_velocity_type), intent(inout) :: Velocity logical :: velocity_truncated integer :: i,j,k integer :: stdoutunit stdoutunit=stdout() velocity_truncated = .false. wrk1 = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec if ( abs(Grd%yu(i,j)) >= truncate_velocity_lat & .and. abs(Velocity%u(i,j,k,n,taup1)) > truncate_velocity_value) then wrk1(i,j,k) = Velocity%u(i,j,k,n,taup1) Velocity%u(i,j,k,n,taup1) = sign(truncate_velocity_value,Velocity%u(i,j,k,n,taup1)) velocity_truncated = .true. endif enddo enddo enddo if(truncate_verbose .and. velocity_truncated) then write(stdoutunit,'(/a)') ' Summary of baroclinic velocity truncation points' do k=1,nk do j=jsc,jec do i=isc,iec if (abs(wrk1(i,j,k)) > 0.0) then write(stdoutunit,'(a,i1,a,i3,a,i3,a,i3,a,e12.3,a,f12.3,a,f12.3,a,i4)') & 'Truncated b/c velocity(' & ,n,': ',i+Dom%ioff,',',j+Dom%joff,',',k,') = ',wrk1(i,j,k) & ,' (m/s) at (lon,lat,klev) = ',Grd%xu(i,j), ',', Grd%yu(i,j), ',', Grd%kmu(i,j) endif enddo enddo enddo endif end subroutine velocity_truncate ! NAME="velocity_truncate" !####################################################################### ! ! ! ! Compute checksum for velocity components ! subroutine ocean_velocity_chksum(Velocity, index, write_advection) type(ocean_velocity_type), intent(in) :: Velocity integer, intent(in) :: index logical, intent(in), optional :: write_advection logical :: writeadvection integer :: stdoutunit stdoutunit=stdout() if(.not. write_a_restart) then write(stdoutunit,'(/a)') '==>Warning from ocean_velocity_mod (ocean_velocity_end): NO restart written.' call mpp_error(WARNING,'==>Warning from ocean_velocity_mod (ocean_velocity_end): NO restart written.') return endif if (PRESENT(write_advection)) then writeadvection = write_advection else writeadvection = .true. endif wrk1_v = 0.0 wrk1_v(isc:iec,jsc:jec,:,1) = Velocity%u(isc:iec,jsc:jec,:,1,index)*Grd%umask(isc:iec,jsc:jec,:) wrk1_v(isc:iec,jsc:jec,:,2) = Velocity%u(isc:iec,jsc:jec,:,2,index)*Grd%umask(isc:iec,jsc:jec,:) write(stdoutunit,*) 'Zonal velocity chksum = ', mpp_chksum(wrk1_v(isc:iec,jsc:jec,:,1)) write(stdoutunit,*) 'Meridional velocity chksum = ', mpp_chksum(wrk1_v(isc:iec,jsc:jec,:,2)) if(tendency==TWO_LEVEL .and. writeadvection) then wrk1_v(isc:iec,jsc:jec,:,1) = Velocity%advection(isc:iec,jsc:jec,:,1,index)*Grd%umask(isc:iec,jsc:jec,:) wrk1_v(isc:iec,jsc:jec,:,2) = Velocity%advection(isc:iec,jsc:jec,:,2,index)*Grd%umask(isc:iec,jsc:jec,:) write(stdoutunit,*) 'Advection of u chksum = ', mpp_chksum(wrk1_v(isc:iec,jsc:jec,:,1)) write(stdoutunit,*) 'Advection of v chksum = ', mpp_chksum(wrk1_v(isc:iec,jsc:jec,:,2)) endif return end subroutine ocean_velocity_chksum ! NAME="ocean_velocity_chksum" !####################################################################### ! ! ! ! Remap in the vertical from s-coordinate to depth and then send to ! diagnostic manager. ! ! This routine is mostly of use for terrain following vertical ! coordinates, which generally deviate a lot from depth or pressure ! coordinates. The zstar and pstar coordinates are very similar ! to z or pressure, so there is no need to do the remapping for ! purposes of visualization. ! ! The routine needs to be made more general and faster. ! It also has been found to be problematic, so it is NOT ! recommended. It remains here as a template for a better algorithm. ! Remapping methods in Ferret are much better. ! ! Use rho_dzu weighting to account for nonBoussinesq. ! ! Author: Stephen.Griffies@noaa.gov ! ! ! subroutine remap_s_to_depth(Thickness, Time, array_in, nvelocity) type(ocean_thickness_type), intent(in) :: Thickness type(ocean_time_type), intent(in) :: Time real, intent(in), dimension(isd:,jsd:,:) :: array_in integer, intent(in) :: nvelocity integer :: i, j, k, kk, tau tau = Time%tau wrk1 = 0.0 wrk2 = 0.0 wrk3 = 0.0 k=1 do kk=1,nk do j=jsd,jed do i=isd,ied if(Thickness%depth_zu(i,j,kk) < Grd%zw(k)) then wrk1(i,j,k) = wrk1(i,j,k) + array_in(i,j,kk)*Thickness%rho_dzu(i,j,kk,tau) wrk2(i,j,k) = wrk2(i,j,k) + Thickness%rho_dzu(i,j,kk,tau) endif enddo enddo enddo do k=1,nk-1 do kk=1,nk do j=jsd,jed do i=isd,ied if(Grd%zw(k) <= Thickness%depth_zu(i,j,kk) .and. Thickness%depth_zu(i,j,kk) < Grd%zw(k+1)) then wrk1(i,j,k+1) = wrk1(i,j,k+1) + array_in(i,j,kk)*Thickness%rho_dzu(i,j,kk,tau) wrk2(i,j,k+1) = wrk2(i,j,k+1) + Thickness%rho_dzu(i,j,kk,tau) endif enddo enddo enddo enddo k=nk do kk=1,nk do j=jsd,jed do i=isd,ied if(Thickness%depth_zu(i,j,kk) > Grd%zw(k) ) then wrk1(i,j,k) = wrk1(i,j,k) + array_in(i,j,kk)*Thickness%rho_dzu(i,j,kk,tau) wrk2(i,j,k) = wrk2(i,j,k) + Thickness%rho_dzu(i,j,kk,tau) endif enddo enddo enddo do k=1,nk do j=jsd,jed do i=isd,ied wrk3(i,j,k) = Grd%umask_depth(i,j,k)*wrk1(i,j,k)/(wrk2(i,j,k)+epsln) enddo enddo enddo used = send_data (id_u_on_depth(nvelocity), wrk3(:,:,:), & Time%model_time,rmask=Grd%umask_depth(:,:,:), & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) end subroutine remap_s_to_depth ! NAME="remap_s_to_depth" end module ocean_velocity_mod