!=======================================================================
!BOP
!
! !MODULE: ice_transport_remap - horizontal transport via incremental remapping
!
! !DESCRIPTION:
!
! Transports quantities using the second-order conservative remapping
! scheme developed by John Dukowicz and John Baumgardner (DB) and modified
! for sea ice by William Lipscomb and Elizabeth Hunke.
!
! References:
!
! Dukowicz, J. K., and J. R. Baumgardner, 2000: Incremental
! remapping as a transport/advection algorithm, J. Comput. Phys.,
! 160, 318-335.
!
! Lipscomb, W. H., and E. C. Hunke, 2004: Modeling sea ice
! transport using incremental remapping, Mon. Wea. Rev., 132,
! 1341-1354.
!
! !REVISION HISTORY:
! SVN:$Id: ice_transport_remap.F 33 2006-11-13 19:51:14Z eclare $
!
! authors William H. Lipscomb, LANL
! John Baumgardner, LANL
!
! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
! 2004-05: Block structure added (WHL)
! 2006: Moved remap driver to ice_transport_driver
! Geometry changes:
! (1) Reconstruct fields in stretched logically rectangular coordinates
! (2) Modify geometry so that the area flux across each edge
! can be specified (following an idea of Mats Bentsen)
!
! !INTERFACE:
!
module ice_transport_remap 5,5
!
! !USES:
!
use ice_kinds_mod
use ice_communicate, only: my_task, master_task, MPI_COMM_ICE
use ice_domain_size
use ice_constants
use ice_fileunits, only: nu_diag
use perf_mod, only: t_startf, t_stopf, t_barrierf
!
!EOP
!
implicit none
save
private
public :: init_remap, horizontal_remap, make_masks
integer (kind=int_kind), parameter :: &
max_ntrace = 2+max_ntrcr+nilyr+nslyr ! hice,hsno,qice,qsno,trcr
integer (kind=int_kind), parameter :: &
ngroups = 6 ,&! number of groups of triangles that
! contribute transports across each edge
nvert = 3 ! number of vertices in a triangle
! for triangle integral formulas
real (kind=dbl_kind), parameter :: &
p5625m = -9._dbl_kind/16._dbl_kind ,&
p52083 = 25._dbl_kind/48._dbl_kind
logical (kind=log_kind), parameter :: bugcheck = .false.
!=======================================================================
! Here is some information about how the incremental remapping scheme
! works in CICE and how it can be adapted for use in other models.
!
! The remapping routine is designed to transport a generic mass-like
! field (in CICE, the ice fractional area) along with an arbitrary number
! of tracers in two dimensions. The velocity components are assumed
! to lie at grid cell corners and the transported scalars at cell centers.
! Incremental remapping has the following desirable properties:
!
! (1) Tracer monotonicity is preserved. That is, no new local
! extrema are produced in fields like ice thickness or internal
! energy.
! (2) The reconstucted mass and tracer fields vary linearly in x and y.
! This means that remapping is 2nd-order accurate in space,
! except where horizontal gradients are limited to preserve
! monotonicity.
! (3) There are economies of scale. Transporting a single field
! is rather expensive, but additional fields have a relatively
! low marginal cost.
!
! The following generic conservation equations may be solved:
!
! dm/dt = del*(u*m) (0)
! d(m*T1)/dt = del*(u*m*T1) (1)
! d(m*T1*T2)/dt = del*(u*m*T1*T2) (2)
! d(m*T1*T2*T3)/dt = del*(u*m*T1*T2*T3) (3)
!
! where d is a partial derivative, del is the 2D divergence operator,
! u is the horizontal velocity, m is the mass density field, and
! T1, T2, and T3 are tracers.
!
! In CICE, these equations have the form
!
! da/dt = del*(u*a) (4)
! dv/dt = d(a*h)/dt = del*(u*a*h) (5)
! de/dt = d(a*h*q)/dt = del*(u*a*h*q) (6)
! d(aT)/dt = del*(u*a*t) (7)
!
! where a = fractional ice area, v = ice/snow volume, h = v/a = thickness,
! e = ice/snow internal energy (J/m^2), q = e/v = internal energy per
! unit volume (J/m^3), and T is a tracer. These equations express
! conservation of ice area, volume, internal energy, and area-weighted
! tracer, respectively.
!
! (Note: In CICE, a, v and e are prognostic quantities from which
! h and q are diagnosed. The remapping routine works with tracers,
! which means that h and q must be derived from a, v, and e before
! calling the remapping routine.)
!
! Earlier versions of CICE assumed fixed ice and snow density.
! Beginning with CICE 4.0, the ice and snow density can be variable.
! In this case, equations (5) and (6) are replaced by
!
! dv/dt = d(a*h)/dt = del*(u*a*h) (8)
! dm/dt = d(a*h*rho)/dt = del*(u*a*h*rho) (9)
! de/dt = d(a*h*rho*qm)/dt = del*(u*a*h*rho*qm) (10)
!
! where rho = density and qm = internal energy per unit mass (J/kg).
! Eq. (9) expresses mass conservation, which in the variable-density
! case is no longer equivalent to volume conservation (8).
!
! Tracers satisfying equations of the form (1) are called "type 1."
! In CICE the paradigmatic type 1 tracers are hi and hs.
!
! Tracers satisfying equations of the form (2) are called "type 2".
! The paradigmatic type 2 tracers are qi and qs (or rhoi and rhos
! in the variable-density case).
!
! Tracers satisfying equations of the form (3) are called "type 3."
! The paradigmatic type 3 tracers are qmi and qms in the variable-density
! case. There are no such tracers in the constant-density case.
!
! The fields a, T1, and T2 are reconstructed in each grid cell with
! 2nd-order accuracy. T3 is reconstructed with 1st-order accuracy
! (i.e., it is transported in upwind fashion) in order to avoid
! additional mathematical complexity.
!
! The mass-like field lives in the array "mm" (shorthand for mean
! mass) and the tracers fields in the array "tm" (mean tracers).
! In order to transport tracers correctly, the remapping routine
! needs to know the tracers types and relationships. This is done
! as follows:
!
! Each field in the "tm" array is assigned an index, 1:max_ntrace.
! (Note: max_ntrace is not the same as max_ntrcr, the number of tracers
! in the trcrn state variable array. For remapping purposes we
! have additional tracers hi, hs, qi and qs.)
! For standard CICE with ntrcr = 1, nilyr = 4, and nslyr = 1, the
! indexing is as follows:
! 1 = hi
! 2 = hs
! 3 = Ts
! 4-7 = qi
! 8 = qs
!
! The tracer types (1,2,3) are contained in the "tracer_type" array.
! For standard CICE:
!
! tracer_type = (1 1 1 2 2 2 2 2)
!
! Type 2 and type 3 tracers are said to depend on type 1 tracers.
! For instance, qi depends on hi, which is to say that
! there is a conservation equation of the form (2) or (6).
! Thus we define a "depend" array. For standard CICE:
!
! depend = (0 0 0 1 1 1 1 2)
!
! which implies that elements 1-3 (hi, hs, Ts) are type 1,
! elements 4-7 (qi) depend on element 1 (hi), and element 8 (qs)
! depends on element 2 (hs).
!
! We also define a logical array "has_dependents". In standard CICE:
!
! has_dependents = (T T F F F F F F),
!
! which means that only elements 1 and 2 (hi and hs) have dependent
! tracers.
!
! For the variable-density case, things are a bit more complicated.
! Suppose we have 4 variable-density ice layers and one variable-
! density snow layer. Then the indexing is as follows:
! 1 = hi
! 2 = hs
! 3 = Ts
! 4-7 = rhoi
! 8 = rhos
! 9-12 = qmi
! 13 = qms
!
! The key arrays are:
!
! tracer_type = (1 1 1 2 2 2 2 2 3 3 3 3 3)
!
! depend = (0 0 0 1 1 1 1 2 4 5 6 7 8)
!
! has_dependents = (T T F T T T T T F F F F F)
!
! which imply that hi and hs are type 1 with dependents rhoi and rhos,
! while rhoi and rhos are type 2 with dependents qmi and qms.
!
! Tracers added to the ntrcr array are handled automatically
! by the remapping with little extra coding. It is necessary
! only to provide the correct type and dependency information.
!
! When using this routine in other models, most of the tracer dependency
! apparatus may be irrelevant. In a layered ocean model, for example,
! the transported fields are the layer thickness h (the mass density
! field) and two or more tracers (T, S, and various trace species).
! Suppose there are just two tracers, T and S. Then the tracer arrays
! have the values:
!
! tracer_type = (1 1)
! depend = (0 0)
! has_dependents = (F F)
!
! which is to say that all tracer transport equations are of the form (1).
!
! The tracer dependency arrays are optional input arguments for the
! main remapping subroutine. If these arrays are not passed in, they
! take on the default values tracer_type(:) = 1, depend(:) = 0, and
! has_dependents(:) = F, which are appropriate for most purposes.
!
! Another optional argument is integral_order. If integral_order = 1,
! then the triangle integrals are exact for linear functions of x and y.
! If integral_order = 2, these integrals are exact for both linear and
! quadratic functions. If integral_order = 3, integrals are exact for
! cubic functions as well. If all tracers are of type 1, then the
! integrals of mass*tracer are quadratic, and integral_order = 2 is
! sufficient. In CICE, where there are type 2 tracers, we integrate
! functions of the form mass*tracer1*tracer2. Thus integral_order = 3
! is required for exactness, though integral_order = 2 may be good enough
! in practice.
!
! Finally, a few words about the edgearea fields:
!
! In earlier versions of this scheme, the divergence of the velocity
! field implied by the remapping was, in general, different from the
! value of del*u computed in the dynamics. For energetic consistency
! (in CICE as well as in layered ocean models such as HYPOP),
! these two values should agree. This can be ensured by setting
! l_fixed_area = T and specifying the area transported across each grid
! cell edge in the arrays edgearea_e and edgearea_n. The departure
! regions are then tweaked, following an idea by Mats Bentsen, such
! that they have the desired area. If l_fixed_area = F, these regions
! are not tweaked, and the edgearea arrays are output variables.
!
!=======================================================================
contains
!=======================================================================
!
!BOP
!
! !IROUTINE: init_remap - initialize grid quantities used for remapping
!
! !INTERFACE:
!
subroutine init_remap 1,5
!
! !DESCRIPTION:
!
! Grid quantities used by the remapping transport scheme
!
! !REVISION HISTORY:
!
! author William H. Lipscomb, LANL
!
! !USES:
!
use ice_boundary
use ice_domain
use ice_blocks
use ice_grid, only: dxt, dyt, &
xav, yav, xxav, xyav, yyav, &
xxxav, xxyav, xyyav, yyyav
use ice_exit
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
integer (kind=int_kind) :: &
i, j, iblk ! standard indices
! Compute grid cell average geometric quantities on the scaled
! rectangular grid with dx = 1, dy = 1.
!
! Note: On a rectangular grid, the integral of any odd function
! of x or y = 0.
!$OMP PARALLEL DO PRIVATE(iblk,i,j)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
xav(i,j,iblk) = c0
yav(i,j,iblk) = c0
!!! These formulas would be used on a rectangular grid
!!! with dimensions (dxt, dyt):
!!! xxav(i,j,iblk) = dxt(i,j,iblk)**2 / c12
!!! yyav(i,j,iblk) = dyt(i,j,iblk)**2 / c12
xxav(i,j,iblk) = c1/c12
yyav(i,j,iblk) = c1/c12
xyav(i,j,iblk) = c0
xxxav(i,j,iblk) = c0
xxyav(i,j,iblk) = c0
xyyav(i,j,iblk) = c0
yyyav(i,j,iblk) = c0
enddo
enddo
enddo
!$OMP END PARALLEL DO
end subroutine init_remap
!=======================================================================
!BOP
!
! !IROUTINE: horizontal_remap - incremental remapping transport scheme
!
! !INTERFACE:
!
subroutine horizontal_remap (dt, ntrace, & 1,38
uvel, vvel, &
mm, tm, &
l_fixed_area, &
edgearea_e, edgearea_n, &
tracer_type_in, depend_in, &
has_dependents_in, &
integral_order_in, &
l_dp_midpt_in)
!
! !DESCRIPTION:
! Solve the transport equations for one timestep using the incremental
! remapping scheme developed by John Dukowicz and John Baumgardner (DB)
! and modified for sea ice by William Lipscomb and Elizabeth Hunke.
!
! This scheme preserves monotonicity of ice area and tracers. That is,
! it does not produce new extrema. It is second-order accurate in space,
! except where gradients are limited to preserve monotonicity.
!
! This version of the remapping allows the user to specify the areal
! flux across each edge, based on an idea developed by Mats Bentsen.
!
! !REVISION HISTORY:
!
! author William H. Lipscomb, LANL
! 2006: Moved driver (subroutine transport_remap) into separate module.
! Geometry changes (logically rectangular coordinates, fixed
! area fluxes)
!
! !USES:
!
use ice_boundary
use ice_global_reductions
use ice_domain
use ice_blocks
use ice_grid, only: HTE, HTN, dxt, dyt, dxu, dyu, &
tarea, tarear, hm, &
xav, yav, xxav, xyav, yyav, &
xxxav, xxyav, xyyav, yyyav
use ice_exit
use ice_calendar, only: istep1
use ice_timers
!
! !INPUT/OUTPUT PARAMETERS:
!
real (kind=dbl_kind), intent(in) :: &
dt ! time step
integer (kind=int_kind), intent(in) :: &
ntrace ! number of tracers in use
real (kind=dbl_kind), intent(in), &
dimension(nx_block,ny_block,max_blocks) :: &
uvel ,&! x-component of velocity (m/s)
vvel ! y-component of velocity (m/s)
real (kind=dbl_kind), intent(inout), &
dimension (nx_block,ny_block,0:ncat,max_blocks) :: &
mm ! mean mass values in each grid cell
real (kind=dbl_kind), intent(inout), &
dimension (nx_block,ny_block,max_ntrace,ncat,max_blocks) :: &
tm ! mean tracer values in each grid cell
!-------------------------------------------------------------------
! If l_fixed_area is true, the area of each departure region is
! computed in advance (e.g., by taking the divergence of the
! velocity field and passed to locate_triangles. The departure
! regions are adjusted to obtain the desired area.
! If false, edgearea is computed in locate_triangles and passed out.
!-------------------------------------------------------------------
logical, intent(in) :: &
l_fixed_area ! if true, edgearea_e and edgearea_n are prescribed
! if false, edgearea is computed here and passed out
real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks), &
intent(inout) :: &
edgearea_e ,&! area of departure regions for east edges
edgearea_n ! area of departure regions for north edges
integer (kind=int_kind), dimension (ntrace), intent(in), &
optional :: &
tracer_type_in ,&! = 1, 2, or 3 (see comments above)
depend_in ! tracer dependencies (see above)
logical (kind=log_kind), dimension (ntrace), intent(in), &
optional :: &
has_dependents_in ! true if a tracer has dependent tracers
integer (kind=int_kind), intent(in), optional :: &
integral_order_in ! polynomial order for triangle integrals
logical (kind=log_kind), intent(in), optional :: &
l_dp_midpt_in ! if true, find departure points using
! corrected midpoint velocity
!
!EOP
!
! local variables
integer (kind=int_kind), dimension (ntrace) :: &
tracer_type ,&! = 1, 2, or 3 (see comments above)
depend ! tracer dependencies (see above)
logical (kind=log_kind), dimension (ntrace) :: &
has_dependents ! true if a tracer has dependent tracers
integer (kind=int_kind) :: &
integral_order ! polynomial order for triangle integrals
logical (kind=log_kind) :: &
l_dp_midpt ! if true, find departure points using
! corrected midpoint velocity
integer (kind=int_kind) :: &
i, j ,&! horizontal indices
iblk ,&! block indices
ilo,ihi,jlo,jhi,&! beginning and end of physical domain
n ! ice category index
integer (kind=int_kind), dimension(0:ncat,max_blocks) :: &
icellsnc ! number of cells with ice
integer (kind=int_kind), &
dimension(nx_block*ny_block,0:ncat,max_blocks) :: &
indxinc, indxjnc ! compressed i/j indices
type (block) :: &
this_block ! block information for current block
real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
dpx ,&! x coordinates of departure points at cell corners
dpy ! y coordinates of departure points at cell corners
real (kind=dbl_kind), dimension(nx_block,ny_block,0:ncat,max_blocks) :: &
mc ,&! mass at geometric center of cell
mx, my ,&! limited derivative of mass wrt x and y
mmask ! = 1. if mass is present, = 0. otherwise
real (kind=dbl_kind), &
dimension (nx_block,ny_block,max_ntrace,ncat,max_blocks) :: &
tc ,&! tracer values at geometric center of cell
tx, ty ,&! limited derivative of tracer wrt x and y
tmask ! = 1. if tracer is present, = 0. otherwise
real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat,max_blocks) :: &
mflxe, mflxn ! mass transports across E and N cell edges
real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrace,ncat,max_blocks) :: &
mtflxe, mtflxn ! mass*tracer transports across E and N cell edges
real (kind=dbl_kind), dimension (nx_block,ny_block,ngroups,max_blocks) :: &
triarea ! area of east-edge departure triangle
real (kind=dbl_kind), dimension (nx_block,ny_block,0:nvert,ngroups,max_blocks) :: &
xp, yp ! x and y coordinates of special triangle points
! (need 4 points for triangle integrals)
integer (kind=int_kind), &
dimension (nx_block,ny_block,ngroups,max_blocks) :: &
iflux ,&! i index of cell contributing transport
jflux ! j index of cell contributing transport
integer (kind=int_kind), dimension(ngroups,max_blocks) :: &
icellsng ! number of cells with ice
integer (kind=int_kind), &
dimension(nx_block*ny_block,ngroups,max_blocks) :: &
indxing, indxjng ! compressed i/j indices
logical (kind=log_kind), dimension(max_blocks) :: &
l_stop ! if true, abort the model
integer (kind=int_kind), dimension(max_blocks) :: &
istop, jstop ! indices of grid cell where model aborts
character (len=char_len), dimension(max_blocks) :: &
edge ! 'north' or 'east'
real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
worka, &
workb, &
workc, &
workd
real (kind=dbl_kind), dimension (nx_block,ny_block,2,max_blocks) :: &
dpwork
real (kind=dbl_kind), dimension(nx_block,ny_block,2,0:ncat,max_blocks) :: &
mwork
! call t_barrierf ('cice_dyn_remap1_BARRIER',MPI_COMM_ICE)
! call t_startf ('cice_dyn_remap1')
l_stop = .false.
istop = 0
jstop = 0
!-------------------------------------------------------------------
! Initialize various remapping arrays and options
! These are either passed in as optional arguments or set to the
! default values.
!-------------------------------------------------------------------
if (present(tracer_type_in)) then
tracer_type(:) = tracer_type_in(:)
else
tracer_type(:) = 1
endif
if (present(depend_in)) then
depend(:) = depend_in(:)
else
depend(:) = 0
endif
if (present(has_dependents_in)) then
has_dependents(:) = has_dependents_in(:)
else
has_dependents(:) = .false.
endif
if (present(integral_order_in)) then
integral_order = integral_order_in
else
integral_order = 2 ! quadratic integrals
endif
if (present(l_dp_midpt_in)) then
l_dp_midpt = l_dp_midpt_in
else
l_dp_midpt = .false.
endif
worka(:,:) = c1
workb(:,:) = c1
workc(:,:) = c1
workd(:,:) = c1
!---!-------------------------------------------------------------------
!---! Remap the ice area and associated tracers.
!---! Remap the open water area (without tracers).
!---!-------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblk,this_block,ilo,ihi,jlo,jhi,n)
do iblk = 1, nblocks
this_block = get_block(blocks_ice(iblk),iblk)
ilo = this_block%ilo
ihi = this_block%ihi
jlo = this_block%jlo
jhi = this_block%jhi
!-------------------------------------------------------------------
! Compute masks and count ice cells.
! Masks are used to prevent tracer values in cells without ice from
! being used to compute tracer gradients.
!-------------------------------------------------------------------
call make_masks (nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
nghost, ntrace, &
has_dependents, icellsnc (:,iblk), &
indxinc(:,:,iblk), indxjnc (:,:,iblk), &
mm (:,:,:,iblk), mmask (:,:,:,iblk), &
tm (:,:,:,:,iblk), tmask(:,:,:,:,iblk))
!-------------------------------------------------------------------
! Construct linear fields, limiting gradients to preserve monotonicity.
! Note: Pass in unit arrays instead of true distances HTE, HTN, etc.
! The resulting gradients are in scaled coordinates.
!-------------------------------------------------------------------
! open water
call construct_fields(nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
nghost, ntrace, &
tracer_type, depend, &
has_dependents, icellsnc (0,iblk), &
indxinc(:,0,iblk), indxjnc(:,0,iblk), &
! HTN (:,:,iblk), HTE (:,:,iblk), &
worka (:,:), workb (:,:), &
hm (:,:,iblk), xav (:,:,iblk), &
yav (:,:,iblk), xxav (:,:,iblk), &
xyav (:,:,iblk), yyav (:,:,iblk), &
xxxav (:,:,iblk), xxyav (:,:,iblk), &
xyyav (:,:,iblk), yyyav (:,:,iblk), &
! dxt (:,:,iblk), dyt (:,:,iblk), &
workc (:,:), workd (:,:), &
mm (:,:,0,iblk), mc (:,:,0,iblk), &
mx (:,:,0,iblk), my (:,:,0,iblk), &
mmask(:,:,0,iblk) )
! ice categories
do n = 1, ncat
call construct_fields(nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
nghost, ntrace, &
tracer_type, depend, &
has_dependents, icellsnc (n,iblk), &
indxinc (:,n,iblk), indxjnc(:,n,iblk), &
! HTN (:,:,iblk), HTE (:,:,iblk), &
worka (:,:), workb (:,:), &
hm (:,:,iblk), xav (:,:,iblk), &
yav (:,:,iblk), xxav (:,:,iblk), &
xyav (:,:,iblk), yyav (:,:,iblk), &
xxxav (:,:,iblk), xxyav (:,:,iblk), &
xyyav (:,:,iblk), yyyav (:,:,iblk), &
! dxt (:,:,iblk), dyt (:,:,iblk), &
workc (:,:), workd (:,:), &
mm (:,:,n,iblk), mc (:,:,n,iblk), &
mx (:,:,n,iblk), my (:,:,n,iblk), &
mmask (:,:,n,iblk), &
tm (:,:,:,n,iblk), tc (:,:,:,n,iblk), &
tx (:,:,:,n,iblk), ty (:,:,:,n,iblk), &
tmask(:,:,:,n,iblk) )
enddo ! n
!-------------------------------------------------------------------
! Given velocity field at cell corners, compute departure points
! of trajectories.
!-------------------------------------------------------------------
call departure_points(nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
nghost, dt, &
uvel(:,:,iblk), vvel(:,:,iblk), &
dxu (:,:,iblk), dyu (:,:,iblk), &
HTN (:,:,iblk), HTE (:,:,iblk), &
dpx (:,:,iblk), dpy (:,:,iblk), &
l_dp_midpt, l_stop (iblk), &
istop (iblk), jstop (iblk))
if (l_stop(iblk)) then
this_block = get_block(blocks_ice(iblk),iblk)
write(nu_diag,*) 'istep1, my_task, iblk =', &
istep1, my_task, iblk
write (nu_diag,*) 'Global block:', this_block%block_id
if (istop(iblk) > 0 .and. jstop(iblk) > 0) &
write(nu_diag,*) 'Global i and j:', &
this_block%i_glob(istop(iblk)), &
this_block%j_glob(jstop(iblk))
call abort_ice('remap transport: bad departure points')
endif
enddo ! iblk
!$OMP END PARALLEL DO
! call t_stopf ('cice_dyn_remap1')
!-------------------------------------------------------------------
! Ghost cell updates
! If nghost >= 2, these calls are not needed
!-------------------------------------------------------------------
if (nghost==1) then
! call t_barrierf ('cice_dyn_remap_update1_BARRIER',MPI_COMM_ICE)
! call t_startf ('cice_dyn_remap_update1')
call ice_timer_start(timer_bound)
! departure points
! load em up
!$OMP PARALLEL DO PRIVATE(iblk,i,j)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
dpwork(i,j,1, iblk) = dpx(i,j,iblk)
dpwork(i,j,2, iblk) = dpy(i,j,iblk)
mwork (i,j,1,:,iblk) = mx (i,j,:,iblk)
mwork (i,j,2,:,iblk) = my (i,j,:,iblk)
enddo
enddo
enddo
!$OMP END PARALLEL DO
!jw call ice_HaloUpdate (dpx, halo_info, &
!jw field_loc_NEcorner, field_type_vector)
!jw call ice_HaloUpdate (dpy, halo_info, &
!jw field_loc_NEcorner, field_type_vector)
call ice_HaloUpdate (dpwork, halo_info, &
field_loc_NEcorner, field_type_vector)
! mass field
call ice_HaloUpdate (mc, halo_info, &
field_loc_center, field_type_scalar)
!jw call ice_HaloUpdate (mx, halo_info, &
!jw field_loc_center, field_type_vector)
!jw call ice_HaloUpdate (my, halo_info, &
!jw field_loc_center, field_type_vector)
call ice_HaloUpdate (mwork, halo_info, &
field_loc_center, field_type_vector)
!$OMP PARALLEL DO PRIVATE(iblk,i,j)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
dpx(i,j, iblk) = dpwork(i,j,1, iblk)
dpy(i,j, iblk) = dpwork(i,j,2, iblk)
mx (i,j,:,iblk) = mwork (i,j,1,:,iblk)
my (i,j,:,iblk) = mwork (i,j,2,:,iblk)
enddo
enddo
enddo
!$OMP END PARALLEL DO
! call t_stopf ('cice_dyn_remap_update1')
! call t_barrierf ('cice_dyn_remap_update2_BARRIER',MPI_COMM_ICE)
! call t_startf ('cice_dyn_remap_update2')
! tracer fields
call ice_HaloUpdate (tc(:,:,1:ntrace,:,:), halo_info, &
field_loc_center, field_type_scalar)
call ice_HaloUpdate (tx(:,:,1:ntrace,:,:), halo_info, &
field_loc_center, field_type_vector)
call ice_HaloUpdate (ty(:,:,1:ntrace,:,:), halo_info, &
field_loc_center, field_type_vector)
call ice_timer_stop(timer_bound)
! call t_stopf ('cice_dyn_remap_update2')
endif ! nghost
! call t_barrierf ('cice_dyn_remap2_BARRIER',MPI_COMM_ICE)
! call t_startf ('cice_dyn_remap2')
!$OMP PARALLEL DO PRIVATE(iblk,this_block,ilo,ihi,jlo,jhi,n)
do iblk = 1, nblocks
this_block = get_block(blocks_ice(iblk),iblk)
ilo = this_block%ilo
ihi = this_block%ihi
jlo = this_block%jlo
jhi = this_block%jhi
!-------------------------------------------------------------------
! Transports for east cell edges.
!-------------------------------------------------------------------
!-------------------------------------------------------------------
! Compute areas and vertices of departure triangles.
!-------------------------------------------------------------------
edge(iblk) = 'east'
call locate_triangles(nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
nghost, edge(iblk), &
icellsng (:,iblk), &
indxing (:,:,iblk), indxjng(:,:,iblk), &
dpx (:,:,iblk), dpy (:,:,iblk), &
dxu (:,:,iblk), dyu (:,:,iblk), &
xp (:,:,:,:,iblk), yp (:,:,:,:,iblk), &
iflux (:,:,:,iblk), jflux(:,:,:,iblk), &
triarea (:,:,:,iblk), l_fixed_area, &
edgearea_e(:,:,iblk))
!-------------------------------------------------------------------
! Given triangle vertices, compute coordinates of triangle points
! needed for transport integrals.
!-------------------------------------------------------------------
call triangle_coordinates (nx_block, ny_block, &
integral_order, icellsng (:,iblk), &
indxing(:,:,iblk), indxjng(:,:,iblk), &
xp (:,:,:,:,iblk), yp (:,:,:,:,iblk))
! Compute the transport across east cell edges by summing contributions
! from each triangle.
!-------------------------------------------------------------------
! open water
call transport_integrals(nx_block, ny_block, &
ntrace, icellsng (:,iblk), &
indxing(:,:,iblk), indxjng (:,:,iblk), &
tracer_type, depend, &
integral_order, triarea(:,:,:,iblk), &
iflux(:,:,:,iblk), jflux (:,:,:,iblk), &
xp (:,:,:,:,iblk), yp (:,:,:,:,iblk), &
mc (:,:,0,iblk), mx (:,:,0,iblk), &
my (:,:,0,iblk), mflxe (:,:,0,iblk))
! ice categories
do n = 1, ncat
call transport_integrals &
(nx_block, ny_block, &
ntrace, icellsng (:,iblk), &
indxing(:,:,iblk), indxjng (:,:,iblk), &
tracer_type, depend, &
integral_order, triarea (:,:,:,iblk), &
iflux(:,:,:,iblk), jflux (:,:,:,iblk), &
xp (:,:,:,:,iblk), yp (:,:,:,:,iblk), &
mc (:,:,n,iblk), mx (:,:,n,iblk), &
my (:,:,n,iblk), mflxe (:,:,n,iblk), &
tc (:,:,:,n,iblk), tx (:,:,:,n,iblk), &
ty (:,:,:,n,iblk), mtflxe(:,:,:,n,iblk))
enddo
!-------------------------------------------------------------------
! Repeat for north edges
!-------------------------------------------------------------------
edge(iblk) = 'north'
call locate_triangles(nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
nghost, edge(iblk), &
icellsng (:,iblk), &
indxing (:,:,iblk), indxjng(:,:,iblk), &
dpx (:,:,iblk), dpy (:,:,iblk), &
dxu (:,:,iblk), dyu (:,:,iblk), &
xp (:,:,:,:,iblk), yp (:,:,:,:,iblk), &
iflux (:,:,:,iblk), jflux(:,:,:,iblk), &
triarea (:,:,:,iblk), l_fixed_area, &
edgearea_n(:,:,iblk))
call triangle_coordinates (nx_block, ny_block, &
integral_order, icellsng (:,iblk), &
indxing(:,:,iblk), indxjng(:,:,iblk), &
xp (:,:,:,:,iblk), yp (:,:,:,:,iblk))
! open water
call transport_integrals(nx_block, ny_block, &
ntrace, icellsng (:,iblk), &
indxing(:,:,iblk), indxjng(:,:,iblk), &
tracer_type, depend, &
integral_order, triarea(:,:,:,iblk), &
iflux(:,:,:,iblk), jflux (:,:,:,iblk), &
xp (:,:,:,:,iblk), yp (:,:,:,:,iblk), &
mc (:,:,0,iblk), mx (:,:,0,iblk), &
my (:,:,0,iblk), mflxn (:,:,0,iblk))
! ice categories
do n = 1, ncat
call transport_integrals &
(nx_block, ny_block, &
ntrace, icellsng (:,iblk), &
indxing(:,:,iblk), indxjng (:,:,iblk), &
tracer_type, depend, &
integral_order, triarea (:,:,:,iblk), &
iflux(:,:,:,iblk), jflux (:,:,:,iblk), &
xp (:,:,:,:,iblk), yp (:,:,:,:,iblk), &
mc (:,:,n,iblk), mx (:,:,n,iblk), &
my (:,:,n,iblk), mflxn (:,:,n,iblk), &
tc (:,:,:,n,iblk), tx (:,:,:,n,iblk), &
ty (:,:,:,n,iblk), mtflxn(:,:,:,n,iblk))
enddo ! n
!-------------------------------------------------------------------
! Update the ice area and tracers.
!-------------------------------------------------------------------
! open water
call update_fields (nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
ntrace, &
tracer_type, depend, &
tarear(:,:,iblk), l_stop(iblk), &
istop(iblk), jstop (iblk), &
mflxe(:,:,0,iblk), mflxn(:,:,0,iblk), &
mm (:,:,0,iblk))
if (l_stop(iblk)) then
this_block = get_block(blocks_ice(iblk),iblk)
write (nu_diag,*) 'istep1, my_task, iblk, cat =', &
istep1, my_task, iblk, '0'
write (nu_diag,*) 'Global block:', this_block%block_id
if (istop(iblk) > 0 .and. jstop(iblk) > 0) &
write(nu_diag,*) 'Global i and j:', &
this_block%i_glob(istop(iblk)), &
this_block%j_glob(jstop(iblk))
call abort_ice ('ice remap_transport: negative area (open water)')
endif
! ice categories
do n = 1, ncat
call update_fields(nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
ntrace, &
tracer_type, depend, &
tarear(:,:,iblk), l_stop(iblk), &
istop(iblk), jstop(iblk), &
mflxe (:,:, n,iblk), mflxn (:,:, n,iblk), &
mm (:,:, n,iblk), &
mtflxe(:,:,:,n,iblk), mtflxn(:,:,:,n,iblk), &
tm (:,:,:,n,iblk))
if (l_stop(iblk)) then
this_block = get_block(blocks_ice(iblk),iblk)
write (nu_diag,*) 'istep1, my_task, iblk, cat =', &
istep1, my_task, iblk, n
write (nu_diag,*) 'Global block:', this_block%block_id
if (istop(iblk) > 0 .and. jstop(iblk) > 0) &
write(nu_diag,*) 'Global i and j:', &
this_block%i_glob(istop(iblk)), &
this_block%j_glob(jstop(iblk))
call abort_ice ('ice remap_transport: negative area (ice)')
endif
enddo ! n
enddo ! iblk
!$OMP END PARALLEL DO
! call t_stopf ('cice_dyn_remap2')
end subroutine horizontal_remap
!=======================================================================
!
!BOP
!
! !IROUTINE: make_masks - make area and tracer masks
!
! !INTERFACE:
!
subroutine make_masks (nx_block, ny_block, & 2
ilo, ihi, jlo, jhi, &
nghost, ntrace, &
has_dependents, &
icells, &
indxi, indxj, &
mm, mmask, &
tm, tmask)
!
! !DESCRIPTION:
!
! Make area and tracer masks.
!
! If an area is masked out (mm < puny), then the values of tracers
! in that grid cell are assumed to have no physical meaning.
!
! Similarly, if a tracer with dependents is masked out
! (abs(tm) < puny), then the values of its dependent tracers in that
! grid cell are assumed to have no physical meaning.
! For example, the enthalpy value has no meaning if the thickness
! is zero.
!
! !REVISION HISTORY:
!
! author William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block ,&! block dimensions
ilo,ihi,jlo,jhi ,&! beginning and end of physical domain
nghost ,&! number of ghost cells
ntrace ! number of tracers in use
logical (kind=log_kind), dimension (ntrace), intent(in) :: &
has_dependents ! true if a tracer has dependent tracers
integer (kind=int_kind), dimension(0:ncat), intent(out) :: &
icells ! number of cells with ice
integer (kind=int_kind), dimension(nx_block*ny_block,0:ncat), &
intent(out) :: &
indxi ,&! compressed i/j indices
indxj
real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat), &
intent(in) :: &
mm ! mean ice area in each grid cell
real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat), &
intent(out) :: &
mmask ! = 1. if ice is present, else = 0.
real (kind=dbl_kind), dimension (nx_block, ny_block, max_ntrace, ncat), &
intent(in), optional :: &
tm ! mean tracer values in each grid cell
real (kind=dbl_kind), dimension (nx_block, ny_block, max_ntrace, ncat), &
intent(out), optional :: &
tmask ! = 1. if tracer is present, else = 0.
!
!EOP
!
integer (kind=int_kind) :: &
i, j, ij ,&! horizontal indices
n ,&! ice category index
nt ! tracer index
do n = 0, ncat
do ij = 1, nx_block*ny_block
indxi(ij,n) = 0
indxj(ij,n) = 0
enddo
enddo
!-------------------------------------------------------------------
! open water mask
!-------------------------------------------------------------------
icells(0) = 0
do j = 1, ny_block
do i = 1, nx_block
if (mm(i,j,0) > puny) then
mmask(i,j,0) = c1
icells(0) = icells(0) + 1
ij = icells(0)
indxi(ij,0) = i
indxj(ij,0) = j
else
mmask(i,j,0) = c0
endif
enddo
enddo
do n = 1, ncat
!-------------------------------------------------------------------
! Find grid cells where ice is present.
!-------------------------------------------------------------------
icells(n) = 0
do j = 1, ny_block
do i = 1, nx_block
if (mm(i,j,n) > puny) then
icells(n) = icells(n) + 1
ij = icells(n)
indxi(ij,n) = i
indxj(ij,n) = j
endif ! mm > puny
enddo
enddo
!-------------------------------------------------------------------
! ice area mask
!-------------------------------------------------------------------
mmask(:,:,n) = c0
do ij = 1, icells(n)
i = indxi(ij,n)
j = indxj(ij,n)
mmask(i,j,n) = c1
enddo
!-------------------------------------------------------------------
! tracer masks
!-------------------------------------------------------------------
if (present(tm)) then
tmask(:,:,:,n) = c0
do nt = 1, ntrace
if (has_dependents(nt)) then
do ij = 1, icells(n)
i = indxi(ij,n)
j = indxj(ij,n)
if (abs(tm(i,j,nt,n)) > puny) then
tmask(i,j,nt,n) = c1
endif
enddo
endif
enddo
endif ! present(tm)
!-------------------------------------------------------------------
! Redefine icells
! For nghost = 1, exclude ghost cells
! For nghost = 2, include one layer of ghost cells
!-------------------------------------------------------------------
icells(n) = 0
do j = jlo-nghost+1, jhi+nghost-1
do i = ilo-nghost+1, ihi+nghost-1
if (mm(i,j,n) > puny) then
icells(n) = icells(n) + 1
ij = icells(n)
indxi(ij,n) = i
indxj(ij,n) = j
endif ! mm > puny
enddo
enddo
enddo ! n
end subroutine make_masks
!=======================================================================
!
!BOP
!
! !IROUTINE: construct_fields - construct fields of ice area and tracers
!
! !INTERFACE:
!
subroutine construct_fields (nx_block, ny_block, & 2,3
ilo, ihi, jlo, jhi, &
nghost, ntrace, &
tracer_type, depend, &
has_dependents, icells, &
indxi, indxj, &
HTN, HTE, &
hm, xav, &
yav, xxav, &
xyav, yyav, &
xxxav, xxyav, &
xyyav, yyyav, &
dxt, dyt, &
mm, mc, &
mx, my, &
mmask, &
tm, tc, &
tx, ty, &
tmask)
!
! !DESCRIPTION:
!
! Construct fields of ice area and tracers.
!
! !REVISION HISTORY:
!
! authors William H. Lipscomb, LANL
! John R. Baumgardner, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block ,&! block dimensions
ilo,ihi,jlo,jhi ,&! beginning and end of physical domain
nghost ,&! number of ghost cells
ntrace ,&! number of tracers in use
icells ! number of cells with mass
integer (kind=int_kind), dimension (ntrace), intent(in) :: &
tracer_type ,&! = 1, 2, or 3 (see comments above)
depend ! tracer dependencies (see above)
logical (kind=log_kind), dimension (ntrace), intent(in) :: &
has_dependents ! true if a tracer has dependent tracers
integer (kind=int_kind), dimension(nx_block*ny_block), intent(in) :: &
indxi ,&! compressed i/j indices
indxj
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(in) :: &
hm ,&! land/boundary mask, thickness (T-cell)
HTN ,&! length of northern edge of T-cell (m)
HTE ,&! length of eastern edge of T-cell (m)
xav, yav ,&! mean T-cell values of x, y
xxav, xyav, yyav ,&! mean T-cell values of xx, xy, yy
xxxav,xxyav,xyyav,yyyav,&! mean T-cell values of , xxy, xyy, yyy
dxt ,&! grid cell width (m)
dyt ! grid cell height (m)
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(in) :: &
mm ,&! mean value of mass field
mmask ! = 1. if ice is present, = 0. otherwise
real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrace), &
intent(in), optional :: &
tm ,&! mean tracer
tmask ! = 1. if tracer is present, = 0. otherwise
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(out) :: &
mc ,&! mass value at geometric center of cell
mx, my ! limited derivative of mass wrt x and y
real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrace), &
intent(out), optional :: &
tc ,&! tracer at geometric center of cell
tx, ty ! limited derivative of tracer wrt x and y
!
!EOP
!
integer (kind=int_kind) :: &
i, j ,&! horizontal indices
nt, nt1 ,&! tracer indices
ij ! combined i/j horizontal index
real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
mxav ,&! x coordinate of center of mass
myav ! y coordinate of center of mass
real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrace) :: &
mtxav ,&! x coordinate of center of mass*tracer
mtyav ! y coordinate of center of mass*tracer
real (kind=dbl_kind) :: &
w1, w2, w3, w4, w5, w6, w7 ! work variables
!-------------------------------------------------------------------
! Compute field values at the geometric center of each grid cell,
! and compute limited gradients in the x and y directions.
!
! For second order accuracy, each state variable is approximated as
! a field varying linearly over x and y within each cell. For each
! category, the integrated value of m(x,y) over the cell must
! equal mm(i,j,n)*tarea(i,j), where tarea(i,j) is the cell area.
! Similarly, the integrated value of m(x,y)*t(x,y) must equal
! the total mass*tracer, mm(i,j,n)*tm(i,j,n)*tarea(i,j).
!
! These integral conditions are satisfied for linear fields if we
! stipulate the following:
! (1) The mean mass, mm, is equal to the mass at the cell centroid.
! (2) The mean value tm1 of type 1 tracers is equal to the value
! at the center of mass.
! (3) The mean value tm2 of type 2 tracers is equal to the value
! at the center of mass*tm1, where tm2 depends on tm1.
! (See comments at the top of the module.)
!
! We want to find the value of each state variable at a standard
! reference point, which we choose to be the geometric center of
! the cell. The geometric center is located at the intersection
! of the line joining the midpoints of the north and south edges
! with the line joining the midpoints of the east and west edges.
! To find the value at the geometric center, we must know the
! location of the cell centroid/center of mass, along with the
! mean value and the gradients with respect to x and y.
!
! The cell gradients are first computed from the difference between
! values in the neighboring cells, then limited by requiring that
! no new extrema are created within the cell.
!
! For rectangular coordinates the centroid and the geometric
! center coincide, which means that some of the equations in this
! subroutine could be simplified. However, the full equations
! are retained for generality.
!-------------------------------------------------------------------
!-------------------------------------------------------------------
! Initialize
!-------------------------------------------------------------------
do j = 1, ny_block
do i = 1, nx_block
mc(i,j) = c0
mx(i,j) = c0
my(i,j) = c0
mxav(i,j) = c0
myav(i,j) = c0
enddo
enddo
if (present(tm)) then
do nt = 1, ntrace
do j = 1, ny_block
do i = 1, nx_block
tc(i,j,nt) = c0
tx(i,j,nt) = c0
ty(i,j,nt) = c0
enddo
enddo
enddo
endif
! limited gradient of mass field in each cell (except masked cells)
! Note: The gradient is computed in scaled coordinates with
! dxt = dyt = hte = htn = 1.
call limited_gradient (nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
nghost, &
mm, hm, &
xav, yav, &
HTN, HTE, &
dxt, dyt, &
mx, my)
do ij = 1,icells ! ice is present
i = indxi(ij)
j = indxj(ij)
! mass field at geometric center
mc(i,j) = mm(i,j) - xav(i,j)*mx(i,j) &
- yav(i,j)*my(i,j)
enddo ! ij
! tracers
if (present(tm)) then
do ij = 1,icells ! cells with mass
i = indxi(ij)
j = indxj(ij)
! center of mass (mxav,myav) for each cell
mxav(i,j) = (mx(i,j)*xxav(i,j) &
+ my(i,j)*xyav(i,j) &
+ mc(i,j)*xav (i,j)) / mm(i,j)
myav(i,j) = (mx(i,j)*xyav(i,j) &
+ my(i,j)*yyav(i,j) &
+ mc(i,j)*yav(i,j)) / mm(i,j)
enddo
do nt = 1, ntrace
if (tracer_type(nt)==1) then ! independent of other tracers
call limited_gradient(nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
nghost, &
tm(:,:,nt), mmask, &
mxav, myav, &
HTN, HTE, &
dxt, dyt, &
tx(:,:,nt), ty(:,:,nt))
if (has_dependents(nt)) then ! need center of area*tracer
do j = 1, ny_block
do i = 1, nx_block
mtxav(i,j,nt) = c0
mtyav(i,j,nt) = c0
enddo
enddo
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells ! Note: no tx or ty in ghost cells
! (bound calls are later)
i = indxi(ij)
j = indxj(ij)
! tracer value at geometric center
tc(i,j,nt) = tm(i,j,nt) - tx(i,j,nt)*mxav(i,j) &
- ty(i,j,nt)*myav(i,j)
if (tmask(i,j,nt) > puny) then
! center of area*tracer
w1 = mc(i,j)*tc(i,j,nt)
w2 = mc(i,j)*tx(i,j,nt) &
+ mx(i,j)*tc(i,j,nt)
w3 = mc(i,j)*ty(i,j,nt) &
+ my(i,j)*tc(i,j,nt)
w4 = mx(i,j)*tx(i,j,nt)
w5 = mx(i,j)*ty(i,j,nt) &
+ my(i,j)*tx(i,j,nt)
w6 = my(i,j)*ty(i,j,nt)
w7 = c1 / (mm(i,j)*tm(i,j,nt))
mtxav(i,j,nt) = (w1*xav (i,j) + w2*xxav (i,j) &
+ w3*xyav (i,j) + w4*xxxav(i,j) &
+ w5*xxyav(i,j) + w6*xyyav(i,j)) &
* w7
mtyav(i,j,nt) = (w1*yav(i,j) + w2*xyav (i,j) &
+ w3*yyav(i,j) + w4*xxyav(i,j) &
+ w5*xyyav(i,j) + w6*yyyav(i,j)) &
* w7
endif ! tmask
enddo ! ij
else ! no dependents
do ij = 1, icells ! mass is present
i = indxi(ij)
j = indxj(ij)
! tracer value at geometric center
tc(i,j,nt) = tm(i,j,nt) - tx(i,j,nt)*mxav(i,j) &
- ty(i,j,nt)*myav(i,j)
enddo ! ij
endif ! has_dependents
elseif (tracer_type(nt)==2) then ! tracer nt depends on nt1
nt1 = depend(nt)
call limited_gradient(nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
nghost, &
tm(:,:,nt), tmask(:,:,nt1), &
mtxav(:,:,nt1), mtyav(:,:,nt1), &
HTN, HTE, &
dxt, dyt, &
tx(:,:,nt), ty(:,:,nt))
do ij = 1, icells ! ice is present
i = indxi(ij)
j = indxj(ij)
tc(i,j,nt) = tm(i,j,nt) &
- tx(i,j,nt) * mtxav(i,j,nt1) &
- ty(i,j,nt) * mtyav(i,j,nt1)
enddo ! ij
elseif (tracer_type(nt)==3) then ! upwind approx; gradient = 0
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
tc(i,j,nt) = tm(i,j,nt)
! tx(i,j,nt) = c0 ! already initialized to 0.
! ty(i,j,nt) = c0
enddo ! ij
endif ! tracer_type
enddo ! ntrace
endif ! present (tm)
end subroutine construct_fields
!=======================================================================
!
!BOP
!
! !IROUTINE: limited_gradient - limited gradient of a scalar field
!
! !INTERFACE:
!
subroutine limited_gradient (nx_block, ny_block, & 3,25
ilo, ihi, jlo, jhi, &
nghost, &
phi, phimask, &
cnx, cny, &
HTN, HTE, &
dxt, dyt, &
gx, gy)
!
! !DESCRIPTION:
!
! Compute a limited gradient of the scalar field phi in scaled coordinates.
! "Limited" means that we do not create new extrema in phi. For
! instance, field values at the cell corners can neither exceed the
! maximum of phi(i,j) in the cell and its eight neighbors, nor fall
! below the minimum.
!
! !REVISION HISTORY:
!
! authors William H. Lipscomb, LANL
! John R. Baumgardner, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block,&! block dimensions
ilo,ihi,jlo,jhi ,&! beginning and end of physical domain
nghost ! number of ghost cells
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent (in) :: &
phi ,&! input tracer field (mean values in each grid cell)
cnx ,&! x-coordinate of phi relative to geometric center of cell
cny ,&! y-coordinate of phi relative to geometric center of cell
dxt ,&! grid cell width (m)
dyt ,&! grid cell height (m)
phimask ,&
! phimask(i,j) = 1 if phi(i,j) has physical meaning, = 0 otherwise.
! For instance, aice has no physical meaning in land cells,
! and hice no physical meaning where aice = 0.
HTN ,&! length of northern edge of T-cell (m)
HTE ! length of eastern edge of T-cell (m)
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(out) :: &
gx ,&! limited x-direction gradient
gy ! limited y-direction gradient
!
!EOP
!
integer (kind=int_kind) :: &
i, j, ij ,&! standard indices
icells ! number of cells to limit
integer (kind=int_kind), dimension(nx_block*ny_block) :: &
indxi, indxj ! combined i/j horizontal indices
real (kind=dbl_kind) :: &
phi_nw, phi_n, phi_ne ,&! values of phi in 8 neighbor cells
phi_w, phi_e ,&
phi_sw, phi_s, phi_se ,&
qmn, qmx ,&! min and max value of phi within grid cell
pmn, pmx ,&! min and max value of phi among neighbor cells
w1, w2, w3, w4 ! work variables
real (kind=dbl_kind) :: &
gxtmp, gytmp ! temporary term for x- and y- limited gradient
gx(:,:) = c0
gy(:,:) = c0
! For nghost = 1, loop over physical cells and update ghost cells later
! For nghost = 2, loop over a layer of ghost cells and skip the update
icells = 0
do j = jlo-nghost+1, jhi+nghost-1
do i = ilo-nghost+1, ihi+nghost-1
if (phimask(i,j) > puny) then
!jw icells = icells + 1
!jw indxi(icells) = i
!jw indxj(icells) = j
!jw endif ! phimask > puny
!jw enddo
!jw enddo
!jw do ij = 1, icells
!jw i = indxi(ij)
!jw j = indxj(ij)
! Store values of phi in the 8 neighbor cells.
! Note: phimask = 1. or 0. If phimask = 1., use the true value;
! if phimask = 0., use the home cell value so that non-physical
! values of phi do not contribute to the gradient.
phi_nw = phimask(i-1,j+1) * phi(i-1,j+1) &
+ (c1-phimask(i-1,j+1))* phi(i,j)
phi_n = phimask(i,j+1) * phi(i,j+1) &
+ (c1-phimask(i,j+1)) * phi(i,j)
phi_ne = phimask(i+1,j+1) * phi(i+1,j+1) &
+ (c1-phimask(i+1,j+1))* phi(i,j)
phi_w = phimask(i-1,j) * phi(i-1,j) &
+ (c1-phimask(i-1,j)) * phi(i,j)
phi_e = phimask(i+1,j) * phi(i+1,j) &
+ (c1-phimask(i+1,j)) * phi(i,j)
phi_sw = phimask(i-1,j-1) * phi(i-1,j-1) &
+ (c1-phimask(i-1,j-1))* phi(i,j)
phi_s = phimask(i,j-1) * phi(i,j-1) &
+ (c1-phimask(i,j-1)) * phi(i,j)
phi_se = phimask(i+1,j-1) * phi(i+1,j-1) &
+ (c1-phimask(i+1,j-1))* phi(i,j)
! unlimited gradient components
! (factors of two cancel out)
gxtmp = (phi_e - phi(i,j)) / (dxt(i,j) + dxt(i+1,j)) &
+ (phi(i,j) - phi_w) / (dxt(i-1,j) + dxt(i,j) )
gytmp = (phi_n - phi(i,j)) / (dyt(i,j) + dyt(i,j+1)) &
+ (phi(i,j) - phi_s) / (dyt(i,j-1) + dyt(i,j) )
! minimum and maximum among the nine local cells
pmn = min (phi_nw, phi_n, phi_ne, phi_w, phi(i,j), &
phi_e, phi_sw, phi_s, phi_se)
pmx = max (phi_nw, phi_n, phi_ne, phi_w, phi(i,j), &
phi_e, phi_sw, phi_s, phi_se)
pmn = pmn - phi(i,j)
pmx = pmx - phi(i,j)
! minimum and maximum deviation of phi within the cell
w1 = (p5*HTN(i,j) - cnx(i,j)) * gxtmp &
+ (p5*HTE(i,j) - cny(i,j)) * gytmp
w2 = (p5*HTN(i,j-1) - cnx(i,j)) * gxtmp &
- (p5*HTE(i,j) + cny(i,j)) * gytmp
w3 = -(p5*HTN(i,j-1) + cnx(i,j)) * gxtmp &
- (p5*HTE(i-1,j) + cny(i,j)) * gytmp
w4 = (p5*HTE(i-1,j) - cny(i,j)) * gytmp &
- (p5*HTN(i,j) + cnx(i,j)) * gxtmp
qmn = min (w1, w2, w3, w4)
qmx = max (w1, w2, w3, w4)
! Watch for underflows here
! the limiting coefficient
if (abs(qmn) > 1.0e-300_dbl_kind) then ! 'abs(qmn) > puny' not sufficient
w1 = max(c0, pmn/qmn)
else
w1 = c1
endif
if (abs(qmx) > 1.0e-300_dbl_kind) then
w2 = max(c0, pmx/qmx)
else
w2 = c1
endif
w1 = min(c1, w1, w2)
! Limit the gradient components
gx(i,j) = w1 * gxtmp
gy(i,j) = w1 * gytmp
!jw enddo ! ij
endif
enddo
enddo
end subroutine limited_gradient
!=======================================================================
!BOP
!
! !IROUTINE: departure_points - compute departure points of trajectories
!
! !INTERFACE:
!
subroutine departure_points (nx_block, ny_block, & 1
ilo, ihi, jlo, jhi, &
nghost, dt, &
uvel, vvel, &
dxu, dyu, &
HTN, HTE, &
dpx, dpy, &
l_dp_midpt, l_stop, &
istop, jstop)
!
! !DESCRIPTION:
!
! Given velocity fields on cell corners, compute departure points
! of back trajectories in nondimensional coordinates.
!
! !REVISION HISTORY:
!
! author William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block,&! block dimensions
ilo,ihi,jlo,jhi, &! beginning and end of physical domain
nghost ! number of ghost cells
real (kind=dbl_kind), intent(in) :: &
dt ! time step (s)
real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
uvel ,&! x-component of velocity (m/s)
vvel ,&! y-component of velocity (m/s)
dxu ,&! E-W dimensions of U-cell (m)
dyu ,&! N-S dimensions of U-cell (m)
HTN ,&! length of north face of T-cell (m)
HTE ! length of east face of T-cell (m)
real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: &
dpx ,&! coordinates of departure points (m)
dpy ! coordinates of departure points (m)
logical (kind=log_kind), intent(in) :: &
l_dp_midpt ! if true, find departure points using
! corrected midpoint velocity
logical (kind=log_kind), intent(inout) :: &
l_stop ! if true, abort on return
integer (kind=int_kind), intent(inout) :: &
istop, jstop ! indices of grid cell where model aborts
!
!EOP
!
integer (kind=int_kind) :: &
i, j, i2, j2 ! horizontal indices
real (kind=dbl_kind) :: &
mpx, mpy ,&! coordinates of midpoint of back trajectory,
! relative to cell corner
mpxt, mpyt ,&! midpoint coordinates relative to cell center
ump, vmp ! corrected velocity at midpoint
!-------------------------------------------------------------------
! Estimate departure points.
! This estimate is 1st-order accurate in time; improve accuracy by
! using midpoint approximation (to add later).
! For nghost = 1, loop over physical cells and update ghost cells later.
! For nghost = 2, loop over a layer of ghost cells and skip update.
!-------------------------------------------------------------------
dpx(:,:) = c0
dpy(:,:) = c0
do j = jlo-nghost+1, jhi+nghost-1
do i = ilo-nghost+1, ihi+nghost-1
dpx(i,j) = -dt*uvel(i,j)
dpy(i,j) = -dt*vvel(i,j)
! Check for values out of bounds (more than one grid cell away)
if (dpx(i,j) < -HTN(i,j) .or. dpx(i,j) > HTN(i+1,j) .or. &
dpy(i,j) < -HTE(i,j) .or. dpy(i,j) > HTE(i,j+1)) then
l_stop = .true.
istop = i
jstop = j
endif
enddo
enddo
if (l_stop) then
i = istop
j = jstop
write (nu_diag,*) ' '
write (nu_diag,*) &
'Warning: Departure points out of bounds in remap'
write (nu_diag,*) 'my_task, i, j =', my_task, i, j
write (nu_diag,*) 'dpx, dpy =', dpx(i,j), dpy(i,j)
write (nu_diag,*) 'HTN(i,j), HTN(i+1,j) =', HTN(i,j), HTN(i+1,j)
write (nu_diag,*) 'HTE(i,j), HTE(i,j+1) =', HTE(i,j), HTE(i,j+1)
return
endif
if (l_dp_midpt) then ! find dep pts using corrected midpt velocity
do j = jlo-nghost+1, jhi+nghost-1
do i = ilo-nghost+1, ihi+nghost-1
if (uvel(i,j)/=c0 .or. vvel(i,j)/=c0) then
!-------------------------------------------------------------------
! Scale departure points to coordinate system in which grid cells
! have sides of unit length.
!-------------------------------------------------------------------
dpx(i,j) = dpx(i,j) / dxu(i,j)
dpy(i,j) = dpy(i,j) / dyu(i,j)
!-------------------------------------------------------------------
! Estimate midpoint of backward trajectory relative to corner (i,j).
!-------------------------------------------------------------------
mpx = p5 * dpx(i,j)
mpy = p5 * dpy(i,j)
!-------------------------------------------------------------------
! Determine the indices (i2,j2) of the cell where the trajectory lies.
! Compute the coordinates of the midpoint of the backward trajectory
! relative to the cell center in a stretch coordinate system
! with vertices at (1/2, 1/2), (1/2, -1/2), etc.
!-------------------------------------------------------------------
if (mpx >= c0 .and. mpy >= c0) then ! cell (i+1,j+1)
i2 = i+1
j2 = j+1
mpxt = mpx - p5
mpyt = mpy - p5
elseif (mpx < c0 .and. mpy < c0) then ! cell (i,j)
i2 = i
j2 = j
mpxt = mpx + p5
mpyt = mpy + p5
elseif (mpx >= c0 .and. mpy < c0) then ! cell (i+1,j)
i2 = i+1
j2 = j
mpxt = mpx - p5
mpyt = mpy + p5
elseif (mpx < c0 .and. mpy >= c0) then ! cell (i,j+1)
i2 = i
j2 = j+1
mpxt = mpx + p5
mpyt = mpy - p5
endif
!-------------------------------------------------------------------
! Using a bilinear approximation, estimate the velocity at the
! trajectory midpoint in the (i2,j2) reference frame.
!-------------------------------------------------------------------
ump = uvel(i2-1,j2-1)*(mpxt-p5)*(mpyt-p5) &
- uvel(i2, j2-1)*(mpxt+p5)*(mpyt-p5) &
+ uvel(i2, j2 )*(mpxt+p5)*(mpyt+p5) &
- uvel(i2-1,j2 )*(mpxt-p5)*(mpyt+p5)
vmp = vvel(i2-1,j2-1)*(mpxt-p5)*(mpyt-p5) &
- vvel(i2, j2-1)*(mpxt+p5)*(mpyt-p5) &
+ vvel(i2, j2 )*(mpxt+p5)*(mpyt+p5) &
- vvel(i2-1,j2 )*(mpxt-p5)*(mpyt+p5)
!-------------------------------------------------------------------
! Use the midpoint velocity to estimate the coordinates of the
! departure point relative to corner (i,j).
!-------------------------------------------------------------------
dpx(i,j) = -dt * ump
dpy(i,j) = -dt * vmp
endif ! nonzero velocity
enddo ! i
enddo ! j
endif ! l_dp_midpt
end subroutine departure_points
!=======================================================================
!
!BOP
!
! !IROUTINE: locate_triangles - triangle info for cell edges
!
! !INTERFACE:
!
subroutine locate_triangles (nx_block, ny_block, & 2
ilo, ihi, jlo, jhi, &
nghost, edge, &
icells, &
indxi, indxj, &
dpx, dpy, &
dxu, dyu, &
xp, yp, &
iflux, jflux, &
triarea, &
l_fixed_area, edgearea)
!
! !DESCRIPTION:
!
! Compute areas and vertices of transport triangles for north or
! east cell edges.
!
! !REVISION HISTORY:
!
! authors William H. Lipscomb, LANL
! John R. Baumgardner, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block,&! block dimensions
ilo,ihi,jlo,jhi ,&! beginning and end of physical domain
nghost ! number of ghost cells
character (len=char_len), intent(in) :: &
edge ! 'north' or 'east'
real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: &
dpx ,&! x coordinates of departure points at cell corners
dpy ,&! y coordinates of departure points at cell corners
dxu ,&! E-W dimension of U-cell (m)
dyu ! N-S dimension of U-cell (m)
real (kind=dbl_kind), dimension (nx_block,ny_block,0:nvert,ngroups), &
intent(out) :: &
xp, yp ! coordinates of triangle vertices
real (kind=dbl_kind), dimension (nx_block,ny_block,ngroups), &
intent(out) :: &
triarea ! area of departure triangle
integer (kind=int_kind), dimension (nx_block,ny_block,ngroups), &
intent(out) :: &
iflux ,&! i index of cell contributing transport
jflux ! j index of cell contributing transport
integer (kind=int_kind), dimension (ngroups), intent(out) :: &
icells ! number of cells where triarea > puny
integer (kind=int_kind), dimension (nx_block*ny_block,ngroups), &
intent(out) :: &
indxi ,&! compressed index in i-direction
indxj ! compressed index in j-direction
logical, intent(in) :: &
l_fixed_area ! if true, the area of each departure region is
! passed in as edgearea
! if false, edgearea if determined internally
! and is passed out
real (kind=dbl_kind), dimension(nx_block,ny_block), intent(inout) :: &
edgearea ! area of departure region for each edge
! edgearea > 0 for eastward/northward flow
!
!EOP
!
integer (kind=int_kind) :: &
i, j, ij, ic ,&! horizontal indices
ib, ie, jb, je ,&! limits for loops over edges
ng, nv ,&! triangle indices
ishift, jshift ! differences between neighbor cells
integer (kind=int_kind) :: &
icellsd ! number of cells where departure area > 0.
integer (kind=int_kind), dimension (nx_block*ny_block) :: &
indxid ,&! compressed index in i-direction
indxjd ! compressed index in j-direction
real (kind=dbl_kind), dimension(nx_block,ny_block) :: &
dx, dy ,&! scaled departure points
areafac_c ,&! area scale factor at center of edge
areafac_l ,&! area scale factor at left corner
areafac_r ! area scale factor at right corner
real (kind=dbl_kind) :: &
xcl, ycl ,&! coordinates of left corner point
! (relative to midpoint of edge)
xdl, ydl ,&! left departure point
xil, yil ,&! left intersection point
xcr, ycr ,&! right corner point
xdr, ydr ,&! right departure point
xir, yir ,&! right intersection point
xic, yic ,&! x-axis intersection point
xicl, yicl ,&! left-hand x-axis intersection point
xicr, yicr ,&! right-hand x-axis intersection point
xdm, ydm ,&! midpoint of segment connecting DL and DR;
! shifted if l_fixed_area = T
dxc ,&! xcr - xcl
dxd ,&! xdr - xdl
md ,&! slope of line connecting DL and DR
mdl ,&! slope of line connecting DL and DM
mdr ,&! slope of line connecting DR and DM
ishift_tl, jshift_tl ,&! i,j indices of TL cell relative to edge
ishift_bl, jshift_bl ,&! i,j indices of BL cell relative to edge
ishift_tr, jshift_tr ,&! i,j indices of TR cell relative to edge
ishift_br, jshift_br ,&! i,j indices of BR cell relative to edge
ishift_tc, jshift_tc ,&! i,j indices of TC cell relative to edge
ishift_bc, jshift_bc ,&! i,j indices of BC cell relative to edge
area1, area2 ,&! temporary triangle areas
area3, area4 ,&!
area_c ,&! center polygon area
w1, w2 ! work variables
real (kind=dbl_kind), dimension (nx_block,ny_block,ngroups) :: &
areafact ! = 1 for positive flux, -1 for negative
real (kind=dbl_kind), dimension(nx_block,ny_block) :: &
areasum ! sum of triangle areas for a given edge
!-------------------------------------------------------------------
! Triangle notation:
! For each edge, there are 20 triangles that can contribute,
! but many of these are mutually exclusive. It turns out that
! at most 5 triangles can contribute to transport integrals at once.
!
! See Figure 3 in DB for pictures of these triangles.
! See Table 1 in DB for logical conditions.
!
! For the north edge, DB refer to these triangles as:
! (1) NW, NW1, W, W2
! (2) NE, NE1, E, E2
! (3) NW2, W1, NE2, E1
! (4) H1a, H1b, N1a, N1b
! (5) H2a, H2b, N2a, N2b
!
! For the east edge, DB refer to these triangles as:
! (1) NE, NE1, N, N2
! (2) SE, SE1, S, S2
! (3) NE2, N1, SE2, S1
! (4) H1a, H1b, E1a, E2b
! (5) H2a, H2b, E2a, E2b
!
! The code below works for either north or east edges.
! The respective triangle labels are:
! (1) TL, TL1, BL, BL2
! (2) TR, TR1, BR, BR2
! (3) TL2, BL1, TR2, BR1
! (4) BC1a, BC1b, TC1a, TC2b
! (5) BC2a, BC2b, TC2a, TC2b
!
! where the cell labels are:
!
! | |
! TL | TC | TR (top left, center, right)
! | |
! ------------------------
! | |
! BL | BC | BR (bottom left, center, right)
! | |
!
! and the transport is across the edge between cells TC and TB.
!
! Departure points are scaled to a local coordinate system
! whose origin is at the midpoint of the edge.
! In this coordinate system, the lefthand corner CL = (-0.5,0)
! and the righthand corner CR = (0.5, 0).
!-------------------------------------------------------------------
!-------------------------------------------------------------------
! Initialize
!-------------------------------------------------------------------
areafac_c(:,:) = c0
areafac_l(:,:) = c0
areafac_r(:,:) = c0
do ng = 1, ngroups
do j = 1, ny_block
do i = 1, nx_block
triarea (i,j,ng) = c0
areafact(i,j,ng) = c0
iflux (i,j,ng) = i
jflux (i,j,ng) = j
enddo
enddo
do nv = 0, nvert
do j = 1, ny_block
do i = 1, nx_block
xp(i,j,nv,ng) = c0
yp(i,j,nv,ng) = c0
enddo
enddo
enddo
enddo
if (trim(edge) == 'north') then
! loop size
ib = ilo
ie = ihi
jb = jlo - nghost ! lowest j index is a ghost cell
je = jhi
! index shifts for neighbor cells
ishift_tl = -1
jshift_tl = 1
ishift_bl = -1
jshift_bl = 0
ishift_tr = 1
jshift_tr = 1
ishift_br = 1
jshift_br = 0
ishift_tc = 0
jshift_tc = 1
ishift_bc = 0
jshift_bc = 0
! area scale factor
do j = jb, je
do i = ib, ie
areafac_l(i,j) = dxu(i-1,j)*dyu(i-1,j)
areafac_r(i,j) = dxu(i,j)*dyu(i,j)
areafac_c(i,j) = p5*(areafac_l(i,j) + areafac_r(i,j))
enddo
enddo
else ! east edge
! loop size
ib = ilo - nghost ! lowest i index is a ghost cell
ie = ihi
jb = jlo
je = jhi
! index shifts for neighbor cells
ishift_tl = 1
jshift_tl = 1
ishift_bl = 0
jshift_bl = 1
ishift_tr = 1
jshift_tr = -1
ishift_br = 0
jshift_br = -1
ishift_tc = 1
jshift_tc = 0
ishift_bc = 0
jshift_bc = 0
! area scale factors
do j = jb, je
do i = ib, ie
areafac_l(i,j) = dxu(i,j)*dyu(i,j)
areafac_r(i,j) = dxu(i,j-1)*dyu(i,j-1)
areafac_c(i,j) = p5 * (areafac_l(i,j) + areafac_r(i,j))
enddo
enddo
endif
!-------------------------------------------------------------------
! Compute mask for edges with nonzero departure areas
!-------------------------------------------------------------------
if (l_fixed_area) then
icellsd = 0
do j = jb, je
do i = ib, ie
if (edgearea(i,j) /= c0) then
icellsd = icellsd + 1
indxid(icellsd) = i
indxjd(icellsd) = j
endif
enddo
enddo
else
icellsd = 0
if (trim(edge) == 'north') then
do j = jb, je
do i = ib, ie
if (dpx(i-1,j)/=c0 .or. dpy(i-1,j)/=c0 &
.or. &
dpx(i,j)/=c0 .or. dpy(i,j)/=c0) then
icellsd = icellsd + 1
indxid(icellsd) = i
indxjd(icellsd) = j
endif
enddo
enddo
else ! east edge
do j = jb, je
do i = ib, ie
if (dpx(i,j-1)/=c0 .or. dpy(i,j-1)/=c0 &
.or. &
dpx(i,j)/=c0 .or. dpy(i,j)/=c0) then
icellsd = icellsd + 1
indxid(icellsd) = i
indxjd(icellsd) = j
endif
enddo
enddo
endif ! edge = north/east
endif ! l_fixed_area
!-------------------------------------------------------------------
! Scale the departure points
!-------------------------------------------------------------------
do j = 1, je
do i = 1, ie
dx(i,j) = dpx(i,j) / dxu(i,j)
dy(i,j) = dpy(i,j) / dyu(i,j)
enddo
enddo
!-------------------------------------------------------------------
! Compute departure regions, divide into triangles, and locate
! vertices of each triangle.
! Work in a nondimensional coordinate system in which lengths are
! scaled by the local metric coefficients (dxu and dyu).
! Note: The do loop includes north faces of the j = 1 ghost cells
! when edge = 'north'. The loop includes east faces of i = 1
! ghost cells when edge = 'east'.
!-------------------------------------------------------------------
do ij = 1, icellsd
i = indxid(ij)
j = indxjd(ij)
xcl = -p5
ycl = c0
xcr = p5
ycr = c0
! Departure points
if (trim(edge) == 'north') then ! north edge
xdl = xcl + dx(i-1,j)
ydl = ycl + dy(i-1,j)
xdr = xcr + dx(i,j)
ydr = ycr + dy(i,j)
else ! east edge; rotate trajectory by pi/2
xdl = xcl - dy(i,j)
ydl = ycl + dx(i,j)
xdr = xcr - dy(i,j-1)
ydr = ycr + dx(i,j-1)
endif
xdm = p5 * (xdr + xdl)
ydm = p5 * (ydr + ydl)
! Intersection points
xil = xcl
yil = (xcl*(ydm-ydl) + xdm*ydl - xdl*ydm) / (xdm - xdl)
xir = xcr
yir = (xcr*(ydr-ydm) - xdm*ydr + xdr*ydm) / (xdr - xdm)
md = (ydr - ydl) / (xdr - xdl)
if (abs(md) > puny) then
xic = xdl - ydl/md
else
xic = c0
endif
yic = c0
xicl = xic
yicl = yic
xicr = xic
yicr = yic
!-------------------------------------------------------------------
! Locate triangles in TL cell (NW for north edge, NE for east edge)
! and BL cell (W for north edge, N for east edge).
!-------------------------------------------------------------------
if (yil > c0 .and. xdl < xcl .and. ydl >= c0) then
! TL (group 1)
ng = 1
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xil
yp (i,j,2,ng) = yil
xp (i,j,3,ng) = xdl
yp (i,j,3,ng) = ydl
iflux (i,j,ng) = i + ishift_tl
jflux (i,j,ng) = j + jshift_tl
areafact(i,j,ng) = -areafac_l(i,j)
elseif (yil < c0 .and. xdl < xcl .and. ydl < c0) then
! BL (group 1)
ng = 1
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xdl
yp (i,j,2,ng) = ydl
xp (i,j,3,ng) = xil
yp (i,j,3,ng) = yil
iflux (i,j,ng) = i + ishift_bl
jflux (i,j,ng) = j + jshift_bl
areafact(i,j,ng) = areafac_l(i,j)
elseif (yil < c0 .and. xdl < xcl .and. ydl >= c0) then
! TL1 (group 1)
ng = 1
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xdl
yp (i,j,2,ng) = ydl
xp (i,j,3,ng) = xic
yp (i,j,3,ng) = yic
iflux (i,j,ng) = i + ishift_tl
jflux (i,j,ng) = j + jshift_tl
areafact(i,j,ng) = areafac_l(i,j)
! BL1 (group 3)
ng = 3
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xic
yp (i,j,2,ng) = yic
xp (i,j,3,ng) = xil
yp (i,j,3,ng) = yil
iflux (i,j,ng) = i + ishift_bl
jflux (i,j,ng) = j + jshift_bl
areafact(i,j,ng) = areafac_l(i,j)
elseif (yil > c0 .and. xdl < xcl .and. ydl < c0) then
! TL2 (group 3)
ng = 3
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xil
yp (i,j,2,ng) = yil
xp (i,j,3,ng) = xic
yp (i,j,3,ng) = yic
iflux (i,j,ng) = i + ishift_tl
jflux (i,j,ng) = j + jshift_tl
areafact(i,j,ng) = -areafac_l(i,j)
! BL2 (group 1)
ng = 1
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xic
yp (i,j,2,ng) = yic
xp (i,j,3,ng) = xdl
yp (i,j,3,ng) = ydl
iflux (i,j,ng) = i + ishift_bl
jflux (i,j,ng) = j + jshift_bl
areafact(i,j,ng) = -areafac_l(i,j)
endif ! TL and BL triangles
!-------------------------------------------------------------------
! Locate triangles in TR cell (NE for north edge, SE for east edge)
! and in BR cell (E for north edge, S for east edge).
!-------------------------------------------------------------------
if (yir > c0 .and. xdr >= xcr .and. ydr >= c0) then
! TR (group 2)
ng = 2
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xdr
yp (i,j,2,ng) = ydr
xp (i,j,3,ng) = xir
yp (i,j,3,ng) = yir
iflux (i,j,ng) = i + ishift_tr
jflux (i,j,ng) = j + jshift_tr
areafact(i,j,ng) = -areafac_r(i,j)
elseif (yir < c0 .and. xdr >= xcr .and. ydr < c0) then
! BR (group 2)
ng = 2
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xir
yp (i,j,2,ng) = yir
xp (i,j,3,ng) = xdr
yp (i,j,3,ng) = ydr
iflux (i,j,ng) = i + ishift_br
jflux (i,j,ng) = j + jshift_br
areafact(i,j,ng) = areafac_r(i,j)
elseif (yir < c0 .and. xdr >= xcr .and. ydr >= c0) then
! TR1 (group 2)
ng = 2
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xic
yp (i,j,2,ng) = yic
xp (i,j,3,ng) = xdr
yp (i,j,3,ng) = ydr
iflux (i,j,ng) = i + ishift_tr
jflux (i,j,ng) = j + jshift_tr
areafact(i,j,ng) = areafac_r(i,j)
! BR1 (group 3)
ng = 3
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xir
yp (i,j,2,ng) = yir
xp (i,j,3,ng) = xic
yp (i,j,3,ng) = yic
iflux (i,j,ng) = i + ishift_br
jflux (i,j,ng) = j + jshift_br
areafact(i,j,ng) = areafac_r(i,j)
elseif (yir > c0 .and. xdr >= xcr .and. ydr < c0) then
! TR2 (group 3)
ng = 3
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xic
yp (i,j,2,ng) = yic
xp (i,j,3,ng) = xir
yp (i,j,3,ng) = yir
iflux (i,j,ng) = i + ishift_tr
jflux (i,j,ng) = j + jshift_tr
areafact(i,j,ng) = -areafac_r(i,j)
! BR2 (group 2)
ng = 2
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xdr
yp (i,j,2,ng) = ydr
xp (i,j,3,ng) = xic
yp (i,j,3,ng) = yic
iflux (i,j,ng) = i + ishift_br
jflux (i,j,ng) = j + jshift_br
areafact(i,j,ng) = -areafac_r(i,j)
endif ! TR and BR triangles
!-------------------------------------------------------------------
! Redefine departure points if not located in central cells (TC or BC)
!-------------------------------------------------------------------
if (xdl < xcl) then
xdl = xil
ydl = yil
endif
if (xdr > xcr) then
xdr = xir
ydr = yir
endif
!-------------------------------------------------------------------
! For l_fixed_area = T, shift the midpoint so that the departure
! region has the prescribed area
!-------------------------------------------------------------------
if (l_fixed_area) then
! Sum the areas of the left and right triangles.
! Note that yp(i,j,1,ng) = 0 for all triangles, so we can
! drop those terms from the area formula.
ng = 1
area1 = p5 * ( (xp(i,j,2,ng)-xp(i,j,1,ng)) * &
yp(i,j,3,ng) &
- yp(i,j,2,ng) * &
(xp(i,j,3,ng)-xp(i,j,1,ng)) ) &
* areafact(i,j,ng)
ng = 2
area2 = p5 * ( (xp(i,j,2,ng)-xp(i,j,1,ng)) * &
yp(i,j,3,ng) &
- yp(i,j,2,ng) * &
(xp(i,j,3,ng)-xp(i,j,1,ng)) ) &
* areafact(i,j,ng)
ng = 3
area3 = p5 * ( (xp(i,j,2,ng)-xp(i,j,1,ng)) * &
yp(i,j,3,ng) &
- yp(i,j,2,ng) * &
(xp(i,j,3,ng)-xp(i,j,1,ng)) ) &
* areafact(i,j,ng)
!-----------------------------------------------------------
! Check whether the central triangles lie in one grid cell or two.
! If all are in one grid cell, then adjust the area of the central
! region so that the sum of all triangle areas is equal to the
! prescribed value.
! If two triangles are in one grid cell and one is in the other,
! then compute the area of the lone triangle using an area factor
! corresponding to the adjacent corner. This is necessary to prevent
! negative masses in some rare cases on curved grids. Then adjust
! the area of the remaining two-triangle region so that the sum of
! all triangle areas has the prescribed value.
!-----------------------------------------------------------
if (ydl*ydr >= c0) then ! Both DPs lie on same side of x-axis
! compute required area of central departure region
area_c = edgearea(i,j) - area1 - area2 - area3
! shift midpoint so that the area of remaining triangles = area_c
w1 = c2*area_c/areafac_c(i,j) &
+ (xdr-xcl)*ydl + (xcr-xdl)*ydr
w2 = (xdr-xdl)**2 + (ydr-ydl)**2
w1 = w1/w2
xdm = xdm + (ydr - ydl) * w1
ydm = ydm - (xdr - xdl) * w1
! compute left and right intersection points
mdl = (ydm - ydl) / (xdm - xdl)
mdr = (ydr - ydm) / (xdr - xdm)
if (abs(mdl) > puny) then
xicl = xdl - ydl/mdl
else
xicl = c0
endif
yicl = c0
if (abs(mdr) > puny) then
xicr = xdr - ydr/mdr
else
xicr = c0
endif
yicr = c0
elseif (xic < c0) then ! fix ICL = IC
xicl = xic
yicl = yic
! compute midpoint between ICL and DR
xdm = p5 * (xdr + xicl)
ydm = p5 * ydr
! compute area of triangle adjacent to left corner
area4 = p5 * (xcl - xic) * ydl * areafac_l(i,j)
area_c = edgearea(i,j) - area1 - area2 - area3 - area4
! shift midpoint so that area of remaining triangles = area_c
w1 = c2*area_c/areafac_c(i,j) + (xcr-xic)*ydr
w2 = (xdr-xic)**2 + ydr**2
w1 = w1/w2
xdm = xdm + ydr*w1
ydm = ydm - (xdr - xic) * w1
! compute ICR
mdr = (ydr - ydm) / (xdr - xdm)
if (abs(mdr) > puny) then
xicr = xdr - ydr/mdr
else
xicr = c0
endif
yicr = c0
elseif (xic >= c0) then ! fix ICR = IR
xicr = xic
yicr = yic
! compute midpoint between ICR and DL
xdm = p5 * (xicr + xdl)
ydm = p5 * ydl
area4 = p5 * (xic - xcr) * ydr * areafac_r(i,j)
area_c = edgearea(i,j) - area1 - area2 - area3 - area4
! shift midpoint so that area of remaining triangles = area_c
w1 = c2*area_c/areafac_c(i,j) + (xic-xcl)*ydl
w2 = (xic-xdl)**2 + ydl**2
w1 = w1/w2
xdm = xdm - ydl*w1
ydm = ydm - (xic - xdl) * w1
! compute ICL
mdl = (ydm - ydl) / (xdm - xdl)
if (abs(mdl) > puny) then
xicl = xdl - ydl/mdl
else
xicl = c0
endif
yicl = c0
endif ! ydl*ydr >= c0
endif ! l_fixed_area
!-------------------------------------------------------------------
! Locate triangles in BC cell (H for both north and east edges)
! and TC cell (N for north edge and E for east edge).
!-------------------------------------------------------------------
! Start with cases where both DPs lie in the same grid cell
if (ydl >= c0 .and. ydr >= c0 .and. ydm >= c0) then
! TC1a (group 4)
ng = 4
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xcr
yp (i,j,2,ng) = ycr
xp (i,j,3,ng) = xdl
yp (i,j,3,ng) = ydl
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)
! TC2a (group 5)
ng = 5
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xdr
yp (i,j,2,ng) = ydr
xp (i,j,3,ng) = xdl
yp (i,j,3,ng) = ydl
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)
! TC3a (group 6)
ng = 6
xp (i,j,1,ng) = xdl
yp (i,j,1,ng) = ydl
xp (i,j,2,ng) = xdr
yp (i,j,2,ng) = ydr
xp (i,j,3,ng) = xdm
yp (i,j,3,ng) = ydm
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)
elseif (ydl >= c0 .and. ydr >= c0 .and. ydm < c0) then ! rare
! TC1b (group 4)
ng = 4
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xicl
yp (i,j,2,ng) = yicl
xp (i,j,3,ng) = xdl
yp (i,j,3,ng) = ydl
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)
! TC2b (group 5)
ng = 5
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xdr
yp (i,j,2,ng) = ydr
xp (i,j,3,ng) = xicr
yp (i,j,3,ng) = yicr
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)
! BC3b (group 6)
ng = 6
xp (i,j,1,ng) = xicr
yp (i,j,1,ng) = yicr
xp (i,j,2,ng) = xicl
yp (i,j,2,ng) = yicl
xp (i,j,3,ng) = xdm
yp (i,j,3,ng) = ydm
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_c(i,j)
elseif (ydl < c0 .and. ydr < c0 .and. ydm < c0) then
! BC1a (group 4)
ng = 4
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xdl
yp (i,j,2,ng) = ydl
xp (i,j,3,ng) = xcr
yp (i,j,3,ng) = ycr
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_c(i,j)
! BC2a (group 5)
ng = 5
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xdl
yp (i,j,2,ng) = ydl
xp (i,j,3,ng) = xdr
yp (i,j,3,ng) = ydr
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_c(i,j)
! BC3a (group 6)
ng = 6
xp (i,j,1,ng) = xdl
yp (i,j,1,ng) = ydl
xp (i,j,2,ng) = xdm
yp (i,j,2,ng) = ydm
xp (i,j,3,ng) = xdr
yp (i,j,3,ng) = ydr
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_c(i,j)
elseif (ydl < c0 .and. ydr < c0 .and. ydm >= c0) then ! rare
! BC1b (group 4)
ng = 4
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xdl
yp (i,j,2,ng) = ydl
xp (i,j,3,ng) = xicl
yp (i,j,3,ng) = yicl
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_c(i,j)
! BC2b (group 5)
ng = 5
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xicr
yp (i,j,2,ng) = yicr
xp (i,j,3,ng) = xdr
yp (i,j,3,ng) = ydr
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_c(i,j)
! TC3b (group 6)
ng = 6
xp (i,j,1,ng) = xicl
yp (i,j,1,ng) = yicl
xp (i,j,2,ng) = xicr
yp (i,j,2,ng) = yicr
xp (i,j,3,ng) = xdm
yp (i,j,3,ng) = ydm
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)
! Now consider cases where the two DPs lie in different grid cells
! For these cases, one triangle is given the area factor associated
! with the adjacent corner, to avoid rare negative masses on curved grids.
elseif (ydl >= c0 .and. ydr < c0 .and. xic >= c0 &
.and. ydm >= c0) then
! TC1b (group 4)
ng = 4
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xicr
yp (i,j,2,ng) = yicr
xp (i,j,3,ng) = xdl
yp (i,j,3,ng) = ydl
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)
! BC2b (group 5)
ng = 5
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xicr
yp (i,j,2,ng) = yicr
xp (i,j,3,ng) = xdr
yp (i,j,3,ng) = ydr
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_r(i,j)
! TC3b (group 6)
ng = 6
xp (i,j,1,ng) = xdl
yp (i,j,1,ng) = ydl
xp (i,j,2,ng) = xicr
yp (i,j,2,ng) = yicr
xp (i,j,3,ng) = xdm
yp (i,j,3,ng) = ydm
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)
elseif (ydl >= c0 .and. ydr < c0 .and. xic >= c0 &
.and. ydm < c0 ) then ! less common
! TC1b (group 4)
ng = 4
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xicl
yp (i,j,2,ng) = yicl
xp (i,j,3,ng) = xdl
yp (i,j,3,ng) = ydl
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)
! BC2b (group 5)
ng = 5
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xicr
yp (i,j,2,ng) = yicr
xp (i,j,3,ng) = xdr
yp (i,j,3,ng) = ydr
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_r(i,j)
! BC3b (group 6)
ng = 6
xp (i,j,1,ng) = xicr
yp (i,j,1,ng) = yicr
xp (i,j,2,ng) = xicl
yp (i,j,2,ng) = yicl
xp (i,j,3,ng) = xdm
yp (i,j,3,ng) = ydm
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_c(i,j)
elseif (ydl >= c0 .and. ydr < c0 .and. xic < c0 &
.and. ydm < c0) then
! TC1b (group 4)
ng = 4
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xicl
yp (i,j,2,ng) = yicl
xp (i,j,3,ng) = xdl
yp (i,j,3,ng) = ydl
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_l(i,j)
! BC2b (group 5)
ng = 5
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xicl
yp (i,j,2,ng) = yicl
xp (i,j,3,ng) = xdr
yp (i,j,3,ng) = ydr
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_c(i,j)
! BC3b (group 6)
ng = 6
xp (i,j,1,ng) = xdr
yp (i,j,1,ng) = ydr
xp (i,j,2,ng) = xicl
yp (i,j,2,ng) = yicl
xp (i,j,3,ng) = xdm
yp (i,j,3,ng) = ydm
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_c(i,j)
elseif (ydl >= c0 .and. ydr < c0 .and. xic < c0 &
.and. ydm >= c0) then ! less common
! TC1b (group 4)
ng = 4
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xicl
yp (i,j,2,ng) = yicl
xp (i,j,3,ng) = xdl
yp (i,j,3,ng) = ydl
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_l(i,j)
! BC2b (group 5)
ng = 5
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xicr
yp (i,j,2,ng) = yicr
xp (i,j,3,ng) = xdr
yp (i,j,3,ng) = ydr
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_c(i,j)
! TC3b (group 6)
ng = 6
xp (i,j,1,ng) = xicl
yp (i,j,1,ng) = yicl
xp (i,j,2,ng) = xicr
yp (i,j,2,ng) = yicr
xp (i,j,3,ng) = xdm
yp (i,j,3,ng) = ydm
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)
elseif (ydl < c0 .and. ydr >= c0 .and. xic < c0 &
.and. ydm >= c0) then
! BC1b (group 4)
ng = 4
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xdl
yp (i,j,2,ng) = ydl
xp (i,j,3,ng) = xicl
yp (i,j,3,ng) = yicl
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_l(i,j)
! TC2b (group 5)
ng = 5
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xdr
yp (i,j,2,ng) = ydr
xp (i,j,3,ng) = xicl
yp (i,j,3,ng) = yicl
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)
! TC3b (group 6)
ng = 6
xp (i,j,1,ng) = xicl
yp (i,j,1,ng) = yicl
xp (i,j,2,ng) = xdr
yp (i,j,2,ng) = ydr
xp (i,j,3,ng) = xdm
yp (i,j,3,ng) = ydm
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)
elseif (ydl < c0 .and. ydr >= c0 .and. xic < c0 &
.and. ydm < c0) then ! less common
! BC1b (group 4)
ng = 4
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xdl
yp (i,j,2,ng) = ydl
xp (i,j,3,ng) = xicl
yp (i,j,3,ng) = yicl
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_l(i,j)
! TC2b (group 5)
ng = 5
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xdr
yp (i,j,2,ng) = ydr
xp (i,j,3,ng) = xicr
yp (i,j,3,ng) = yicr
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)
! BC3b (group 6)
ng = 6
xp (i,j,1,ng) = xicr
yp (i,j,1,ng) = yicr
xp (i,j,2,ng) = xicl
yp (i,j,2,ng) = yicl
xp (i,j,3,ng) = xdm
yp (i,j,3,ng) = ydm
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_c(i,j)
elseif (ydl < c0 .and. ydr >= c0 .and. xic >= c0 &
.and. ydm < c0) then
! BC1b (group 4)
ng = 4
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xdl
yp (i,j,2,ng) = ydl
xp (i,j,3,ng) = xicr
yp (i,j,3,ng) = yicr
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_c(i,j)
! TC2b (group 5)
ng = 5
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xdr
yp (i,j,2,ng) = ydr
xp (i,j,3,ng) = xicr
yp (i,j,3,ng) = yicr
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_r(i,j)
! BC3b (group 6)
ng = 6
xp (i,j,1,ng) = xicr
yp (i,j,1,ng) = yicr
xp (i,j,2,ng) = xdl
yp (i,j,2,ng) = ydl
xp (i,j,3,ng) = xdm
yp (i,j,3,ng) = ydm
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_c(i,j)
elseif (ydl < c0 .and. ydr >= c0 .and. xic >= c0 &
.and. ydm >= c0) then ! less common
! BC1b (group 4)
ng = 4
xp (i,j,1,ng) = xcl
yp (i,j,1,ng) = ycl
xp (i,j,2,ng) = xdl
yp (i,j,2,ng) = ydl
xp (i,j,3,ng) = xicl
yp (i,j,3,ng) = yicl
iflux (i,j,ng) = i + ishift_bc
jflux (i,j,ng) = j + jshift_bc
areafact(i,j,ng) = areafac_c(i,j)
! TC2b (group 5)
ng = 5
xp (i,j,1,ng) = xcr
yp (i,j,1,ng) = ycr
xp (i,j,2,ng) = xdr
yp (i,j,2,ng) = ydr
xp (i,j,3,ng) = xicr
yp (i,j,3,ng) = yicr
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_r(i,j)
! TC3b (group 6)
ng = 6
xp (i,j,1,ng) = xicl
yp (i,j,1,ng) = yicl
xp (i,j,2,ng) = xicr
yp (i,j,2,ng) = yicr
xp (i,j,3,ng) = xdm
yp (i,j,3,ng) = ydm
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)
endif ! TC and BC triangles
enddo ! ij
!-------------------------------------------------------------------
! Compute triangle areas with appropriate sign.
! These are found by computing the area in scaled coordinates and
! multiplying by a scale factor (areafact).
! Note that the scale factor is positive for fluxes out of the cell
! and negative for fluxes into the cell.
!
! Note: The triangle area formula below gives A >=0 iff the triangle
! points x1, x2, and x3 are taken in counterclockwise order.
! These points are defined above in such a way that the
! order is nearly always CCW.
! In rare cases, we may compute A < 0. In this case,
! the quadrilateral departure area is equal to the
! difference of two triangle areas instead of the sum.
! The fluxes work out correctly in the end.
!
! Also compute the cumulative area transported across each edge.
! If l_fixed_area = T, this area is compared to edgearea as a bug check.
! If l_fixed_area = F, this area is passed as an output array.
!-------------------------------------------------------------------
areasum(:,:) = c0
do ng = 1, ngroups
icells(ng) = 0
do ij = 1, icellsd
i = indxid(ij)
j = indxjd(ij)
triarea(i,j,ng) = p5 * ( (xp(i,j,2,ng)-xp(i,j,1,ng)) * &
(yp(i,j,3,ng)-yp(i,j,1,ng)) &
- (yp(i,j,2,ng)-yp(i,j,1,ng)) * &
(xp(i,j,3,ng)-xp(i,j,1,ng)) ) &
* areafact(i,j,ng)
if (abs(triarea(i,j,ng)) < eps16*areafac_c(i,j)) then
triarea(i,j,ng) = c0
else
icells(ng) = icells(ng) + 1
ic = icells(ng)
indxi(ic,ng) = i
indxj(ic,ng) = j
endif
areasum(i,j) = areasum(i,j) + triarea(i,j,ng)
enddo ! ij
enddo ! ng
if (l_fixed_area) then
if (bugcheck) then ! set bugcheck = F to speed up code
do ij = 1, icellsd
i = indxid(ij)
j = indxjd(ij)
if (abs(areasum(i,j) - edgearea(i,j)) > eps13*areafac_c(i,j)) then
print*, ''
print*, 'Areas do not add up: m, i, j, edge =', &
my_task, i, j, trim(edge)
print*, 'edgearea =', edgearea(i,j)
print*, 'areasum =', areasum(i,j)
print*, 'areafac_c =', areafac_c(i,j)
print*, ''
print*, 'Triangle areas:'
do ng = 1, ngroups ! not vector friendly
if (abs(triarea(i,j,ng)) > eps16*abs(areafact(i,j,ng))) then
print*, ng, triarea(i,j,ng)
endif
enddo
endif
enddo
endif ! bugcheck
else ! l_fixed_area = F
do ij = 1, icellsd
i = indxid(ij)
j = indxjd(ij)
edgearea(i,j) = areasum(i,j)
enddo
endif ! l_fixed_area
!-------------------------------------------------------------------
! Transform triangle vertices to a scaled coordinate system centered
! in the cell containing the triangle.
!-------------------------------------------------------------------
if (trim(edge) == 'north') then
do ng = 1, ngroups
do nv = 1, nvert
do ij = 1, icells(ng)
i = indxi(ij,ng)
j = indxj(ij,ng)
ishift = iflux(i,j,ng) - i
jshift = jflux(i,j,ng) - j
xp(i,j,nv,ng) = xp(i,j,nv,ng) - c1*ishift
yp(i,j,nv,ng) = yp(i,j,nv,ng) + p5 - c1*jshift
enddo ! ij
enddo ! nv
enddo ! ng
else ! east edge
do ng = 1, ngroups
do nv = 1, nvert
do ij = 1, icells(ng)
i = indxi(ij,ng)
j = indxj(ij,ng)
ishift = iflux(i,j,ng) - i
jshift = jflux(i,j,ng) - j
! Note rotation of pi/2 here
w1 = xp(i,j,nv,ng)
xp(i,j,nv,ng) = yp(i,j,nv,ng) + p5 - c1*ishift
yp(i,j,nv,ng) = -w1 - c1*jshift
enddo ! ij
enddo ! nv
enddo ! ng
endif
if (bugcheck) then
do ng = 1, ngroups
do nv = 1, nvert
do j = jb, je
do i = ib, ie
if (abs(triarea(i,j,ng)) > puny) then
if (abs(xp(i,j,nv,ng)) > p5+puny) then
print*, ''
print*, 'WARNING: xp =', xp(i,j,nv,ng)
print*, 'm, i, j, ng, nv =', my_task, i, j, ng, nv
! print*, 'yil,xdl,xcl,ydl=',yil,xdl,xcl,ydl
! print*, 'yir,xdr,xcr,ydr=',yir,xdr,xcr,ydr
! print*, 'ydm=',ydm
! stop
endif
if (abs(yp(i,j,nv,ng)) > p5+puny) then
print*, ''
print*, 'WARNING: yp =', yp(i,j,nv,ng)
print*, 'm, i, j, ng, nv =', my_task, i, j, ng, nv
endif
endif ! triarea
enddo
enddo
enddo
enddo
endif ! bugcheck
end subroutine locate_triangles
!=======================================================================
!
!BOP
! !IROUTINE: triangle_coordinates - find coordinates of quadrature points
!
! !INTERFACE:
!
subroutine triangle_coordinates (nx_block, ny_block, & 2
integral_order, icells, &
indxi, indxj, &
xp, yp)
!
! !DESCRIPTION:
!
! For each triangle, find the coordinates of the quadrature points needed
! to compute integrals of linear, quadratic, or cubic polynomials,
! using formulas from A.H. Stroud, Approximate Calculation of Multiple
! Integrals, Prentice-Hall, 1971. (Section 8.8, formula 3.1.)
! Linear functions can be integrated exactly by evaluating the function
! at just one point (the midpoint). Quadratic functions require
! 3 points, and cubics require 4 points.
! The default is cubic, but the code can be sped up slightly using
! linear or quadratic integrals, usually with little loss of accuracy.
!
! The formulas are as follows:
!
! I1 = integral of f(x,y)*dA
! = A * f(x0,y0)
! where A is the traingle area and (x0,y0) is the midpoint.
!
! I2 = A * (f(x1,y1) + f(x2,y2) + f(x3,y3))
! where these three points are located halfway between the midpoint
! and the three vertics of the triangle.
!
! I3 = A * [ -9/16 * f(x0,y0)
! + 25/48 * (f(x1,y1) + f(x2,y2) + f(x3,y3))]
! where (x0,y0) is the midpoint, and the other three points are
! located 2/5 of the way from the midpoint to the three vertices.
!
! !REVISION HISTORY:
!
! author William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block,&! block dimensions
integral_order ! polynomial order for quadrature integrals
integer (kind=int_kind), dimension (ngroups), intent(in) :: &
icells ! number of cells where triarea > puny
integer (kind=int_kind), dimension (nx_block*ny_block,ngroups), &
intent(in) :: &
indxi ,&! compressed index in i-direction
indxj ! compressed index in j-direction
real (kind=dbl_kind), intent(inout), &
dimension (nx_block, ny_block, 0:nvert, ngroups) :: &
xp, yp ! coordinates of triangle points
!
!EOP
!
integer (kind=int_kind) :: &
i, j, ij ,&! horizontal indices
ng ! triangle index
if (integral_order == 1) then ! linear (1-point formula)
do ng = 1, ngroups
do ij = 1, icells(ng)
i = indxi(ij,ng)
j = indxj(ij,ng)
! coordinates of midpoint
xp(i,j,0,ng) = p333 &
* (xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng))
yp(i,j,0,ng) = p333 &
* (yp(i,j,1,ng) + yp(i,j,2,ng) + yp(i,j,3,ng))
enddo ! ij
enddo ! ng
elseif (integral_order == 2) then ! quadratic (3-point formula)
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ng = 1, ngroups
do ij = 1, icells(ng)
i = indxi(ij,ng)
j = indxj(ij,ng)
! coordinates of midpoint
xp(i,j,0,ng) = p333 &
* (xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng))
yp(i,j,0,ng) = p333 &
* (yp(i,j,1,ng) + yp(i,j,2,ng) + yp(i,j,3,ng))
! coordinates of the 3 points needed for integrals
xp(i,j,1,ng) = p5*xp(i,j,1,ng) + p5*xp(i,j,0,ng)
yp(i,j,1,ng) = p5*yp(i,j,1,ng) + p5*yp(i,j,0,ng)
xp(i,j,2,ng) = p5*xp(i,j,2,ng) + p5*xp(i,j,0,ng)
yp(i,j,2,ng) = p5*yp(i,j,2,ng) + p5*yp(i,j,0,ng)
xp(i,j,3,ng) = p5*xp(i,j,3,ng) + p5*xp(i,j,0,ng)
yp(i,j,3,ng) = p5*yp(i,j,3,ng) + p5*yp(i,j,0,ng)
enddo ! ij
enddo ! ng
else ! cubic (4-point formula)
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ng = 1, ngroups
do ij = 1, icells(ng)
i = indxi(ij,ng)
j = indxj(ij,ng)
! coordinates of midpoint
xp(i,j,0,ng) = p333 &
* (xp(i,j,1,ng) + xp(i,j,2,ng) + xp(i,j,3,ng))
yp(i,j,0,ng) = p333 &
* (yp(i,j,1,ng) + yp(i,j,2,ng) + yp(i,j,3,ng))
! coordinates of the other 3 points needed for integrals
xp(i,j,1,ng) = p4*xp(i,j,1,ng) + p6*xp(i,j,0,ng)
yp(i,j,1,ng) = p4*yp(i,j,1,ng) + p6*yp(i,j,0,ng)
xp(i,j,2,ng) = p4*xp(i,j,2,ng) + p6*xp(i,j,0,ng)
yp(i,j,2,ng) = p4*yp(i,j,2,ng) + p6*yp(i,j,0,ng)
xp(i,j,3,ng) = p4*xp(i,j,3,ng) + p6*xp(i,j,0,ng)
yp(i,j,3,ng) = p4*yp(i,j,3,ng) + p6*yp(i,j,0,ng)
enddo ! ij
enddo ! ng
endif
end subroutine triangle_coordinates
!=======================================================================
!
!BOP
!
! !IROUTINE: transport_integrals - compute transports across each edge
!
! !INTERFACE:
!
subroutine transport_integrals (nx_block, ny_block, & 4
ntrace, icells, &
indxi, indxj, &
tracer_type, depend, &
integral_order, triarea, &
iflux, jflux, &
xp, yp, &
mc, mx, &
my, mflx, &
tc, tx, &
ty, mtflx)
!
! !DESCRIPTION:
!
! Compute the transports across each edge by integrating the mass
! and tracers over each departure triangle.
! Input variables have the same meanings as in the main subroutine.
! Repeated use of certain sums makes the calculation more efficient.
! Integral formulas are described in triangle_coordinates subroutine.
!
! !REVISION HISTORY:
!
! author William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block ,&! block dimensions
ntrace ,&! number of tracers in use
integral_order ! polynomial order for quadrature integrals
integer (kind=int_kind), dimension (ntrace), intent(in) :: &
tracer_type ,&! = 1, 2, or 3 (see comments above)
depend ! tracer dependencies (see above)
integer (kind=int_kind), dimension (ngroups), intent(in) :: &
icells ! number of cells where triarea > puny
integer (kind=int_kind), dimension (nx_block*ny_block,ngroups), &
intent(in) :: &
indxi ,&! compressed index in i-direction
indxj ! compressed index in j-direction
real (kind=dbl_kind), intent(in), &
dimension (nx_block, ny_block, 0:nvert, ngroups) :: &
xp, yp ! coordinates of triangle points
real (kind=dbl_kind), intent(in), &
dimension (nx_block, ny_block, ngroups) :: &
triarea ! triangle area
integer (kind=int_kind), intent(in), &
dimension (nx_block, ny_block, ngroups) :: &
iflux ,&
jflux
real (kind=dbl_kind), intent(in), &
dimension (nx_block, ny_block) :: &
mc, mx, my
real (kind=dbl_kind), intent(out), &
dimension (nx_block, ny_block) :: &
mflx
real (kind=dbl_kind), intent(in), &
dimension (nx_block, ny_block, max_ntrace), optional :: &
tc, tx, ty
real (kind=dbl_kind), intent(out), &
dimension (nx_block, ny_block, max_ntrace), optional :: &
mtflx
!
!EOP
!
integer (kind=int_kind) :: &
i, j, ij ,&! horizontal indices of edge
i2, j2 ,&! horizontal indices of cell contributing transport
ng ,&! triangle index
nt, nt1 ,&! tracer indices
ilo,ihi,jlo,jhi ! beginning and end of physical domain
real (kind=dbl_kind) :: &
m0, m1, m2, m3 ,&! mass field at internal points
w0, w1, w2, w3 ! work variables
real (kind=dbl_kind), dimension (nx_block, ny_block) :: &
msum, mxsum, mysum ,&! sum of mass, mass*x, and mass*y
mxxsum, mxysum, myysum ! sum of mass*x*x, mass*x*y, mass*y*y
real (kind=dbl_kind), dimension (nx_block, ny_block, max_ntrace) :: &
mtsum ,&! sum of mass*tracer
mtxsum ,&! sum of mass*tracer*x
mtysum ! sum of mass*tracer*y
!-------------------------------------------------------------------
! Initialize
!-------------------------------------------------------------------
mflx(:,:) = c0
if (present(mtflx)) then
do nt = 1, ntrace
mtflx(:,:,nt) = c0
enddo
endif
!-------------------------------------------------------------------
! Main loop
!-------------------------------------------------------------------
do ng = 1, ngroups
if (integral_order == 1) then ! linear (1-point formula)
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells(ng)
i = indxi(ij,ng)
j = indxj(ij,ng)
i2 = iflux(i,j,ng)
j2 = jflux(i,j,ng)
! mass transports
m0 = mc(i2,j2) + xp(i,j,0,ng)*mx(i2,j2) &
+ yp(i,j,0,ng)*my(i2,j2)
msum(i,j) = m0
mflx(i,j) = mflx(i,j) + triarea(i,j,ng)*msum(i,j)
! quantities needed for tracer transports
mxsum(i,j) = m0*xp(i,j,0,ng)
mxxsum(i,j) = mxsum(i,j)*xp(i,j,0,ng)
mxysum(i,j) = mxsum(i,j)*yp(i,j,0,ng)
mysum(i,j) = m0*yp(i,j,0,ng)
myysum(i,j) = mysum(i,j)*yp(i,j,0,ng)
enddo ! ij
elseif (integral_order == 2) then ! quadratic (3-point formula)
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells(ng)
i = indxi(ij,ng)
j = indxj(ij,ng)
i2 = iflux(i,j,ng)
j2 = jflux(i,j,ng)
! mass transports
! Weighting factor of 1/3 is incorporated into the ice
! area terms m1, m2, and m3.
m1 = p333 * (mc(i2,j2) + xp(i,j,1,ng)*mx(i2,j2) &
+ yp(i,j,1,ng)*my(i2,j2))
m2 = p333 * (mc(i2,j2) + xp(i,j,2,ng)*mx(i2,j2) &
+ yp(i,j,2,ng)*my(i2,j2))
m3 = p333 * (mc(i2,j2) + xp(i,j,3,ng)*mx(i2,j2) &
+ yp(i,j,3,ng)*my(i2,j2))
msum(i,j) = m1 + m2 + m3
mflx(i,j) = mflx(i,j) + triarea(i,j,ng)*msum(i,j)
! quantities needed for mass_tracer transports
w1 = m1 * xp(i,j,1,ng)
w2 = m2 * xp(i,j,2,ng)
w3 = m3 * xp(i,j,3,ng)
mxsum(i,j) = w1 + w2 + w3
mxxsum(i,j) = w1*xp(i,j,1,ng) + w2*xp(i,j,2,ng) &
+ w3*xp(i,j,3,ng)
mxysum(i,j) = w1*yp(i,j,1,ng) + w2*yp(i,j,2,ng) &
+ w3*yp(i,j,3,ng)
w1 = m1 * yp(i,j,1,ng)
w2 = m2 * yp(i,j,2,ng)
w3 = m3 * yp(i,j,3,ng)
mysum(i,j) = w1 + w2 + w3
myysum(i,j) = w1*yp(i,j,1,ng) + w2*yp(i,j,2,ng) &
+ w3*yp(i,j,3,ng)
enddo ! ij
else ! cubic (4-point formula)
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells(ng)
i = indxi(ij,ng)
j = indxj(ij,ng)
i2 = iflux(i,j,ng)
j2 = jflux(i,j,ng)
! mass transports
! Weighting factors are incorporated into the
! terms m0, m1, m2, and m3.
m0 = p5625m * (mc(i2,j2) + xp(i,j,0,ng)*mx(i2,j2) &
+ yp(i,j,0,ng)*my(i2,j2))
m1 = p52083 * (mc(i2,j2) + xp(i,j,1,ng)*mx(i2,j2) &
+ yp(i,j,1,ng)*my(i2,j2))
m2 = p52083 * (mc(i2,j2) + xp(i,j,2,ng)*mx(i2,j2) &
+ yp(i,j,2,ng)*my(i2,j2))
m3 = p52083 * (mc(i2,j2) + xp(i,j,3,ng)*mx(i2,j2) &
+ yp(i,j,3,ng)*my(i2,j2))
msum(i,j) = m0 + m1 + m2 + m3
mflx(i,j) = mflx(i,j) + triarea(i,j,ng)*msum(i,j)
! quantities needed for tracer transports
w0 = m0 * xp(i,j,0,ng)
w1 = m1 * xp(i,j,1,ng)
w2 = m2 * xp(i,j,2,ng)
w3 = m3 * xp(i,j,3,ng)
mxsum(i,j) = w0 + w1 + w2 + w3
mxxsum(i,j) = w0*xp(i,j,0,ng) + w1*xp(i,j,1,ng) &
+ w2*xp(i,j,2,ng) + w3*xp(i,j,3,ng)
mxysum(i,j) = w0*yp(i,j,0,ng) + w1*yp(i,j,1,ng) &
+ w2*yp(i,j,2,ng) + w3*yp(i,j,3,ng)
w0 = m0 * yp(i,j,0,ng)
w1 = m1 * yp(i,j,1,ng)
w2 = m2 * yp(i,j,2,ng)
w3 = m3 * yp(i,j,3,ng)
mysum(i,j) = w0 + w1 + w2 + w3
myysum(i,j) = w0*yp(i,j,0,ng) + w1*yp(i,j,1,ng) &
+ w2*yp(i,j,2,ng) + w3*yp(i,j,3,ng)
enddo ! ij
endif ! integral_order
! mass * tracer transports
if (present(mtflx)) then
do nt = 1, ntrace
if (tracer_type(nt)==1) then ! does not depend on another tracer
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells(ng)
i = indxi(ij,ng)
j = indxj(ij,ng)
i2 = iflux(i,j,ng)
j2 = jflux(i,j,ng)
mtsum(i,j,nt) = msum(i,j) * tc(i2,j2,nt) &
+ mxsum(i,j) * tx(i2,j2,nt) &
+ mysum(i,j) * ty(i2,j2,nt)
mtflx(i,j,nt) = mtflx(i,j,nt) &
+ triarea(i,j,ng) * mtsum(i,j,nt)
! quantities needed for dependent tracers
mtxsum(i,j,nt) = mxsum(i,j) * tc(i2,j2,nt) &
+ mxxsum(i,j) * tx(i2,j2,nt) &
+ mxysum(i,j) * ty(i2,j2,nt)
mtysum(i,j,nt) = mysum(i,j) * tc(i2,j2,nt) &
+ mxysum(i,j) * tx(i2,j2,nt) &
+ myysum(i,j) * ty(i2,j2,nt)
enddo ! ij
elseif (tracer_type(nt)==2) then ! depends on another tracer
nt1 = depend(nt)
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells(ng)
i = indxi(ij,ng)
j = indxj(ij,ng)
i2 = iflux(i,j,ng)
j2 = jflux(i,j,ng)
mtsum(i,j,nt) = mtsum(i,j,nt1) * tc(i2,j2,nt) &
+ mtxsum(i,j,nt1) * tx(i2,j2,nt) &
+ mtysum(i,j,nt1) * ty(i2,j2,nt)
mtflx(i,j,nt) = mtflx(i,j,nt) &
+ triarea(i,j,ng) * mtsum(i,j,nt)
enddo ! ij
elseif (tracer_type(nt)==3) then ! depends on two tracers
nt1 = depend(nt)
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells(ng)
i = indxi(ij,ng)
j = indxj(ij,ng)
i2 = iflux(i,j,ng)
j2 = jflux(i,j,ng)
! upwind approx (tx=ty=0) for type 3 tracers
mtsum(i,j,nt) = mtsum(i,j,nt1) * tc(i2,j2,nt)
mtflx(i,j,nt) = mtflx(i,j,nt) &
+ triarea(i,j,ng) * mtsum(i,j,nt)
enddo ! ij
endif ! tracer type
enddo ! ntrace
endif ! present(mtflx)
enddo ! ng
end subroutine transport_integrals
!=======================================================================
!
!BOP
!
! !IROUTINE: update_fields - compute new area and tracers
!
! !INTERFACE:
!
subroutine update_fields (nx_block, ny_block, & 2
ilo, ihi, jlo, jhi, &
ntrace, &
tracer_type, depend, &
tarear, l_stop, &
istop, jstop, &
mflxe, mflxn, &
mm, &
mtflxe, mtflxn, &
tm)
!
! !DESCRIPTION:
!
! Given transports through cell edges, compute new area and tracers.
!
! !REVISION HISTORY:
!
! author William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block,&! block dimensions
ilo,ihi,jlo,jhi ,&! beginning and end of physical domain
ntrace ! number of tracers in use
integer (kind=int_kind), dimension (ntrace), intent(in) :: &
tracer_type ,&! = 1, 2, or 3 (see comments above)
depend ! tracer dependencies (see above)
real (kind=dbl_kind), dimension (nx_block, ny_block), &
intent(in) :: &
mflxe, mflxn ,&! mass transport across east and north cell edges
tarear ! 1/tarea
real (kind=dbl_kind), dimension (nx_block, ny_block), &
intent(inout) :: &
mm ! mass field (mean)
real (kind=dbl_kind), dimension (nx_block, ny_block, max_ntrace), &
intent(in), optional :: &
mtflxe, mtflxn ! mass*tracer transport across E and N cell edges
real (kind=dbl_kind), dimension (nx_block, ny_block, max_ntrace), &
intent(inout), optional :: &
tm ! tracer fields
logical (kind=log_kind), intent(inout) :: &
l_stop ! if true, abort on return
integer (kind=int_kind), intent(inout) :: &
istop, jstop ! indices of grid cell where model aborts
!
!EOP
!
integer (kind=int_kind) :: &
i, j ,&! horizontal indices
nt, nt1, nt2 ! tracer indices
real (kind=dbl_kind), dimension(nx_block,ny_block,max_ntrace) :: &
mtold ! old mass*tracer
real (kind=dbl_kind) :: &
w1, w2 ! work variables
integer (kind=int_kind), dimension(nx_block*ny_block) :: &
indxi ,&! compressed indices in i and j directions
indxj
integer (kind=int_kind) :: &
icells ,&! number of cells with mm > 0.
ij ! combined i/j horizontal index
!-------------------------------------------------------------------
! Save starting values of mass*tracer
!-------------------------------------------------------------------
if (present(tm)) then
do nt = 1, ntrace
if (tracer_type(nt)==1) then ! does not depend on other tracers
do j = jlo, jhi
do i = ilo, ihi
mtold(i,j,nt) = mm(i,j) * tm(i,j,nt)
enddo ! i
enddo ! j
elseif (tracer_type(nt)==2) then ! depends on another tracer
nt1 = depend(nt)
do j = jlo, jhi
do i = ilo, ihi
mtold(i,j,nt) = mm(i,j) * tm(i,j,nt1) * tm(i,j,nt)
enddo ! i
enddo ! j
elseif (tracer_type(nt)==3) then ! depends on two tracers
nt1 = depend(nt)
nt2 = depend(nt1)
do j = jlo, jhi
do i = ilo, ihi
mtold(i,j,nt) = mm(i,j) &
* tm(i,j,nt2) * tm(i,j,nt1) * tm(i,j,nt)
enddo ! i
enddo ! j
endif ! depend(nt) = 0
enddo ! nt
endif ! present(tm)
!-------------------------------------------------------------------
! Update mass field
!-------------------------------------------------------------------
do j = jlo, jhi
do i = ilo, ihi
w1 = mflxe(i,j) - mflxe(i-1,j) &
+ mflxn(i,j) - mflxn(i,j-1)
mm(i,j) = mm(i,j) - w1*tarear(i,j)
if (mm(i,j) < -puny) then ! abort with negative value
l_stop = .true.
istop = i
jstop = j
elseif (mm(i,j) < c0) then ! set to zero
mm(i,j) = c0
endif
enddo
enddo
if (l_stop) then
i = istop
j = jstop
w1 = mflxe(i,j) - mflxe(i-1,j) &
+ mflxn(i,j) - mflxn(i,j-1)
write (nu_diag,*) ' '
write (nu_diag,*) 'New mass < 0, i, j =', i, j
write (nu_diag,*) 'Old mass =', mm(i,j) + w1*tarear(i,j)
write (nu_diag,*) 'New mass =', mm(i,j)
write (nu_diag,*) 'Net transport =', -w1*tarear(i,j)
return
endif
!-------------------------------------------------------------------
! Update tracers
!-------------------------------------------------------------------
if (present(tm)) then
icells = 0
do j = jlo, jhi
do i = ilo, ihi
if (mm(i,j) > c0) then ! grid cells with positive areas
icells = icells + 1
indxi(icells) = i
indxj(icells) = j
endif
enddo ! i
enddo ! j
do nt = 1, ntrace
do j = jlo, jhi
do i = ilo, ihi
tm(i,j,nt) = c0
enddo
enddo
if (tracer_type(nt)==1) then ! does not depend on other tracers
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
w1 = mtflxe(i,j,nt) - mtflxe(i-1,j,nt) &
+ mtflxn(i,j,nt) - mtflxn(i,j-1,nt)
tm(i,j,nt) = (mtold(i,j,nt) - w1*tarear(i,j)) &
/ mm(i,j)
enddo ! ij
elseif (tracer_type(nt)==2) then ! depends on another tracer
nt1 = depend(nt)
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (abs(tm(i,j,nt1)) > c0) then
w1 = mtflxe(i,j,nt) - mtflxe(i-1,j,nt) &
+ mtflxn(i,j,nt) - mtflxn(i,j-1,nt)
tm(i,j,nt) = (mtold(i,j,nt) - w1*tarear(i,j)) &
/ (mm(i,j) * tm(i,j,nt1))
endif
enddo ! ij
elseif (tracer_type(nt)==3) then ! depends on two tracers
nt1 = depend(nt)
nt2 = depend(nt1)
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (abs(tm(i,j,nt1)) > c0 .and. &
abs(tm(i,j,nt2)) > c0) then
w1 = mtflxe(i,j,nt) - mtflxe(i-1,j,nt) &
+ mtflxn(i,j,nt) - mtflxn(i,j-1,nt)
tm(i,j,nt) = (mtold(i,j,nt) - w1*tarear(i,j)) &
/ (mm(i,j) * tm(i,j,nt2) * tm(i,j,nt1))
endif
enddo ! ij
endif ! tracer_type
enddo ! nt
endif ! present(tm)
end subroutine update_fields
!=======================================================================
end module ice_transport_remap
!=======================================================================