!=======================================================================
!BOP
!
! !MODULE: ice_transport_driver - drivers for ice transport
!
! !DESCRIPTION:
!
! Drivers for remapping and upwind ice transport
!
! !REVISION HISTORY:
! SVN:$Id: ice_transport_upwind.F 28 2006-11-03 20:32:53Z eclare $
!
! authors: Elizabeth C. Hunke and William H. Lipscomb, LANL
!
! 2004: Revised by William Lipscomb from ice_transport_mpdata.
! Stripped out mpdata, retained upwind, and added block structure.
! 2006: Incorporated remap transport driver and renamed from
! ice_transport_upwind.
!
! !INTERFACE:
module ice_transport_driver 4,5
!
! !USES:
use ice_kinds_mod
use ice_communicate, only: my_task, master_task
use ice_domain_size
use ice_constants
use ice_fileunits
!
!EOP
!
implicit none
save
character (len=char_len) :: &
advection ! type of advection scheme used
! 'upwind' => 1st order donor cell scheme
! 'remap' => remapping scheme
logical, parameter :: & ! if true, prescribe area flux across each edge
l_fixed_area = .false.
! NOTE: For remapping, hice, hsno, qice, and qsno are considered tracers.
! max_ntrace is not equal to max_ntrcr!
! ntrace is not equal to ntrcr!
integer (kind=int_kind), parameter :: &
max_ntrace = 2+max_ntrcr+nilyr+nslyr ! hice,hsno,qice,qsno,trcr
integer (kind=int_kind) :: &
ntrace ! number of tracers in use
integer (kind=int_kind), dimension (:), allocatable :: &
tracer_type ,&! = 1, 2, or 3 (see comments below)
depend ! tracer dependencies (see below)
logical (kind=log_kind), dimension (:), allocatable :: &
has_dependents ! true if a tracer has dependent tracers
integer (kind=int_kind), parameter :: &
integral_order = 3 ! polynomial order of quadrature integrals
! linear=1, quadratic=2, cubic=3
logical (kind=log_kind), parameter :: &
l_dp_midpt = .true. ! if true, find departure points using
! corrected midpoint velocity
!=======================================================================
contains
!=======================================================================
!BOP
!
! !IROUTINE: init_transport - initializations for horizontal transport
!
! !INTERFACE:
!
subroutine init_transport 1,8
!
! !DESCRIPTION:
!
! This subroutine is a wrapper for init_remap, which initializes the
! remapping transport scheme. If the model is run with upwind
! transport, no initializations are necessary.
!
! !REVISION HISTORY:
!
! authors William H. Lipscomb, LANL
!
! !USES:
!
use ice_state, only: trcr_depend, ntrcr
use ice_exit
use ice_timers
use ice_transport_remap, only: init_remap
!
!EOP
!
integer (kind=int_kind) :: &
k, nt, nt1 ! tracer indices
call ice_timer_start(timer_advect) ! advection
ntrace = 2+ntrcr+nilyr+nslyr ! hice,hsno,qice,qsno,trcr
allocate (tracer_type (ntrace), &
depend (ntrace), &
has_dependents(ntrace))
if (trim(advection)=='remap') then
!lipscomb - two branches for now; consolidate later
! define tracer dependency arrays
! see comments in remapping routine
depend(1:2) = 0 ! hice, hsno
tracer_type(1:2) = 1 ! no dependency
k = 2
do nt = 1, ntrcr
depend(k+nt) = trcr_depend(nt) ! 0 for ice area tracers
! 1 for ice volume tracers
! 2 for snow volume tracers
if (trcr_depend(nt) == 0) then
tracer_type(k+nt) = 1
else ! trcr_depend = 1 or 2
tracer_type(k+nt) = 2
endif
enddo
k = k + ntrcr
depend(k+1:k+nilyr) = 1 ! qice depends on hice
tracer_type(k+1:k+nilyr) = 2
k = k + nilyr
depend(k+1:k+nslyr) = 2 ! qsno depends on hsno
tracer_type(k+1:k+nslyr) = 2
has_dependents = .false.
do nt = 1, ntrace
if (depend(nt) > 0) then
nt1 = depend(nt)
has_dependents(nt1) = .true.
if (nt1 > nt) then
write(nu_diag,*) &
'Tracer nt2 =',nt,' depends on tracer nt1 =',nt1
call abort_ice &
('ice: remap transport: Must have nt2 > nt1')
endif
endif
enddo ! ntrace
call init_remap ! grid quantities
else ! upwind
continue
endif
call ice_timer_stop(timer_advect) ! advection
end subroutine init_transport
!=======================================================================
!BOP
!
! !IROUTINE: transport_remap - wrapper for remapping transport scheme
!
! !INTERFACE:
!
subroutine transport_remap (dt) 1,37
!
! !DESCRIPTION:
!
! This subroutine solves the transport equations for one timestep
! using the conservative 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.
!
! !REVISION HISTORY:
!
! authors William H. Lipscomb, LANL
!
! !USES:
!
use ice_boundary
use ice_global_reductions
use ice_domain
use ice_blocks
use ice_state
use ice_grid, only: tarea, HTE, HTN
use ice_exit
use ice_calendar, only: istep1, istep, diagfreq
use ice_timers
use ice_transport_remap, only: horizontal_remap, make_masks
!
! !INPUT/OUTPUT PARAMETERS:
!
real (kind=dbl_kind), intent(in) :: &
dt ! time step
!
!EOP
!
! local variables
integer (kind=int_kind) :: &
i, j ,&! horizontal indices
iblk ,&! block index
ilo,ihi,jlo,jhi,&! beginning and end of physical domain
n ,&! ice category index
nt, nt1, nt2 ! tracer indices
real (kind=dbl_kind), &
dimension (nx_block,ny_block,0:ncat,max_blocks) :: &
aim ,&! mean ice category areas in each grid cell
aimask ! = 1. if ice is present, = 0. otherwise
real (kind=dbl_kind), &
dimension (nx_block,ny_block,max_ntrace,ncat,max_blocks) :: &
trm ,&! mean tracer values in each grid cell
trmask ! = 1. if tracer is present, = 0. otherwise
logical (kind=log_kind) :: &
l_stop ! if true, abort the model
integer (kind=int_kind) :: &
istop, jstop ! indices of grid cell where model aborts
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
!-------------------------------------------------------------------
! 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.
!-------------------------------------------------------------------
real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: &
edgearea_e ,&! area of departure regions for east edges
edgearea_n ! area of departure regions for north edges
! variables related to optional bug checks
logical (kind=log_kind), parameter :: &
l_conservation_check = .true. ,&! if true, check conservation
l_monotonicity_check = .true. ! if true, check monotonicity
real (kind=dbl_kind), dimension(0:ncat) :: &
asum_init ,&! initial global ice area
asum_final ! final global ice area
real (kind=dbl_kind), dimension(max_ntrace,ncat) :: &
atsum_init ,&! initial global ice area*tracer
atsum_final ! final global ice area*tracer
real (kind=dbl_kind), dimension (:,:,:,:,:), allocatable :: &
tmin ,&! local min tracer
tmax ! local max tracer
integer (kind=int_kind) :: alloc_error
real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: &
work1
call ice_timer_start(timer_advect) ! advection
!---!-------------------------------------------------------------------
!---! Prepare for remapping.
!---! Initialize, update ghost cells, fill tracer arrays.
!---!-------------------------------------------------------------------
l_stop = .false.
istop = 0
jstop = 0
!-------------------------------------------------------------------
! Compute open water area in each grid cell.
! Note: An aggregate_area call is needed only if the open
! water area has changed since the previous call.
! Here we assume that aice0 is up to date.
!-------------------------------------------------------------------
! do iblk = 1, nblocks
! call aggregate_area (nx_block, ny_block,
! iblk, &
! aicen(:,:,:,iblk), &
! aice (:,:, iblk), &
! aice0(:,:, iblk))
! enddo
!-------------------------------------------------------------------
! Ghost cell updates for state variables.
! Commented out because ghost cells are updated after cleanup_itd.
!-------------------------------------------------------------------
! call ice_timer_start(timer_bound)
! call ice_HaloUpdate (aice0, halo_info, &
! field_loc_center, field_type_scalar)
! call bound_state (aicen, trcrn, &
! vicen, vsnon, &
! eicen, esnon)
! call ice_timer_stop(timer_bound)
!-------------------------------------------------------------------
! Ghost cell updates for ice velocity.
! Commented out because ghost cell velocities are computed
! in ice_dyn_evp.
!-------------------------------------------------------------------
! call ice_timer_start(timer_bound)
! call ice_HaloUpdate (uvel, halo_info, &
! field_loc_NEcorner, field_type_vector)
! call ice_HaloUpdate (vvel, halo_info, &
! field_loc_NEcorner, field_type_vector)
! call ice_timer_stop(timer_bound)
!$OMP PARALLEL DO PRIVATE(iblk)
do iblk = 1, nblocks
!-------------------------------------------------------------------
! Fill arrays with fields to be remapped.
!-------------------------------------------------------------------
call state_to_tracers(nx_block, ny_block, &
ntrcr, ntrace, &
aice0(:,:, iblk), &
aicen(:,:,:,iblk), trcrn(:,:,:,:,iblk), &
vicen(:,:,:,iblk), vsnon(:,:, :,iblk), &
eicen(:,:,:,iblk), esnon(:,:, :,iblk), &
aim (:,:,:,iblk), trm (:,:,:,:,iblk))
enddo
!$OMP END PARALLEL DO
!---!-------------------------------------------------------------------
!---! Optional conservation and monotonicity checks.
!---!-------------------------------------------------------------------
if (l_conservation_check .and. mod(istep,diagfreq) == 0) then
!-------------------------------------------------------------------
! Compute initial values of globally conserved quantities.
!-------------------------------------------------------------------
do n = 0, ncat
asum_init(n) = global_sum(aim(:,:,n,:), distrb_info, &
field_loc_center, tarea)
enddo
do n = 1, ncat
do nt = 1, ntrace
if (tracer_type(nt)==1) then ! does not depend on another tracer
atsum_init(nt,n) = &
global_sum_prod(trm(:,:,nt,n,:), aim(:,:,n,:), &
distrb_info, field_loc_center, &
tarea)
elseif (tracer_type(nt)==2) then ! depends on another tracer
nt1 = depend(nt)
do iblk = 1, nblocks
do j= 1,ny_block
do i = 1,nx_block
work1(i,j,iblk) = trm(i,j,nt,n,iblk)*trm(i,j,nt1,n,iblk)
end do
end do
end do
atsum_init(nt,n) = &
global_sum_prod(work1(:,:,:), aim(:,:,n,:), &
distrb_info, field_loc_center, &
tarea)
elseif (tracer_type(nt)==3) then ! depends on two tracers
nt1 = depend(nt)
nt2 = depend(nt1)
do iblk = 1, nblocks
do j= 1,ny_block
do i = 1,nx_block
work1(i,j,iblk) = trm(i,j,nt,n,iblk)*trm(i,j,nt1,n,iblk) &
*trm(i,j,nt2,n,iblk)
end do
end do
end do
atsum_init(nt,n) = &
global_sum_prod(work1(:,:,:), aim(:,:,n,:), &
distrb_info, field_loc_center, &
tarea)
endif ! tracer_type
enddo ! nt
enddo ! n
endif ! l_conservation_check
if (l_monotonicity_check .and. mod(istep,diagfreq) == 0) then
allocate(tmin(nx_block,ny_block,ntrace,ncat,max_blocks), &
tmax(nx_block,ny_block,ntrace,ncat,max_blocks), &
STAT=alloc_error)
if (alloc_error /= 0) &
call abort_ice ('ice: allocation error')
tmin(:,:,:,:,:) = c0
tmax(:,:,:,:,:) = c0
!$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.
! Masks are used to prevent tracer values in cells without ice
! from being used in the monotonicity check.
!-------------------------------------------------------------------
call make_masks (nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
nghost, ntrace, &
has_dependents, icellsnc(:,iblk), &
indxinc(:,:,iblk), indxjnc(:,:,iblk), &
aim(:,:,:,iblk), aimask(:,:,:,iblk), &
trm(:,:,:,:,iblk), trmask(:,:,:,:,iblk))
!-------------------------------------------------------------------
! Compute local max and min of tracer fields.
!-------------------------------------------------------------------
do n = 1, ncat
call local_max_min &
(nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
trm (:,:,:,n,iblk), &
tmin(:,:,:,n,iblk), tmax (:,:,:,n,iblk), &
aimask(:,:,n,iblk), trmask(:,:,:,n,iblk))
enddo
enddo
!$OMP END PARALLEL DO
call ice_timer_start(timer_bound)
call ice_HaloUpdate (tmin, halo_info, &
field_loc_center, field_type_scalar)
call ice_HaloUpdate (tmax, halo_info, &
field_loc_center, field_type_scalar)
call ice_timer_stop(timer_bound)
!$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
do n = 1, ncat
call quasilocal_max_min (nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
tmin(:,:,:,n,iblk), &
tmax(:,:,:,n,iblk))
enddo
enddo
!$OMP END PARALLEL DO
endif ! l_monotonicity_check
!-------------------------------------------------------------------
! Main remapping routine: Step ice area and tracers forward in time.
!-------------------------------------------------------------------
! If l_fixed_area is true, compute edgearea by taking the divergence
! of the velocity field. Otherwise, initialize edgearea.
!-------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblk,i,j)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
edgearea_e(i,j,iblk) = c0
edgearea_n(i,j,iblk) = c0
enddo
enddo
enddo
!$OMP END PARALLEL DO
if (l_fixed_area) then
!$OMP PARALLEL DO PRIVATE(iblk,this_block,i,j,ilo,ihi,jlo,jhi)
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
do j = jlo, jhi
do i = ilo-1, ihi
edgearea_e(i,j,iblk) = (uvel(i,j,iblk) + uvel(i,j-1,iblk)) &
* p5 * HTE(i,j,iblk) * dt
enddo
enddo
do j = jlo-1, jhi
do i = ilo, ihi
edgearea_n(i,j,iblk) = (vvel(i,j,iblk) + vvel(i-1,j,iblk)) &
* p5 * HTN(i,j,iblk) * dt
enddo
enddo
enddo ! iblk
!$OMP END PARALLEL DO
endif
call horizontal_remap (dt, ntrace, &
uvel (:,:,:), vvel (:,:,:), &
aim (:,:,:,:), trm (:,:,:,:,:), &
l_fixed_area, &
edgearea_e(:,:,:), edgearea_n(:,:,:), &
tracer_type, depend, &
has_dependents, integral_order, &
l_dp_midpt)
!-------------------------------------------------------------------
! Given new fields, recompute state variables.
!-------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblk)
do iblk = 1, nblocks
call tracers_to_state (nx_block, ny_block, &
ntrcr, ntrace, &
aim (:,:,:,iblk), trm (:,:,:,:,iblk), &
aice0(:,:, iblk), &
aicen(:,:,:,iblk), trcrn(:,:,:,:,iblk), &
vicen(:,:,:,iblk), vsnon(:,:, :,iblk), &
eicen(:,:,:,iblk), esnon(:,:, :,iblk))
enddo ! iblk
!$OMP END PARALLEL DO
!-------------------------------------------------------------------
! Ghost cell updates for state variables.
!-------------------------------------------------------------------
call ice_timer_start(timer_bound)
call bound_state (aicen, trcrn, &
vicen, vsnon, &
eicen, esnon)
call ice_timer_stop(timer_bound)
!---!-------------------------------------------------------------------
!---! Optional conservation and monotonicity checks
!---!-------------------------------------------------------------------
!-------------------------------------------------------------------
! Compute final values of globally conserved quantities.
! Check global conservation of area and area*tracers. (Optional)
!-------------------------------------------------------------------
if (l_conservation_check .and. mod(istep,diagfreq) == 0) then
do n = 0, ncat
asum_final(n) = global_sum(aim(:,:,n,:), distrb_info, &
field_loc_center, tarea)
enddo
do n = 1, ncat
do nt = 1, ntrace
if (tracer_type(nt)==1) then ! does not depend on another tracer
atsum_final(nt,n) = &
global_sum_prod(trm(:,:,nt,n,:), aim(:,:,n,:), &
distrb_info, field_loc_center, &
tarea)
elseif (tracer_type(nt)==2) then ! depends on another tracer
nt1 = depend(nt)
do iblk = 1, nblocks
do j= 1,ny_block
do i = 1,nx_block
work1(i,j,iblk) = trm(i,j,nt,n,iblk)*trm(i,j,nt1,n,iblk)
end do
end do
end do
atsum_final(nt,n) = &
global_sum_prod(work1(:,:,:), aim(:,:,n,:), &
distrb_info, field_loc_center, &
tarea)
elseif (tracer_type(nt)==3) then ! depends on two tracers
nt1 = depend(nt)
nt2 = depend(nt1)
do iblk = 1, nblocks
do j= 1,ny_block
do i = 1,nx_block
work1(i,j,iblk) = trm(i,j,nt,n,iblk)*trm(i,j,nt1,n,iblk) &
*trm(i,j,nt2,n,iblk)
end do
end do
end do
atsum_final(nt,n) = &
global_sum_prod(work1(:,:,:), aim(:,:,n,:), &
distrb_info, field_loc_center, &
tarea)
endif ! tracer_type
enddo ! nt
enddo ! n
if (my_task == master_task) then
call global_conservation (l_stop, &
asum_init(0), asum_final(0))
if (l_stop) then
write (nu_diag,*) 'istep1 =', istep1
write (nu_diag,*) 'transport: conservation error, cat 0'
l_stop = .false.
call abort_ice('ice remap transport: conservation error')
endif
do n = 1, ncat
call global_conservation &
(l_stop, &
asum_init(n), asum_final(n), &
atsum_init(:,n), atsum_final(:,n))
if (l_stop) then
write (nu_diag,*) 'istep1, cat =', &
istep1, n
write (nu_diag,*) 'transport: conservation error, cat ',n
l_stop = .false.
call abort_ice &
('ice remap transport: conservation error')
endif
enddo ! n
endif ! my_task = master_task
endif ! l_conservation_check
!-------------------------------------------------------------------
! Check tracer monotonicity. (Optional)
!-------------------------------------------------------------------
if (l_monotonicity_check .and. mod(istep,diagfreq) == 0) then
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
do n = 1, ncat
call check_monotonicity &
(nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
iblk, &
tmin(:,:,:,n,iblk), tmax(:,:,:,n,iblk), &
aim (:,:, n,iblk), trm (:,:,:,n,iblk), &
l_stop, &
istop, jstop)
if (l_stop) then
write (nu_diag,*) 'istep1, my_task, iblk, cat =', &
istep1, my_task, iblk, n
write (nu_diag,*) 'i_glob, j_glob',this_block%i_glob(istop), &
this_block%j_glob(jstop)
call abort_ice('ice remap transport: monotonicity error')
endif
enddo ! n
enddo ! iblk
deallocate(tmin, tmax, STAT=alloc_error)
if (alloc_error /= 0) call abort_ice ('deallocation error')
endif ! l_monotonicity_check
call ice_timer_stop(timer_advect) ! advection
end subroutine transport_remap
!=======================================================================
!BOP
!
! !IROUTINE: transport_upwind - upwind transport
!
! !INTERFACE:
!
subroutine transport_upwind (dt) 1,22
!
! !DESCRIPTION:
!
! Computes the transport equations for one timestep using upwind. Sets
! several fields into a work array and passes it to upwind routine.
!
! !REVISION HISTORY:
!
! same as module
!
! !USES:
!
use ice_boundary
use ice_blocks
use ice_domain
use ice_state
use ice_grid, only: HTE, HTN, tarea
use ice_timers
!
! !INPUT/OUTPUT PARAMETERS:
!
real (kind=dbl_kind), intent(in) :: &
dt ! time step
!
!EOP
!
integer (kind=int_kind) :: &
narr ! number of state variable arrays
! not including eicen, esnon
integer (kind=int_kind) :: &
i, j, iblk ,&! horizontal indices
ilo,ihi,jlo,jhi ! beginning and end of physical domain
real (kind=dbl_kind), dimension (nx_block,ny_block,nblocks) :: &
uee, vnn ! cell edge velocities
real (kind=dbl_kind), &
dimension (:,:,:,:), allocatable :: &
works ! work array
type (block) :: &
this_block ! block information for current block
call ice_timer_start(timer_advect) ! advection
narr = 1 + ncat*(3+ntrcr) ! max number of state variable arrays
! not including eicen, esnon
allocate (works(nx_block,ny_block,narr,max_blocks))
!-------------------------------------------------------------------
! Get ghost cell values of state variables.
! (Assume velocities are already known for ghost cells, also.)
!-------------------------------------------------------------------
! call bound_state (aicen, trcrn, &
! vicen, vsnon, &
! eicen, esnon)
uee(:,:,:) = c0
vnn(:,:,:) = c0
!-------------------------------------------------------------------
! Average corner velocities to edges.
!-------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblk,this_block,i,j,ilo,ihi,jlo,jhi)
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
do j = jlo, jhi
do i = ilo, ihi
uee(i,j,iblk) = p5*(uvel(i,j,iblk) + uvel(i,j-1,iblk))
vnn(i,j,iblk) = p5*(vvel(i,j,iblk) + vvel(i-1,j,iblk))
enddo
enddo
enddo
!$OMP END PARALLEL DO
call ice_timer_start(timer_bound)
call ice_HaloUpdate (uee, halo_info, &
field_loc_Eface, field_type_vector)
call ice_HaloUpdate (vnn, halo_info, &
field_loc_Nface, field_type_vector)
call ice_timer_stop(timer_bound)
!$OMP PARALLEL DO PRIVATE(iblk,this_block,ilo,ihi,jlo,jhi)
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
!-----------------------------------------------------------------
! fill work arrays with fields to be advected
!-----------------------------------------------------------------
call state_to_work (nx_block, ny_block, &
ntrcr, &
narr, trcr_depend, &
aicen (:,:, :,iblk), trcrn (:,:,:,:,iblk), &
vicen (:,:, :,iblk), vsnon (:,:, :,iblk), &
aice0 (:,:, iblk), works (:,:, :,iblk))
!-----------------------------------------------------------------
! advect
!-----------------------------------------------------------------
call upwind_field (nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
dt, &
narr, works(:,:,:,iblk), &
uee(:,:,iblk), vnn (:,:,iblk), &
HTE(:,:,iblk), HTN (:,:,iblk), &
tarea(:,:,iblk))
call upwind_field (nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
dt, &
ntilyr, eicen(:,:,:,iblk), &
uee(:,:,iblk), vnn (:,:,iblk), &
HTE(:,:,iblk), HTN (:,:,iblk), &
tarea(:,:,iblk))
call upwind_field (nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
dt, &
ntslyr, esnon(:,:,:,iblk), &
uee(:,:,iblk), vnn (:,:,iblk), &
HTE(:,:,iblk), HTN (:,:,iblk), &
tarea(:,:,iblk))
!-----------------------------------------------------------------
! convert work arrays back to state variables
!-----------------------------------------------------------------
call work_to_state (nx_block, ny_block, &
ntrcr, &
narr, trcr_depend, &
aicen(:,:, :,iblk), trcrn (:,:,:,:,iblk), &
vicen(:,:, :,iblk), vsnon (:,:, :,iblk), &
aice0(:,:, iblk), works (:,:, :,iblk))
enddo ! iblk
!$OMP END PARALLEL DO
!-------------------------------------------------------------------
! Ghost cell updates for state variables.
!-------------------------------------------------------------------
call ice_timer_start(timer_bound)
call bound_state (aicen, trcrn, &
vicen, vsnon, &
eicen, esnon)
call ice_timer_stop(timer_bound)
call ice_timer_stop(timer_advect) ! advection
deallocate(works)
end subroutine transport_upwind
!=======================================================================
! The next few subroutines (through check_monotonicity) are called
! by transport_remap.
!=======================================================================
!
!BOP
!
! !IROUTINE: state_to_tracers -fill ice area and tracer arrays
!
! !INTERFACE:
!
subroutine state_to_tracers (nx_block, ny_block, & 1,1
ntrcr, ntrace, &
aice0, &
aicen, trcrn, &
vicen, vsnon, &
eicen, esnon, &
aim, trm)
!
! !DESCRIPTION:
!
! Fill ice area and tracer arrays.
! Assume that the advected tracers are hicen, hsnon, trcrn,
! qicen(1:nilyr), and qsnon(1:nslyr).
! This subroutine must be modified if a different set of tracers
! is to be transported. The rule for ordering tracers
! is that a dependent tracer (such as qice) must have a larger
! tracer index than the tracer it depends on (i.e., hice).
!
! !REVISION HISTORY:
!
! author William H. Lipscomb, LANL
!
! !USES:
!
use ice_itd, only: ilyr1, slyr1
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
ntrcr , & ! number of tracers in use
ntrace ! number of tracers in use incl. hi, hs
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(in) :: &
aice0 ! fractional open water area
real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), &
intent(in) :: &
aicen ,&! fractional ice area
vicen ,&! volume per unit area of ice (m)
vsnon ! volume per unit area of snow (m)
real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrcr,ncat), &
intent(in) :: &
trcrn ! ice area tracers
real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr), &
intent(in) :: &
eicen ! energy of melting for each ice layer (J/m^2)
real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr), &
intent(in) :: &
esnon ! energy of melting for each snow layer (J/m^2)
real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat), &
intent(out):: &
aim ! mean ice area in each grid cell
real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrace,ncat), &
intent(out) :: &
trm ! mean tracer values in each grid cell
real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
worka, &
workb
!
!EOP
!
integer (kind=int_kind) :: &
i, j, k, n ,&! standard indices
it, kt ,&! tracer indices
ij ! combined i/j index
real (kind=dbl_kind) :: &
w1 ! work variable
integer (kind=int_kind), dimension(nx_block*ny_block,0:ncat) :: &
indxi ,&! compressed i/j indices
indxj
integer (kind=int_kind), dimension(0:ncat) :: &
icells ! number of cells with ice
worka(:,:) = c0
workb(:,:) = c0
aim(:,:,0) = aice0(:,:)
do n = 1, ncat
trm(:,:,:,n) = c0
!-------------------------------------------------------------------
! Find grid cells where ice is present and fill area array.
!-------------------------------------------------------------------
icells(n) = 0
do j = 1, ny_block
do i = 1, nx_block
aim(i,j,n) = aicen(i,j,n)
if (aim(i,j,n) > puny) then
icells(n) = icells(n) + 1
ij = icells(n)
indxi(ij,n) = i
indxj(ij,n) = j
endif ! aim > puny
enddo
enddo
!-------------------------------------------------------------------
! Fill tracer array
! Note: If aice > 0, then hice > 0, but we can have hsno = 0.
! Alse note: We transport qice*nilyr rather than qice, so as to
! avoid extra operations here and in tracers_to_state.
!-------------------------------------------------------------------
do ij = 1, icells(n)
i = indxi(ij,n)
j = indxj(ij,n)
w1 = c1 / aim(i,j,n)
worka(i,j) = c1 / vicen(i,j,n)
trm(i,j,1,n) = vicen(i,j,n) * w1 ! hice
trm(i,j,2,n) = vsnon(i,j,n) * w1 ! hsno
if (trm(i,j,2,n) > puny) then
workb(i,j) = c1 / vsnon(i,j,n)
else
workb(i,j) = c0
endif
enddo
kt = 2
do it = 1, ntrcr
do ij = 1, icells(n)
i = indxi(ij,n)
j = indxj(ij,n)
trm(i,j,kt+it,n) = trcrn(i,j,it,n) ! ice area tracers
enddo
enddo
kt = kt + ntrcr
do k =1, nilyr
do ij = 1, icells(n)
i = indxi(ij,n)
j = indxj(ij,n)
trm(i,j,kt+k,n) = eicen(i,j,ilyr1(n)+k-1)*worka(i,j) ! qice
enddo ! ij
enddo ! ilyr
kt = kt + nilyr
do k = 1, nslyr
do ij = 1, icells(n)
i = indxi(ij,n)
j = indxj(ij,n)
if (trm(i,j,2,n) > puny) & ! hsno > puny
trm(i,j,kt+k,n) = esnon(i,j,slyr1(n)+k-1)*workb(i,j) & ! qsno
+ rhos*Lfresh
enddo ! ij
enddo ! nslyr
enddo ! ncat
end subroutine state_to_tracers
!=======================================================================
!BOP
!
! !IROUTINE: tracers_to_state - convert tracer array to state variables
!
! !INTERFACE:
!
subroutine tracers_to_state (nx_block, ny_block, & 1,1
ntrcr, ntrace, &
aim, trm, &
aice0, &
aicen, trcrn, &
vicen, vsnon, &
eicen, esnon)
!
! !DESCRIPTION:
!
! Convert area and tracer arrays back to state variables.
!
! !REVISION HISTORY:
!
! author William H. Lipscomb, LANL
!
! !USES:
!
use ice_itd, only: ilyr1, slyr1
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
ntrcr , & ! number of tracers in use
ntrace ! number of tracers in use incl. hi, hs
real (kind=dbl_kind), dimension (nx_block,ny_block,0:ncat), &
intent(in) :: &
aim ! fractional ice area
real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrace,ncat), &
intent(in) :: &
trm ! mean tracer values in each grid cell
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(inout) :: &
aice0 ! fractional ice area
real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), &
intent(inout) :: &
aicen ,&! fractional ice area
vicen ,&! volume per unit area of ice (m)
vsnon ! volume per unit area of snow (m)
real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrcr,ncat), &
intent(inout) :: &
trcrn ! tracers
real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr), &
intent(inout) :: &
eicen ! energy of melting for each ice layer (J/m^2)
real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr), &
intent(inout) :: &
esnon ! energy of melting for each snow layer (J/m^2)
!
!EOP
!
integer (kind=int_kind) :: &
i, j, k, n ,&! standard indices
it, kt ,&! tracer indices
icells ,&! number of cells with ice
ij
integer (kind=int_kind), dimension (nx_block*ny_block) :: &
indxi, indxj ! compressed indices
aice0(:,:) = aim(:,:,0)
do n = 1, ncat
icells = 0
do j = 1, ny_block
do i = 1, nx_block
if (aim(i,j,n) > c0) then
icells = icells + 1
indxi(icells) = i
indxj(icells) = j
endif
enddo
enddo
!-------------------------------------------------------------------
! Compute state variables.
!-------------------------------------------------------------------
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
aicen(i,j,n) = aim(i,j,n)
vicen(i,j,n) = aim(i,j,n)*trm(i,j,1,n) ! aice*hice
vsnon(i,j,n) = aim(i,j,n)*trm(i,j,2,n) ! aice*hsno
enddo ! ij
kt = 2
do it = 1, ntrcr
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
trcrn(i,j,it,n) = trm(i,j,kt+it,n) ! ice tracers
enddo ! ij
enddo ! ntrcr
kt = kt + ntrcr
do k = 1, nilyr
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
eicen(i,j,ilyr1(n)+k-1) = vicen(i,j,n)*trm(i,j,kt+k,n)
enddo ! ij
enddo ! nilyr
kt = kt + nilyr
do k = 1, nslyr
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
esnon(i,j,slyr1(n)+k-1) = (trm(i,j,kt+k,n) - rhos*Lfresh) &
* vsnon(i,j,n)
enddo ! ij
enddo ! nslyr
enddo ! ncat
end subroutine tracers_to_state
!=======================================================================
!
!BOP
!
! !IROUTINE: global_conservation - check for changes in conserved quantities
!
! !INTERFACE:
!
subroutine global_conservation (l_stop, & 2
asum_init, asum_final, &
atsum_init, atsum_final)
!
! !DESCRIPTION:
!
! Check whether values of conserved quantities have changed.
! An error probably means that ghost cells are treated incorrectly.
!
! !REVISION HISTORY:
!
! author William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
real (kind=dbl_kind), intent(in) :: &
asum_init ,&! initial global ice area
asum_final ! final global ice area
real (kind=dbl_kind), dimension(max_ntrace), intent(in), optional :: &
atsum_init ,&! initial global ice area*tracer
atsum_final ! final global ice area*tracer
logical (kind=log_kind), intent(inout) :: &
l_stop ! if true, abort on return
!
!EOP
!
integer (kind=int_kind) :: &
nt ! tracer index
real (kind=dbl_kind) :: &
diff ! difference between initial and final values
if (asum_init > puny) then
diff = asum_final - asum_init
if (abs(diff/asum_init) > puny) then
l_stop = .true.
write (nu_diag,*)
write (nu_diag,*) 'Ice area conserv error'
write (nu_diag,*) 'Initial global area =', asum_init
write (nu_diag,*) 'Final global area =', asum_final
write (nu_diag,*) 'Fractional error =', abs(diff)/asum_init
write (nu_diag,*) 'asum_final-asum_init =', diff
endif
endif
if (present(atsum_init)) then
do nt = 1, ntrace
if (abs(atsum_init(nt)) > puny) then
diff = atsum_final(nt) - atsum_init(nt)
if (abs(diff/atsum_init(nt)) > puny) then
l_stop = .true.
write (nu_diag,*)
write (nu_diag,*) 'area*tracer conserv error'
write (nu_diag,*) 'tracer index =', nt
write (nu_diag,*) 'Initial global area*tracer =', &
atsum_init(nt)
write (nu_diag,*) 'Final global area*tracer =', &
atsum_final(nt)
write (nu_diag,*) 'Fractional error =', &
abs(diff)/atsum_init(nt)
write (nu_diag,*) 'atsum_final-atsum_init =', diff
endif
endif
enddo
endif ! present(atsum_init)
end subroutine global_conservation
!=======================================================================
!BOP
!
! !IROUTINE: local_max_min - compute local max and min of a scalar field
!
! !INTERFACE:
!
subroutine local_max_min (nx_block, ny_block, & 1
ilo, ihi, jlo, jhi, &
trm, &
tmin, tmax, &
aimask, trmask)
!
! !DESCRIPTION:
!
! At each grid point, compute the local max and min of a scalar
! field phi: i.e., the max and min values in the nine-cell region
! consisting of the home cell and its eight neighbors.
!
! To extend to the neighbors of the neighbors (25 cells in all),
! follow this call with a call to quasilocal_max_min.
!
! !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
real (kind=dbl_kind), intent(in), &
dimension(nx_block,ny_block) :: &
aimask ! ice area mask
real (kind=dbl_kind), intent(in), &
dimension (nx_block,ny_block,max_ntrace) :: &
trm ,&! tracer fields
trmask ! tracer mask
real (kind=dbl_kind), intent(out), &
dimension (nx_block,ny_block,max_ntrace) :: &
tmin ,&! local min tracer
tmax ! local max tracer
!
!EOP
!
integer (kind=int_kind) :: &
i, j ,&! horizontal indices
nt, nt1 ! tracer indices
real (kind=dbl_kind), dimension(nx_block,ny_block) :: &
phimask ! aimask or trmask, as appropriate
real (kind=dbl_kind) :: &
phi_nw, phi_n, phi_ne ,&! field values in 8 neighbor cells
phi_w, phi_e ,&
phi_sw, phi_s, phi_se
do nt = 1, ntrace
if (tracer_type(nt)==1) then ! does not depend on another tracer
do j = 1, ny_block
do i = 1, nx_block
phimask(i,j) = aimask(i,j)
enddo
enddo
else ! depends on another tracer
nt1 = depend(nt)
do j = 1, ny_block
do i = 1, nx_block
phimask(i,j) = trmask(i,j,nt1)
enddo
enddo
endif
!-----------------------------------------------------------------------
! Store values of trm in the 8 neighbor cells.
! If aimask = 1, use the true value; otherwise use the home cell value
! so that non-physical values of phi do not contribute to the gradient.
!-----------------------------------------------------------------------
do j = jlo, jhi
do i = ilo, ihi
phi_nw = phimask(i-1,j+1) * trm(i-1,j+1,nt) &
+ (c1-phimask(i-1,j+1))* trm(i, j, nt)
phi_n = phimask(i, j+1) * trm(i, j+1,nt) &
+ (c1-phimask(i, j+1))* trm(i, j, nt)
phi_ne = phimask(i+1,j+1) * trm(i+1,j+1,nt) &
+ (c1-phimask(i+1,j+1))* trm(i, j, nt)
phi_w = phimask(i-1,j) * trm(i-1,j, nt) &
+ (c1-phimask(i-1,j)) * trm(i, j, nt)
phi_e = phimask(i+1,j) * trm(i+1,j, nt) &
+ (c1-phimask(i+1,j)) * trm(i, j, nt)
phi_sw = phimask(i-1,j-1) * trm(i-1,j-1,nt) &
+ (c1-phimask(i-1,j-1))* trm(i, j, nt)
phi_s = phimask(i, j-1) * trm(i, j-1,nt) &
+ (c1-phimask(i, j-1))* trm(i, j, nt)
phi_se = phimask(i+1,j-1) * trm(i+1,j-1,nt) &
+ (c1-phimask(i+1,j-1))* trm(i, j, nt)
!-----------------------------------------------------------------------
! Compute the minimum and maximum among the nine local cells.
!-----------------------------------------------------------------------
tmax(i,j,nt) = max (phi_nw, phi_n, phi_ne, phi_w, &
trm(i,j,nt), phi_e, phi_sw, phi_s, phi_se)
tmin(i,j,nt) = min (phi_nw, phi_n, phi_ne, phi_w, &
trm(i,j,nt), phi_e, phi_sw, phi_s, phi_se)
enddo ! i
enddo ! j
enddo ! nt
end subroutine local_max_min
!=======================================================================
!BOP
!
! !IROUTINE: quasilocal_max_min - look one grid cell farther away
!
! !INTERFACE:
!
subroutine quasilocal_max_min (nx_block, ny_block, & 1
ilo, ihi, jlo, jhi, &
tmin, tmax)
!
! !DESCRIPTION:
!
! Extend the local max and min by one grid cell in each direction.
! Incremental remapping is monotone for the "quasilocal" max and min,
! but in rare cases may violate monotonicity for the local max and min.
!
! !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
real (kind=dbl_kind), intent(inout), &
dimension (nx_block,ny_block,ntrace) :: &
tmin ,&! local min tracer
tmax ! local max tracer
!
!EOP
!
integer (kind=int_kind) :: &
i, j ,&! horizontal indices
nt ! tracer index
do nt = 1, ntrace
do j = jlo, jhi
do i = ilo, ihi
tmax(i,j,nt) = &
max (tmax(i-1,j+1,nt), tmax(i,j+1,nt), tmax(i+1,j+1,nt), &
tmax(i-1,j, nt), tmax(i,j, nt), tmax(i+1,j, nt), &
tmax(i-1,j-1,nt), tmax(i,j-1,nt), tmax(i+1,j-1,nt))
tmin(i,j,nt) = &
min (tmin(i-1,j+1,nt), tmin(i,j+1,nt), tmin(i+1,j+1,nt), &
tmin(i-1,j, nt), tmin(i,j, nt), tmin(i+1,j, nt), &
tmin(i-1,j-1,nt), tmin(i,j-1,nt), tmin(i+1,j-1,nt))
enddo ! i
enddo ! j
enddo
end subroutine quasilocal_max_min
!======================================================================
!
!BOP
!
! !IROUTINE: check_monotonicity - check bounds on new tracer values
!
! !INTERFACE:
!
subroutine check_monotonicity (nx_block, ny_block, & 1
ilo, ihi, jlo, jhi, &
iblk, &
tmin, tmax, &
aim, trm, &
l_stop, &
istop, jstop)
!
! !DESCRIPTION:
!
! At each grid point, make sure that the new tracer values
! fall between the local max and min values before transport.
!
! !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
iblk ! block index (diagnostic only)
real (kind=dbl_kind), intent(in), &
dimension (nx_block,ny_block) :: &
aim ! new ice area
real (kind=dbl_kind), intent(in), &
dimension (nx_block,ny_block,max_ntrace) :: &
trm ! new tracers
real (kind=dbl_kind), intent(in), &
dimension (nx_block,ny_block,ntrace) :: &
tmin ,&! local min tracer
tmax ! local max tracer
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) :: &
w1, w2 ! work variables
logical (kind=log_kind), dimension (nx_block, ny_block) :: &
l_check ! if true, check monotonicity
do nt = 1, ntrace
!-------------------------------------------------------------------
! Load logical array to identify tracers that need checking.
!-------------------------------------------------------------------
if (tracer_type(nt)==1) then ! does not depend on another tracer
do j = jlo, jhi
do i = ilo, ihi
if (aim(i,j) > puny) then
l_check(i,j) = .true.
else
l_check(i,j) = .false.
endif
enddo
enddo
elseif (tracer_type(nt)==2) then ! depends on another tracer
nt1 = depend(nt)
do j = jlo, jhi
do i = ilo, ihi
if (abs(trm(i,j,nt1)) > puny .and. aim(i,j) > puny) then
l_check(i,j) = .true.
else
l_check(i,j) = .false.
endif
enddo
enddo
elseif (tracer_type(nt)==3) then ! depends on two tracers
nt1 = depend(nt)
nt2 = depend(nt1)
do j = jlo, jhi
do i = ilo, ihi
if (abs(trm(i,j,nt1)) > puny .and. &
abs(trm(i,j,nt2)) > puny .and. aim(i,j) > puny) then
l_check(i,j) = .true.
else
l_check(i,j) = .false.
endif
enddo
enddo
endif
!-------------------------------------------------------------------
! Make sure new values lie between tmin and tmax
!-------------------------------------------------------------------
do j = jlo, jhi
do i = ilo, ihi
if (l_check(i,j)) then
! w1 and w2 allow for roundoff error when abs(trm) is big
w1 = max(c1, abs(tmin(i,j,nt)))
w2 = max(c1, abs(tmax(i,j,nt)))
if (trm(i,j,nt) < tmin(i,j,nt)-w1*puny) then
l_stop = .true.
istop = i
jstop = j
write (nu_diag,*) ' '
write (nu_diag,*) 'new tracer < tmin'
write (nu_diag,*) 'i, j, nt =', i, j, nt
write (nu_diag,*) 'new tracer =', trm (i,j,nt)
write (nu_diag,*) 'tmin =' , tmin(i,j,nt)
write (nu_diag,*) 'ice area =' , aim(i,j)
elseif (trm(i,j,nt) > tmax(i,j,nt)+w2*puny) then
l_stop = .true.
istop = i
jstop = j
write (nu_diag,*) ' '
write (nu_diag,*) 'new tracer > tmax'
write (nu_diag,*) 'i, j, nt =', i, j, nt
write (nu_diag,*) 'new tracer =', trm (i,j,nt)
write (nu_diag,*) 'tmax =' , tmax(i,j,nt)
write (nu_diag,*) 'ice area =' , aim(i,j)
endif
endif
enddo ! i
enddo ! j
enddo ! nt
end subroutine check_monotonicity
!=======================================================================
! The remaining subroutines are called by transport_upwind.
!=======================================================================
!BOP
!
! !IROUTINE: state_to_work - fill work arrays with state variables
!
! !INTERFACE:
!
subroutine state_to_work (nx_block, ny_block, & 1,1
ntrcr, &
narr, trcr_depend, &
aicen, trcrn, &
vicen, vsnon, &
aice0, works)
!
! !DESCRIPTION:
!
! Fill work array with state variables in preparation for upwind transport
!
! !REVISION HISTORY:
!
! same as module
!
! !USES:
!
use ice_itd, only: ilyr1, slyr1
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block ,&! block dimensions
ntrcr , & ! number of tracers in use
narr ! number of 2D state variable arrays in works array
integer (kind=int_kind), dimension (max_ntrcr), intent(in) :: &
trcr_depend ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), &
intent(in) :: &
aicen ,&! concentration of ice
vicen ,&! volume per unit area of ice (m)
vsnon ! volume per unit area of snow (m)
real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrcr,ncat), &
intent(in) :: &
trcrn ! ice tracers
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(in) :: &
aice0 ! concentration of open water
real (kind=dbl_kind), dimension(nx_block,ny_block,narr), &
intent (out) :: &
works ! work array
!
!EOP
!
integer (kind=int_kind) :: &
i, j, k, n, it ,&! counting indices
narrays ! counter for number of state variable arrays
!-----------------------------------------------------------------
! This array is used for performance (balance memory/cache vs
! number of bound calls); a different number of arrays may perform
! better depending on the machine used, number of processors, etc.
! --tested on SGI R2000, using 4 pes for the ice model under MPI
!-----------------------------------------------------------------
do j = 1, ny_block
do i = 1, nx_block
works(i,j,1) = aice0(i,j)
enddo
enddo
narrays = 1
do n=1, ncat
do j = 1, ny_block
do i = 1, nx_block
works(i,j,narrays+1) = aicen(i,j,n)
works(i,j,narrays+2) = vicen(i,j,n)
works(i,j,narrays+3) = vsnon(i,j,n)
enddo ! i
enddo ! j
narrays = narrays + 3
do it = 1, ntrcr
if (trcr_depend(it) == 0) then
do j = 1, ny_block
do i = 1, nx_block
works(i,j,narrays+it) = aicen(i,j,n)*trcrn(i,j,it,n)
enddo
enddo
elseif (trcr_depend(it) == 1) then
do j = 1, ny_block
do i = 1, nx_block
works(i,j,narrays+it) = vicen(i,j,n)*trcrn(i,j,it,n)
enddo
enddo
elseif (trcr_depend(it) == 2) then
do j = 1, ny_block
do i = 1, nx_block
works(i,j,narrays+it) = vsnon(i,j,n)*trcrn(i,j,it,n)
enddo
enddo
endif
enddo
narrays = narrays + ntrcr
enddo ! n
if (narr /= narrays) write(nu_diag,*) &
"Wrong number of arrays in transport bound call"
end subroutine state_to_work
!=======================================================================
!BOP
!
! !IROUTINE: work_to_state - convert work arrays back to state variables
!
! !INTERFACE:
!
subroutine work_to_state (nx_block, ny_block, & 1,2
ntrcr, &
narr, trcr_depend, &
aicen, trcrn, &
vicen, vsnon, &
aice0, works)
!
! !DESCRIPTION:
!
! Convert work array back to state variables
!
! !REVISION HISTORY:
!
! same as module
!
! !USES:
!
use ice_itd, only: ilyr1, slyr1, compute_tracers
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent (in) :: &
nx_block, ny_block, & ! block dimensions
ntrcr , & ! number of tracers in use
narr ! number of 2D state variable arrays in works array
integer (kind=int_kind), dimension (max_ntrcr), intent(in) :: &
trcr_depend ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
real (kind=dbl_kind), intent (in) :: &
works (nx_block,ny_block,narr)
real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), &
intent(out) :: &
aicen ,&! concentration of ice
vicen ,&! volume per unit area of ice (m)
vsnon ! volume per unit area of snow (m)
real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrcr,ncat), &
intent(out) :: &
trcrn ! ice tracers
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(out) :: &
aice0 ! concentration of open water
!
!EOP
!
integer (kind=int_kind) :: &
i, j, k, n , it,&! counting indices
narrays ,&! counter for number of state variable arrays
icells ! number of ocean/ice cells
integer (kind=int_kind), dimension (nx_block*ny_block) :: &
indxi, indxj
real (kind=dbl_kind), dimension (nx_block*ny_block,narr) :: &
work
! for call to compute_tracers
icells = 0
do j = 1, ny_block
do i = 1, nx_block
icells = icells + 1
indxi(icells) = i
indxj(icells) = j
work (icells,:) = works(i,j,:)
enddo
enddo
do j=1,ny_block
do i=1,nx_block
aice0(i,j) = works(i,j,1)
enddo
enddo
narrays = 1 ! aice0 is first array
do n=1,ncat
do j=1,ny_block
do i=1,nx_block
aicen(i,j,n) = works(i,j,narrays+1)
vicen(i,j,n) = works(i,j,narrays+2)
vsnon(i,j,n) = works(i,j,narrays+3)
enddo
enddo
narrays = narrays + 3
call compute_tracers (nx_block, ny_block, &
icells, indxi, indxj, &
ntrcr, trcr_depend, &
work (:,narrays+1:narrays+ntrcr), &
aicen(:,:,n), &
vicen(:,:,n), vsnon(:,:,n), &
trcrn(:,:,:,n))
narrays = narrays + ntrcr
enddo ! ncat
end subroutine work_to_state
!=======================================================================
!BOP
!
! !IROUTINE: upwind_field - advection according to upwind
!
! !INTERFACE:
!
subroutine upwind_field (nx_block, ny_block, & 3,6
ilo, ihi, jlo, jhi, &
dt, &
narrays, phi, &
uee, vnn, &
HTE, HTN, &
tarea)
!
! !DESCRIPTION:
!
! upwind transport algorithm
!
! !REVISION HISTORY:
!
! same as module
!
! !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
narrays ! number of 2D arrays to be transported
real (kind=dbl_kind), intent(in) :: &
dt ! time step
real (kind=dbl_kind), dimension(nx_block,ny_block,narrays), &
intent(inout) :: &
phi ! scalar field
real (kind=dbl_kind), dimension(nx_block,ny_block), &
intent(in):: &
uee, vnn ! cell edge velocities
real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: &
HTE ,&! length of east cell edge
HTN ,&! length of north cell edge
tarea ! grid cell area
!
!EOP
!
integer (kind=int_kind) :: &
i, j, k, n ! standard indices
real (kind=dbl_kind) :: &
upwind, y1, y2, a, h ! function
real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
worka, &
workb
!-------------------------------------------------------------------
! Define upwind function
!-------------------------------------------------------------------
upwind(y1,y2,a,h) = p5*dt*h*((a+abs(a))*y1+(a-abs(a))*y2)
!-------------------------------------------------------------------
! upwind transport
!-------------------------------------------------------------------
worka(:,:) = c0
workb(:,:) = c0
do n = 1, narrays
do j = 1, jhi
do i = 1, ihi
worka(i,j)= &
upwind(phi(i,j,n),phi(i+1,j,n),uee(i,j),HTE(i,j))
workb(i,j)= &
upwind(phi(i,j,n),phi(i,j+1,n),vnn(i,j),HTN(i,j))
enddo
enddo
do j = jlo, jhi
do i = ilo, ihi
phi(i,j,n) = phi(i,j,n) - ( worka(i,j)-worka(i-1,j) &
+ workb(i,j)-workb(i,j-1) ) &
/ tarea(i,j)
enddo
enddo
enddo ! narrays
end subroutine upwind_field
!=======================================================================
end module ice_transport_driver
!=======================================================================