!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! 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 !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #include module transport_matrix_mod !{ ! ! Samar Khatiwala ! ! ! Rick Slater ! ! Stephen M. Griffies ! ! Jennifer Simeon ! ! ! ! Transport Matrix Method ! ! ! ! Transport Matrix Method, for use in finding an approximate steady ! state of tracers. ! ! Ported to MOM4p0d by Samar Khatiwala spk@ldeo.columbia.edu ! June-July 2007 ! ! Ported to MOM4p1 by Stephen.Griffies@noaa.gov ! July 2007 ! ! Some code clean-up; increased conformity to tracer module standard ! jes.13JUN08 ! ! further code cleanup by Richard.Slater@noaa.gov Jan2009. ! ! **Preliminary testing only. ! ! This module saves out the components of the explicit and implicit ! transport matrix. spk uses the output from the transport matrix module ! to assemble the actual transport matrix with his matlab code. ! ! In conjunction with the Newton-Krylov (NK) solver, the transport matrix ! is used to do an accelerated forward model in an iterative manner ! which allows the NK solver to reach solution convergence with greater ! efficiency. ! ! To create the transport matrix building blocks, passive, prognostic ! "tracers" are simulated. These "tracers" are, in practice, dye ! tracers initialized to 1 at one or more grid points and initialized ! to 0 elsewhere. Since the "area of influence" of ! a tracer at a particular grid point is determined by the advection scheme, ! it becomes more practical to introduce a number of dye tracers at other ! grid points, such that the "tracers" are staggered and will have ! non-overlapping areas of influence. These non-overlapping regions are ! called by spk, "tiles". The total number of dye tracers ! needed is independent of the horizontal grid resolution, but ! rather dependent upon the advection scheme's area of influence, i.e. ! the total number of grid points contained in a "tile". For a ! simple linear advection scheme, the number of tracers needed is about ! 10 x (number of vertical levels). ! spk has some special matlab code that generates the initial condition ! for the "tracers", given the grid_spec of the model. ! ! The behavior of the tracers' advection and diffusion is averaged over ! time (which the user may specify) and is saved out as a component ! of the later-to-be-assembled transport matrix. ! ! To run the transport_matrix module, an initial condition ! must be created by spk. ! ! Include the transport_matrix field_table and diag_table entries ! for the xml. Examples of these follow. ! ! Sample field table. ! ------------------------------------------------- ! "tracer_packages","ocean_mod","transport_matrix" ! ! names = '01', '02', '03' ! horizontal-advection-scheme = mdfl_sweby ! vertical-advection-scheme = mdfl_sweby ! / ! ! -------------------------------------------------- ! ! Sample diag table entry: enter as many "tracers" as you need ! with the naming convention exp_tm_# and imp_tm_# ! where # is a string as given abaove in "names" ! -------------------------------------------------- !"transport_matrix","exp_tm_01", "exp_tm_01" ,"ocean_transport_matrix","all",.false.,"none",1 !"transport_matrix","imp_tm_01", "imp_tm_01" ,"ocean_transport_matrix","all",.false.,"none",1 !"transport_matrix","exp_tm_02", "exp_tm_02" ,"ocean_transport_matrix","all",.false.,"none",1 !"transport_matrix","imp_tm_02", "imp_tm_02" ,"ocean_transport_matrix","all",.false.,"none",1 !"transport_matrix","exp_tm_03", "exp_tm_03" ,"ocean_transport_matrix","all",.false.,"none",1 !"transport_matrix","imp_tm_03", "imp_tm_03" ,"ocean_transport_matrix","all",.false.,"none",1 ! -------------------------------------------------- ! ! ! SPK NOTES: ! 1) The calling sequence is as follows: ! Top level driver (e.g., ocean_solo) ! -> S/R ocean_model_init(Ocean, Time_init, Time_in, Time_step_ocean, ensemble_ocean) ! Time%init = Time_in .eq. Time_init ! Time%Time_init = Time_init ! Time%Time_step = Time_step_ocean ! Time%model_time = Time_in ! Time%itt = 0 ! -> S/R ocean_prog_tracer_init ! -> S/R ocean_tpm_init ! -> S/R transport_matrix_init ! -> S/R ocean_tpm_start ! -> S/R transport_matrix_start ! Start time stepping loop ! do nc=1, num_cpld_calls ! do no=1, num_ocean_calls ! ocean_seg_start = ( no .eq. 1 ) ! ocean_seg_end = ( no .eq. num_ocean_calls ) ! -> S/R update_ocean_model(Ice_ocean_boundary, Ocean_sfc, & ! ocean_seg_start, ocean_seg_end, num_ocean_calls) ! Time%model_time = Time%model_time + Time%Time_step ! Time%itt = Time%itt+1 ! -> S/R ocean_tracer, S/R update_ocean_tracer ! do explicit transport ! -> S/R transport_matrix_store_explicit ! (accumulate explicit matrix and reset tracer field to initial condition) ! do implicit transport ! -> S/R ocean_tpm_tracer ! -> S/R transport_matrix_store_implicit ! (accumulate implicit matrix and reset tracer field to initial condition) ! -> S/R transport_matrix_write(.FALSE.) (time average and write matrices) ! Time = Time + Time_step_ocean ! enddo ! enddo ! -> S/R ocean_model_end ! -> S/R ocean_tpm_end ! -> S/R transport_matrix_end ! -> S/R transport_matrix_write(.TRUE.) ! (time average and write matrices for multi year runs) ! ! 2) Time counters are incremented BEFORE calling S/R update_ocean_tracer, so the first time ! transport_matrix_store_explicit is called, itt (and hence myIter) will be 1. ! The namelist parameter matrixStoreStartIter indicating the iteration number to begin ! accumulating matrices at should be RELATIVE to the current model run start (unlike ! in the MIT GCM where it refers to an absolute counter. ! I am not entirely certain this is handled correctly below. Things might be off by ! 1 time step. ! ! ! ! ! ! ! ! Khatiwala, S., M. Visbeck, M.A. Cane, 2005. ! Accelerated simulation of passive tracers in ocean circulation models. ! Ocean Modelling, 9, 51-69. ! ! ! ! $Id: transport_matrix.F90,v 16.0.2.1.32.1.52.1 2009/10/10 00:43:09 nnz Exp $ ! use field_manager_mod, only: fm_string_len, fm_path_name_len, fm_field_name_len use field_manager_mod, only: fm_get_length, fm_get_value use mpp_mod, only: stdout, mpp_error, FATAL use diag_manager_mod, only: send_data, register_diag_field use ocean_tpm_util_mod, only: otpm_set_tracer_package, otpm_set_prog_tracer use ocean_types_mod, only: ocean_time_type, ocean_prog_tracer_type ! ! force all variables to be "typed" ! implicit none ! ! Set all variables to be private by default private ! ! Public routines ! public :: transport_matrix_init public :: transport_matrix_start public :: transport_matrix_store_explicit public :: transport_matrix_store_implicit ! ! Private routines ! ! ! Public parameters ! ! ! Private parameters ! character(len=fm_field_name_len), parameter :: package_name = 'transport_matrix' character(len=48), parameter :: diag_name = 'transport_matrix' character(len=48), parameter :: mod_name = 'transport_matrix_mod' character(len=fm_string_len), parameter :: default_restart_file = 'transport_matrix.res.nc' ! ! Public variables ! logical, public :: do_transport_matrix ! ! Private types ! type instance_type character(len=fm_field_name_len) :: name real, dimension(:,:,:), _ALLOCATABLE :: initial_value _NULL integer :: index integer :: id_extm integer :: id_imtm end type instance_type ! ! Private variables ! type(instance_type), dimension(:), allocatable :: instance integer :: instances integer :: package_index contains !####################################################################### ! ! ! ! Set up any extra fields needed by tracer package manager ! ! subroutine transport_matrix_init !{ ! ! local parameters ! character(len=64), parameter :: sub_name = 'transport_matrix_init' character(len=256), parameter :: error_header = & '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' character(len=256), parameter :: warn_header = & '==>Warning from ' // trim(mod_name) // '(' // trim(sub_name) // '):' character(len=256), parameter :: note_header = & '==>Note from ' // trim(mod_name) // '(' // trim(sub_name) // '):' ! !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- ! integer :: n character(len=fm_field_name_len) :: name character(len=fm_path_name_len) :: path_to_names character(len=fm_field_name_len+1) :: suffix character(len=fm_field_name_len+3) :: long_suffix integer :: stdoutunit stdoutunit=stdout() ! ! Initialize the ocean transport matrix package ! package_index = otpm_set_tracer_package(package_name, & caller = trim(mod_name) // '(' // trim(sub_name) // ')', & restart_file = default_restart_file) ! ! Check whether to use this package ! path_to_names = '/ocean_mod/tracer_packages/' // trim(package_name) // '/names' instances = fm_get_length(path_to_names) if (instances .lt. 0) then !{ call mpp_error(FATAL, trim(error_header) // ' Could not get number of instances') endif !} ! ! Check some things ! write (stdoutunit,*) if (instances .eq. 0) then !{ write (stdoutunit,*) & trim(note_header), ' No instances' do_transport_matrix = .false. else !}{ write (stdoutunit,*) trim(note_header), ' ', instances, ' instances' do_transport_matrix = .true. endif !} ! ! Return if we don't want to use this package ! if (.not. do_transport_matrix) then !{ return endif !} ! ! allocate the instance structure ! allocate ( instance(instances) ) ! ! loop over the names, saving them into the instance array do n = 1, instances !{ if (fm_get_value(path_to_names, name, index = n)) then !{ instance(n)%name = name else !}{ write (name,*) n call mpp_error(FATAL, trim(error_header) // & ' Bad field name for index ' // trim(name)) endif !} enddo !} n do n = 1, instances !{ !---------------------------------------------------------------------------------------- ! Register the prognostic tracers !---------------------------------------------------------------------------------------- ! Set up the instance name from field table name = instance(n)%name if (name(1:1) .eq. '_') then !{ suffix = ' ' long_suffix = ' ' else !}{ suffix = '_' // name long_suffix = ' (' // trim(name) // ')' endif !} ! determine the tracer names for this instance and assign indices instance(n)%index = otpm_set_prog_tracer('ptr' // trim(suffix), & package_name, longname = 'Transport matrix passive tracer' // trim(long_suffix), & units = 'unit/kg', flux_units = 'unit/m^2/s', & caller = trim(mod_name)//'('//trim(sub_name)//')') enddo !} n return end subroutine transport_matrix_init !} ! NAME="transort_matrix_init" !####################################################################### ! ! ! ! ! ! subroutine transport_matrix_start(Time, T_prog, isd, ied, jsd, jed, nk, & grid_tracer_axes) !{ !----------------------------------------------------------------------- ! modules (have to come first) !----------------------------------------------------------------------- ! !----------------------------------------------------------------------- ! Arguments !----------------------------------------------------------------------- ! type(ocean_time_type), intent(in) :: Time type(ocean_prog_tracer_type), dimension(:), intent(in) :: T_prog integer, intent(in) :: isd integer, intent(in) :: ied integer, intent(in) :: jsd integer, intent(in) :: jed integer, intent(in) :: nk integer, dimension(3), intent(in) :: grid_tracer_axes ! ! local parameters ! character(len=64), parameter :: sub_name = 'transport_matrix_start' character(len=256), parameter :: error_header = & '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' character(len=256), parameter :: warn_header = & '==>Warning from ' // trim(mod_name) // '(' // trim(sub_name) // '):' character(len=256), parameter :: note_header = & '==>Note from ' // trim(mod_name) // '(' // trim(sub_name) // '):' ! ! local variables ! integer :: n character(len=fm_field_name_len) :: name character(len=fm_field_name_len+3) :: long_suffix character(len=fm_field_name_len+1) :: suffix ! ! save initial values ! do n = 1, instances !{ allocate( instance(n)%initial_value(isd:ied,jsd:jed,nk) ) instance(n)%initial_value(:,:,:) = T_prog(instance(n)%index)%field(:,:,:,Time%tau) enddo !} n !----------------------------------------------------------------------- ! Set up analyses !----------------------------------------------------------------------- ! ! ! register the diag fields, i.e. the TMs ! do n = 1, instances !{ ! Set up the instance name from field table name = instance(n)%name if (name(1:1) .eq. '_') then !{ suffix = ' ' long_suffix = ' ' else !}{ suffix = '_' // name long_suffix = ' (' // trim(name) // ')' endif !} ! determine the tracer names for this instance and register instance(n)%id_extm = register_diag_field(diag_name, & 'exp_tm' // trim(suffix), grid_tracer_axes(1:3), Time%model_time, & 'Explicit transport matrix for passive tracer ' // trim(long_suffix), ' ', & missing_value = -1.0e+10) instance(n)%id_imtm = register_diag_field(diag_name, & 'imp_tm' // trim(suffix), grid_tracer_axes(1:3), Time%model_time, & 'Implicit transport matrix for passive tracer ' // trim(long_suffix), ' ', & missing_value = -1.0e+10) enddo !} n return end subroutine transport_matrix_start !} ! NAME="transport_matrix_start" !####################################################################### ! ! ! ! For the time explicit tendencies. ! ! subroutine transport_matrix_store_explicit(Time, T_prog, isd, ied, jsd, jed, nk, & isc, iec, jsc, jec, dtts, grid_tmask) type(ocean_time_type), intent(in) :: Time type(ocean_prog_tracer_type), dimension(:), intent(inout) :: T_prog integer, intent(in) :: isd integer, intent(in) :: ied integer, intent(in) :: jsd integer, intent(in) :: jed integer, intent(in) :: nk integer, intent(in) :: isc integer, intent(in) :: iec integer, intent(in) :: jsc integer, intent(in) :: jec real, intent(in) :: dtts real, dimension(isd:ied,jsd:jed,nk), intent(in) :: grid_tmask ! ! local parameters ! character(len=64), parameter :: sub_name = 'transport_matrix_store_explicit' character(len=256), parameter :: error_header = & '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' character(len=256), parameter :: warn_header = & '==>Warning from ' // trim(mod_name) // '(' // trim(sub_name) // '):' character(len=256), parameter :: note_header = & '==>Note from ' // trim(mod_name) // '(' // trim(sub_name) // '):' ! ! local variables ! integer :: n real, dimension(isd:ied,jsd:jed,nk) :: temp logical :: used if(.not. do_transport_matrix) return do n = 1, instances !{ if (instance(n)%id_extm .gt. 0) then !{ temp(:,:,:) = (T_prog(instance(n)%index)%field(:,:,:,Time%taup1) - instance(n)%initial_value(:,:,:)) / dtts used = send_data(instance(n)%id_extm, temp, Time%model_time, rmask = grid_tmask, & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) endif !} T_prog(instance(n)%index)%field(:,:,:,Time%taup1) = instance(n)%initial_value(:,:,:) enddo !} n return end subroutine transport_matrix_store_explicit !} ! NAME="transport_matrix_store_explicit" !####################################################################### ! ! ! ! For the time implicit tendencies. ! ! subroutine transport_matrix_store_implicit(Time, T_prog, isd, ied, jsd, jed, nk, & isc, iec, jsc, jec, grid_tmask) !{ !--------------------------------------------------------------- ! Modules first !--------------------------------------------------------------- !--------------------------------------------------------------- ! Arguments !--------------------------------------------------------------- type(ocean_time_type), intent(in) :: Time type(ocean_prog_tracer_type), dimension(:), intent(inout) :: T_prog integer, intent(in) :: isd integer, intent(in) :: ied integer, intent(in) :: jsd integer, intent(in) :: jed integer, intent(in) :: nk integer, intent(in) :: isc integer, intent(in) :: iec integer, intent(in) :: jsc integer, intent(in) :: jec real, dimension(isd:ied,jsd:jed,nk), intent(in) :: grid_tmask ! ! local parameters ! character(len=64), parameter :: sub_name = 'transport_matrix_store_implicit' character(len=256), parameter :: error_header = & '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' character(len=256), parameter :: warn_header = & '==>Warning from ' // trim(mod_name) // '(' // trim(sub_name) // '):' character(len=256), parameter :: note_header = & '==>Note from ' // trim(mod_name) // '(' // trim(sub_name) // '):' ! ! local variables ! integer :: n logical :: used do n = 1, instances !{ if (instance(n)%id_imtm .gt. 0) then !{ used = send_data(instance(n)%id_imtm, T_prog(instance(n)%index)%field(:,:,:,Time%taup1), & Time%model_time, rmask = grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) endif !} T_prog(instance(n)%index)%field(:,:,:,Time%taup1) = instance(n)%initial_value(:,:,:) enddo !} n return end subroutine transport_matrix_store_implicit !} ! NAME="transport_matrix_store_implicit" end module transport_matrix_mod !}