diff --git a/src/mom5/ocean_param/sources/ocean_sponges_eta.F90 b/src/mom5/ocean_param/sources/ocean_sponges_eta.F90 index ad01f5881e..78911d00c9 100644 --- a/src/mom5/ocean_param/sources/ocean_sponges_eta.F90 +++ b/src/mom5/ocean_param/sources/ocean_sponges_eta.F90 @@ -1,395 +1,400 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! !! +!! 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_sponges_eta_mod -! -! Paul Sandery -! -! -! -! Weighted eta tendency [meter*meter/sec] from sponges. -! -! -! -! -! This module applies sponge to eta. The sponges -! can occur at any location and with any distribution in the domain, and -! with any time step and damping rate. Sponges occur where positive -! inverse restore times occur in the field passed to sponge_init. An -! array of eta tendencies due to the sponges is augmented through a -! call to sponge_eta_source. The array of eta tendencies must be -! reset to zero between calls. -! -! Different damping rates can be specified by making -! calls to register_sponge_rate - no sponges are applied to fields for -! which uniformly zero inverse damping rates are set with a call to -! register_sponge_rate. The value towards which a field is damped is -! set with calls to register_sponge_field; successive calls are used to -! set up linear interpolation of this restore rate. -! -! Sponge data and damping coefficients are 2 dimensional. -! -! The user is responsible for providing (and registering) the data on -! the model grid of values towards which the etas are being driven. -! -! -! -! -! -! -! For using this module. Default use_this_module=.false. -! -! -! -! For case when damping coefficients are full 3d field of values. -! Default damp_coeff_3d=.false., which means damping coeffs are -! 2d horizontal array. -! -! -! -! -use diag_manager_mod, only: register_diag_field -use fms_mod, only: write_version_number, open_namelist_file, close_file -use fms_mod, only: file_exist -use fms_mod, only: open_namelist_file, check_nml_error, close_file -use fms_mod, only: read_data, lowercase, FATAL, WARNING, stdout, stdlog -use mpp_mod, only: input_nml_file, mpp_sum, mpp_error, mpp_max -use mpp_domains_mod, only: mpp_global_sum -use time_interp_external_mod, only: init_external_field, time_interp_external -use time_manager_mod, only: time_type, set_date, get_time -use time_manager_mod, only: operator( + ), operator( - ), operator( // ) -use time_manager_mod, only: operator( > ), operator( == ), operator( <= ) - -use ocean_domains_mod, only: get_local_indices -use ocean_parameters_mod, only: missing_value, rho0 -use ocean_types_mod, only: ocean_domain_type, ocean_grid_type, ocean_options_type -use ocean_types_mod, only: ocean_time_type, ocean_External_mode_type -use ocean_workspace_mod, only: wrk1_2d, wrk2_2d -use ocean_util_mod, only: diagnose_2d - -implicit none - -private + ! + ! Paul Sandery + ! + ! + ! + ! Weighted eta tendency [meter*meter/sec] from sponges. + ! + ! + ! + ! + ! This module applies sponge to eta. The sponges + ! can occur at any location and with any distribution in the domain, and + ! with any time step and damping rate. Sponges occur where positive + ! inverse restore times occur in the field passed to sponge_init. An + ! array of eta tendencies due to the sponges is augmented through a + ! call to sponge_eta_source. The array of eta tendencies must be + ! reset to zero between calls. + ! + ! Different damping rates can be specified by making + ! calls to register_sponge_rate - no sponges are applied to fields for + ! which uniformly zero inverse damping rates are set with a call to + ! register_sponge_rate. The value towards which a field is damped is + ! set with calls to register_sponge_field; successive calls are used to + ! set up linear interpolation of this restore rate. + ! + ! Sponge data and damping coefficients are 2 dimensional. + ! + ! The user is responsible for providing (and registering) the data on + ! the model grid of values towards which the etas are being driven. + ! + ! RASF: Added Paul's adaptive restoring scheme + ! + ! + ! + ! + ! + ! + ! For using this module. Default use_this_module=.false. + ! + ! + ! + ! For case when damping coefficients are full 3d field of values. + ! Default damp_coeff_3d=.false., which means damping coeffs are + ! 2d horizontal array. + ! + ! + ! + ! + ! + ! + ! + ! Use adaptive restoring scheme. + ! + ! + ! Normalise???? + ! + ! + ! Force surface height to a fixed value on first time step. + ! + ! + ! ... + ! + ! + ! Shortest restoring timescale? + ! + ! + ! ... + ! + ! + ! Power law for restoring. RMSE is raised to this power. + ! + ! + use diag_manager_mod, only: register_diag_field, send_data + use fms_mod, only: write_version_number, open_namelist_file, close_file + use fms_mod, only: file_exist + use fms_mod, only: open_namelist_file, check_nml_error, close_file + use fms_mod, only: read_data, lowercase, FATAL, WARNING, stdout, stdlog + use mpp_mod, only: mpp_sum, mpp_error, mpp_max + use mpp_domains_mod, only: mpp_global_sum, mpp_update_domains + use time_interp_external_mod, only: init_external_field, time_interp_external + use time_manager_mod, only: time_type, set_date, get_time + use time_manager_mod, only: operator( + ), operator( - ), operator( // ) + use time_manager_mod, only: operator( > ), operator( == ), operator( <= ) + + use ocean_domains_mod, only: get_local_indices + use ocean_parameters_mod, only: missing_value, rho0 + use ocean_types_mod, only: ocean_domain_type, ocean_grid_type, ocean_options_type + use ocean_types_mod, only: ocean_time_type, ocean_External_mode_type + use ocean_workspace_mod, only: wrk1_2d, wrk2_2d + + implicit none + + private #include -type ocean_sponge_type - integer :: id ! time_interp_external index - character(len=32) :: name ! eta name corresponding to sponge - real, dimension(:,:) , pointer :: damp_coeff2d => NULL() ! 2d inverse damping rate (eta units/ sec) -end type ocean_sponge_type - -type(ocean_sponge_type), allocatable, dimension(:) :: Sponge -type(ocean_domain_type), pointer :: Dom => NULL() -type(ocean_grid_type), pointer :: Grd => NULL() - -public ocean_sponges_eta_init -public sponge_eta_source - -character(len=126) :: version = '$Id: ocean_sponges_eta.F90,v 20.0 2013/12/14 00:16:22 fms Exp $' -character (len=128) :: tagname = '$Name: tikal $' - -! for diagnostics -logical :: used -integer, dimension(:), allocatable :: id_sponge_tend -logical :: module_is_initialized = .false. -logical :: use_this_module = .false. - -! Adaptive restoring -real :: athresh = 0.5 -real :: npower = 1.0 -real :: sdiffo = 0.5 -real :: lambda = 0.0083 -real :: taumin = 720 -logical :: use_adaptive_restore = .false. -logical :: use_sponge_after_init = .false. -logical :: use_normalising = .false. -logical :: use_hard_thump = .false. -integer :: days_to_restore = 1 -integer :: secs_to_restore = 0 - -integer :: secs_restore -integer :: initial_day, initial_secs - -namelist /ocean_sponges_eta_nml/ use_this_module -namelist /ocean_sponges_eta_ofam_nml/ & - use_adaptive_restore, use_sponge_after_init, use_normalising, & - use_hard_thump, athresh, taumin, lambda, npower, days_to_restore, & - secs_to_restore + type ocean_sponge_type + integer :: id ! time_interp_external index + character(len=32) :: name ! eta name corresponding to sponge + real, dimension(:,:) , pointer :: damp_coeff2d => NULL() ! 2d inverse damping rate (eta units/ sec) + end type ocean_sponge_type + + type(ocean_sponge_type), allocatable, dimension(:) :: Sponge + type(ocean_domain_type), pointer :: Dom => NULL() + type(ocean_grid_type), pointer :: Grd => NULL() + + public ocean_sponges_eta_init + public sponge_eta_source + + character(len=126) :: version = '$Id: ocean_sponges_eta.F90,v 1.1.2.4.42.1 2009/10/10 00:42:52 nnz Exp $' + character (len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' + + ! for diagnostics + logical :: used + integer, dimension(:), allocatable :: id_sponge_tend + logical :: module_is_initialized = .false. + logical :: use_this_module = .false. + + ! Adaptive restoring + real :: athresh = 0.5 + real :: npower = 1.0 + real :: sdiffo = 0.5 + real :: lambda = 0.0083 + real :: taumin = 3600 + logical :: use_adaptive_restore = .false. + logical :: use_sponge_after_init = .false. + logical :: use_normalising = .false. + logical :: use_hard_thump = .false. + integer :: secs_to_restore = 0 + integer :: days_to_restore = 0 + integer :: secs_restore + integer :: initial_day, initial_secs + + namelist /ocean_sponges_eta_nml/ use_this_module + namelist /ocean_sponges_eta_OFAM_nml/ use_adaptive_restore,use_sponge_after_init, use_normalising, use_hard_thump, & + athresh, taumin,lambda,npower,days_to_restore, secs_to_restore contains - -!####################################################################### -! -! -! -! This subroutine is intended to be used to initialize the sponges. -! Everything in this subroutine is a user prototype, and should be replacable. -! -! -subroutine ocean_sponges_eta_init(Grid, Domain, Time, dtime, Ocean_options) - - type(ocean_grid_type), intent(in), target :: Grid - type(ocean_domain_type), intent(in), target :: Domain - type(ocean_time_type), intent(in) :: Time - real, intent(in) :: dtime - type(ocean_options_type), intent(inout) :: Ocean_options - - integer :: i, j - integer :: ioun, io_status, ierr - integer :: secs, days - real :: dtimer - - character(len=128) :: name - - integer :: stdoutunit,stdlogunit - stdoutunit=stdout();stdlogunit=stdlog() - - if ( module_is_initialized ) then - call mpp_error(FATAL, & - '==>Error in ocean_sponges_eta_mod (ocean_sponges_init): module already initialized') - endif - - module_is_initialized = .TRUE. - - allocate( Sponge(1) ) - - call write_version_number(version, tagname) - - ! provide for namelist over-ride of default values -#ifdef INTERNAL_FILE_NML - read(input_nml_file, nml=ocean_sponges_eta_nml, iostat=io_status) - ierr = check_nml_error(io_status, 'ocean_sponges_eta_nml') - read(input_nml_file, nml=ocean_sponges_eta_ofam_nml, iostat=io_status) - ierr = check_nml_error(io_status, 'ocean_sponges_eta_ofam_nml') -#else - ioun = open_namelist_file() - read(ioun, ocean_sponges_eta_nml, IOSTAT=io_status) - ierr = check_nml_error(io_status, 'ocean_sponges_eta_nml') - rewind(ioun) - read(ioun, ocean_sponges_eta_ofam_nml, IOSTAT=io_status) - ierr = check_nml_error(io_status, 'ocean_sponges_eta_nml') - call close_file (ioun) -#endif - write (stdlogunit, ocean_sponges_eta_nml) - write (stdoutunit, '(/)') - write (stdoutunit, ocean_sponges_eta_nml) - - write (stdlogunit, ocean_sponges_eta_ofam_nml) - write (stdoutunit, '(/)') - write (stdoutunit, ocean_sponges_eta_ofam_nml) - - Dom => Domain - Grd => Grid - -#ifndef MOM_STATIC_ARRAYS - call get_local_indices(Domain, isd, ied, jsd, jed, isc, iec, jsc, jec) - nk = Grid%nk + !####################################################################### + ! + ! + ! + ! This subroutine is intended to be used to initialize the sponges. + ! Everything in this subroutine is a user prototype, and should be replacable. + ! + ! + subroutine ocean_sponges_eta_init(Grid, Domain, Time, dtime, Ocean_options) + + type(ocean_grid_type), intent(in), target :: Grid + type(ocean_domain_type), intent(in), target :: Domain + type(ocean_time_type), intent(in) :: Time + real, intent(in) :: dtime + type(ocean_options_type), intent(inout) :: Ocean_options + + integer :: i, j, k, n + integer :: ioun, io_status, ierr + integer :: secs, days + real :: dtimer + + character(len=128) :: name + + integer :: stdoutunit,stdlogunit + stdoutunit=stdout();stdlogunit=stdlog() + + if ( module_is_initialized ) then + call mpp_error(FATAL, & + '==>Error in ocean_sponges_eta_mod (ocean_sponges_init): module already initialized') + endif + + module_is_initialized = .TRUE. + + allocate( Sponge(1) ) + + call write_version_number() + + ! provide for namelist over-ride of default values + ioun = open_namelist_file() + read (ioun,ocean_sponges_eta_nml,IOSTAT=io_status) + write (stdlogunit,ocean_sponges_eta_nml) + write (stdoutunit,'(/)') + write (stdoutunit,ocean_sponges_eta_nml) + ierr = check_nml_error(io_status,'ocean_sponges_eta_nml') + call close_file (ioun) + + ! provide for namelist over-ride of default values + ioun = open_namelist_file() + read (ioun,ocean_sponges_eta_OFAM_nml,IOSTAT=io_status) + write (stdlogunit,ocean_sponges_eta_OFAM_nml) + write (stdoutunit,'(/)') + write (stdoutunit,ocean_sponges_eta_OFAM_nml) + ierr = check_nml_error(io_status,'ocean_sponges_eta_OFAM_nml') + call close_file (ioun) + + Dom => Domain + Grd => Grid + +#ifndef MOM4_STATIC_ARRAYS + call get_local_indices(Domain, isd, ied, jsd, jed, isc, iec, jsc, jec) + nk = Grid%nk #endif - Sponge(1)%id = -1 - dtimer = 1.0/dtime - - if(use_this_module) then - write(stdoutunit,*)'==>Note from ocean_sponges_eta_mod: Using this module.' - Ocean_options%ocean_sponges_eta= 'Used ocean eta sponges.' - else - write(stdoutunit,*)'==>Note from ocean_sponges_eta_mod: NOT using this module.' - Ocean_options%ocean_sponges_eta= 'Did NOT use ocean eta sponges.' - return - endif - - ! Set up some initial times - call get_time(Time%model_time, secs, days) - initial_day = days - initial_secs = secs - - !read damping rates for eta - if (use_adaptive_restore) then - allocate(Sponge(1)%damp_coeff2d(isd:ied, jsd:jed)) - Sponge(1)%damp_coeff2d(:,:) = 0.0 - else - name = 'INPUT/eta_sponge_coeff.nc' - if (file_exist(name)) then - write(stdoutunit,*) '==> Using sponge damping times specified from file ',name - allocate(Sponge(1)%damp_coeff2d(isd:ied,jsd:jed)) - Sponge(1)%damp_coeff2d(:,:) = 0.0 - call read_data(name,'coeff',Sponge(1)%damp_coeff2d,domain=Dom%domain2d,timelevel=1) - - do j=jsc,jec - do i=isc,iec - if(Grd%tmask(i,j,1) == 0.0) then - Sponge(1)%damp_coeff2d(i,j) = 0.0 - endif - enddo - enddo - - - ! modify damping rates to allow restoring to be solved implicitly - ! note: test values between zero and 4.0e-8 revert to damping rates defined above - do j=jsc,jec - do i=isc,iec - if (dtime*Sponge(1)%damp_coeff2d(i,j) > 4.0e-8) then - if (dtime*Sponge(1)%damp_coeff2d(i,j) > 37.0) then - Sponge(1)%damp_coeff2d(i,j) = dtimer - else - Sponge(1)%damp_coeff2d(i,j) = (1.0 - exp(-dtime*Sponge(1)%damp_coeff2d(i,j))) * dtimer - endif - else if (dtime*Sponge(1)%damp_coeff2d(i,j) <= 0.0) then - Sponge(1)%damp_coeff2d(i,j) = 0.0 - endif - enddo - enddo - - endif ! endif for fileexist 'INPUT/eta_sponge_coeff.nc' - endif - - - name = 'INPUT/eta_sponge.nc' - if(file_exist(trim(name))) then - Sponge(1)%id = init_external_field(name,'eta_t',domain=Domain%domain2d) - if (Sponge(1)%id < 1) then - if ( use_adaptive_restore ) then - call mpp_error(FATAL,& - '==>Error: in ocean_sponges_eta_mod: adaptive restoring is specified but sponge values are not') - else + Sponge(1)%id = -1 + dtimer = 1.0/dtime + + if(use_this_module) then + write(stdoutunit,*)'==>Note from ocean_sponges_eta_mod: Using this module.' + Ocean_options%ocean_sponges_eta= 'Used ocean eta sponges.' + else + write(stdoutunit,*)'==>Note from ocean_sponges_eta_mod: NOT using this module.' + Ocean_options%ocean_sponges_eta= 'Did NOT use ocean eta sponges.' + return + endif + + !set up some initial times + call get_time(Time%model_time,secs,days) + initial_day = days + initial_secs = secs + + !read damping rates for eta + if (.not. use_hard_thump .and. .not. use_adaptive_restore) then + allocate(Sponge(1)%damp_coeff2d(isd:ied,jsd:jed)) + Sponge(1)%damp_coeff2d(:,:) = 0.0 + name = 'INPUT/eta_sponge_coeff.nc' + if (file_exist(name)) then + write(stdoutunit,*) '==> Using sponge damping times specified from file ',name + allocate(Sponge(1)%damp_coeff2d(isd:ied,jsd:jed)) + Sponge(1)%damp_coeff2d(:,:) = 0.0 + call read_data(name,'coeff',Sponge(1)%damp_coeff2d,domain=Dom%domain2d,timelevel=1) + do j=jsc,jec + do i=isc,iec + if(Grd%tmask(i,j,1) == 0.0) then + Sponge(1)%damp_coeff2d(i,j) = 0.0 + endif + enddo + enddo + !modify damping rates to allow restoring to be solved implicitly + !note: test values between zero and 4.0e-8 revert to damping rates defined above + do j=jsc,jec + do i=isc,iec + if (dtime*Sponge(1)%damp_coeff2d(i,j) > 4.0e-8) then + if (dtime*Sponge(1)%damp_coeff2d(i,j) > 37.0) then + Sponge(1)%damp_coeff2d(i,j) = dtimer + else + Sponge(1)%damp_coeff2d(i,j) = (1.0 - exp(-dtime*Sponge(1)%damp_coeff2d(i,j))) * dtimer + endif + else if (dtime*Sponge(1)%damp_coeff2d(i,j) <= 0.0) then + Sponge(1)%damp_coeff2d(i,j) = 0.0 + endif + enddo + enddo + endif !(file_exist(name)) + endif !.not. use_hard_thump .and. .not. use_adaptive_restore + + name = 'INPUT/eta_sponge.nc' + if(file_exist(trim(name))) then + Sponge(1)%id = init_external_field(name,'eta_t',domain=Domain%domain2d) + if (Sponge(1)%id < 1) then call mpp_error(FATAL,& - '==>Error: in ocean_sponges_eta_mod: damping rates are specified but sponge values are not') - endif - endif - write(stdoutunit,*) '==> Using increment data specified from file '//trim(name) - else - write(stdoutunit,*) '==> '//trim(name)//' not found. Increment not being applied ' - endif - - allocate (id_sponge_tend(1)) - id_sponge_tend = -1 - - id_sponge_tend = register_diag_field ('ocean_model', 'eta_t_sponge_tend', & - Grd%tracer_axes(1:2), Time%model_time, 'eta_t tendency due to sponge', & - 'm/s', missing_value=missing_value, range=(/-1.e10,1.e10/)) - - -end subroutine ocean_sponges_eta_init -! NAME="ocean_sponges_init" - - -!####################################################################### -! -! -! -! This subroutine calculates thickness weighted and density weighted -! time tendencies due to damping by sponges or damping through adaptive -! restoring. -! -! -subroutine sponge_eta_source(Time, Ext_mode) - - type(ocean_time_type), intent(in) :: Time - type(ocean_external_mode_type), intent(inout) :: Ext_mode - - integer :: taum1, tau - integer :: i, j - - real :: sdiff, sum_val, adaptive_coeff - integer :: secs, days, numsecs - - logical :: do_adaptive_restore = .false. ! Only restore in specified time period. - logical :: do_normal_restore = .false. ! Only restore in specified time period. - logical, save :: first_pass = .true. - real, save :: num_surface_wet_cells - - if(.not. use_this_module) return - - taum1 = Time%taum1 - tau = Time%tau - wrk1_2d = 0.0 - wrk2_2d = 0.0 - - if (use_adaptive_restore) then - if (first_pass) then - num_surface_wet_cells = mpp_global_sum(Dom%domain2d, Grd%tmask(:,:,1)) - end if - - secs_restore = days_to_restore * 86400 + secs_to_restore - sdiff = 0.0 - call get_time(Time%model_time, secs, days) - numsecs = (days - initial_day) * 86400 + secs - initial_secs - - do_adaptive_restore = (numsecs < secs_restore) - do_normal_restore = (.not. do_adaptive_restore) .and. use_sponge_after_init - else - do_normal_restore = .true. - end if - - if (Sponge(1)%id > 0) then - - call time_interp_external(Sponge(1)%id, Time%model_time, wrk1_2d) - - if (do_adaptive_restore) then - if (first_pass) then - if (use_hard_thump) then - do j = jsd, jed - do i = isd, ied - Ext_mode%eta_t(i, j, taum1) = wrk1_2d(i, j) - Ext_mode%eta_t(i, j, tau) = wrk1_2d(i, j) - end do - end do - end if - - wrk2_2d = abs(Ext_mode%eta_t(:,:,taum1) - wrk1_2d(:,:)) & - * Grd%tmask(:,:,1) - sum_val = sum(wrk2_2d(isc:iec, jsc:jec)) - - call mpp_sum(sum_val) - sdiffo = sum_val / num_surface_wet_cells - write (stdout(), *) 'mean eta_diff', sdiffo - - if (.not. use_normalising) then - sdiffo = 1.0 - else - sdiffo = maxval(wrk2_2d) + '==>Error: in ocean_sponges_eta_mod: sponge values are required but not specified') + endif + write(stdoutunit,*) '==> Using sponge data specified from file '//trim(name) + else + write(stdoutunit,*) '==> '//trim(name)//' not found. Sponge not being applied ' + endif + + allocate (id_sponge_tend(1)) + id_sponge_tend = -1 ! needed? + id_sponge_tend = register_diag_field ('ocean_model', 'eta_t_sponge_tend', & + Grd%tracer_axes(1:2), Time%model_time, 'eta_t tendency due to sponge', & + 'm/s', missing_value=missing_value, range=(/-1.e-2,1.e-2/)) + + end subroutine ocean_sponges_eta_init + ! NAME="ocean_sponges_init" + + + !####################################################################### + ! + ! + ! + ! This subroutine calculates thickness weighted and density weighted + ! time tendencies due to damping by sponges or damping through adaptive + ! restoring. + ! + ! + subroutine sponge_eta_source(Time, Ext_mode) + type(ocean_time_type), intent(in) :: Time + type(ocean_external_mode_type), intent(inout) :: Ext_mode + integer :: taum1, tau + integer :: i, j + real :: sdiff=0.0, sum_val, adaptive_coeff + integer :: secs, days, numsecs + logical :: do_adaptive_restore = .false. ! Only restore in specified time period. + logical,save :: first_pass = .true. + real, save :: num_surface_wet_cells + + if(.not. use_this_module) return + + taum1 = Time%taum1 + tau = Time%tau + wrk1_2d = 0.0 + wrk2_2d = 0.0 + call get_time(Time%model_time,secs,days) + !----------------------------------------------------------------------------------------------------------------- + if (Sponge(1)%id > 0) then + call time_interp_external(Sponge(1)%id, Time%model_time, wrk1_2d) + if (first_pass) then + if (use_hard_thump) then + do j=jsd,jed + do i=isd,ied + Ext_mode%eta_t(i,j,taum1)=wrk1_2d(i,j) ! 6: + Ext_mode%eta_t(i,j,tau)=wrk1_2d(i,j) + enddo + enddo + call mpp_update_domains (Ext_mode%eta_t(:,:,taum1) ,Dom%domain2d) + call mpp_update_domains (Ext_mode%eta_t(:,:,tau) ,Dom%domain2d) + wrk2_2d=abs(Ext_mode%eta_t(:,:,taum1)-wrk1_2d(:,:))*Grd%tmask(:,:,1) ! 8: 9/10 only get called if H=0 so 8 must mean H=1 + endif !use_hard_thump + sum_val=sum(wrk2_2d(isc:iec,jsc:jec)) + call mpp_sum(sum_val) + num_surface_wet_cells = mpp_global_sum(Dom%domain2d,Grd%tmask(:,:,1)) ! 1: + sdiffo = sum_val/num_surface_wet_cells + write (stdout(),*)'mean absolute deviation eta ', sdiffo + if ( use_normalising == .false.) then + sdiffo=1.0 + else + sdiffo = maxval(wrk2_2d) call mpp_max(sdiffo) - sdiffo = 1.01 * sdiffo - end if - - first_pass = .false. - end if - - ! Calculate idealised forcing tendency timescale for each timestep and - ! array element - do j = jsd, jed - do i = isd, ied - sdiff = abs(Ext_mode%eta_t(i, j, taum1) - wrk1_2d(i, j)) - sdiff = max(sdiff, 1e-9) !minimum tolerance - adaptive_coeff = 1. / (taumin - (lambda * real(secs_restore)) & - / (((sdiff / sdiffo)**npower) & - * log(1 - athresh))) - wrk2_2d(i,j) = adaptive_coeff & - * (wrk1_2d(i, j) - Ext_mode%eta_t(i, j, taum1)) - end do - end do - - else if (do_normal_restore) then - do j = jsc, jec - do i = isc, iec - wrk2_2d(i,j) = Sponge(1)%damp_coeff2d(i, j) & - * (wrk1_2d(i, j) - Ext_mode%eta_t(i, j, taum1)) - end do - end do - end if - - do j = jsc, jec - do i = isc, iec - Ext_mode%source(i,j) = Ext_mode%source(i,j) + rho0 * wrk2_2d(i,j) - end do - end do - end if - - call diagnose_2d(Time, Grd, id_sponge_tend(1), wrk2_2d(:,:)) - - return - -end subroutine sponge_eta_source -! NAME="sponge_eta_source" - - + sdiffo = 1.01*sdiffo + endif + endif !first_pass + !-------------------------------------------------------------------------------------------------------------- + if (use_hard_thump == .false.) then + if (use_adaptive_restore) then + secs_restore = days_to_restore*86400 + secs_to_restore + numsecs = (days-initial_day)*86400 + secs - initial_secs + do_adaptive_restore = ( numsecs < secs_restore ) + endif + if (do_adaptive_restore) then + ! calculate idealised forcing tendency timescale for each timestep and array element + do j=jsd,jed + do i=isd,ied + sdiff=abs(Ext_mode%eta_t(i,j,taum1)-wrk1_2d(i,j)) + sdiff = max(sdiff,1e-9) ! minimum tolerance + adaptive_coeff =1./(taumin - (lambda*real(secs_restore))/(((sdiff/sdiffo)**npower)*log(1-athresh))) + wrk2_2d(i,j) = adaptive_coeff * (wrk1_2d(i,j) - Ext_mode%eta_t(i,j,taum1)) + enddo + enddo + else !normal restore + do j=jsd,jed + do i=isd,ied + wrk2_2d(i,j) = Sponge(1)%damp_coeff2d(i,j)*(wrk1_2d(i,j) - Ext_mode%eta_t(i,j,taum1)) + enddo + enddo + endif + endif !use_hard_thump + !-------------------------------------------------------------------------------------------------------------- + do j=jsc,jec + do i=isc,iec + Ext_mode%source(i,j) = Ext_mode%source(i,j) + rho0*wrk2_2d(i,j) + enddo + enddo + endif !Sponge(1)%id + !----------------------------------------------------------------------------------------------------------------- + if (id_sponge_tend(1) > 0) used = send_data(id_sponge_tend(1), & + wrk2_2d(:,:), Time%model_time, rmask=Grd%tmask(:,:,1), & + is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) + + first_pass = .false. + + return + + end subroutine sponge_eta_source + ! NAME="sponge_eta_source" end module ocean_sponges_eta_mod diff --git a/src/mom5/ocean_param/sources/ocean_sponges_tracer.F90 b/src/mom5/ocean_param/sources/ocean_sponges_tracer.F90 index bf98f924c5..ea53910a29 100644 --- a/src/mom5/ocean_param/sources/ocean_sponges_tracer.F90 +++ b/src/mom5/ocean_param/sources/ocean_sponges_tracer.F90 @@ -1,496 +1,562 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! !! +!! 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_sponges_tracer_mod -! -! Bonnie Samuels -! -! -! R.W. Hallberg -! -! -! M.J. Harrison -! -! -! P. S. Swathi -! -! -! S. M. Griffies -! -! -! -! Thickness weighted tracer tendency [tracer*meter/sec] from sponges. -! -! -! -! This module applies sponges to tracers. The sponges -! can occur at any location and with any distribution in the domain, and -! with any time step and damping rate. Sponges occur where positive -! inverse restore times occur in the field passed to sponge_init. An -! array of tracer tendencies due to the sponges is augmented through a -! call to sponge_tracer_source. The array of tracer tendencies must be -! reset to zero between calls. -! -! Different damping rates can be specified for each field by making -! calls to register_sponge_rate - no sponges are applied to fields for -! which uniformly zero inverse damping rates are set with a call to -! register_sponge_rate. The value towards which a field is damped is -! set with calls to register_sponge_field; successive calls are used to -! set up linear interpolation of this restore rate. -! -! Sponge data and damping coefficients are generally 3 dimensional. -! -! The user is responsible for providing (and registering) the data on -! the model grid of values towards which the tracers are being driven. -! -! -! -! -! For using this module. Default use_this_module=.false. -! -! -! For case when damping coefficients are full 3d field of values. -! Default damp_coeff_3d=.false., which means damping coeffs are -! 2d horizontal array. -! -! -! -use diag_manager_mod, only: register_diag_field -use fms_mod, only: write_version_number, open_namelist_file, close_file -use fms_mod, only: file_exist -use fms_mod, only: open_namelist_file, check_nml_error, close_file -use fms_mod, only: read_data, lowercase, FATAL, WARNING, stdout, stdlog -use mpp_mod, only: input_nml_file, mpp_sum, mpp_error, mpp_max -use time_interp_external_mod, only: init_external_field, time_interp_external -use time_manager_mod, only: time_type, set_date, get_time -use time_manager_mod, only: operator( + ), operator( - ), operator( // ) -use time_manager_mod, only: operator( > ), operator( == ), operator( <= ) - -use ocean_domains_mod, only: get_local_indices -use ocean_parameters_mod, only: missing_value -use ocean_types_mod, only: ocean_domain_type, ocean_grid_type, ocean_thickness_type -use ocean_types_mod, only: ocean_prog_tracer_type, ocean_options_type, ocean_time_type -use ocean_workspace_mod, only: wrk1, wrk2 -use ocean_util_mod, only: diagnose_3d - -implicit none - -private + ! + ! Bonnie Samuels + ! + ! + ! R.W. Hallberg + ! + ! + ! M.J. Harrison + ! + ! + ! P. S. Swathi + ! + ! + ! R. C. Pacanowski + ! + ! + ! S. M. Griffies + ! + ! + ! + ! Thickness weighted tracer tendency [tracer*meter/sec] from sponges. + ! + ! + ! + ! This module applies sponges to tracers. The sponges + ! can occur at any location and with any distribution in the domain, and + ! with any time step and damping rate. Sponges occur where positive + ! inverse restore times occur in the field passed to sponge_init. An + ! array of tracer tendencies due to the sponges is augmented through a + ! call to sponge_tracer_source. The array of tracer tendencies must be + ! reset to zero between calls. + ! + ! Different damping rates can be specified for each field by making + ! calls to register_sponge_rate - no sponges are applied to fields for + ! which uniformly zero inverse damping rates are set with a call to + ! register_sponge_rate. The value towards which a field is damped is + ! set with calls to register_sponge_field; successive calls are used to + ! set up linear interpolation of this restore rate. + ! + ! Sponge data and damping coefficients are generally 3 dimensional. + ! + ! The user is responsible for providing (and registering) the data on + ! the model grid of values towards which the tracers are being driven. + ! + ! OFAM MODIFICATIONS: + ! + ! Adaptive nudiging is permitted. + ! We can limit the length of time for which nudging is applied. + ! There is an option to restore to a climatology in the bottom 4 layers + ! We can limit the temp and salt to fix values. + + ! + ! + ! + ! + ! For using this module. Default use_this_module=.false. + ! + ! + ! For case when damping coefficients are full 3d field of values. + ! Default damp_coeff_3d=.false., which means damping coeffs are + ! 2d horizontal array. + ! + ! + ! + use diag_manager_mod, only: register_diag_field, send_data + use fms_mod, only: write_version_number, open_namelist_file, close_file + use fms_mod, only: file_exist + use fms_mod, only: open_namelist_file, check_nml_error, close_file + use fms_mod, only: read_data, lowercase, FATAL, WARNING, stdout, stdlog + use mpp_mod, only: mpp_sum, mpp_error, mpp_max, mpp_min,mpp_pe + use time_interp_external_mod, only: init_external_field, time_interp_external + use time_manager_mod, only: time_type, set_date, get_time + use time_manager_mod, only: operator( + ), operator( - ), operator( // ) + use time_manager_mod, only: operator( > ), operator( == ), operator( <= ) + use mpp_domains_mod, only: mpp_update_domains + use ocean_domains_mod, only: get_local_indices + use ocean_parameters_mod, only: missing_value + use ocean_types_mod, only: ocean_domain_type, ocean_grid_type, ocean_thickness_type + use ocean_types_mod, only: ocean_prog_tracer_type, ocean_options_type, ocean_time_type + use ocean_workspace_mod, only: wrk1, wrk2 + + implicit none + + private #include -type ocean_sponge_type - integer :: id ! time_interp_external index - character(len=32) :: name ! tracer name corresponding to sponge - real, dimension(:,:,:), pointer :: damp_coeff => NULL() ! 3d inverse damping rate (tracer units/ sec) - real, dimension(:,:) , pointer :: damp_coeff2d => NULL() ! 2d inverse damping rate (tracer units/ sec) -end type ocean_sponge_type - - -type(ocean_sponge_type), allocatable, dimension(:) :: Sponge -type(ocean_domain_type), pointer :: Dom => NULL() -type(ocean_grid_type), pointer :: Grd => NULL() - -public ocean_sponges_tracer_init -public sponge_tracer_source - -character(len=126) :: version = '$Id: ocean_sponges_tracer.F90,v 20.0 2013/12/14 00:16:24 fms Exp $' -character (len=128) :: tagname = '$Name: tikal $' - -! for diagnostics -logical :: used -integer, dimension(:), allocatable :: id_sponge_tend - -integer :: num_prog_tracers = 0 -logical :: module_is_initialized = .FALSE. -logical :: damp_coeff_3d = .false. -logical :: use_this_module = .false. - -! Adaptive restoring -real :: athresh = 0.5 -real :: npower = 1.0 -real :: lambda = 0.0083 -real :: taumin = 720 -logical :: use_adaptive_restore = .false. -logical :: use_sponge_after_init = .false. -logical :: use_normalising = .false. -logical :: use_hard_thump = .false. -logical :: deflate = .false. -real :: deflate_fraction = 0.6 -integer :: secs_to_restore = 0 -integer :: days_to_restore = 1 - -logical :: limit_temp = .false. -real :: limit_temp_min = -1.8 -real :: limit_temp_restore = 10800. - -logical :: limit_salt = .false. -real :: limit_salt_min = 0.01 -real :: limit_salt_restore = 3600. - -integer :: secs_restore -integer :: initial_day -integer :: initial_secs -real, allocatable :: sdiffo(:) - -namelist /ocean_sponges_tracer_nml/ use_this_module, damp_coeff_3d -namelist /ocean_sponges_tracer_ofam_nml/ & - use_adaptive_restore, use_sponge_after_init, use_normalising, & - use_hard_thump, athresh, taumin, lambda, npower, days_to_restore, & - secs_to_restore, deflate, deflate_fraction, & - limit_temp, limit_temp_min, limit_temp_restore, & - limit_salt, limit_salt_min, limit_salt_restore + type ocean_sponge_type + integer :: id ! time_interp_external index + character(len=32) :: name ! tracer name corresponding to sponge + real, dimension(:,:,:), pointer :: damp_coeff => NULL() ! 3d inverse damping rate (tracer units/ sec) + real, dimension(:,:) , pointer :: damp_coeff2d => NULL() ! 2d inverse damping rate (tracer units/ sec) + end type ocean_sponge_type + + type ocean_cars_type + integer :: id ! time_interp_external index + character(len=32) :: name ! tracer name corresponding to sponge + real, dimension(:,:,:) , pointer :: field => NULL() ! field to be damped + end type ocean_cars_type + + + + type(ocean_sponge_type), allocatable, dimension(:) :: Sponge + type(ocean_cars_type), allocatable, dimension(:) :: Cars + type(ocean_domain_type), pointer :: Dom => NULL() + type(ocean_grid_type), pointer :: Grd => NULL() + + public ocean_sponges_tracer_init + public sponge_tracer_source + + character(len=126) :: version = '$Id: ocean_sponges_tracer.F90,v 1.1.2.4.42.1 2009/10/10 00:42:53 nnz Exp $' + character (len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' + + ! for diagnostics + logical :: used + integer, dimension(:), allocatable :: id_sponge_tend + + integer :: num_prog_tracers = 0 + logical :: module_is_initialized = .FALSE. + logical :: damp_coeff_3d = .false. + logical :: use_this_module = .false. + + ! Adaptive restoring + + real :: athresh = 0.5 + real :: npower = 1.0 + real :: lambda = 0.0083 + real :: taumin = 3600 + logical :: use_adaptive_restore = .false. + logical :: use_sponge_after_init = .false. + logical :: use_normalising = .false. + logical :: use_hard_thump = .false. + logical :: use_increment = .false. + integer :: secs_to_restore = 0 + integer :: days_to_restore = 0 + integer :: secs_restore + integer :: initial_day + integer :: initial_secs + real,allocatable :: sdiffo(:) + logical :: restore_cars_deep = .false. + logical :: limit_temp = .false. + logical :: limit_salt = .false. + real :: limit_temp_max = 40.0 + real :: limit_temp_min = -1.8 + real :: limit_salt_max = 40.0 + real :: limit_salt_min = 0.1 + + namelist /ocean_sponges_tracer_nml/ use_this_module, damp_coeff_3d + namelist /ocean_sponges_tracer_OFAM_nml/ use_adaptive_restore,use_sponge_after_init, use_normalising, use_hard_thump, & + athresh, taumin,lambda,npower,days_to_restore, secs_to_restore, & + restore_cars_deep, limit_temp, limit_salt, & + limit_temp_max, limit_temp_min, limit_salt_max, limit_salt_min contains + !####################################################################### + ! + ! + ! + ! This subroutine is intended to be used to initialize the tracer sponges. + ! Everything in this subroutine is a user prototype, and should be replacable. + ! + ! + subroutine ocean_sponges_tracer_init(Grid, Domain, Time, T_prog, dtime, Ocean_options) + + type(ocean_grid_type), intent(in), target :: Grid + type(ocean_domain_type), intent(in), target :: Domain + type(ocean_time_type), intent(in) :: Time + type(ocean_prog_tracer_type), intent(in) :: T_prog(:) + real, intent(in) :: dtime + type(ocean_options_type), intent(inout) :: Ocean_options + + integer :: i, j, k, n + integer :: ioun, io_status, ierr + integer :: index_temp + integer :: secs, days + real :: dtimer + + character(len=128) :: name + + integer :: stdoutunit,stdlogunit + stdoutunit=stdout();stdlogunit=stdlog() + + if ( module_is_initialized ) then + call mpp_error(FATAL, & + '==>Error in ocean_sponges_tracer_mod (ocean_sponges_tracer_init): module already initialized') + endif -!####################################################################### -! -! -! -! This subroutine is intended to be used to initialize the tracer sponges. -! Everything in this subroutine is a user prototype, and should be replacable. -! -! -subroutine ocean_sponges_tracer_init(Grid, Domain, Time, T_prog, dtime, Ocean_options) - - type(ocean_grid_type), intent(in), target :: Grid - type(ocean_domain_type), intent(in), target :: Domain - type(ocean_time_type), intent(in) :: Time - type(ocean_prog_tracer_type), intent(in) :: T_prog(:) - real, intent(in) :: dtime - type(ocean_options_type), intent(inout) :: Ocean_options - - integer :: i, j, k, n - integer :: ioun, io_status, ierr - integer :: index_temp - integer :: secs, days - real :: dtimer - - character(len=128) :: name - - integer :: stdoutunit,stdlogunit - stdoutunit=stdout();stdlogunit=stdlog() - - if ( module_is_initialized ) then - call mpp_error(FATAL, & - '==>Error in ocean_sponges_tracer_mod (ocean_sponges_tracer_init): module already initialized') - endif - - module_is_initialized = .TRUE. - - num_prog_tracers = size(T_prog(:)) - do n=1,num_prog_tracers - if (trim(T_prog(n)%name) == 'temp') index_temp = n - enddo - - allocate( Sponge(num_prog_tracers) ) - - call write_version_number(version, tagname) - - ! provide for namelist over-ride of default values -#ifdef INTERNAL_FILE_NML - read(input_nml_file, nml=ocean_sponges_tracer_nml, iostat=io_status) - ierr = check_nml_error(io_status, 'ocean_sponges_tracer_nml') - read(input_nml_file, nml=ocean_sponges_tracer_ofam_nml, iostat=io_status) - ierr = check_nml_error(io_status, 'ocean_sponges_tracer_ofam_nml') -#else - ioun = open_namelist_file() - read(ioun, ocean_sponges_tracer_nml, IOSTAT=io_status) - ierr = check_nml_error(io_status, 'ocean_sponges_tracer_nml') - rewind(ioun) - read(ioun, ocean_sponges_tracer_ofam_nml, IOSTAT=io_status) - ierr = check_nml_error(io_status, 'ocean_sponges_tracer_ofam_nml') - call close_file(ioun) -#endif - write(stdlogunit, ocean_sponges_tracer_nml) - write(stdoutunit, '(/)') - write(stdoutunit, ocean_sponges_tracer_nml) - - write(stdlogunit, ocean_sponges_tracer_ofam_nml) - write(stdoutunit, '(/)') - write(stdoutunit, ocean_sponges_tracer_ofam_nml) - - Dom => Domain - Grd => Grid - -#ifndef MOM_STATIC_ARRAYS - call get_local_indices(Domain, isd, ied, jsd, jed, isc, iec, jsc, jec) - nk = Grid%nk + module_is_initialized = .TRUE. + + num_prog_tracers = size(T_prog(:)) + do n=1,num_prog_tracers + if (trim(T_prog(n)%name) == 'temp') index_temp = n + enddo + + allocate( Sponge(num_prog_tracers) ) + allocate( Cars(num_prog_tracers) ) + + call write_version_number() + + ! provide for namelist over-ride of default values + ioun = open_namelist_file() + read (ioun,ocean_sponges_tracer_nml,IOSTAT=io_status) + write (stdlogunit,ocean_sponges_tracer_nml) + write (stdoutunit,'(/)') + write (stdoutunit,ocean_sponges_tracer_nml) + ierr = check_nml_error(io_status,'ocean_sponges_tracer_nml') + call close_file (ioun) + + ! provide for namelist over-ride of default values + ioun = open_namelist_file() + read (ioun,ocean_sponges_tracer_OFAM_nml,IOSTAT=io_status) + write (stdlogunit,ocean_sponges_tracer_OFAM_nml) + write (stdoutunit,'(/)') + write (stdoutunit,ocean_sponges_tracer_OFAM_nml) + ierr = check_nml_error(io_status,'ocean_sponges_tracer_OFAM_nml') + call close_file (ioun) + + Dom => Domain + Grd => Grid + +#ifndef MOM4_STATIC_ARRAYS + call get_local_indices(Domain, isd, ied, jsd, jed, isc, iec, jsc, jec) + nk = Grid%nk #endif - do n=1,num_prog_tracers - Sponge(n)%id = -1 - enddo - - dtimer = 1.0/dtime - - if(use_this_module) then - write(stdoutunit,*)'==>Note from ocean_sponges_tracer_mod: Using this module.' - Ocean_options%ocean_sponges_tracer= 'Used ocean tracer sponges.' - else - write(stdoutunit,*)'==>Note from ocean_sponges_tracer_mod: NOT using ocean tracer sponges.' - Ocean_options%ocean_sponges_tracer= 'Did NOT use ocean tracer sponges.' - return - endif - - if (use_adaptive_restore) then + do n=1,num_prog_tracers + Sponge(n)%id = -1 + enddo allocate(sdiffo(num_prog_tracers)) - ! Set up some initial times - call get_time(Time%model_time, secs, days) - initial_day = days - initial_secs = secs - end if + dtimer = 1.0/dtime - do n = 1, num_prog_tracers - if (use_adaptive_restore) then - allocate(Sponge(n)%damp_coeff(isd:ied,jsd:jed,nk)) - Sponge(n)%damp_coeff=0.0 + if(use_this_module) then + write(stdoutunit,*)'==>Note from ocean_sponges_tracer_mod: Using this module.' + Ocean_options%ocean_sponges_tracer= 'Used ocean tracer sponges.' else - ! read damping rates (1/sec) - - name = 'INPUT/'//trim(T_prog(n)%name)//'_sponge_coeff.nc' - if (file_exist(trim(name))) then - write(stdoutunit,*) '==> Using sponge damping times specified from file '//trim(name) - allocate(Sponge(n)%damp_coeff(isd:ied,jsd:jed,nk)) - allocate(Sponge(n)%damp_coeff2d(isd:ied,jsd:jed)) - Sponge(n)%damp_coeff(:,:,:) = 0.0 - Sponge(n)%damp_coeff2d(:,:) = 0.0 - - if(damp_coeff_3d) then - call read_data(name,'coeff',Sponge(n)%damp_coeff,domain=Domain%domain2d,timelevel=1) - else - call read_data(name,'coeff',Sponge(n)%damp_coeff2d,domain=Domain%domain2d,timelevel=1) - do k=1,nk - do j=jsc,jec - do i=isc,iec - Sponge(n)%damp_coeff(i,j,k) = Sponge(n)%damp_coeff2d(i,j) - enddo - enddo - enddo - endif - - do k=1,nk - do j=jsc,jec - do i=isc,iec - if(Grd%tmask(i,j,k) == 0.0) then - Sponge(n)%damp_coeff(i,j,k) = 0.0 - endif - enddo - enddo - enddo - - ! modify damping rates to allow restoring to be solved implicitly - ! note: test values between zero and 4.0e-8 revert to damping rates defined above - - do k=1,nk - do j=jsc,jec - do i=isc,iec - if (dtime*Sponge(n)%damp_coeff(i,j,k) > 4.0e-8) then - if (dtime*Sponge(n)%damp_coeff(i,j,k) > 37.0) then - Sponge(n)%damp_coeff(i,j,k) = dtimer - else - Sponge(n)%damp_coeff(i,j,k) = (1.0 - exp(-dtime*Sponge(n)%damp_coeff(i,j,k))) * dtimer - endif - else if (dtime*Sponge(n)%damp_coeff(i,j,k) <= 0.0) then - Sponge(n)%damp_coeff(i,j,k) = 0.0 - endif - enddo - enddo - enddo - else - write(stdoutunit,*) '==> '//trim(name)//' not found. Sponge not being applied ' - cycle - endif + write(stdoutunit,*)'==>Note from ocean_sponges_tracer_mod: NOT using ocean tracer sponges.' + Ocean_options%ocean_sponges_tracer= 'Did NOT use ocean tracer sponges.' + return endif - ! Read restoring data + !set up some initial times + call get_time(Time%model_time,secs,days) + initial_day = days + initial_secs = secs - name = 'INPUT/'//trim(T_prog(n)%name)//'_sponge.nc' - Sponge(n)%id = init_external_field(name,T_prog(n)%name,domain=Domain%domain2d) - if (Sponge(n)%id < 1) then - call mpp_error(FATAL,& - '==>Error: in ocean_sponges_tracer_mod: damping rates are specified but sponge values are not') + if ( limit_temp ) then + write(stdoutunit,*) '==> Preventing freezing by restoring to -1.8C where applicable' + endif + if ( limit_salt ) then + write(stdoutunit,*) '==> Preventing excessive freshening by restoring salinity to 0.1 where applicable' + write(stdoutunit,*) '==> Preventing excessive saltiness by restoring salinity to 42.5 where applicable' endif - write(stdoutunit,*) '==> Using sponge data specified from file '//trim(name) - enddo - - ! register diagnostic output - allocate (id_sponge_tend(num_prog_tracers)) - id_sponge_tend = -1 - - do n=1,num_prog_tracers - if(n==index_temp) then - id_sponge_tend(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_sponge_tend',& - Grd%tracer_axes(1:3), Time%model_time, 'rho*dzt*cp*heating due to sponge', & - trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e10,1.e10/)) - else - id_sponge_tend(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_sponge_tend',& - Grd%tracer_axes(1:3), Time%model_time, 'rho*dzt*tendency due to sponge', & - trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e10,1.e10/)) - endif - enddo - - -end subroutine ocean_sponges_tracer_init -! NAME="ocean_sponges_tracer_init" - - -!####################################################################### -! -! -! -! This subroutine calculates thickness weighted and density weighted -! time tendencies of tracers due to damping by sponges. -! -! -subroutine sponge_tracer_source(Time, Thickness, T_prog) - - type(ocean_time_type), intent(in) :: Time - type(ocean_thickness_type), intent(in) :: Thickness - type(ocean_prog_tracer_type), intent(inout) :: T_prog(:) - - integer :: taum1, tau - integer :: i, j, k, n - - real :: sdiff, sum_val, adaptive_coeff - integer :: secs, days, numsecs - logical :: do_adaptive_restore = .false. ! Only restore in specified time period. - logical :: do_normal_restore = .false. ! Only restore in specified time period. - logical, save :: first_pass = .true. - - if(.not. use_this_module) return - - taum1 = Time%taum1 - tau = Time%tau - wrk1 = 0.0 - - if (use_adaptive_restore) then - secs_restore = days_to_restore * 86400 + secs_to_restore - sdiff = 0.0 - call get_time(Time%model_time, secs, days) - numsecs = (days - initial_day) * 86400 + secs - initial_secs - - do_adaptive_restore = (numsecs < secs_restore) - do_normal_restore = (.not. do_adaptive_restore) .and. use_sponge_after_init - else - do_normal_restore = .true. - endif - - do n=1,size(T_prog(:)) - - ! Need to reinitialise wrk2 here due to limiting of temperature. - wrk2 = 0.0 - - if (Sponge(n)%id > 0) then - - ! get sponge value for current time - call time_interp_external(Sponge(n)%id, Time%model_time, wrk1) - if (do_adaptive_restore) then - if (first_pass) then - if (use_hard_thump) then - do k = 1, nk - do j = jsd, jed - do i = isd, ied - T_prog(n)%field(i,j,k,taum1) = wrk1(i,j,k) - T_prog(n)%field(i,j,k,tau) = wrk1(i,j,k) - end do - end do - end do - end if - wrk2 = abs(T_prog(n)%field(:,:,:,taum1) - wrk1(:,:,:)) & - * Grd%tmask(:,:,:) - sum_val = sum(wrk2(isc:iec,jsc:jec,:)) - call mpp_sum(sum_val) - sdiffo(n) = sum_val / Grd%wet_t_points - if (.not. use_normalising) then - sdiffo(n) = 1.0 - else - sdiffo(n) = maxval(wrk2) + + do n=1,num_prog_tracers + if ( restore_cars_deep ) then + name = 'INPUT/'//trim(T_prog(n)%name)//'_cars.nc' + if (file_exist(trim(name))) then + write(stdoutunit,*) '==> damping '//trim(T_prog(n)%name)//' at depth from file '//trim(name) + allocate(Cars(n)%field(isd:ied,jsd:jed,nk)) + call read_data(name,T_prog(n)%name,Cars(n)%field,domain=Domain%domain2d,timelevel=1) + Cars(n)%id = 1 + else + Cars(n)%id = 0 + endif + endif + + ! read damping rates (1/sec) + if (.not. use_hard_thump .and. .not. use_adaptive_restore) then + allocate(Sponge(n)%damp_coeff(isd:ied,jsd:jed,nk)) + Sponge(n)%damp_coeff=0.0 + name = 'INPUT/'//trim(T_prog(n)%name)//'_sponge_coeff.nc' + if (file_exist(trim(name))) then + write(stdoutunit,*) '==> Using sponge damping times specified from file '//trim(name) + allocate(Sponge(n)%damp_coeff(isd:ied,jsd:jed,nk)) + allocate(Sponge(n)%damp_coeff2d(isd:ied,jsd:jed)) + Sponge(n)%damp_coeff(:,:,:) = 0.0 + Sponge(n)%damp_coeff2d(:,:) = 0.0 + if(damp_coeff_3d) then + call read_data(name,'coeff',Sponge(n)%damp_coeff,domain=Domain%domain2d,timelevel=1) + else + call read_data(name,'coeff',Sponge(n)%damp_coeff2d,domain=Domain%domain2d,timelevel=1) + do k=1,nk + do j=jsc,jec + do i=isc,iec + Sponge(n)%damp_coeff(i,j,k) = Sponge(n)%damp_coeff2d(i,j) + enddo + enddo + enddo + endif + do k=1,nk + do j=jsc,jec + do i=isc,iec + if(Grd%tmask(i,j,k) == 0.0) then + Sponge(n)%damp_coeff(i,j,k) = 0.0 + endif + enddo + enddo + enddo + ! modify damping rates to allow restoring to be solved implicitly + ! note: test values between zero and 4.0e-8 revert to damping rates defined above + do k=1,nk + do j=jsc,jec + do i=isc,iec + if (dtime*Sponge(n)%damp_coeff(i,j,k) > 4.0e-8) then + if (dtime*Sponge(n)%damp_coeff(i,j,k) > 37.0) then + Sponge(n)%damp_coeff(i,j,k) = dtimer + else + Sponge(n)%damp_coeff(i,j,k) = (1.0 - exp(-dtime*Sponge(n)%damp_coeff(i,j,k))) * dtimer + endif + else if (dtime*Sponge(n)%damp_coeff(i,j,k) <= 0.0) then + Sponge(n)%damp_coeff(i,j,k) = 0.0 + endif + enddo + enddo + enddo + else + write(stdoutunit,*) '==> '//trim(name)//' not found. Sponge damping coefficients not being applied ' + endif !(file_exist(trim(name))) + endif !.not. use_hard_thump .and. .not. use_adaptive_restore + + ! read restoring data + name = 'INPUT/'//trim(T_prog(n)%name)//'_sponge.nc' + if (file_exist(trim(name)) ) then + Sponge(n)%id = init_external_field(name,T_prog(n)%name,domain=Domain%domain2d) + if ( Sponge(n)%id < 1) then + call mpp_error(FATAL,& + '==>Error: in ocean_sponges_tracer_mod: Sponge values are required but not specified') + endif + write(stdoutunit,*) '==> Using restoring data specified from file '//trim(name) + else + write(stdoutunit,*) '==> '//trim(name)//' not found. Sponge not being applied ' + endif + enddo !n=1,num_prog_tracers + + ! register diagnostic output + allocate (id_sponge_tend(num_prog_tracers)) + id_sponge_tend = -1 + + do n=1,num_prog_tracers + if(n==index_temp) then + id_sponge_tend(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_sponge_tend',& + Grd%tracer_axes(1:3), Time%model_time, 'rho*dzt*cp*heating due to sponge', & + trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e6,1.e6/)) + else + id_sponge_tend(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_sponge_tend',& + Grd%tracer_axes(1:3), Time%model_time, 'rho*dzt*tendency due to sponge', & + trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e-1,1.e-1/)) + endif + enddo + + end subroutine ocean_sponges_tracer_init + ! NAME="ocean_sponges_tracer_init" + + !####################################################################### + ! + ! + ! + ! This subroutine calculates thickness weighted and density weighted + ! time tendencies of tracers due to damping by sponges. + ! + ! + subroutine sponge_tracer_source(Time, Thickness, T_prog) + + type(ocean_time_type), intent(in) :: Time + type(ocean_thickness_type), intent(in) :: Thickness + type(ocean_prog_tracer_type), intent(inout) :: T_prog(:) + + integer :: taum1, tau + integer :: i, j, k, n + + real :: sdiff, sum_val, adaptive_coeff + integer :: secs, days, numsecs + logical :: do_adaptive_restore = .false. ! Only restore in specified time period. + logical,save :: first_pass = .true. + real, parameter :: one_over_year = 1./(365.0*86400.0) + + if(.not. use_this_module) return + + taum1 = Time%taum1 + tau = Time%tau + wrk1 = 0.0 + + + sdiff=0.0 + call get_time(Time%model_time,secs,days) + + do n=1,size(T_prog(:)) + ! Need to reinitialise wrk2 here due to limiting of temperature. + wrk2 = 0.0 + !----------------------------------------------------------------------------------------------------------------- + if (Sponge(n)%id > 0) then + ! get sponge value for current time + call time_interp_external(Sponge(n)%id, Time%model_time, wrk1) + !-------------------------------------------------------------------------------------------------------------- + if ( first_pass ) then + if ( use_hard_thump ) then + do k = 1,nk + do j=jsd,jed + do i=isd,ied + T_prog(n)%field(i,j,k,taum1)=wrk1(i,j,k) + T_prog(n)%field(i,j,k,tau)=wrk1(i,j,k) + enddo + enddo + enddo + call mpp_update_domains ( T_prog(n)%field(:,:,:,tau),Dom%domain2d) + call mpp_update_domains ( T_prog(n)%field(:,:,:,taum1),Dom%domain2d) + wrk2=abs(T_prog(n)%field(:,:,:,taum1)-wrk1(:,:,:))*Grd%tmask(:,:,:) + endif !use_hard_thump + sum_val=sum(wrk2(isc:iec,jsc:jec,:)) + call mpp_sum(sum_val) + sdiffo(n) = sum_val/Grd%wet_t_points + write (stdout(),*)'mean absolute deviation ',trim(T_prog(n)%name),' = ',sdiffo(n) + if ( .not. use_normalising ) then + sdiffo(n)=1.0 + else + sdiffo(n) = maxval(wrk2) call mpp_max(sdiffo(n)) - sdiffo(n) = 1.01 * sdiffo(n) - endif - endif - - do k = 1, nk - do j = jsd, jed - do i = isd, ied - sdiff = abs(T_prog(n)%field(i,j,k,taum1) - wrk1(i,j,k)) - sdiff = max(sdiff, 1e-9) ! Minimum tolerance - adaptive_coeff = 1.0 / (taumin - & - (lambda * real(secs_restore)) & - / (((sdiff / sdiffo(n))**npower) & - * log(1 - athresh))) - wrk2(i,j,k) = Thickness%rho_dzt(i,j,k,tau) * adaptive_coeff & - * (wrk1(i,j,k) - T_prog(n)%field(i,j,k,taum1)) - enddo + sdiffo(n) = 1.01*sdiffo(n) + endif + endif !first_pass + !-------------------------------------------------------------------------------------------------------------- + if (use_hard_thump == .false.) then + if (use_adaptive_restore) then + secs_restore = days_to_restore*86400 + secs_to_restore + numsecs = (days-initial_day)*86400 + secs - initial_secs + do_adaptive_restore = ( numsecs < secs_restore ) + endif + if ( do_adaptive_restore ) then + do k=1,nk + do j=jsd,jed + do i=isd,ied + sdiff=abs(T_prog(n)%field(i,j,k,taum1)-wrk1(i,j,k)) + sdiff = max(sdiff,1e-9) !minimum tolerance + adaptive_coeff = 1.0/(taumin - (lambda*real(secs_restore))/(((sdiff/sdiffo(n))**npower)*log(1-athresh))) + wrk2(i,j,k) = Thickness%rho_dzt(i,j,k,tau)* adaptive_coeff & + *(wrk1(i,j,k)-T_prog(n)%field(i,j,k,taum1)) + enddo + enddo + enddo + else !normal restore + do k=1,nk + do j=jsd,jed + do i=isd,ied + wrk2(i,j,k) = Thickness%rho_dzt(i,j,k,tau)*Sponge(n)%damp_coeff(i,j,k) & + *(wrk1(i,j,k)-T_prog(n)%field(i,j,k,taum1)) + enddo + enddo + enddo + endif + endif !use_hard_thump + !-------------------------------------------------------------------------------------------------------------- + endif !Sponge(1)%id + !----------------------------------------------------------------------------------------------------------------- + if ( restore_cars_deep .and. Cars(n)%id == 1 ) then + do k=nk-4,nk + do j=jsd,jed + do i=isd,ied + wrk2(i,j,k) = Thickness%rho_dzt(i,j,k,tau)*one_over_year & + *(Cars(n)%field(i,j,k)-T_prog(n)%field(i,j,k,taum1)) + enddo enddo - enddo - else if (do_normal_restore) then - do k = 1, nk - do j = jsc, jec - do i = isc, iec - wrk2(i,j,k) = Thickness%rho_dzt(i,j,k,tau) & - * Sponge(n)%damp_coeff(i,j,k) & - * (wrk1(i,j,k) - T_prog(n)%field(i,j,k,taum1)) - end do - end do - end do - end if - - end if - - ! These tracer relaxation schemes were in the original OFAM and AusCOM - ! source code, and may be needed by other OFAM-based experiments. - ! - ! For now, let's keep them in the same place, but they are now turned off - ! by default. They may need to be moved or deleted in the future. - - ! Limit to -1.8 deg with 3hour restoring (got other limits here for initialising. Not sure if needed (namelist use?). - if (limit_temp .and. trim(T_prog(n)%name) == 'temp') then - do k=1,nk - do j=jsc,jec - do i=isc,iec - wrk2(i,j,k) = wrk2(i,j,k) + Thickness%rho_dzt(i,j,k,tau) & - * max(limit_temp_min - T_prog(n)%field(i,j,k,taum1), 0.0) & - / limit_temp_restore -! wrk2(i,j,k) = wrk2(i,j,k)+Thickness%rho_dzt(i,j,k,tau)*min(43.0-T_prog(n)%field(i,j,k,taum1),0.0)/3600.0 + enddo + endif + !-------------------------------------------------------------------------------------------------------------- + !this limiting scheme restores to sponge data rather than hard limits + if (trim(T_prog(n)%name) == 'temp' .and. limit_temp) then + do k=1,nk + do j=jsd,jed + do i=isd,ied + if ( Grd%tmask(i,j,k)==1.0 ) then + if ( T_prog(n)%field(i,j,k,taum1).ge.limit_temp_max ) then + sdiff=abs(limit_temp_max-T_prog(n)%field(i,j,k,taum1)) + sdiff = max(sdiff,1e-9) !minimum tolerance + adaptive_coeff = 1.0/(taumin-(lambda*real(86400))/(((sdiff/sdiffo(n))**npower)*log(1-athresh))) + wrk2(i,j,k) = Grd%tmask(i,j,k)*Thickness%rho_dzt(i,j,k,tau)*adaptive_coeff & + *(limit_temp_max-T_prog(n)%field(i,j,k,taum1)) + endif + if ( T_prog(n)%field(i,j,k,taum1).le.limit_temp_min ) then + sdiff=abs(limit_temp_min-T_prog(n)%field(i,j,k,taum1)) + sdiff = max(sdiff,1e-9) !minimum tolerance + adaptive_coeff = 1.0/(taumin-(lambda*real(86400))/(((sdiff/sdiffo(n))**npower)*log(1-athresh))) + wrk2(i,j,k) = Grd%tmask(i,j,k)*Thickness%rho_dzt(i,j,k,tau)*adaptive_coeff & + *(limit_temp_min-T_prog(n)%field(i,j,k,taum1)) + endif + endif + enddo enddo enddo - enddo - end if - - if (limit_salt .and. trim(T_prog(n)%name) == 'salt') then + endif !trim + !-------------------------------------------------------------------------------------------------------------- + if (trim(T_prog(n)%name) == 'salt' .and. limit_salt) then + do k=1,nk + do j=jsd,jed + do i=isd,ied + if ( Grd%tmask(i,j,k)==real(1) ) then + if ( T_prog(n)%field(i,j,k,taum1).ge.limit_salt_max ) then + sdiff=abs(limit_salt_max-T_prog(n)%field(i,j,k,taum1)) + sdiff = max(sdiff,1e-9) !minimum tolerance + adaptive_coeff = 1.0/(taumin-(lambda*real(86400))/(((sdiff/sdiffo(n))**npower)*log(1-athresh))) + wrk2(i,j,k) = Grd%tmask(i,j,k)*Thickness%rho_dzt(i,j,k,tau)*adaptive_coeff & + *(limit_salt_max-T_prog(n)%field(i,j,k,taum1)) + endif + if ( T_prog(n)%field(i,j,k,taum1).le.limit_salt_min ) then + sdiff=abs(limit_salt_min-T_prog(n)%field(i,j,k,taum1)) + sdiff = max(sdiff,1e-9) !minimum tolerance + adaptive_coeff = 1.0/(taumin-(lambda*real(86400))/(((sdiff/sdiffo(n))**npower)*log(1-athresh))) + wrk2(i,j,k) = Grd%tmask(i,j,k)*Thickness%rho_dzt(i,j,k,tau)*adaptive_coeff & + *(limit_salt_min-T_prog(n)%field(i,j,k,taum1)) + endif + endif + enddo + enddo + enddo + endif !trim + !-------------------------------------------------------------------------------------------------------------- do k=1,nk do j=jsc,jec - do i=isc,iec - wrk2(i,j,k) = wrk2(i,j,k) + Thickness%rho_dzt(i,j,k,tau) & - * max(limit_salt_min - T_prog(n)%field(i,j,k,taum1), 0.0) & - / limit_salt_restore + do i=isc,iec + T_prog(n)%th_tendency(i,j,k) = T_prog(n)%th_tendency(i,j,k) + wrk2(i,j,k) enddo enddo enddo - end if - - ! Update tendency - do k = 1, nk - do j = jsc, jec - do i = isc, iec - T_prog(n)%th_tendency(i,j,k) = T_prog(n)%th_tendency(i,j,k) & - + wrk2(i,j,k) - end do - end do - end do - - if (id_sponge_tend(n) > 0) call diagnose_3d(Time, Grd, id_sponge_tend(n), & - T_prog(n)%conversion*wrk2(:,:,:)) - - enddo - - first_pass = .false. - - return + if (id_sponge_tend(n) > 0) used = send_data(id_sponge_tend(n), & + T_prog(n)%conversion*wrk2(:,:,:), Time%model_time, rmask=Grd%tmask(:,:,:), & + is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) + !-------------------------------------------------------------------------------------------------------------- + enddo !do n=1,size(T_prog(:)) -end subroutine sponge_tracer_source -! NAME="sponge_tracer_source" + first_pass = .false. + return + end subroutine sponge_tracer_source + ! NAME="sponge_tracer_source" end module ocean_sponges_tracer_mod diff --git a/src/mom5/ocean_param/sources/ocean_sponges_velocity.F90 b/src/mom5/ocean_param/sources/ocean_sponges_velocity.F90 index 3b344ee693..c72bf5075b 100644 --- a/src/mom5/ocean_param/sources/ocean_sponges_velocity.F90 +++ b/src/mom5/ocean_param/sources/ocean_sponges_velocity.F90 @@ -1,564 +1,603 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! !! +!! 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_sponges_velocity_mod -! -! Paul Sandery -! -! -! -! Thickness weighted velocity tendency [meter*meter/sec*sec] from sponges. -! -! -! -! -! This module applies sponges to currents. The sponges -! can occur at any location and with any distribution in the domain, and -! with any time step and damping rate. Sponges occur where positive -! inverse restore times occur in the field passed to sponge_init. An -! array of tracer tendencies due to the sponges is augmented through a -! call to sponge_tracer_source. The array of tracer tendencies must be -! reset to zero between calls. -! -! Different damping rates can be specified for each field by making -! calls to register_sponge_rate - no sponges are applied to fields for -! which uniformly zero inverse damping rates are set with a call to -! register_sponge_rate. The value towards which a field is damped is -! set with calls to register_sponge_field; successive calls are used to -! set up linear interpolation of this restore rate. -! -! Sponge data and damping coefficients are generally 3 dimensional. -! -! The user is responsible for providing (and registering) the data on -! the model grid of values towards which the currents are being driven. -! -! -! -! -! -! -! For using this module. Default use_this_module=.false. -! -! -! -! For case when damping coefficients are full 3d field of values. -! Default damp_coeff_3d=.false., which means damping coeffs are -! 2d horizontal array. -! -! -! -! -use diag_manager_mod, only: register_diag_field -use fms_mod, only: write_version_number, open_namelist_file, close_file -use fms_mod, only: file_exist -use fms_mod, only: open_namelist_file, check_nml_error, close_file -use fms_mod, only: read_data, lowercase, FATAL, WARNING, stdout, stdlog -use mpp_mod, only: input_nml_file, mpp_sum, mpp_error, mpp_max -use time_interp_external_mod, only: init_external_field, time_interp_external -use time_manager_mod, only: time_type, set_date, get_time -use time_manager_mod, only: operator( + ), operator( - ), operator( // ) -use time_manager_mod, only: operator( > ), operator( == ), operator( <= ) - -use ocean_domains_mod, only: get_local_indices -use ocean_parameters_mod, only: missing_value, rho0 -use ocean_types_mod, only: ocean_domain_type, ocean_grid_type, ocean_thickness_type -use ocean_types_mod, only: ocean_time_type, ocean_velocity_type, ocean_options_type -use ocean_workspace_mod, only: wrk1, wrk2 -use ocean_util_mod, only: diagnose_3d_u - -implicit none - -private + ! + ! Paul Sandery + ! + ! + ! + ! Thickness weighted velocity tendency [meter*meter/sec*sec] from sponges. + ! + ! + ! + ! + ! This module applies sponges to currents. The sponges + ! can occur at any location and with any distribution in the domain, and + ! with any time step and damping rate. Sponges occur where positive + ! inverse restore times occur in the field passed to sponge_init. An + ! array of tracer tendencies due to the sponges is augmented through a + ! call to sponge_tracer_source. The array of tracer tendencies must be + ! reset to zero between calls. + ! + ! Different damping rates can be specified for each field by making + ! calls to register_sponge_rate - no sponges are applied to fields for + ! which uniformly zero inverse damping rates are set with a call to + ! register_sponge_rate. The value towards which a field is damped is + ! set with calls to register_sponge_field; successive calls are used to + ! set up linear interpolation of this restore rate. + ! + ! Sponge data and damping coefficients are generally 3 dimensional. + ! + ! The user is responsible for providing (and registering) the data on + ! the model grid of values towards which the currents are being driven. + ! + ! We now allow for the possibility of restoring adaptively according to + ! Sandery et al. (2011). The user may specify parameters as defined in + ! in that paper. + ! + ! + ! + ! + ! + ! + ! For using this module. Default use_this_module=.false. + ! + ! + ! + ! For case when damping coefficients are full 3d field of values. + ! Default damp_coeff_3d=.false., which means damping coeffs are + ! 2d horizontal array. + ! + ! + ! + ! + use diag_manager_mod, only: register_diag_field, send_data + use fms_mod, only: write_version_number, open_namelist_file, close_file + use fms_mod, only: file_exist + use fms_mod, only: open_namelist_file, check_nml_error, close_file + use fms_mod, only: read_data, lowercase, FATAL, WARNING, stdout, stdlog + use mpp_mod, only: mpp_sum, mpp_error, mpp_max + use time_interp_external_mod, only: init_external_field, time_interp_external + use time_manager_mod, only: time_type, set_date, get_time + use time_manager_mod, only: operator( + ), operator( - ), operator( // ) + use time_manager_mod, only: operator( > ), operator( == ), operator( <= ) + use mpp_domains_mod, only: mpp_update_domains + use ocean_domains_mod, only: get_local_indices + use ocean_parameters_mod, only: missing_value, rho0 + use ocean_types_mod, only: ocean_domain_type, ocean_grid_type, ocean_thickness_type + use ocean_types_mod, only: ocean_time_type, ocean_velocity_type, ocean_options_type + use ocean_workspace_mod, only: wrk1, wrk2 + + implicit none + + private #include -type ocean_sponge_type - integer :: id ! time_interp_external index - character(len=32) :: name ! tracer name corresponding to sponge - real, dimension(:,:), pointer :: damp_coeff_u2d => NULL() ! 3d inverse damping rate (tracer units/ sec) - real, dimension(:,:), pointer :: damp_coeff_v2d => NULL() ! 2d inverse damping rate (tracer units/ sec) - real, dimension(:,:,:), pointer :: damp_coeff_u => NULL() ! 3d inverse damping rate for u current (m/ sec^2) - real, dimension(:,:,:), pointer :: damp_coeff_v => NULL() ! 3d inverse damping rate for v current (m/ sec^2) -end type ocean_sponge_type - -type(ocean_sponge_type), allocatable, dimension(:) :: Sponge_u -type(ocean_sponge_type), allocatable, dimension(:) :: Sponge_v - -type(ocean_domain_type), pointer :: Dom => NULL() -type(ocean_grid_type), pointer :: Grd => NULL() - - -public ocean_sponges_velocity_init -public sponge_velocity_source - -character(len=126) :: version = '$Id: ocean_sponges_velocity.F90,v 20.0 2013/12/14 00:16:26 fms Exp $' -character (len=128) :: tagname = '$Name: tikal $' - -! for diagnostics -logical :: used -integer, dimension(:), allocatable :: id_sponge_tend -logical :: module_is_initialized = .false. -logical :: damp_coeff_3d = .false. -logical :: use_this_module = .false. - -! Adaptive restoring -real :: athresh = 0.5 -real :: npower = 1.0 -real :: sdiffo_u = 0.5 -real :: sdiffo_v = 0.5 -real :: lambda = 0.0083 -real :: taumin = 720 -logical :: use_adaptive_restore = .false. -logical :: use_sponge_after_init = .false. -logical :: use_normalising = .false. -logical :: use_hard_thump = .false. -integer :: secs_to_restore = 0 -integer :: days_to_restore = 1 - -integer :: secs_restore -integer :: days_end_restore -integer :: secs_end_restore -integer :: initial_day, initial_secs - -namelist /ocean_sponges_velocity_nml/ use_this_module, damp_coeff_3d -namelist /ocean_sponges_velocity_ofam_nml/ & - use_adaptive_restore, use_sponge_after_init, use_normalising, & - use_hard_thump, athresh, taumin, lambda, npower, days_to_restore, & - secs_to_restore + type ocean_sponge_type + integer :: id ! time_interp_external index + character(len=32) :: name ! tracer name corresponding to sponge + real, dimension(:,:), pointer :: damp_coeff_u2d => NULL() ! 3d inverse damping rate (tracer units/ sec) + real, dimension(:,:), pointer :: damp_coeff_v2d => NULL() ! 2d inverse damping rate (tracer units/ sec) + real, dimension(:,:,:), pointer :: damp_coeff_u => NULL() ! 3d inverse damping rate for u current (m/ sec^2) + real, dimension(:,:,:), pointer :: damp_coeff_v => NULL() ! 3d inverse damping rate for v current (m/ sec^2) + end type ocean_sponge_type + + type(ocean_sponge_type), allocatable, dimension(:) :: Sponge_u + type(ocean_sponge_type), allocatable, dimension(:) :: Sponge_v + + type(ocean_domain_type), pointer :: Dom => NULL() + type(ocean_grid_type), pointer :: Grd => NULL() + + public ocean_sponges_velocity_init + public sponge_velocity_source + + character(len=126) :: version = '$Id: ocean_sponges_velocity.F90,v 1.1.2.3.42.1 2009/10/10 00:42:54 nnz Exp $' + character (len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' + + ! for diagnostics + logical :: used + integer, dimension(:), allocatable :: id_sponge_tend + logical :: module_is_initialized = .false. + logical :: damp_coeff_3d = .false. + logical :: use_this_module = .false. + + ! Adaptive restoring + + real :: athresh = 0.5 + real :: npower = 1.0 + real :: sdiffo_u = 0.5 + real :: sdiffo_v = 0.5 + real :: lambda = 0.0083 + real :: taumin = 3600 + logical :: use_adaptive_restore = .false. + logical :: use_sponge_after_init = .false. + logical :: use_normalising = .false. + logical :: use_hard_thump = .false. + logical :: use_increment = .false. + integer :: secs_to_restore = 0 + integer :: days_to_restore = 0 + integer :: secs_restore + integer :: days_end_restore + integer :: secs_end_restore + integer :: initial_day, initial_secs + logical :: limit_u = .false. + logical :: limit_v = .false. + real :: limit_u_max = 3.0 + real :: limit_u_min = -3.0 + real :: limit_v_max = 3.0 + real :: limit_v_min = -3.0 + + namelist /ocean_sponges_velocity_nml/ use_this_module, damp_coeff_3d + namelist /ocean_sponges_velocity_OFAM_nml/ use_adaptive_restore,use_sponge_after_init, use_normalising,use_hard_thump,& + athresh, taumin,lambda,npower,days_to_restore, secs_to_restore, limit_u, limit_v contains -!####################################################################### -! -! -! -! This subroutine is intended to be used to initialize the sponges. -! Everything in this subroutine is a user prototype, and should be replacable. -! -! -subroutine ocean_sponges_velocity_init(Grid, Domain, Time, dtime, Ocean_options) - - type(ocean_grid_type), intent(in), target :: Grid - type(ocean_domain_type), intent(in), target :: Domain - type(ocean_time_type), intent(in) :: Time - real, intent(in) :: dtime - type(ocean_options_type), intent(inout) :: Ocean_options - - integer :: i, j, k - integer :: ioun, io_status, ierr - integer :: secs, days - real :: dtimer - - character(len=128) :: name - - integer :: stdoutunit,stdlogunit - stdoutunit=stdout();stdlogunit=stdlog() - - if ( module_is_initialized ) then - call mpp_error(FATAL, & - '==>Error in ocean_sponges_velocity_mod (ocean_sponges_velocity_init): module already initialized') - endif - - module_is_initialized = .TRUE. - - allocate( Sponge_u(1) ) - allocate( Sponge_v(1) ) - - call write_version_number(version, tagname) - - ! provide for namelist over-ride of default values -#ifdef INTERNAL_FILE_NML - read(input_nml_file, nml=ocean_sponges_velocity_nml, iostat=io_status) - ierr = check_nml_error(io_status, 'ocean_sponges_velocity_nml') - read(input_nml_file, nml=ocean_sponges_velocity_ofam_nml, iostat=io_status) - ierr = check_nml_error(io_status, 'ocean_sponges_velocity_ofam_nml') -#else - ioun = open_namelist_file() - read(ioun, ocean_sponges_velocity_nml, IOSTAT=io_status) - ierr = check_nml_error(io_status,'ocean_sponges_velocity_nml') - rewind(ioun) - read(ioun, ocean_sponges_velocity_ofam_nml, IOSTAT=io_status) - ierr = check_nml_error(io_status,'ocean_sponges_velocity_ofam_nml') - call close_file(ioun) + !####################################################################### + ! + ! + ! + ! This subroutine is intended to be used to initialize the sponges. + ! Everything in this subroutine is a user prototype, and should be replacable. + ! + ! + subroutine ocean_sponges_velocity_init(Grid, Domain, Time, dtime, Ocean_options) + + type(ocean_grid_type), intent(in), target :: Grid + type(ocean_domain_type), intent(in), target :: Domain + type(ocean_time_type), intent(in) :: Time + real, intent(in) :: dtime + type(ocean_options_type), intent(inout) :: Ocean_options + + integer :: i, j, k, n + integer :: ioun, io_status, ierr + integer :: secs, days + real :: dtimer + + character(len=128) :: name + + integer :: stdoutunit,stdlogunit + stdoutunit=stdout();stdlogunit=stdlog() + + if ( module_is_initialized ) then + call mpp_error(FATAL, & + '==>Error in ocean_sponges_velocity_mod (ocean_sponges_velocity_init): module already initialized') + endif + + module_is_initialized = .TRUE. + + allocate( Sponge_u(1) ) + allocate( Sponge_v(1) ) + + call write_version_number() + + ! provide for namelist over-ride of default values + ioun = open_namelist_file() + read (ioun,ocean_sponges_velocity_nml,IOSTAT=io_status) + write (stdlogunit,ocean_sponges_velocity_nml) + write (stdoutunit,'(/)') + write (stdoutunit,ocean_sponges_velocity_nml) + ierr = check_nml_error(io_status,'ocean_sponges_velocity_nml') + call close_file (ioun) + + ! provide for namelist over-ride of default values + ioun = open_namelist_file() + read (ioun,ocean_sponges_velocity_OFAM_nml,IOSTAT=io_status) + write (stdlogunit,ocean_sponges_velocity_OFAM_nml) + write (stdoutunit,'(/)') + write (stdoutunit,ocean_sponges_velocity_OFAM_nml) + ierr = check_nml_error(io_status,'ocean_sponges_velocity_OFAM_nml') + call close_file (ioun) + + Dom => Domain + Grd => Grid + +#ifndef MOM4_STATIC_ARRAYS + call get_local_indices(Domain, isd, ied, jsd, jed, isc, iec, jsc, jec) + nk = Grid%nk #endif - write(stdlogunit, ocean_sponges_velocity_nml) - write(stdoutunit, '(/)') - write(stdoutunit, ocean_sponges_velocity_nml) - write(stdlogunit, ocean_sponges_velocity_ofam_nml) - write(stdoutunit, '(/)') - write(stdoutunit, ocean_sponges_velocity_ofam_nml) - - Dom => Domain - Grd => Grid - -#ifndef MOM_STATIC_ARRAYS - call get_local_indices(Domain, isd, ied, jsd, jed, isc, iec, jsc, jec) - nk = Grid%nk -#endif - - Sponge_u(1)%id = -1 - Sponge_v(1)%id = -1 - dtimer = 1.0/dtime - - if(use_this_module) then - write(stdoutunit,*)'==>Note from ocean_sponges_velocity_mod: Using this module.' - Ocean_options%ocean_sponges_velocity= 'Used ocean velocity sponges.' - else - write(stdoutunit,*)'==>Note from ocean_sponges_velocity_mod: NOT using this module: no velocity sponges.' - Ocean_options%ocean_sponges_velocity= 'Did NOT use ocean velocity sponges.' - return - endif - - if (use_adaptive_restore) then - ! Set up some initial times - call get_time(Time%model_time, secs, days) - days_end_restore = days + days_to_restore - secs_end_restore = secs + secs_to_restore + Sponge_u(1)%id = -1 + Sponge_v(1)%id = -1 + dtimer = 1.0/dtime + + if(use_this_module) then + write(stdoutunit,*)'==>Note from ocean_sponges_velocity_mod: Using this module.' + Ocean_options%ocean_sponges_velocity= 'Used ocean velocity sponges.' + else + write(stdoutunit,*)'==>Note from ocean_sponges_velocity_mod: NOT using this module: no velocity sponges.' + Ocean_options%ocean_sponges_velocity= 'Did NOT use ocean velocity sponges.' + return + endif + + !set up some initial times + call get_time(Time%model_time,secs,days) + days_end_restore=days + days_to_restore + secs_end_restore=secs + secs_to_restore initial_day = days initial_secs = secs - allocate(Sponge_u(1)%damp_coeff_u(isd:ied, jsd:jed, nk)) - Sponge_u(1)%damp_coeff_u(:,:,:) = 0.0 - - allocate(Sponge_v(1)%damp_coeff_v(isd:ied,jsd:jed,nk)) - Sponge_v(1)%damp_coeff_v(:,:,:) = 0.0 - else - ! Read damping rates for u - name = 'INPUT/u_sponge_coeff.nc' - if (file_exist(name)) then - write(stdoutunit,*) '==> Using sponge damping times specified from file ',name - allocate(Sponge_u(1)%damp_coeff_u(isd:ied,jsd:jed,nk)) - allocate(Sponge_u(1)%damp_coeff_u2d(isd:ied,jsd:jed)) - Sponge_u(1)%damp_coeff_u(:,:,:) = 0.0 - Sponge_u(1)%damp_coeff_u2d(:,:) = 0.0 - - if(damp_coeff_3d) then - call read_data(name,'coeff',Sponge_u(1)%damp_coeff_u,domain=Domain%domain2d,timelevel=1) - else - call read_data(name,'coeff',Sponge_u(1)%damp_coeff_u2d,domain=Domain%domain2d,timelevel=1) - do k=1,nk - do j=jsc,jec - do i=isc,iec - Sponge_u(1)%damp_coeff_u(i,j,k) = Sponge_u(1)%damp_coeff_u2d(i,j) - enddo - enddo - enddo - endif - - do k=1,nk - do j=jsc,jec - do i=isc,iec - if(Grd%umask(i,j,k) == 0.0) then - Sponge_u(1)%damp_coeff_u(i,j,k) = 0.0 - endif - enddo - enddo - enddo - - - ! modify damping rates to allow restoring to be solved implicitly - ! note: test values between zero and 4.0e-8 revert to damping rates defined above - do k=1,nk - do j=jsc,jec - do i=isc,iec - if (dtime*Sponge_u(1)%damp_coeff_u(i,j,k) > 4.0e-8) then - if (dtime*Sponge_u(1)%damp_coeff_u(i,j,k) > 37.0) then + !read damping rates for u + if (.not. use_hard_thump .and. .not. use_adaptive_restore) then + name = 'INPUT/u_sponge_coeff.nc' + if (file_exist(name)) then + write(stdoutunit,*) '==> Using sponge damping times specified from file ',name + allocate(Sponge_u(1)%damp_coeff_u(isd:ied,jsd:jed,nk)) + allocate(Sponge_u(1)%damp_coeff_u2d(isd:ied,jsd:jed)) + Sponge_u(1)%damp_coeff_u(:,:,:) = 0.0 + Sponge_u(1)%damp_coeff_u2d(:,:) = 0.0 + if(damp_coeff_3d) then + call read_data(name,'coeff',Sponge_u(1)%damp_coeff_u,domain=Domain%domain2d,timelevel=1) + else + call read_data(name,'coeff',Sponge_u(1)%damp_coeff_u2d,domain=Domain%domain2d,timelevel=1) + do k=1,nk + do j=jsc,jec + do i=isc,iec + Sponge_u(1)%damp_coeff_u(i,j,k) = Sponge_u(1)%damp_coeff_u2d(i,j) + enddo + enddo + enddo + endif !damp_coeff_3d + do k=1,nk + do j=jsc,jec + do i=isc,iec + if(Grd%umask(i,j,k) == 0.0) then + Sponge_u(1)%damp_coeff_u(i,j,k) = 0.0 + endif + enddo + enddo + enddo + ! modify damping rates to allow restoring to be solved implicitly + ! note: test values between zero and 4.0e-8 revert to damping rates defined above + do k=1,nk + do j=jsc,jec + do i=isc,iec + if (dtime*Sponge_u(1)%damp_coeff_u(i,j,k) > 4.0e-8) then + if (dtime*Sponge_u(1)%damp_coeff_u(i,j,k) > 37.0) then Sponge_u(1)%damp_coeff_u(i,j,k) = dtimer - else + else Sponge_u(1)%damp_coeff_u(i,j,k) = (1.0 - exp(-dtime*Sponge_u(1)%damp_coeff_u(i,j,k))) * dtimer - endif - else if (dtime*Sponge_u(1)%damp_coeff_u(i,j,k) <= 0.0) then - Sponge_u(1)%damp_coeff_u(i,j,k) = 0.0 - endif - enddo - enddo - enddo - endif ! endif for if fileexist INPUT/u_sponge_coeff.nc - - ! Read damping rates for v-sponge - name = 'INPUT/v_sponge_coeff.nc' - if (file_exist(name)) then - write(stdoutunit,*) '==> Using sponge damping times specified from file ',name - allocate(Sponge_v(1)%damp_coeff_v(isd:ied,jsd:jed,nk)) - allocate(Sponge_v(1)%damp_coeff_v2d(isd:ied,jsd:jed)) - Sponge_v(1)%damp_coeff_v(:,:,:) = 0.0 - Sponge_v(1)%damp_coeff_v2d(:,:) = 0.0 - - if(damp_coeff_3d) then - call read_data(name,'coeff',Sponge_v(1)%damp_coeff_v,domain=Domain%domain2d,timelevel=1) - else - call read_data(name,'coeff',Sponge_v(1)%damp_coeff_v2d,domain=Domain%domain2d,timelevel=1) - do k=1,nk - do j=jsc,jec - do i=isc,iec - Sponge_v(1)%damp_coeff_v(i,j,k) = Sponge_v(1)%damp_coeff_v2d(i,j) - enddo - enddo - enddo - endif - - do k=1,nk - do j=jsc,jec - do i=isc,iec - if(Grd%umask(i,j,k) == 0.0) then - Sponge_v(1)%damp_coeff_v(i,j,k) = 0.0 - endif - enddo - enddo - enddo - - ! modify damping rates to allow restoring to be solved implicitly - ! note: test values between zero and 4.0e-8 revert to damping rates defined above - do k=1,nk - do j=jsc,jec - do i=isc,iec - if (dtime*Sponge_v(1)%damp_coeff_v(i,j,k) > 4.0e-8) then - if (dtime*Sponge_v(1)%damp_coeff_v(i,j,k) > 37.0) then + endif + else if (dtime*Sponge_u(1)%damp_coeff_u(i,j,k) <= 0.0) then + Sponge_u(1)%damp_coeff_u(i,j,k) = 0.0 + endif + enddo + enddo + enddo + endif !if fileexist INPUT/u_sponge_coeff.nc + endif !.not. use_hard_thump .and. .not. use_adaptive_restore + + !read damping rates for v-sponge + if (.not. use_hard_thump .and. .not. use_adaptive_restore) then + allocate(Sponge_v(1)%damp_coeff_v(isd:ied,jsd:jed,nk)) + Sponge_v(1)%damp_coeff_v(:,:,:) = 0.0 + name = 'INPUT/v_sponge_coeff.nc' + if (file_exist(name)) then + write(stdoutunit,*) '==> Using sponge damping times specified from file ',name + allocate(Sponge_v(1)%damp_coeff_v(isd:ied,jsd:jed,nk)) + allocate(Sponge_v(1)%damp_coeff_v2d(isd:ied,jsd:jed)) + Sponge_v(1)%damp_coeff_v(:,:,:) = 0.0 + Sponge_v(1)%damp_coeff_v2d(:,:) = 0.0 + if(damp_coeff_3d) then + call read_data(name,'coeff',Sponge_v(1)%damp_coeff_v,domain=Domain%domain2d,timelevel=1) + else + call read_data(name,'coeff',Sponge_v(1)%damp_coeff_v2d,domain=Domain%domain2d,timelevel=1) + do k=1,nk + do j=jsc,jec + do i=isc,iec + Sponge_v(1)%damp_coeff_v(i,j,k) = Sponge_v(1)%damp_coeff_v2d(i,j) + enddo + enddo + enddo + endif! damp_coeff_3d + do k=1,nk + do j=jsc,jec + do i=isc,iec + if(Grd%umask(i,j,k) == 0.0) then + Sponge_v(1)%damp_coeff_v(i,j,k) = 0.0 + endif + enddo + enddo + enddo + ! modify damping rates to allow restoring to be solved implicitly + ! note: test values between zero and 4.0e-8 revert to damping rates defined above + do k=1,nk + do j=jsc,jec + do i=isc,iec + if (dtime*Sponge_v(1)%damp_coeff_v(i,j,k) > 4.0e-8) then + if (dtime*Sponge_v(1)%damp_coeff_v(i,j,k) > 37.0) then Sponge_v(1)%damp_coeff_v(i,j,k) = dtimer - else + else Sponge_v(1)%damp_coeff_v(i,j,k) = (1.0 - exp(-dtime*Sponge_v(1)%damp_coeff_v(i,j,k))) * dtimer - endif - else if (dtime*Sponge_v(1)%damp_coeff_v(i,j,k) <= 0.0) then - Sponge_v(1)%damp_coeff_v(i,j,k) = 0.0 - endif - enddo - enddo - enddo - - endif ! endif for if fileexist INPUT/v_sponge_coeff.nc - end if - - ! read forcing data for u-sponge - name = 'INPUT/u_sponge.nc' - if (file_exist(trim(name)) ) then - Sponge_u(1)%id = init_external_field(name,'u',domain=Domain%domain2d) - if (Sponge_u(1)%id < 1) then - if ( use_adaptive_restore ) then - call mpp_error(FATAL,& - '==>Error: in ocean_sponges_velocity_mod: adaptive forcing is specified but restoring values are not') - else - call mpp_error(FATAL,& - '==>Error: in ocean_sponges_velocity_mod: forcing rates are specified but restoring values are not') - endif - endif - write(stdoutunit,*) '==> Using restoring data specified from file '//trim(name) - else - write(stdoutunit,*) '==> '//trim(name)//' not found. Increment not being applied ' - endif - - - ! read forcing data for v-sponge - name = 'INPUT/v_sponge.nc' - if (file_exist(trim(name)) ) then - Sponge_v(1)%id = init_external_field(name,'v',domain=Domain%domain2d) - if (Sponge_v(1)%id < 1) then - if ( use_adaptive_restore ) then - call mpp_error(FATAL,& - '==>Error: in ocean_sponges_velocity_mod: adaptive forcing is specified but restoring values are not') - else - call mpp_error(FATAL,& - '==>Error: in ocean_sponges_velocity_mod: forcing rates are specified but restoring values are not') - endif - endif - write(stdoutunit,*) '==> Using restoring data specified from file '//trim(name) - else - write(stdoutunit,*) '==> '//trim(name)//' not found. Increment not being applied ' - endif - - - ! register diagnostic output - allocate (id_sponge_tend(2)) - id_sponge_tend = -1 - - id_sponge_tend(1) = register_diag_field ('ocean_model', 'u_sponge_tend', & - Grd%vel_axes_uv(1:3), Time%model_time, 'rho*dzt*u_tendency due to sponge',& - '(kg/m^3)*(m^2/s^2)', missing_value=missing_value, range=(/-1.e10,1.e10/)) - id_sponge_tend(2) = register_diag_field ('ocean_model', 'v_sponge_tend', & - Grd%vel_axes_uv(1:3), Time%model_time, 'rho*dzt*v_tendency due to sponge',& - '(kg/m^3)*(m^2/s^2)', missing_value=missing_value, range=(/-1.e10,1.e10/)) - - -end subroutine ocean_sponges_velocity_init -! NAME="ocean_sponges_velocity_init" - - - -!####################################################################### -! -! -! -! This subroutine calculates thickness weighted and density weighted -! time tendencies due to damping by sponges. -! -! -subroutine sponge_velocity_source(Time, Thickness, Velocity) - - type(ocean_time_type), intent(in) :: Time - type(ocean_thickness_type), intent(in) :: Thickness - type(ocean_velocity_type), intent(inout) :: Velocity - - integer :: secs, days - integer :: taum1, tau - integer :: i, j, k - - real :: sdiff, sum_val, adaptive_coeff - integer :: numsecs - - logical :: do_adaptive_restore = .false. ! Only restore in specified time period. - logical :: do_normal_restore = .false. ! - logical, save :: first_pass = .true. - - if (.not. use_this_module) return - - taum1 = Time%taum1 - tau = Time%tau - - if (use_adaptive_restore) then - secs_restore = days_to_restore * 86400 + secs_to_restore - sdiff = 0.0 - call get_time(Time%model_time, secs, days) - numsecs = (days - initial_day) * 86400 + secs - initial_secs - - do_adaptive_restore = (numsecs < secs_restore) - do_normal_restore = (.not. do_adaptive_restore) .and. use_sponge_after_init - - write(stdout(),*) 'dar,dnr', do_adaptive_restore, do_normal_restore - else - do_normal_restore = .true. - endif - - wrk1 = 0.0 - wrk2 = 0.0 - if (Sponge_u(1)%id > 0) then - - call time_interp_external(Sponge_u(1)%id, Time%model_time, wrk1) - - if ( do_adaptive_restore ) then - if ( first_pass ) then - if ( use_hard_thump ) then - do k = 1,nk - do j=jsd,jed - do i=isd,ied - Velocity%u(i,j,k,1,taum1)=wrk1(i,j,k) - Velocity%u(i,j,k,1,tau)=wrk1(i,j,k) - enddo - enddo - enddo - endif - wrk2=abs(Velocity%u(:,:,:,1,taum1)-wrk1(:,:,:))*Grd%umask(:,:,:) - sum_val=sum(wrk2(isc:iec,jsc:jec,:)) - - call mpp_sum(sum_val) - sdiffo_u = sum_val/Grd%wet_u_points - !if ( mpp_pe() == mpp_root_pe() ) write (stdout(),*)'mean U_diff', sdiffo - write (stdout(),*)'mean U_diff', sdiffo_u - if ( .not. use_normalising ) then - sdiffo_u=1.0 - else - ! Work out max difference + threshhold - sdiffo_u = maxval(wrk2) - call mpp_max(sdiffo_u) - sdiffo_u = 1.01*sdiffo_u - endif + endif + else if (dtime*Sponge_v(1)%damp_coeff_v(i,j,k) <= 0.0) then + Sponge_v(1)%damp_coeff_v(i,j,k) = 0.0 + endif + enddo + enddo + enddo + endif !fileexist INPUT/v_sponge_coeff.nc + endif !.not. use_hard_thump .and. .not. use_adaptive_restore + + ! read forcing data for u-sponge + name = 'INPUT/u_sponge.nc' + if (file_exist(trim(name)) ) then + Sponge_u(1)%id = init_external_field(name,'u',domain=Domain%domain2d) + if (Sponge_u(1)%id < 1) then + call mpp_error(FATAL,& + '==>Error: in ocean_sponges_velocity_mod: Sponge values for u are required but not specified') + endif + write(stdoutunit,*) '==> Using restoring data specified from file '//trim(name) + else + write(stdoutunit,*) '==> '//trim(name)//' not found. Sponge not being applied ' + endif + + ! read forcing data for v-sponge + name = 'INPUT/v_sponge.nc' + if (file_exist(trim(name)) ) then + Sponge_v(1)%id = init_external_field(name,'v',domain=Domain%domain2d) + if (Sponge_v(1)%id < 1) then + call mpp_error(FATAL,& + '==>Error: in ocean_sponges_velocity_mod: Sponge values for v are required but not specified') + endif + write(stdoutunit,*) '==> Using restoring data specified from file '//trim(name) + else + write(stdoutunit,*) '==> '//trim(name)//' not found. Sponge not being applied ' + endif + + ! register diagnostic output + allocate (id_sponge_tend(2)) + id_sponge_tend = -1 + + id_sponge_tend(1) = register_diag_field ('ocean_model', 'u_sponge_tend', & + Grd%vel_axes_uv(1:3), Time%model_time, 'u_tendency due to sponge',& + '(m/s^2)', missing_value=missing_value, range=(/-1.e-2,1.e-2/)) + id_sponge_tend(2) = register_diag_field ('ocean_model', 'v_sponge_tend', & + Grd%vel_axes_uv(1:3), Time%model_time, 'v_tendency due to sponge',& + '(m/s^2)', missing_value=missing_value, range=(/-1.e-2,1.e-2/)) + + + end subroutine ocean_sponges_velocity_init + ! NAME="ocean_sponges_velocity_init" + + !####################################################################### + ! + ! + ! + ! This subroutine calculates thickness weighted and density weighted + ! time tendencies due to damping by sponges. + ! + ! + subroutine sponge_velocity_source(Time, Thickness, Velocity) + + type(ocean_time_type), intent(in) :: Time + type(ocean_thickness_type), intent(in) :: Thickness + type(ocean_velocity_type), intent(inout) :: Velocity + + integer :: secs, days + integer :: taum1, tau + integer :: i, j, k, n + + real :: sdiff, sum_val, adaptive_coeff + integer :: numsecs + + logical :: do_adaptive_restore = .false. ! Only restore in specified time period. + logical,save :: first_pass = .true. + + if (.not. use_this_module) return + + taum1 = Time%taum1 + tau = Time%tau + + wrk1 = 0.0 + wrk2 = 0.0 + sdiff=0.0 + + call get_time(Time%model_time,secs,days) + !----------------------------------------------------------------------------------------------------------------- + !----------------------------------------------------------------------------------------------------------------- + if (Sponge_u(1)%id > 0) then + call time_interp_external(Sponge_u(1)%id, Time%model_time, wrk1) + ! filter + do k=1,nk + do j=jsd,jed + do i=isd,ied + if ( wrk1(i,j,k).ge.limit_u_max ) then + wrk1(i,j,k)=limit_u_max + endif + if ( wrk1(i,j,k).le.limit_u_min ) then + wrk1(i,j,k)=limit_u_min + endif + enddo + enddo + enddo + !----------------------------------------------------------------------------------------------------------------- + if ( first_pass ) then + if ( use_hard_thump ) then + do k = 1,nk + do j=jsd,jed + do i=isd,ied + Velocity%u(i,j,k,1,taum1)=wrk1(i,j,k) + Velocity%u(i,j,k,1,tau)=wrk1(i,j,k) + enddo + enddo + enddo + call mpp_update_domains (Velocity%u(:,:,:,1,taum1) ,Dom%domain2d) + call mpp_update_domains (Velocity%u(:,:,:,1,tau) ,Dom%domain2d) + wrk2=abs(Velocity%u(:,:,:,1,taum1)-wrk1(:,:,:))*Grd%umask(:,:,:) + endif !use_hard_thump + sum_val=sum(wrk2(isc:iec,jsc:jec,:)) + call mpp_sum(sum_val) + sdiffo_u = sum_val/Grd%wet_u_points + write (stdout(),*)'mean absolute deviation u = ', sdiffo_u + if ( .not. use_normalising ) then + sdiffo_u=1.0 + else + !work out max difference + threshhold + sdiffo_u = maxval(wrk2) + call mpp_max(sdiffo_u) + sdiffo_u = 1.01*sdiffo_u endif - - ! Calculate idealised forcing tendency timescale for each timestep and array element - do k=1,nk - do j=jsd,jed - do i=isd,ied - sdiff=abs(Velocity%u(i,j,k,1,taum1)-wrk1(i,j,k)) - sdiff = max(sdiff,1e-9) !minimum tolerance - adaptive_coeff=1.0/(taumin - (lambda*real(secs_restore))/(((sdiff/sdiffo_v)**npower)*log(1-athresh))) - wrk2(i,j,k) = adaptive_coeff*(wrk1(i,j,k) - Velocity%u(i,j,k,1,taum1)) - Velocity%accel(i,j,k,1) = Velocity%accel(i,j,k,1) + Thickness%rho_dzu(i,j,k,tau)*wrk2(i,j,k) - enddo - enddo - enddo - - elseif ( do_normal_restore ) then - do k=1,nk - do j=jsd,jed - do i=isd,ied - wrk2(i,j,k) = Sponge_u(1)%damp_coeff_u(i,j,k)*(wrk1(i,j,k) - & - Velocity%u(i,j,k,1,taum1)) - Velocity%accel(i,j,k,1) = Velocity%accel(i,j,k,1) + & - Thickness%rho_dzu(i,j,k,tau)*wrk2(i,j,k) - enddo - enddo - enddo - endif - - endif - - call diagnose_3d_u(Time, Grd, id_sponge_tend(1), wrk2(:,:,:)) - - wrk1 = 0.0 - wrk2 = 0.0 - if (Sponge_v(1)%id > 0) then - - call time_interp_external(Sponge_v(1)%id, Time%model_time, wrk1) - - if ( do_adaptive_restore ) then - if ( first_pass ) then - if ( use_hard_thump ) then - do k = 1,nk - do j=jsd,jed - do i=isd,ied - Velocity%u(i,j,k,2,taum1)=wrk1(i,j,k) - Velocity%u(i,j,k,2,tau)=wrk1(i,j,k) - enddo - enddo - enddo - endif - wrk2=abs(Velocity%u(:,:,:,2,taum1)-wrk1(:,:,:))*Grd%umask(:,:,:) - sum_val=sum(wrk2(isc:iec,jsc:jec,:)) - - call mpp_sum(sum_val) - sdiffo_v = sum_val/Grd%wet_u_points - - if ( .not. use_normalising ) then - sdiffo_v=1.0 - else - sdiffo_v = maxval(wrk2) - call mpp_max(sdiffo_v) - sdiffo_v = 1.01*sdiffo_v - endif + endif !first_pass + !-------------------------------------------------------------------------------------------------------------- + if (use_hard_thump == .false.) then + if (use_adaptive_restore) then + secs_restore = days_to_restore*86400 + secs_to_restore + numsecs = (days-initial_day)*86400 + secs - initial_secs + do_adaptive_restore = ( numsecs < secs_restore ) endif - - - do k=1,nk - do j=jsd,jed - do i=isd,ied - sdiff=abs(Velocity%u(i,j,k,2,taum1)-wrk1(i,j,k)) - sdiff = max(sdiff,1e-9) !minimum tolerance - adaptive_coeff=1.0/(taumin - (lambda*real(secs_restore))/(((sdiff/sdiffo_v)**npower)*log(1-athresh))) - wrk2(i,j,k) = adaptive_coeff*(wrk1(i,j,k) - Velocity%u(i,j,k,2,taum1)) - Velocity%accel(i,j,k,2) = Velocity%accel(i,j,k,2) + Thickness%rho_dzu(i,j,k,tau)*wrk2(i,j,k) - enddo - enddo - enddo - endif - - if ( do_normal_restore ) then - do k=1,nk - do j=jsc,jec - do i=isc,iec - wrk2(i,j,k) = Sponge_v(1)%damp_coeff_v(i,j,k)*(wrk1(i,j,k) - & - Velocity%u(i,j,k,2,taum1)) - Velocity%accel(i,j,k,2) = Velocity%accel(i,j,k,2) + & - Thickness%rho_dzu(i,j,k,tau)*wrk2(i,j,k) - enddo - enddo - enddo - endif - endif - - call diagnose_3d_u(Time, Grd, id_sponge_tend(2), wrk2(:,:,:)) - - first_pass = .false. - - return - -end subroutine sponge_velocity_source -! NAME="sponge_velocity_source" - + if ( do_adaptive_restore ) then + !calculate idealised forcing tendency timescale for each timestep and array element + do k=1,nk + do j=jsd,jed + do i=isd,ied + sdiff=abs(Velocity%u(i,j,k,1,taum1)-wrk1(i,j,k)) + sdiff = max(sdiff,1e-9) !minimum tolerance + adaptive_coeff=1.0/(taumin - (lambda*real(secs_restore))/(((sdiff/sdiffo_u)**npower)*log(1-athresh))) + wrk2(i,j,k) = adaptive_coeff*(wrk1(i,j,k) - Velocity%u(i,j,k,1,taum1)) + Velocity%accel(i,j,k,1) = Velocity%accel(i,j,k,1) + Thickness%rho_dzu(i,j,k,tau)*wrk2(i,j,k) + enddo + enddo + enddo + else !normal restore + do k=1,nk + do j=jsd,jed + do i=isd,ied + wrk2(i,j,k) = Sponge_u(1)%damp_coeff_u(i,j,k)*(wrk1(i,j,k) - & + Velocity%u(i,j,k,1,taum1)) + Velocity%accel(i,j,k,1) = Velocity%accel(i,j,k,1) + & + Thickness%rho_dzu(i,j,k,tau)*wrk2(i,j,k) + enddo + enddo + enddo + endif !do_adaptive_restore + endif !use_hard_thump + !----------------------------------------------------------------------------------------------------------------- + endif !Sponge_u + !-------------------------------------------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------------------------------------------- + if (id_sponge_tend(1) > 0) then + used = send_data(id_sponge_tend(1), & + wrk2(:,:,:), 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 + !-------------------------------------------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------------------------------------------- + wrk1 = 0.0 + wrk2 = 0.0 + if (Sponge_v(1)%id > 0) then + call time_interp_external(Sponge_v(1)%id, Time%model_time, wrk1) + ! filter + do k=1,nk + do j=jsd,jed + do i=isd,ied + if ( wrk1(i,j,k).ge.limit_v_max ) then + wrk1(i,j,k)=limit_v_max + endif + if ( wrk1(i,j,k).le.limit_v_min ) then + wrk1(i,j,k)=limit_v_min + endif + enddo + enddo + enddo + !-------------------------------------------------------------------------------------------------------------- + if ( first_pass ) then + if ( use_hard_thump ) then + do k = 1,nk + do j=jsd,jed + do i=isd,ied + Velocity%u(i,j,k,2,taum1)=wrk1(i,j,k) + Velocity%u(i,j,k,2,tau)=wrk1(i,j,k) + enddo + enddo + enddo + call mpp_update_domains (Velocity%u(:,:,:,2,taum1) ,Dom%domain2d) + call mpp_update_domains (Velocity%u(:,:,:,2,tau) ,Dom%domain2d) + wrk2=abs(Velocity%u(:,:,:,2,taum1)-wrk1(:,:,:))*Grd%umask(:,:,:) + endif + sum_val=sum(wrk2(isc:iec,jsc:jec,:)) + call mpp_sum(sum_val) + sdiffo_v = sum_val/Grd%wet_u_points + write (stdout(),*)'mean absolute deviation v = ', sdiffo_v + if ( .not. use_normalising ) then + sdiffo_v=1.0 + else + sdiffo_v = maxval(wrk2) + call mpp_max(sdiffo_v) + sdiffo_v = 1.01*sdiffo_v + endif + endif !first_pass + !-------------------------------------------------------------------------------------------------------------- + if (use_hard_thump == .false.) then + if (use_adaptive_restore) then + secs_restore = days_to_restore*86400 + secs_to_restore + numsecs = (days-initial_day)*86400 + secs - initial_secs + do_adaptive_restore = ( numsecs < secs_restore ) + endif + if ( do_adaptive_restore ) then + do k=1,nk + do j=jsd,jed + do i=isd,ied + sdiff=abs(Velocity%u(i,j,k,2,taum1)-wrk1(i,j,k)) + sdiff = max(sdiff,1e-9) !minimum tolerance + adaptive_coeff=1.0/(taumin - (lambda*real(secs_restore))/(((sdiff/sdiffo_v)**npower)*log(1-athresh))) + wrk2(i,j,k) = adaptive_coeff*(wrk1(i,j,k) - Velocity%u(i,j,k,2,taum1)) + Velocity%accel(i,j,k,2) = Velocity%accel(i,j,k,2) + Thickness%rho_dzu(i,j,k,tau)*wrk2(i,j,k) + enddo + enddo + enddo + else !normal restore + do k=1,nk + do j=jsd,jed + do i=isd,ied + wrk2(i,j,k) = Sponge_v(1)%damp_coeff_v(i,j,k)*(wrk1(i,j,k) - & + Velocity%u(i,j,k,2,taum1)) + Velocity%accel(i,j,k,2) = Velocity%accel(i,j,k,2) + & + Thickness%rho_dzu(i,j,k,tau)*wrk2(i,j,k) + enddo + enddo + enddo + endif !do_adaptive_restore + endif !use_hard_thump + !----------------------------------------------------------------------------------------------------------------- + endif !Sponge_v + !-------------------------------------------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------------------------------------------- + if (id_sponge_tend(2) > 0) then + used = send_data(id_sponge_tend(2), & + wrk2(:,:,:), 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 + !-------------------------------------------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------------------------------------------- + first_pass = .false. + + return + + end subroutine sponge_velocity_source + ! NAME="sponge_velocity_source" end module ocean_sponges_velocity_mod