module ocean_xlandinsert_mod
!
! S.M. Griffies
!
!
! K.W. Dixon
!
!
! M.J. Harrison
!
!
!
! Tracer and mass source from cross-land insertion.
!
!
!
! Compute thickness and density weighted tendency [rho*tracer*meter/sec]
! of tracer associated with cross-land insertion. Also compute
! mass source from cross-land insertion.
!
!
!
!
!
! S.M. Griffies, M.J. Harrison, R. C. Pacanowski, and A. Rosati
! A Technical Guide to MOM4 (2003)
! NOAA/Geophysical Fluid Dynamics Laboratory
!
!
!
! Algorithm ensures both total tracer and total mass is conserved.
! Algorithm sets up the insertion points as in xlandmix, and
! transports between the columns as in rivermix.
!
!
!
! 2D domain decomposition is implemented according to following
! notions. If the two points in an xlandinsert pair live within halo
! of each other (i.e., within same local_domain),
! then no added mpp communication required. However, nore generally
! the two points live further away. In this case, xland_domain
! is defined so that its halos incorporate the maximum separation
! of xlandinsert points. New tracer, eta, and grid arrays
! are defined over this extended xland_domain. This added domain
! size will come at some computational cost, so it is advantageous
! to limit the separation between points within an xlandinsert pair.
!
!
!
! The current implementation of xlandinsert has not been generalized
! to allow for communication between points separated by the tripolar fold.
! The problem is in the logic used to compute xland_domain.
! There is nothing fundamental limiting a generalization of the
! logic used to generate xland_domain.
!
!
!
! Many of the user specified values given in USER INPUT are
! model dependent since unresolved straits can become resolved
! in finer mesh models.
!
!
!
! BUG FIX 30 May 2011.
! Trying to mix across the ends of the global domain results in incorrect indexing.
! Local indices need to be incremented/decremented by the global extent.
! Russell.Fiedler@csiro.au
!
!
!
!
!
!
! Needs to be true in order to use this scheme. Default is false.
!
!
! For verbose initialization information.
!
!
! For debugging.
!
!
use constants_mod, only: epsln
use diag_manager_mod, only: register_diag_field
use field_manager_mod, only: MODEL_OCEAN, parse, find_field_index
use field_manager_mod, only: get_field_methods, method_type, get_field_info
use fms_mod, only: stdout, stdlog, FATAL, NOTE, WARNING
use fms_mod, only: write_version_number, open_namelist_file, check_nml_error, close_file
use mpp_domains_mod, only: mpp_update_domains
use mpp_domains_mod, only: BITWISE_EXACT_SUM
use mpp_mod, only: input_nml_file, mpp_error
use ocean_domains_mod, only: get_local_indices, get_global_indices, set_ocean_domain
use ocean_parameters_mod, only: rho0, rho0r, missing_value
use ocean_types_mod, only: ocean_domain_type, ocean_grid_type, ocean_time_type
use ocean_types_mod, only: ocean_thickness_type, ocean_external_mode_type
use ocean_types_mod, only: ocean_prog_tracer_type, ocean_options_type, ocean_density_type
use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk5, wrk1_v
use ocean_util_mod, only: diagnose_2d, diagnose_3d, diagnose_sum
use ocean_tracer_util_mod,only: diagnose_3d_rho
implicit none
private
public ocean_xlandinsert_init
public xlandinsert
private xlandinsert_mass
private watermass_diag_init
private watermass_diag
type(ocean_domain_type), pointer :: Dom => NULL()
type(ocean_grid_type), pointer :: Grd => NULL()
! Variables set inside subroutine ocean_xlandinsert_init
! See comments in ocean_xlandinsert_init for description
integer, dimension(:,:), allocatable :: ixland
integer, dimension(:,:), allocatable :: jxland
integer, dimension(:,:), allocatable :: kxland
real, dimension(:), allocatable :: tauxland
real, dimension(:), allocatable :: tauxlandr
real, dimension(:,:), allocatable :: water_rate ! rate (kg/m^3)*(m/s) of water transferred between columns
! number of xland mixing pairs
integer :: nxland=0
integer :: xland_halo ! halo size needed to perform xlandinsert
integer :: kxland_min ! min of the k-levels for the nxland pairs
integer :: kxland_max ! max of the k-levels for the nxland pairs
type(ocean_domain_type), save :: Xland_domain ! domain for xlandinsert
! variables defined with halos given by xland_halo
real, dimension(:,:), allocatable :: xland_rho_dzt ! mass source on (x,y) space from xlandinsert (kg/m^3)*(m/s)
real, dimension(:,:), allocatable :: rho_dzt_x ! surface value for rho_dzt (kg/m^3)*(m)
real, dimension(:,:), allocatable :: xtx ! zonal position of tracer points (degrees)
real, dimension(:,:), allocatable :: ytx ! meridional position of tracer points (degrees)
real, dimension(:,:), allocatable :: datx ! area of tracer cell (m^2)
real :: dtime ! time step
integer :: unit=6 ! processor zero writes to unit 6
logical, dimension(:), allocatable :: error_xland ! for checking that all xland points are OK.
! work array for eta_tend
real, dimension(:,:,:), allocatable :: theta_tend
real, dimension(:,:,:), allocatable :: salt_tend
! for global area normalization
real :: cellarea_r
! internally set for computing watermass diagnostics
logical :: compute_watermass_diag = .false.
! for diagnostics
logical :: used
integer, dimension(:), allocatable :: id_xland
integer :: id_xland_rho_dzt=-1
integer :: id_neut_rho_xinsert =-1
integer :: id_neut_rho_xinsert_on_nrho =-1
integer :: id_wdian_rho_xinsert =-1
integer :: id_wdian_rho_xinsert_on_nrho =-1
integer :: id_tform_rho_xinsert =-1
integer :: id_tform_rho_xinsert_on_nrho =-1
integer :: id_eta_tend_xinsert =-1
integer :: id_eta_tend_xinsert_glob=-1
integer :: id_neut_temp_xinsert =-1
integer :: id_neut_temp_xinsert_on_nrho =-1
integer :: id_wdian_temp_xinsert =-1
integer :: id_wdian_temp_xinsert_on_nrho =-1
integer :: id_tform_temp_xinsert =-1
integer :: id_tform_temp_xinsert_on_nrho =-1
integer :: id_neut_salt_xinsert =-1
integer :: id_neut_salt_xinsert_on_nrho =-1
integer :: id_wdian_salt_xinsert =-1
integer :: id_wdian_salt_xinsert_on_nrho =-1
integer :: id_tform_salt_xinsert =-1
integer :: id_tform_salt_xinsert_on_nrho =-1
integer :: num_prog_tracers=0
integer :: index_temp=-1
integer :: index_salt=-1
character(len=128) :: version=&
'$Id: ocean_xlandinsert.F90,v 20.0 2013/12/14 00:16:28 fms Exp $'
character (len=128) :: tagname = &
'$Name: tikal $'
integer :: isd, ied, jsd, jed, isc, iec, jsc, jec, nk
integer :: isg, ieg, jsg, jeg
integer :: ixl, jxl, ixl2, jxl2
logical :: verbose_init = .true.
logical :: use_this_module = .false.
logical :: module_is_initialized = .FALSE.
logical :: debug_this_module = .false.
namelist /ocean_xlandinsert_nml/ verbose_init, use_this_module, debug_this_module
contains
!#######################################################################
!
!
!
! Initial set up for crossland insertion of tracers and eta
!
! Checks are performed to ensure that the crossland mixing
! grid locations are valid according to model configuration.
!
! A summary of the locations of crossland points is written out.
!
! User specified inputs in "USER INPUT" section:
!
! ixland and jxland = user specified nxland pairs of i,j grid locations
!
! kxland = user specified uppermost (ktop=kxland(n,1)) and
! deepest (kbot=kxland(n,2)) levels over which crossland
! will be done for each crossland pair. Note the for
! xlandinsert, kxland(nxl,1)=1 is required, since the
! aim is to move excess volume from one column to another.
!
! tauxland = user specified time constant (seconds) setting the rate
! of transport via upwind advection.
!
! NOTE: for ixland, jxland, and kxland pairs, values of the
! (n,1) element should be < the corresponding (n,2) value.
!
!
subroutine ocean_xlandinsert_init(Grid, Domain, Time, Dens, T_prog, Ocean_options, d_time)
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_density_type), intent(in) :: Dens
type(ocean_prog_tracer_type), intent(in) :: T_prog(:)
type(ocean_options_type), intent(inout) :: Ocean_options
real, intent(in) :: d_time
type(method_type), allocatable, dimension(:) :: xland_methods
character(len=32) :: fld_type, fld_name
integer :: n, nt, model, imax, imin
integer :: nxl, lx, xland_halo_test,i
integer :: parse_ok
integer :: io_status, ioun, ierr
real :: ztop
integer :: stdoutunit,stdlogunit
stdoutunit=stdout();stdlogunit=stdlog()
if (module_is_initialized ) then
call mpp_error(FATAL, &
'==>Error in ocean_xlandinsert_mod (ocean_xlandinsert_init): module already initialized')
endif
module_is_initialized = .TRUE.
call write_version_number(version, tagname)
call get_local_indices(Domain,isd,ied,jsd,jed,isc,iec,jsc,jec)
call get_global_indices(Domain,isg,ieg,jsg,jeg)
nk = Grid%nk
Dom => Domain
Grd => Grid
dtime = d_time
cellarea_r = 1.0/(epsln + Grd%tcellsurf)
num_prog_tracers = size(T_prog(:))
do nt=1,num_prog_tracers
if (T_prog(nt)%name == 'temp') index_temp = nt
if (T_prog(nt)%name == 'salt') index_salt = nt
enddo
n = find_field_index(MODEL_OCEAN,'xland_insert')
if (n < 1) then
write(stdoutunit,*)' '
write(stdoutunit,*) &
'==>Warning: ocean_xlandinsert_init found n < 1 for xland_insert table. Will NOT use ocean_xlandinsert.'
write(stdoutunit,*)' '
Ocean_options%xlandinsert = 'Did NOT use xlandinsert.'
return
endif
call get_field_info(n,fld_type,fld_name,model,nxland)
! provide for namelist over-ride of default values
#ifdef INTERNAL_FILE_NML
read (input_nml_file, nml=ocean_xlandinsert_nml, iostat=io_status)
ierr = check_nml_error(io_status,'ocean_xlandinsert_nml')
#else
ioun = open_namelist_file()
read (ioun,ocean_xlandinsert_nml,IOSTAT=io_status)
ierr = check_nml_error(io_status,'ocean_xlandinsert_nml')
call close_file (ioun)
#endif
write (stdlogunit,ocean_xlandinsert_nml)
write (stdoutunit,'(/)')
write (stdoutunit,ocean_xlandinsert_nml)
if(use_this_module) then
if(nxland == 0) then
write(stdoutunit, *) &
'==>Warning: use_this_module=.true. but no mixing points chosen. Is this what you wish?'
call mpp_error(WARNING,&
'==>NOTE from ocean_xlandinsert_mod: No cross-land mixing points chosen')
Ocean_options%xlandinsert = 'Used xlandinsert, but set zero cross-land points, so module does nothing.'
return
else
call mpp_error(NOTE, &
'==>Note from ocean_xlandinsert_mod (ocean_xlandinsert_init): USING this module')
write(stdoutunit, *) &
'==>NOte: Using xlandinsert to connect tracers and mass between non-local ocean cells.'
Ocean_options%xlandinsert = 'Used xlandinsert to transfer tracer between inland seas and open ocean.'
endif
else
call mpp_error(NOTE, '==>Note from ocean_xlandinsert_mod: NOT USING this module')
if(nxland > 0) then
call mpp_error(WARNING, &
'==>Warning: nxland > 0, but use_this_module=.false. NOT using xlandinsert. Is this what you wish?')
endif
Ocean_options%xlandinsert = 'Did NOT use xlandinsert.'
return
endif
write(stdoutunit,'(/a,f10.2)')'==> Note from ocean_xlandinsert_mod: using time step of (secs)', dtime
allocate(xland_methods(nxland))
allocate (ixland(nxland,2))
ixland(:,:) = 0
allocate (jxland(nxland,2))
jxland(:,:) = 0
allocate (kxland(nxland,2))
kxland(:,:) = 0
allocate (tauxland(nxland))
tauxland(:) = 0.0
allocate (tauxlandr(nxland))
tauxlandr(:) = 0.0
allocate (error_xland(nxland))
error_xland(:) = .false.
allocate (water_rate(nxland,2))
water_rate(:,:) = 0.0
call get_field_methods(n,xland_methods)
do i=1, nxland
parse_ok = parse(xland_methods(i)%method_control,'ixland_1',ixland(i,1))
if (parse_ok == 0) call mpp_error(FATAL,'==>Error ocean_xlandinsert_mod: xland table entry "ixland_1" error')
parse_ok = parse(xland_methods(i)%method_control,'ixland_2',ixland(i,2))
if (parse_ok == 0) call mpp_error(FATAL,'==>Error ocean_xlandinsert_mod: xland table entry "ixland_2" error')
parse_ok = parse(xland_methods(i)%method_control,'jxland_1',jxland(i,1))
if (parse_ok == 0) call mpp_error(FATAL,'==>Error ocean_xlandinsert_mod: xland table entry "jxland_1" error')
parse_ok = parse(xland_methods(i)%method_control,'jxland_2',jxland(i,2))
if (parse_ok == 0) call mpp_error(FATAL,'==>Error ocean_xlandinsert_mod: xland table entry "jxland_2" error')
parse_ok = parse(xland_methods(i)%method_control,'kxland_1',kxland(i,1))
if (parse_ok == 0) call mpp_error(FATAL,'==>Error ocean_xlandinsert_mod: xland table entry "kxland_1" error')
parse_ok = parse(xland_methods(i)%method_control,'kxland_2',kxland(i,2))
if (parse_ok == 0) call mpp_error(FATAL,'==>Error ocean_xlandinsert_mod: xland table entry "kxland_2" error')
parse_ok = parse(xland_methods(i)%method_control,'tauxland',tauxland(i))
if (parse_ok == 0) call mpp_error(FATAL,'==>Error ocean_xlandinsert_mod: xland table entry "tauxland" error')
enddo
if(Grid%tripolar) then
write(stdoutunit,'(a)') &
'==>Warning: xlandinsert has not been implemented to connect points across the tripolar fold.'
endif
write(stdoutunit, *) &
'Using xlandinsert to connect tracers and surface height between non-local ocean cells.'
xland_halo = min(isc-isd,ied-iec,jsc-jsd,jed-jec)
kxland_min = nk
kxland_max = 1
do nxl=1,nxland
! check for invalid crossland insertion grid locations
if(kxland(nxl,1) > 1) then
call mpp_error(FATAL,&
'==>Error in ocean_xlandinsert_init: kxland set wrong--must start from top cell w/ kxland(nxl,1)=1.')
endif
if (kxland(nxl,1) > kxland(nxl,2)) then
write (unit,994) nxl, kxland(nxl,1), nxl, kxland(nxl,2)
error_xland(nxl) = .true.
endif
if (kxland(nxl,1) < kxland_min ) kxland_min = kxland(nxl,1)
if (kxland(nxl,2) > kxland_max ) kxland_max = kxland(nxl,2)
do lx = 1,2
if (ixland(nxl,lx) < 1 .or. ixland(nxl,lx) > Grd%ni) then
write(unit,991) nxl, lx, ixland(nxl,lx)
error_xland(nxl) = .true.
endif
if (jxland(nxl,lx) < 1 .or. jxland(nxl,lx) > Grd%nj ) then
write(unit,992) nxl, lx, jxland(nxl,lx)
error_xland(nxl) = .true.
endif
if (kxland(nxl,lx) < 1 .or. kxland(nxl,lx) > Grd%nk ) then
write(unit,993) nxl, lx, kxland(nxl,lx)
error_xland(nxl) = .true.
endif
enddo
! determine size for xland_halo
if(Grd%cyclic_x) then
imax = max(ixland(nxl,1),ixland(nxl,2))
imin = min(ixland(nxl,1),ixland(nxl,2))
xland_halo_test = min( abs(imax-imin), abs(imax-Grd%ni-imin))
else
xland_halo_test = abs(ixland(nxl,1)-ixland(nxl,2))
endif
if (xland_halo_test > xland_halo ) xland_halo = xland_halo_test
xland_halo_test = abs(jxland(nxl,1)-jxland(nxl,2))
if (xland_halo_test > xland_halo ) xland_halo = xland_halo_test
! set tauxlandr as inverse of tauxland, except when tauxland=0.
if(tauxland(nxl) > 0.0) tauxlandr(nxl) = 1.0/tauxland(nxl)
enddo ! end of do-loop nxl=1,nxland
if(xland_halo == 0) then
call mpp_error(FATAL,&
'==>Error in ocean_xlandinsert_init: xland_halo=0 => problems defining xlandinsert points.')
endif
if(xland_halo > Dom%xhalo .or. xland_halo > Dom%yhalo) then
write(stdoutunit,*) &
'Defining extra tracer arrays for xland_domain over k-levels ',kxland_min,' through ',kxland_max
endif
! define arrays and xland_domain
if(xland_halo > Dom%xhalo .or. xland_halo > Dom%yhalo) then
write(stdoutunit,*)'The model local computational domain has a x,y halo = ',Dom%xhalo, Dom%yhalo
write(stdoutunit,*)'This is smaller than the halo required for xlandinsert.'
write(stdoutunit,*)'For xlandinsert, define a new domain type with halo = ', xland_halo
if(xland_halo > (Dom%xhalo+6)) then
write(stdoutunit,*)'=>Warning: ocean_xlandinsert_init determined xland_halo = ', xland_halo
write(stdoutunit,*)' Are you sure you wish to connect such distant locations?'
write(stdoutunit,*)' Are you instead trying to connect points across a tripolar fold?'
write(stdoutunit,*)' If so, then presently this has yet to be coded into xlandinsert.'
write(stdoutunit,*)' Alternative xlandinsert points must be taken, or new code implemented.'
endif
endif
! define xlandinsert domain
call set_ocean_domain(Xland_domain,Grd,xhalo=xland_halo,yhalo=xland_halo,&
name='xlandinsert',maskmap=Dom%maskmap)
! time independent arrays defined over the xland_domain
allocate (datx(isc-xland_halo:iec+xland_halo,jsc-xland_halo:jec+xland_halo))
allocate (xtx(isc-xland_halo:iec+xland_halo,jsc-xland_halo:jec+xland_halo))
allocate (ytx(isc-xland_halo:iec+xland_halo,jsc-xland_halo:jec+xland_halo))
datx(:,:) = 0.0
xtx(:,:) = 0.0
ytx(:,:) = 0.0
datx(isc:iec,jsc:jec) = Grd%dat(isc:iec,jsc:jec)
xtx(isc:iec,jsc:jec) = Grd%xt(isc:iec,jsc:jec)
ytx(isc:iec,jsc:jec) = Grd%yt(isc:iec,jsc:jec)
call mpp_update_domains (datx(:,:), Xland_domain%domain2d)
call mpp_update_domains (xtx(:,:), Xland_domain%domain2d)
call mpp_update_domains (ytx(:,:), Xland_domain%domain2d)
! time dependent arrays defined over the xland_domain
allocate (xland_rho_dzt(isc:iec,jsc:jec))
allocate (rho_dzt_x(isc-xland_halo:iec+xland_halo,jsc-xland_halo:jec+xland_halo))
xland_rho_dzt(:,:) = 0.0
rho_dzt_x(:,:) = 0.0
! Mixing across Straits of Gibraltar in CSIRO MK3.5 model resulted in
! exceeding bounds with the Dec2009 verison of MOM4p1.
! ixl needs to be changed to reflect cyclic BCs.
! Earlier checks were made to see if ixl was negative or too
! large wrt the global domain. However we want ixl in the current
! expanded data domain of xtx.
if(Grd%cyclic_x) then
do nxl=1,nxland
if( at_least_one_in_comp_domain(nxl) ) then
do lx = 1,2
if(.not. on_processor(nxl, lx)) then ! Not in domain. This is WRONG!
ixl = ixland(nxl,lx)-Dom%ioff
if(ixl < lbound(xtx,1) ) then ! Gone outside lower index
ixl2 = ixland(nxl,lx) + ieg
write(unit,*) 'xlandinsert init: ixland outside local domain i=',lbound(xtx,1),':',ubound(xtx,1), &
' Changing index from ', ixland(nxl,lx) ,' to ', ixl2
ixland(nxl,lx) = ixl2
else if ( ixl > ubound(xtx,1) ) then ! Gone outside upper index
ixl2 = ixland(nxl,lx) - ieg
write(unit,*) 'xlandinsert init: ixland outside local domain i=',lbound(xtx,1),':',ubound(xtx,1), &
' Changing index from ', ixland(nxl,lx) ,' to ', ixl2
ixland(nxl,lx) = ixl2
endif
if(.not. on_processor(nxl, lx)) then
write(unit,*) 'xlandinsert init?: ixl =', ixl, ' Bounds ',isc+Dom%ioff-xland_halo,':',iec+Dom%ioff+xland_halo
error_xland(nxl)=.true.
endif
endif
enddo
endif
enddo
endif
do nxl=1,nxland
! check for attempts to connect land rather than sea
do lx = 1,2
if(on_comp_domain(nxl, lx)) then
ixl = ixland(nxl,lx)-Dom%ioff
jxl = jxland(nxl,lx)-Dom%joff
if (kxland(nxl,2) > Grid%kmt(ixl,jxl)) then
write (unit,192) nxl, kxland(nxl,2), ixland(nxl,lx),jxland(nxl,lx),Grid%kmt(ixl,jxl)
error_xland(nxl)=.true.
endif
endif
enddo
! write out summary information for this pair of crossland points
if (kxland(nxl,1) == 1) then
ztop = 0.0
else
ztop = Grid%zw(kxland(nxl,1)-1)
endif
if(at_least_one_in_comp_domain(nxl)) then
ixl = ixland(nxl,1)-Dom%ioff
jxl = jxland(nxl,1)-Dom%joff
ixl2 = ixland(nxl,2)-Dom%ioff
jxl2 = jxland(nxl,2)-Dom%joff
if(verbose_init) then
write(unit,191) nxl, ixland(nxl,1), jxland(nxl,1) &
,xtx(ixl,jxl), ytx(ixl,jxl) &
,ixland(nxl,2), jxland(nxl,2) &
,xtx(ixl2,jxl2), ytx(ixl2,jxl2) &
,kxland(nxl,1), kxland(nxl,2), ztop, Grid%zw(kxland(nxl,2)), &
tauxland(nxl)
endif
endif
enddo
! Bring model down here if there are problems with xlandinsert points.
do nxl=1,nxland
if(error_xland(nxl)) then
call mpp_error(FATAL, &
'==>Error detected in ocean_xlandinsert_mod: Use "grep -i error" to find ALL errors.')
endif
enddo
! register for diag_manager
allocate (id_xland(num_prog_tracers))
id_xland = -1
do n=1,num_prog_tracers
if(T_prog(n)%name == 'temp') then
id_xland(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_xlandinsert', &
Grd%tracer_axes(1:3), Time%model_time, 'xlandinsert*cp*rho*dzt*temp', &
'Watt/m^2', missing_value=missing_value, range=(/-1.e10,1.e10/))
else
id_xland(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_xlandinsert', &
Grd%tracer_axes(1:3), Time%model_time, &
'xlandinsert*rho*dzt*tracer for '//trim(T_prog(n)%name), &
trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e10,1.e10/))
endif
enddo
call watermass_diag_init(Time, Dens)
191 format(/' ===== from ocean_xlandinsert_init ====='/ &
' for crossland sea connection pair number',i3,/ &
' mix (i,j) gridpoint (',i4,',',i4,') [long=',f8.3,' lat=',f8.3,']',/ &
' with (i,j) gridpoint (',i4,',',i4,') [long=',f8.3,' lat=',f8.3,']',/ &
' from level',i3,' to',i3,' [depths of ',f10.3,' to ',f10.3,'m]', / &
' Time scale for insertion via upwind advection is (sec) ',f12.6 /)
192 format(/'=>Error: Problem with crossland insertion: improper k-levels requested.'/&
'kxland(',i2,',2) =',i4, ' exceeds kmt(',i4,',',i4,') = ',i4)
991 format(/'=>Error: problem with crossland tracer insertion:',/, &
' out of bounds i grid location requested'/ &
' ixland(',i4,',',i1,') was set incorrectly set to = ',i8,/)
992 format(/'=>Error: problem with crossland tracer insertion:',/, &
' out of bounds j-row grid location requested'/ &
' jxland(',i4,',',i1,') was set incorrectly set to = ',i8,/)
993 format(/'=>Error: problem with crossland tracer insertion:',/, &
' out of bounds k-level grid location requested'/ &
' kxland(',i4,',',i1,') was set incorrectly set to = ',i8,/)
994 format(/'=>Error: problem with crossland tracer insertion:',/, &
' improper k-values requested'/' kxland(',i3,',1)=',i5, &
' and is not less than or equal to kxland(',i3,',2)=',i5,/)
end subroutine ocean_xlandinsert_init
! NAME="ocean_xlandinsert_init"
!#######################################################################
!
!
!
! Compute thickness and density weighted tracer source
! [(kg/m^3)*tracer*m/s] and mass source (kg/m^3)*(m/s)
! associated with discharge of tracer from surface of a
! thick column into a set of boxes within a thin column.
!
! NOTE:
! Compute rho_dzt_x as if using a geopotential model, since the
! xlandinsert rates are typically first tuned with geopotential models.
! also, we wish to have the transport proportional to differences in
! eta_t, which for zstar and pstar are generally must larger than
! differences in dzt.
!
! Modified by Keith.Dixon (January 2004)
!
! Modified by Stephen.Griffies for general
! vertical coordinates January 2005 and August 2006.
!
!
!
subroutine xlandinsert (Time, Ext_mode, Dens, Thickness, T_prog)
type(ocean_time_type), intent(in) :: Time
type(ocean_external_mode_type), intent(in) :: Ext_mode
type(ocean_density_type), intent(in) :: Dens
type(ocean_thickness_type), intent(inout) :: Thickness
type(ocean_prog_tracer_type), intent(inout) :: T_prog(:)
! tracer concentration on stack
real, dimension(isc-xland_halo:iec+xland_halo,jsc-xland_halo:jec+xland_halo,kxland_min:kxland_max,num_prog_tracers) :: tracerx
real, dimension(nxland,0:2) :: top_rho_dzt
real, dimension(nxland,2) :: area
real, dimension(nxland,2) :: mass
real, dimension(nxland) :: rho_dzt_new_large
real, dimension(nxland) :: rho_dzt_new_small
real :: thkocean ! thickness of column over which we insert seawater
real :: delta(nk) ! fraction inserted into a cell
real :: tracernew(nk)
real :: blanda, tracerextra, zextra, zinsert
real :: tottracer_thin, tottracer_thick
real :: totrate_thin, totrate_thick
integer :: i, j, k, nxl, nz
integer :: tau
integer :: n
integer :: lx
integer :: lx_thin(max(1,nxland))
integer :: lx_thick(max(1,nxland))
if (.not. module_is_initialized ) then
call mpp_error(FATAL, &
'==>Error in ocean_xlandinsert_mod (xlandinsert): module must be initialized')
endif
if (nxland < 1) return
if (.not. use_this_module) return
tau = Time%tau
! tracer independent computations
xland_rho_dzt(:,:) = 0.0
rho_dzt_x(:,:) = 0.0
tracerx(:,:,:,:) = 0.0
if(xland_halo > Dom%xhalo .or. xland_halo > Dom%yhalo) then
rho_dzt_x(isc:iec,jsc:jec) = rho0*(Grd%dzt(1)+ext_mode%eta_t(isc:iec,jsc:jec,tau))
call mpp_update_domains (rho_dzt_x(:,:), Xland_domain%domain2d)
do n=1,num_prog_tracers
do k=1,kxland_max
do j=jsc,jec
do i=isc,iec
tracerx(i,j,k,n) = T_prog(n)%field(i,j,k,tau)
enddo
enddo
enddo
enddo
call mpp_update_domains (tracerx(:,:,:,:), Xland_domain%domain2d)
else
rho_dzt_x(:,:) = rho0*(Grd%dzt(1)+ext_mode%eta_t(:,:,tau))
do n=1,num_prog_tracers
do k=1,kxland_max
tracerx(:,:,k,n) = T_prog(n)%field(:,:,k,tau)
enddo
enddo
endif
do nxl=1,nxland
! calculate only if at least one point in the pair is on the comp-domain
if(at_least_one_in_comp_domain(nxl)) then
! geometric factors
do lx=1,2
ixl = ixland(nxl,lx)-Dom%ioff
jxl = jxland(nxl,lx)-Dom%joff
top_rho_dzt(nxl,lx) = rho_dzt_x(ixl,jxl)
area(nxl,lx) = datx(ixl,jxl)
mass(nxl,lx) = top_rho_dzt(nxl,lx)*area(nxl,lx)
enddo
if(top_rho_dzt(nxl,1) > top_rho_dzt(nxl,2)) then
lx_thick(nxl) = 1
lx_thin(nxl) = 2
else
lx_thick(nxl) = 2
lx_thin(nxl) = 1
endif
! transfer rate
top_rho_dzt(nxl,0) = (mass(nxl,1) + mass(nxl,2))/(area(nxl,1) + area(nxl,2))
water_rate(nxl,lx_thick(nxl)) = (top_rho_dzt(nxl,lx_thick(nxl))-top_rho_dzt(nxl,0))*tauxlandr(nxl)
water_rate(nxl,lx_thin(nxl)) = (top_rho_dzt(nxl,0)-top_rho_dzt(nxl,lx_thin(nxl))) *tauxlandr(nxl)
if(debug_this_module) then
write(*,*)' '
write(*,*)'For xlandinsert pair = ',nxl
write(*,*)'top_rho_dzt(',lx_thick(nxl),') (kg/m^2) = ',top_rho_dzt(nxl,lx_thick(nxl)),' (nxl=',nxl,')'
write(*,*)'top_rho_dzt(',lx_thin(nxl),') (kg/m^2) = ',top_rho_dzt(nxl,lx_thin(nxl)),' (nxl=',nxl,')'
write(*,*)'top_rho_dzt( 0 ) (kg/m^2) = ',top_rho_dzt(nxl,0),' (nxl=',nxl,')'
write(*,*)'water_rate_thin (kg/m^3)*(m/s) = ',water_rate(nxl,lx_thin(nxl)),' (nxl=',nxl,')'
write(*,*)'water_rate_thick (kg/m^3)*(m/s) = ',water_rate(nxl,lx_thick(nxl)),' (nxl=',nxl,')'
write(*,*)'area_thin (m^2) = ',area(nxl,lx_thin(nxl)),' (nxl=',nxl,')'
write(*,*)'area_thick (m^2) = ',area(nxl,lx_thick(nxl)),' (nxl=',nxl,')'
write(*,*) &
'{[area*water_rate](thick) - [area*water_rate](thin)}/(area(thin)+area(thick)) (kg/m^3)*(m/s) = ', &
( water_rate(nxl,lx_thin(nxl))*area(nxl,lx_thin(nxl)) &
-water_rate(nxl,lx_thick(nxl))*area(nxl,lx_thick(nxl))) &
/(area(nxl,lx_thin(nxl)) + area(nxl,lx_thick(nxl))),' (nxl=',nxl,')'
write(*,*) &
' >> check mass rates (thin) ', water_rate(nxl,lx_thin(nxl))*area(nxl,lx_thin(nxl)),' (nxl=',nxl,')'
write(*,*) &
' >> check mass rates (thick) ', water_rate(nxl,lx_thick(nxl))*area(nxl,lx_thick(nxl)),' (nxl=',nxl,')'
write(*,*)' Insert ',water_rate(nxl,lx_thin(nxl))*dtime*area(nxl,lx_thin(nxl)),'kg', &
' of water over ',dtime,' sec timestep (nxl=',nxl,')'
if(on_comp_domain(nxl,lx_thick(nxl))) then
write(*,*)' the thick point is on the computational domain (nxl=',nxl,')'
else
write(*,*)' the thick point is NOT on the computational domain (nxl=',nxl,')'
endif
if(on_comp_domain(nxl,lx_thin(nxl))) then
write(*,*)' the thin point is on the computational domain (nxl=',nxl,')'
else
write(*,*)' the thin point is NOT on the computational domain (nxl=',nxl,')'
endif
endif
! compute tendency to remove water from region of more mass to
! region with less mass.
!
! xland_rho_dzt is added to Thickness%mass_source
!
! rho_dzt_new_large is the estimated tau+1 value of rho_dzt if
! the xland_insert adjustment was the only process contributing
! to the tendency
!
! note that the more massive side can be associated with more that one xland pair,
! so the tendency is accumulated. Also, this is only calculated if the massive
! side is a computational point.
lx=lx_thick(nxl)
if(on_comp_domain(nxl,lx)) then
ixl = ixland(nxl,lx)-Dom%ioff
jxl = jxland(nxl,lx)-Dom%joff
xland_rho_dzt(ixl,jxl) = xland_rho_dzt(ixl,jxl) - water_rate(nxl,lx)
rho_dzt_new_large(nxl) = rho_dzt_x(ixl,jxl) - dtime*water_rate(nxl,lx)
if(debug_this_module) then
write(*,*) ' BEFORE: (thick) rho_dzt=',rho_dzt_x(ixl,jxl),' (nxl=',nxl,')'
write(*,*) ' BEFORE: (thick) sfc mass=',(area(nxl,lx)*top_rho_dzt(nxl,lx)),' (nxl=',nxl,')'
write(*,*) ' >> check sfc mass above =', mass(nxl,lx),' (nxl=',nxl,')'
write(*,*) ' PREDICT: (thick) mass shift per timestep=',(-water_rate(nxl,lx)*area(nxl,lx)*dtime), &
' (nxl=',nxl,')'
write(*,*) ' PREDICT: rho_dzt_new_large=', rho_dzt_new_large(nxl), &
' (nxl=',nxl,')'
write(*,*) ' PREDICT: (thick) sfc mass(tau+1)=',( area(nxl,lx)* &
(top_rho_dzt(nxl,lx)-(water_rate(nxl,lx)*dtime)) ), ' (nxl=',nxl,')'
endif
endif
! compute tendency to add water to region of less mass from region of more mass.
!
! xland_rho_dzt is added to Thickness%mass_source
!
! rho_dzt_new_small is the estimated tau+1 value of rho_dzt if
! the xland_insert adjustment was the only process contributing
! to the tendency
!
! note that the less massive side can be associated with more that one xland pair,
! so the tendency is accumulated. Also, this is only calculated if the less massive
! side is a computational point.
lx = lx_thin(nxl)
if(on_comp_domain(nxl,lx)) then
ixl = ixland(nxl,lx)-Dom%ioff
jxl = jxland(nxl,lx)-Dom%joff
xland_rho_dzt(ixl,jxl) = xland_rho_dzt(ixl,jxl) + water_rate(nxl,lx)
rho_dzt_new_small(nxl) = rho_dzt_x(ixl,jxl) + dtime*water_rate(nxl,lx_thin(nxl))
if(debug_this_module) then
write(*,*) ' BEFORE: (thin) rho_dzt=',rho_dzt_x(ixl,jxl),' (nxl=',nxl,')'
write(*,*) ' BEFORE: (thin) sfc mass=',(area(nxl,lx)*top_rho_dzt(nxl,lx)),' (nxl=',nxl,')'
write(*,*) ' PREDICT: (thin) mass shift per timestep= ',(water_rate(nxl,lx)*area(nxl,lx)*dtime), &
' (nxl=',nxl,')'
write(*,*) ' PREDICT: rho_dzt_new_small =', rho_dzt_new_small(nxl), &
' (nxl=',nxl,')'
write(*,*) ' PREDICT: (thin) sfc mass(tau+1)=',(area(nxl,lx)* &
(top_rho_dzt(nxl,lx)+(water_rate(nxl,lx)*dtime))), ' (nxl=',nxl,')'
endif
endif
endif ! end of at_least_one_in_comp_domain(nxl) if-test
enddo ! end of the nxland do-loop
! do the tracers
do n=1,num_prog_tracers
wrk1(:,:,:) = 0.0
! loop over nxland pairs
do nxl=1,nxland
! perform calculations only if at least one point in the pair is on comp-domain
if(at_least_one_in_comp_domain(nxl)) then
! remove tracer from top of the thick column via upwind advection
lx=lx_thick(nxl)
k =1 ! kxland(nxl,1)=1 for xlandinsert
blanda = 0.0
if(on_comp_domain(nxl,lx)) then
ixl = ixland(nxl,lx)-Dom%ioff
jxl = jxland(nxl,lx)-Dom%joff
blanda = -1.0*water_rate(nxl,lx)*tracerx(ixl,jxl,k,n)
wrk1(ixl,jxl,k) = wrk1(ixl,jxl,k) + blanda
if(debug_this_module) then
! let totrate_thick keep track of the mass weighted tendency of the tracer content and
! tottracer_thick keep track of the total tracer content
tottracer_thick = tracerx(ixl,jxl,k,n)*top_rho_dzt(nxl,lx)*area(nxl,lx)
totrate_thick = blanda*area(nxl,lx)
write(*,*)' '
write(*,*) 'BEFORE: tracer value = ',tracerx(ixl,jxl,k,n), &
'(nxl=',nxl,' tracer#=',n,' k=',k,' thick)'
write(*,*) 'BEFORE: thick mass = ',(top_rho_dzt(nxl,lx)*area(nxl,lx)),' = ',&
top_rho_dzt(nxl,lx),'*',area(nxl,lx), &
'(nxl=',nxl,' tracer#=',n,' k=',k,' thick)'
write(*,*) 'BEFORE: tottracer_thick = ',tottracer_thick, &
'(nxl=',nxl,' tracer#=',n,' k=',k,')'
write(*,*) 'xland_insert totrate_thick = ',totrate_thick, &
'(nxl=',nxl,' tracer#=',n,' k=',k,' thick)'
write(*,*) 'PREDICT AFTER: thick mass = ', area(nxl,lx) * &
( top_rho_dzt(nxl,lx) - (water_rate(nxl,lx)*dtime) ),' (nxl=',nxl,' tracer#=',n,')'
write(*,*) 'PREDICT AFTER: tottracer_thick = ',( tracerx(ixl,jxl,k,n) * area(nxl,lx) * &
( top_rho_dzt(nxl,lx) - (water_rate(nxl,lx)*dtime) ) ),' (nxl=',nxl,' tracer#=',n,')'
endif
endif
lx = lx_thin(nxl)
if(on_comp_domain(nxl,lx)) then
nz = kxland(nxl,2) ! kxland(nxl,2) is the deepest insertion level
ixl = ixland(nxl,lx) -Dom%ioff
jxl = jxland(nxl,lx) -Dom%joff
ixl2 = ixland(nxl,3-lx)-Dom%ioff
jxl2 = jxland(nxl,3-lx)-Dom%joff
! compute fraction inserted into each thin grid cell
thkocean = 0.0
do k=1,nz
thkocean = thkocean + Thickness%rho_dzt(ixl,jxl,k,tau)
enddo
do k=1,nz
delta(k) = Thickness%rho_dzt(ixl,jxl,k,tau)/(epsln+thkocean)
enddo
if(debug_this_module) then
tottracer_thin = 0.0
do k=1,nz
if (k.eq. 1) then
write(*,*) k, tracerx(ixl,jxl,k,n), tracerx(ixl2,jxl2,k,n), ' (nxl=',nxl,' tracer#=',n,')' &
,' thin Thickness%rho_dzt =',Thickness%rho_dzt(ixl,jxl,k,tau)
else
write(*,*) k, tracerx(ixl,jxl,k,n), ' (nxl=',nxl,' tracer#=',n,')' &
,' thin Thickness%rho_dzt =',Thickness%rho_dzt(ixl,jxl,k,tau)
endif
! for tottracer_thin, add in tracer content of destination (thin) column k-levels
tottracer_thin = tottracer_thin + &
( tracerx(ixl,jxl,k,n)*Thickness%rho_dzt(ixl,jxl,k,tau)*area(nxl,lx) )
enddo
write(*,*)' Does Thickness%rho_dzt(ixl,jxl,1,tau) = ',Thickness%rho_dzt(ixl,jxl,1,tau), &
' = ',top_rho_dzt(nxl,lx),' = top_rho_dzt(xnl,lx) ?', &
'(nxl=',nxl,' tracer#=',n,')'
write(*,*)' BEFORE xland_insert: tottracer_thin = ',tottracer_thin, &
' (nxl=',nxl,' tracer#=',n,')'
write(*,*)' '
write(*,*)' tracernew(k) = ((tracerextra*zextra) +', &
' (tracerx(ixl,jxl,k,n)*Thickness%rho_dzt(ixl,jxl,k,tau)) +', &
' (tracerx(ixl2,jxl2,1,n)*zinsert) )', &
' / (zextra + Thickness%rho_dzt(ixl,jxl,k,tau) + zinsert)'
endif
! Insert tracer into thin column by doing an xland_insert adjustment process.
! and then diagnose the effective tendencies from the tau+1 - tau timelevel values.
! Start at the bottom (k=nz) of the destination (thin) column and compute what the
! tracer value would be at the end of the timestep if xland_insert was the only
! process contributing to the tracer tendency. Then move up the thin column, one
! k-level at a time, and at each k-level from k=nz-1 up to k=1 compute what the updated
! tracer value would be accounting for both the insertion of the water from the top of
! the thick column +and+ the bubbling up of the excess mass of water from the level
! below. When the process is done, compare the updated tracer values with the current
! tracer values and convert back into a tendency.
! zextra is the vertical thickness of the excess water bubbling up from the
! layer below due to the accumulated effects of xland_insert
! tracerextra is the tracer value at the tau+1 time step of the water bubbling
! up from the layer below
! zinsert is the surface height change that will occur due to xland_insert
! at the thin (destination) point over the course of this tracer timestep.
! tracernew would be the tau+1 time level tracer value value if the xland_insert
! adjustment was the only process contributing to the tendencies
zextra = 0.0
do k=nz,1,-1
tracernew(k) = 0.0
if (k .EQ. nz) then
tracerextra = 0.0
else
tracerextra = tracernew(k+1)
endif
zinsert = water_rate(nxl,lx)*dtime*delta(k)
tracernew(k) = ((tracerextra*zextra) + &
(tracerx(ixl,jxl,k,n)*Thickness%rho_dzt(ixl,jxl,k,tau)) + &
(tracerx(ixl2,jxl2,1,n)*zinsert) ) &
/ (zextra + Thickness%rho_dzt(ixl,jxl,k,tau) + zinsert)
if(debug_this_module) then
write(*,*)' k=',k, tracernew(k),' (nxl=',nxl,' tracer#=',n,')'
write(*,*)'[nxl=',nxl,'] ',tracernew(k),' = ((',tracerextra,'*',zextra,') +'
write(*,*)'[nxl=',nxl,'] (',tracerx(ixl,jxl,k,n),'*',Thickness%rho_dzt(ixl,jxl,k,tau),') +'
write(*,*)'[nxl=',nxl,'] (',tracerx(ixl2,jxl2,1,n),'*',zinsert,') )'
write(*,*)'[nxl=',nxl,'] / (',zextra,' + ',Thickness%rho_dzt(ixl,jxl,k,tau),' + ',zinsert,')'
endif
zextra = zextra + zinsert
enddo
if(debug_this_module) then
write(*,*)' '
write(*,*) ' CHECK: (thin) sfc mass(tau+1)=',area(nxl,lx)*rho_dzt_new_small(nxl), &
' (nxl=',nxl,')'
write(*,*) ' CHECK: (thin) Does zextra =',zextra,' = ', &
'(rho_dzt_new_small(nxl) - rho_dzt_x(ixl,jxl)) = ',(rho_dzt_new_small(nxl) - rho_dzt_x(ixl,jxl)), &
'(dtime*water_rate(nxl,lx_thin(nxl)))=',(dtime*water_rate(nxl,lx_thin(nxl)))
! accumulate tottracer_thin at tau+1
tottracer_thin = 0.0
tottracer_thin = tracernew(1)*area(nxl,lx)* &
(Thickness%rho_dzt(ixl,jxl,1,tau)-rho_dzt_x(ixl,jxl)+rho_dzt_new_small(nxl))
if(debug_this_module) then
write(*,*) ' >nxl=',nxl,' tracer#=',n,') tottracer_thin = ', tottracer_thin,' = '
write(*,*) ' >nxl=',nxl,' ',tracernew(1),'*',area(nxl,lx),'*'
write(*,*) ' >nxl=',nxl,' (',Thickness%rho_dzt(ixl,jxl,1,tau),'-',rho_dzt_x(ixl,jxl),&
'+',rho_dzt_new_small(nxl),')'
write(*,*) ' >nxl=',nxl,' which is [',(Thickness%rho_dzt(ixl,jxl,1,tau)-rho_dzt_x(ixl,jxl)&
+rho_dzt_new_small(nxl)),']'
endif
do k=1,nz
write(*,*) k, tracernew(k), tracerx(ixl2,jxl2,k,n)
if (k .NE. 1) then
tottracer_thin = tottracer_thin + &
(tracernew(k)*area(nxl,lx)*Thickness%rho_dzt(ixl,jxl,k,tau))
if(debug_this_module) then
write(*,*) ' >nxl=',nxl,' tottracer_thin = ', tottracer_thin,' = tottracer_thin'
write(*,*) ' >nxl=',nxl,' (',tracernew(k),'*',area(nxl,lx),'*',Thickness%rho_dzt(ixl,jxl,k,tau),')'
endif
endif
enddo
write(*,*)' AFTER xland_insert: tottracer_thin = ',tottracer_thin, &
' (nxl=',nxl,' tracer#=',n,')'
endif
! convert change in tracer amounts between tau a tau+1 timesteps
! into tendencies and store into the wrk1 array
blanda = 0.0
do k=1,nz
if (k .eq. 1) then
blanda = ( ( tracernew(1)*(Thickness%rho_dzt(ixl,jxl,1,tau)-rho_dzt_x(ixl,jxl)+rho_dzt_new_small(nxl)) ) &
- ( tracerx(ixl,jxl,1,n)*Thickness%rho_dzt(ixl,jxl,1,tau) ) )/dtime
else
blanda = Thickness%rho_dzt(ixl,jxl,k,tau)*(tracernew(k) - tracerx(ixl,jxl,k,n))/dtime
endif
wrk1(ixl,jxl,k) = wrk1(ixl,jxl,k) + blanda
if (debug_this_module) then
if (k .eq. 1) then
totrate_thin = 0.0
endif
totrate_thin = totrate_thin + (blanda*area(nxl,lx))
write(*,*) 'xland_insert tendency = ',blanda,'(nxl=',nxl,' tracer#=',n,' k=',k,' thin)'
endif
end do
if (debug_this_module) then
write(*,*) ' '
write(*,*) 'xland_insert totrate_thin = ',totrate_thin , &
'(nxl=',nxl,' tracer#=',n,' thin )'
endif
endif ! end of on_comp_domain
endif ! end of at_least_one_in_comp_domain(nxl)
enddo ! end of nxland
! fill tendency array
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) + wrk1(i,j,k)
enddo
enddo
enddo
! send diagnostics
if(id_xland(n) > 0) then
call diagnose_3d(Time, Grd, id_xland(n), T_prog(n)%conversion*wrk1(:,:,:))
endif
! for watermass diagnostics
if(n==index_temp) then
theta_tend(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
theta_tend(i,j,k) = wrk1(i,j,k)
enddo
enddo
enddo
endif
if(n==index_salt) then
salt_tend(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
salt_tend(i,j,k) = wrk1(i,j,k)
enddo
enddo
enddo
endif
enddo ! end of num_prog_tracers do-loop
! fill mass source array
call xlandinsert_mass(Time, Thickness)
! call some diagnostics
call watermass_diag(Time, Dens)
end subroutine xlandinsert
! NAME="xlandinsert"
!#######################################################################
!
!
!
! Compute the mass source (kg/m^3)*meter/sec.
!
! Note that xlandinsert has already been called, so xland_rho_dzt
! has been filled.
!
!
subroutine xlandinsert_mass (Time, Thickness)
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(inout) :: Thickness
integer :: i,j
if (.not. module_is_initialized ) then
call mpp_error(FATAL, &
'==>Error in ocean_xlandinsert_mod (xlandinsert_mass): module must be initialized')
endif
do j=jsc,jec
do i=isc,iec
Thickness%mass_source(i,j,1) = Thickness%mass_source(i,j,1) + xland_rho_dzt(i,j)
enddo
enddo
call mpp_update_domains (Thickness%mass_source(:,:,1), Dom%domain2d)
call diagnose_2d(Time, Grd, id_xland_rho_dzt, xland_rho_dzt(:,:))
end subroutine xlandinsert_mass
! NAME="xlandinsert_mass"
!#######################################################################
!
!
!
! Function to see if at least one of the two points in a crossland pair
! is within the computational domain for the processor.
!
!
!
! Integer labeling the particular xlandinsert pair
!
!
function at_least_one_in_comp_domain(nxl)
integer, intent(in) :: nxl
integer :: lx
logical :: in_domain(2), at_least_one_in_comp_domain
do lx=1,2
if(isc+Dom%ioff <= ixland(nxl,lx) .and. ixland(nxl,lx) <= iec+Dom%ioff .and. &
jsc+Dom%joff <= jxland(nxl,lx) .and. jxland(nxl,lx) <= jec+Dom%joff) then
in_domain(lx) = .true.
else
in_domain(lx) = .false.
endif
enddo
at_least_one_in_comp_domain = .false.
if(in_domain(1) .and. in_domain(2) ) at_least_one_in_comp_domain = .true.
if(in_domain(1) .and. .not.in_domain(2)) at_least_one_in_comp_domain = .true.
if(.not.in_domain(1) .and. in_domain(2)) at_least_one_in_comp_domain = .true.
end function at_least_one_in_comp_domain
! NAME="at_least_one_in_comp_domain"
!#######################################################################
!
!
!
! Determine if the point is in comp-domain for the processor
!
!
!
! Integer labeling the particular xlandinsert pair
!
!
! lx=1,2 labels the point within an xlandinsert pair
!
!
function on_comp_domain(nxl, lx)
integer, intent(in) :: nxl, lx
logical :: on_comp_domain
if(isc+Dom%ioff <= ixland(nxl,lx) .and. ixland(nxl,lx) <= iec+Dom%ioff .and. &
jsc+Dom%joff <= jxland(nxl,lx) .and. jxland(nxl,lx) <= jec+Dom%joff) then
on_comp_domain = .true.
else
on_comp_domain = .false.
endif
end function on_comp_domain
! NAME="on_comp_domain"
!#######################################################################
!
!
!
! Determine if the point is on processor
!
!
function on_processor(nxl, lx)
integer, intent(in) :: nxl, lx
logical :: on_processor
if(isc+Dom%ioff-xland_halo <= ixland(nxl,lx) .and. ixland(nxl,lx) <= iec+Dom%ioff+xland_halo .and. &
jsd+Dom%joff-xland_halo <= jxland(nxl,lx) .and. jxland(nxl,lx) <= jed+Dom%joff+xland_halo) then
on_processor = .true.
else
on_processor = .false.
endif
end function on_processor
! NAME="on_processor"
!#######################################################################
!
!
!
! Initialization of watermass diagnostic output files.
!
!
subroutine watermass_diag_init(Time, Dens)
type(ocean_time_type), intent(in) :: Time
type(ocean_density_type), intent(in) :: Dens
integer :: stdoutunit
stdoutunit=stdout()
compute_watermass_diag = .false.
id_neut_rho_xinsert = register_diag_field ('ocean_model', 'neut_rho_xinsert',&
Grd%tracer_axes(1:3), Time%model_time, &
'update of locally ref potrho from xlandinsert', &
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_rho_xinsert > 0) compute_watermass_diag = .true.
id_neut_rho_xinsert_on_nrho = register_diag_field ('ocean_model', &
'neut_rho_xinsert_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'update of locally ref potrho from xlandinsert as binned to neutral density layers',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_rho_xinsert_on_nrho > 0) compute_watermass_diag = .true.
id_wdian_rho_xinsert = register_diag_field ('ocean_model', 'wdian_rho_xinsert',&
Grd%tracer_axes(1:3), Time%model_time, &
'dianeutral mass transport due to xlandinsert', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_rho_xinsert > 0) compute_watermass_diag = .true.
id_wdian_rho_xinsert_on_nrho = register_diag_field ('ocean_model', &
'wdian_rho_xinsert_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'dianeutral mass transport from xlandinsert as binned to neutral density layers',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_rho_xinsert_on_nrho > 0) compute_watermass_diag = .true.
id_tform_rho_xinsert = register_diag_field ('ocean_model', &
'tform_rho_xinsert', Grd%tracer_axes(1:3), Time%model_time,&
'watermass transform from xlandinsert', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_rho_xinsert_on_nrho > 0) compute_watermass_diag = .true.
id_tform_rho_xinsert_on_nrho = register_diag_field ('ocean_model', &
'tform_rho_xinsert_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'watermass transform from xlandinsert as binned to neutral density layers',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_rho_xinsert_on_nrho > 0) compute_watermass_diag = .true.
id_eta_tend_xinsert= -1
id_eta_tend_xinsert= register_diag_field ('ocean_model','eta_tend_xinsert',&
Grd%tracer_axes(1:2), Time%model_time, &
'non-Bouss steric sea level tendency from xinsert tendency', 'm/s', &
missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_eta_tend_xinsert > 0) compute_watermass_diag=.true.
id_eta_tend_xinsert_glob= -1
id_eta_tend_xinsert_glob= register_diag_field ('ocean_model', 'eta_tend_xinsert_glob',&
Time%model_time, &
'global mean non-bouss steric sea level tendency from xinsert tendency', &
'm/s', missing_value=missing_value, range=(/-1e10,1.e10/))
if(id_eta_tend_xinsert_glob > 0) compute_watermass_diag=.true.
! temperature effects
id_neut_temp_xinsert = register_diag_field ('ocean_model', 'neut_temp_xinsert',&
Grd%tracer_axes(1:3), Time%model_time, &
'temp related update of locally ref potrho from xlandinsert', &
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_temp_xinsert > 0) compute_watermass_diag = .true.
id_neut_temp_xinsert_on_nrho = register_diag_field ('ocean_model', &
'neut_temp_xinsert_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related update of locally ref potrho from xlandinsert as binned to neutral density layers',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_temp_xinsert_on_nrho > 0) compute_watermass_diag = .true.
id_wdian_temp_xinsert = register_diag_field ('ocean_model', 'wdian_temp_xinsert',&
Grd%tracer_axes(1:3), Time%model_time, &
'temp related dianeutral mass transport due to xlandinsert', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_temp_xinsert > 0) compute_watermass_diag = .true.
id_wdian_temp_xinsert_on_nrho = register_diag_field ('ocean_model', &
'wdian_temp_xinsert_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related dianeutral mass transport from xlandinsert as binned to neutral density layers',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_temp_xinsert_on_nrho > 0) compute_watermass_diag = .true.
id_tform_temp_xinsert = register_diag_field ('ocean_model', &
'tform_temp_xinsert', Grd%tracer_axes(1:3), Time%model_time,&
'temp related watermass transform from xlandinsert', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_temp_xinsert_on_nrho > 0) compute_watermass_diag = .true.
id_tform_temp_xinsert_on_nrho = register_diag_field ('ocean_model', &
'tform_temp_xinsert_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'temp related watermass transform from xlandinsert as binned to neutral density layers',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_temp_xinsert_on_nrho > 0) compute_watermass_diag = .true.
! salinity effects
id_neut_salt_xinsert = register_diag_field ('ocean_model', 'neut_salt_xinsert',&
Grd%tracer_axes(1:3), Time%model_time, &
'salt related update of locally ref potrho from xlandinsert', &
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_salt_xinsert > 0) compute_watermass_diag = .true.
id_neut_salt_xinsert_on_nrho = register_diag_field ('ocean_model', &
'neut_salt_xinsert_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related update of locally ref potrho from xlandinsert as binned to neutral density layers',&
'(kg/m^3)/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_neut_salt_xinsert_on_nrho > 0) compute_watermass_diag = .true.
id_wdian_salt_xinsert = register_diag_field ('ocean_model', 'wdian_salt_xinsert',&
Grd%tracer_axes(1:3), Time%model_time, &
'salt related dianeutral mass transport due to xlandinsert', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_salt_xinsert > 0) compute_watermass_diag = .true.
id_wdian_salt_xinsert_on_nrho = register_diag_field ('ocean_model', &
'wdian_salt_xinsert_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related dianeutral mass transport from xlandinsert as binned to neutral density layers',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_wdian_salt_xinsert_on_nrho > 0) compute_watermass_diag = .true.
id_tform_salt_xinsert = register_diag_field ('ocean_model', &
'tform_salt_xinsert', Grd%tracer_axes(1:3), Time%model_time,&
'salt related watermass transform from xlandinsert', &
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_salt_xinsert_on_nrho > 0) compute_watermass_diag = .true.
id_tform_salt_xinsert_on_nrho = register_diag_field ('ocean_model', &
'tform_salt_xinsert_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, &
'salt related watermass transform from xlandinsert as binned to neutral density layers',&
'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/))
if(id_tform_salt_xinsert_on_nrho > 0) compute_watermass_diag = .true.
allocate(theta_tend(isd:ied,jsd:jed,nk))
allocate(salt_tend(isd:ied,jsd:jed,nk))
theta_tend(:,:,:) = 0.0
salt_tend(:,:,:) = 0.0
if(compute_watermass_diag) then
write(stdoutunit,'(/a/)') &
'==>Note: running ocean_xlandinsert_mod w/ compute_watermass_diag=.true.'
endif
end subroutine watermass_diag_init
! NAME="watermass_diag_init"
!#######################################################################
!
!
!
! Diagnose effects from xlandinsert on watermass transformation.
!
!
subroutine watermass_diag(Time, Dens)
type(ocean_time_type), intent(in) :: Time
type(ocean_density_type), intent(in) :: Dens
integer :: i,j,k,tau
real, dimension(isd:ied,jsd:jed) :: eta_tend
if (.not.module_is_initialized) then
call mpp_error(FATAL, &
'==>Error from ocean_xlandinsert (watermass_diag): module needs initialization ')
endif
if(.not. compute_watermass_diag) return
tau = Time%tau
! full effects from temperature and salinity
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
wrk5(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Dens%drhodT(i,j,k)*theta_tend(i,j,k) &
+Dens%drhodS(i,j,k)*salt_tend(i,j,k)
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
wrk5(i,j,k) =-wrk1(i,j,k)/(epsln+Dens%rho(i,j,k,tau)**2) ! for eta_tend
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_rho_xinsert, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_rho_xinsert, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_rho_xinsert, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_rho_xinsert_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_rho_xinsert_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_rho_xinsert_on_nrho, wrk4)
if(id_eta_tend_xinsert > 0 .or. id_eta_tend_xinsert_glob > 0) then
eta_tend(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
eta_tend(i,j) = eta_tend(i,j) + wrk5(i,j,k)
enddo
enddo
enddo
call diagnose_2d(Time, Grd, id_eta_tend_xinsert, eta_tend(:,:))
call diagnose_sum(Time, Grd, Dom, id_eta_tend_xinsert_glob, eta_tend, cellarea_r)
endif
! effects from temperature
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Dens%drhodT(i,j,k)*theta_tend(i,j,k)
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_temp_xinsert, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_temp_xinsert, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_temp_xinsert, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_temp_xinsert_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_temp_xinsert_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_temp_xinsert_on_nrho, wrk4)
! effects from salinity
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = Dens%drhodS(i,j,k)*salt_tend(i,j,k)
wrk2(i,j,k) = wrk1(i,j,k)*Dens%rho_dztr_tau(i,j,k)
wrk3(i,j,k) = wrk2(i,j,k)*Dens%stratification_factor(i,j,k)
wrk4(i,j,k) = wrk1(i,j,k)*Grd%dat(i,j)*Dens%watermass_factor(i,j,k)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_neut_salt_xinsert, wrk2(:,:,:))
call diagnose_3d(Time, Grd, id_wdian_salt_xinsert, wrk3(:,:,:))
call diagnose_3d(Time, Grd, id_tform_salt_xinsert, wrk4(:,:,:))
call diagnose_3d_rho(Time, Dens, id_neut_salt_xinsert_on_nrho, wrk2)
call diagnose_3d_rho(Time, Dens, id_wdian_salt_xinsert_on_nrho, wrk3)
call diagnose_3d_rho(Time, Dens, id_tform_salt_xinsert_on_nrho, wrk4)
end subroutine watermass_diag
! NAME="watermass_diag"
end module ocean_xlandinsert_mod