!=======================================================================
!BOP
!
! !MODULE: ice_itd - initialize and redistribute ice in the ITD
!
! !DESCRIPTION:
!
! Routines to initialize the ice thickness distribution and
! utilities to redistribute ice among categories. These routines
! are not specific to a particular numerical implementation.
!
! See Bitz, C.M., and W.H. Lipscomb, 1999:
! An energy-conserving thermodynamic model of sea ice,
! J. Geophys. Res., 104, 15,669--15,677.
!
! See Bitz, C.M., M.M. Holland, A.J. Weaver, M. Eby, 2001:
! Simulating the ice-thickness distribution in a climate model,
! J. Geophys. Res., 106, 2441--2464.
!
! !REVISION HISTORY:
! SVN:$Id: ice_itd.F90 138 2008-07-08 20:39:37Z eclare $
!
! authors: C. M. Bitz, UW
! William H. Lipscomb and Elizabeth C. Hunke, LANL
!
! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
!
! 2004 WHL: Added multiple snow layers, block structure, cleanup_itd
! 2006 ECH: Added WMO standard ice thickness categories as kcatbound=2
! Streamlined for efficiency
! Converted to free source form (F90)
!
! !INTERFACE:
!
module ice_itd 26,6
!
! !USES:
!
use ice_kinds_mod
use ice_communicate, only: my_task, master_task
use ice_domain_size
use ice_constants
use ice_fileunits
use ice_exit
!
!EOP
!
implicit none
save
integer (kind=int_kind) :: &
kitd , & ! type of itd conversions
! 0 = delta function
! 1 = linear remap
kcatbound , & ! 0 = old category boundary formula
! 1 = new formula giving round numbers
! 2 = WMO standard
ilyr1 (ncat), & ! array position of top ice layer in each cat
ilyrn (ncat), & ! array position of bottom ice layer in each cat
slyr1 (ncat), & ! array position of top snow layer in each cat
slyrn (ncat) ! array position of bottom snow layer in each cat
real (kind=dbl_kind), parameter :: &
hi_min = p01 ! minimum ice thickness allowed (m)
real (kind=dbl_kind) :: &
hin_max(0:ncat) ! category limits (m)
character (len=35) :: c_hi_range(ncat)
!-------------------------------------------------------------------
! a note regarding hi_min and hin_max(0):
! both represent a minimum ice thickness. hin_max(0) is
! intended to be used for particular numerical implementations
! of category conversions in the ice thickness distribution.
! hi_min is a more general purpose parameter, but is specifically
! for maintaining stability in the thermodynamics.
! hin_max(0) = 0.1 m for the delta function itd
! hin_max(0) = 0.0 m for linear remapping
!
! Also note that the upper limit on the thickest category
! is only used for the linear remapping scheme
! and it is not a true upper limit on the thickness
!-------------------------------------------------------------------
!=======================================================================
contains
!=======================================================================
!BOP
!
! !IROUTINE: init_itd - initalize area fraction and thickness boundaries for ITD
!
! !INTERFACE:
!
subroutine init_itd 1
!
! !DESCRIPTION:
!
! Initialize area fraction and thickness boundaries for the itd model
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL
! C. M. Bitz, UW
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
integer (kind=int_kind) :: &
n ! thickness category index
real (kind=dbl_kind) :: &
cc1, cc2, cc3, & ! parameters for kcatbound = 0
x1 , &
rn , & ! real(n)
rncat , & ! real(ncat)
d1 , & ! parameters for kcatbound = 1 (m)
d2
real (kind=dbl_kind), dimension(5) :: wmo5 ! data for wmo itd
real (kind=dbl_kind), dimension(6) :: wmo6 ! data for wmo itd
real (kind=dbl_kind), dimension(7) :: wmo7 ! data for wmo itd
character(len=8) :: c_hinmax1,c_hinmax2
character(len=2) :: c_nc
rncat = real(ncat, kind=dbl_kind)
d1 = 3.0_dbl_kind / rncat
d2 = 0.5_dbl_kind / rncat
!-----------------------------------------------------------------
! Choose category boundaries based on one of three options.
!
! The first formula (kcatbound = 0) was used in Lipscomb (2001)
! and in CICE versions 3.0 and 3.1.
!
! The second formula is more user-friendly in the sense that it
! is easy to obtain round numbers for category boundaries:
!
! H(n) = n * [d1 + d2*(n-1)]
!
! Default values are d1 = 300/ncat, d2 = 50/ncat.
! For ncat = 5, boundaries in cm are 60, 140, 240, 360, which are
! close to the standard values given by the first formula.
! For ncat = 10, boundaries in cm are 30, 70, 120, 180, 250, 330,
! 420, 520, 630.
!
! The third option provides support for World Meteorological
! Organization classification based on thickness. The full
! WMO thickness distribution is used if ncat = 7; if ncat=5
! or ncat = 6, some of the thinner categories are combined.
! For ncat = 5, boundaries are 30, 70, 120, 200, >200 cm.
! For ncat = 6, boundaries are 15, 30, 70, 120, 200, >200 cm.
! For ncat = 7, boundaries are 10, 15, 30, 70, 120, 200, >200 cm.
!-----------------------------------------------------------------
if (kcatbound == 0) then ! original scheme
if (kitd == 1) then
! linear remapping itd category limits
cc1 = c3/rncat
cc2 = c15*cc1
cc3 = c3
hin_max(0) = c0 ! minimum ice thickness, m
else
! delta function itd category limits
cc1 = max(1.1_dbl_kind/rncat,c1*hi_min)
cc2 = c25*cc1
cc3 = 2.25_dbl_kind
! hin_max(0) should not be zero
! use some caution in making it less than 0.10
hin_max(0) = hi_min ! minimum ice thickness, m
endif ! kitd
do n = 1, ncat
x1 = real(n-1,kind=dbl_kind) / rncat
hin_max(n) = hin_max(n-1) &
+ cc1 + cc2*(c1 + tanh(cc3*(x1-c1)))
enddo
elseif (kcatbound == 1) then ! new scheme
hin_max(0) = c0
do n = 1, ncat
rn = real(n, kind=dbl_kind)
hin_max(n) = rn * (d1 + (rn-c1)*d2)
enddo
elseif (kcatbound == 2) then ! WMO standard
if (ncat == 5) then
! thinnest 3 categories combined
data wmo5 / 0.30_dbl_kind, 0.70_dbl_kind, &
1.20_dbl_kind, 2.00_dbl_kind, &
999._dbl_kind /
hin_max(0) = c0
do n = 1, ncat
hin_max(n) = wmo5(n)
enddo
elseif (ncat == 6) then
! thinnest 2 categories combined
data wmo6 / 0.15_dbl_kind, &
0.30_dbl_kind, 0.70_dbl_kind, &
1.20_dbl_kind, 2.00_dbl_kind, &
999._dbl_kind /
hin_max(0) = c0
do n = 1, ncat
hin_max(n) = wmo6(n)
enddo
elseif (ncat == 7) then
! all thickness categories
data wmo7 / 0.10_dbl_kind, 0.15_dbl_kind, &
0.30_dbl_kind, 0.70_dbl_kind, &
1.20_dbl_kind, 2.00_dbl_kind, &
999._dbl_kind /
hin_max(0) = c0
do n = 1, ncat
hin_max(n) = wmo7(n)
enddo
else
write (nu_diag,*) 'kcatbound=3 (WMO) must have ncat=5, 6 or 7'
stop
endif
endif ! kcatbound
if (my_task == master_task) then
write (nu_diag,*) ' '
write (nu_diag,*) 'hin_max(n-1) < Cat n < hin_max(n)'
do n = 1, ncat
write (nu_diag,*) hin_max(n-1),' < Cat ',n, ' < ',hin_max(n)
! Write integer n to character string
write (c_nc, '(i2)') n
! Write hin_max to character string
write (c_hinmax1, '(f5.3)') hin_max(n-1)
write (c_hinmax2, '(f5.3)') hin_max(n)
! Save character string to write to history file
c_hi_range(n)=c_hinmax1//'m < hi Cat '//c_nc//' < '// &
c_hinmax2//'m'
enddo
write (nu_diag,*) ' '
endif
!-----------------------------------------------------------------
! vectors identifying first and last layer in each category
!-----------------------------------------------------------------
ilyr1(1) = 1 ! if nilyr = 4
ilyrn(1) = nilyr ! ilyr1 = { 1,5,9 }
do n = 2, ncat ! ilyrn = { 4,8,12} etc
ilyr1(n) = ilyrn(n-1) + 1
ilyrn(n) = ilyrn(n-1) + nilyr
enddo
slyr1(1) = 1
slyrn(1) = nslyr
do n = 2, ncat
slyr1(n) = slyrn(n-1) + 1
slyrn(n) = slyrn(n-1) + nslyr
enddo
end subroutine init_itd
!=======================================================================
!BOP
!
! !IROUTINE: aggregate - aggregate ice state variables
!
! !INTERFACE:
!
subroutine aggregate (nx_block, ny_block, & 6,1
aicen, trcrn, &
vicen, vsnon, &
eicen, esnon, &
aice, trcr, &
vice, vsno, &
eice, esno, &
aice0, tmask, &
ntrcr, trcr_depend)
!
! !DESCRIPTION:
!
! Aggregate ice state variables over thickness categories.
!
! !REVISION HISTORY:
!
! authors: C. M. Bitz, UW
! W. H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
ntrcr ! number of tracers in use
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,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)
logical (kind=log_kind), dimension (nx_block,ny_block), &
intent(in) :: &
tmask ! land/boundary mask, thickness (T-cell)
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), &
intent(out) :: &
aice , & ! concentration of ice
vice , & ! volume per unit area of ice (m)
vsno , & ! volume per unit area of snow (m)
eice , & ! energy of melt. of ice (J/m^2)
esno , & ! energy of melt. of snow layer (J/m^2)
aice0 ! concentration of open water
real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrcr), &
intent(out) :: &
trcr ! ice tracers
!
!EOP
!
integer (kind=int_kind) :: &
icells ! number of ocean/ice cells
integer (kind=int_kind), dimension (nx_block*ny_block) :: &
indxi, & ! compressed indices in i/j directions
indxj
integer (kind=int_kind) :: &
i, j, k, n, it, &
ij ! combined i/j horizontal index
real (kind=dbl_kind), dimension (:,:), allocatable :: &
atrcr ! sum of aicen*trcrn or vicen*trcrn or vsnon*trcrn
!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------
icells = 0
do j = 1, ny_block
do i = 1, nx_block
if (tmask(i,j)) then
icells = icells + 1
indxi(icells) = i
indxj(icells) = j
endif ! tmask
aice0(i,j) = c1
aice (i,j) = c0
vice (i,j) = c0
vsno (i,j) = c0
eice (i,j) = c0
esno (i,j) = c0
enddo
enddo
allocate (atrcr(icells,ntrcr))
!-----------------------------------------------------------------
! Aggregate
!-----------------------------------------------------------------
atrcr(:,:) = c0
do n = 1, ncat
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
aice(i,j) = aice(i,j) + aicen(i,j,n)
vice(i,j) = vice(i,j) + vicen(i,j,n)
vsno(i,j) = vsno(i,j) + vsnon(i,j,n)
enddo ! ij
do it = 1, ntrcr
if (trcr_depend(it) == 0) then ! ice area tracer
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
atrcr(ij,it) = atrcr(ij,it) &
+ trcrn(i,j,it,n)*aicen(i,j,n)
enddo ! ij
elseif (trcr_depend(it) == 1) then ! ice volume tracer
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
atrcr(ij,it) = atrcr(ij,it) &
+ trcrn(i,j,it,n)*vicen(i,j,n)
enddo ! ij
elseif (trcr_depend(it) ==2) then ! snow volume tracer
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
atrcr(ij,it) = atrcr(ij,it) &
+ trcrn(i,j,it,n)*vsnon(i,j,n)
enddo ! ij
endif ! trcr_depend
enddo ! ntrcr
do k = 1, nilyr
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
eice(i,j) = eice(i,j) + eicen(i,j,ilyr1(n)+k-1)
enddo
enddo ! nilyr
do k = 1, nslyr
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
esno(i,j) = esno(i,j) + esnon(i,j,slyr1(n)+k-1)
enddo
enddo ! nslyr
enddo ! ncat
! Open water fraction
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
aice0(i,j) = max (c1 - aice(i,j), c0)
enddo ! ij
! Tracers
call compute_tracers (nx_block, ny_block, &
icells, indxi, indxj, &
ntrcr, trcr_depend, &
atrcr(:,:), aice(:,:), &
vice (:,:), vsno(:,:), &
trcr(:,:,:))
deallocate (atrcr)
end subroutine aggregate
!=======================================================================
!BOP
!
! !IROUTINE: aggregate_area - aggregate ice area
!
! !INTERFACE:
!
subroutine aggregate_area (nx_block, ny_block, & 3
aicen, aice, aice0)
!
! !DESCRIPTION:
!
! Aggregate ice area (but not other state variables) over thickness
! categories.
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb, LANL
! modified Jan 2004 by Clifford Chen, Fujitsu
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block ! block dimensions
real (kind=dbl_kind), dimension (:,:,:), intent(in) :: &
aicen ! concentration of ice
real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
aice, & ! concentration of ice
aice0 ! concentration of open water
!
!EOP
!
integer (kind=int_kind) :: i, j, n
!-----------------------------------------------------------------
! Aggregate
!-----------------------------------------------------------------
aice(:,:) = c0
do n = 1, ncat
do j = 1, ny_block
do i = 1, nx_block
aice(i,j) = aice(i,j) + aicen(i,j,n)
enddo ! i
enddo ! j
enddo ! n
do j = 1, ny_block
do i = 1, nx_block
! open water fraction
aice0(i,j) = max (c1 - aice(i,j), c0)
enddo ! i
enddo ! j
end subroutine aggregate_area
!=======================================================================
!BOP
!
! !IROUTINE: rebin - rebins thicknesses into defined categories
!
! !INTERFACE:
!
subroutine rebin (nx_block, ny_block, & 3,2
icells, indxi, indxj, &
ntrcr, trcr_depend, &
aicen, trcrn, &
vicen, vsnon, &
eicen, esnon, &
l_stop, &
istop, jstop)
!
! !DESCRIPTION:
!
! Rebins thicknesses into defined categories
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
icells , & ! number of grid cells with ice
ntrcr ! number of tracers in use
integer (kind=int_kind), dimension (nx_block*ny_block), &
intent(in) :: &
indxi, indxj ! compressed i/j indices
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(inout) :: &
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(inout) :: &
trcrn ! ice 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)
logical (kind=log_kind), intent(out) :: &
l_stop ! if true, abort on return
integer (kind=int_kind), intent(out) :: &
istop, jstop ! indices of grid cell where model aborts
!
!EOP
!
integer (kind=int_kind) :: &
i,j , & ! horizontal indices
n , & ! category index
ij ! combined horizontal index
logical (kind=log_kind) :: &
shiftflag ! = .true. if ice must be shifted
integer (kind=int_kind), dimension (icells,ncat) :: &
donor ! donor category index
real (kind=dbl_kind), dimension (icells,ncat) :: &
daice , & ! ice area transferred
dvice , & ! ice volume transferred
hicen ! ice thickness for each cat (m)
!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------
l_stop = .false.
istop = 0
jstop = 0
do n = 1, ncat
do ij = 1, icells ! aice(i,j) > puny
i = indxi(ij)
j = indxj(ij)
donor(ij,n) = 0
daice(ij,n) = c0
dvice(ij,n) = c0
!-----------------------------------------------------------------
! Compute ice thickness.
!-----------------------------------------------------------------
if (aicen(i,j,n) > puny) then
hicen(ij,n) = vicen(i,j,n) / aicen(i,j,n)
else
hicen(ij,n) = c0
endif
enddo ! ij
enddo ! n
!-----------------------------------------------------------------
! make sure thickness of cat 1 is at least hin_max(0)
!-----------------------------------------------------------------
do ij = 1, icells ! aice(i,j) > puny
i = indxi(ij)
j = indxj(ij)
if (aicen(i,j,1) > puny) then
if (hicen(ij,1) <= hin_max(0) .and. hin_max(0) > c0 ) then
aicen(i,j,1) = vicen(i,j,1) / hin_max(0)
hicen(ij,1) = hin_max(0)
endif
endif
enddo ! ij
!-----------------------------------------------------------------
! If a category thickness is not in bounds, shift the
! entire area, volume, and energy to the neighboring category
!-----------------------------------------------------------------
!-----------------------------------------------------------------
! Move thin categories up
!-----------------------------------------------------------------
do n = 1, ncat-1 ! loop over category boundaries
!-----------------------------------------------------------------
! identify thicknesses that are too big
!-----------------------------------------------------------------
shiftflag = .false.
do ij = 1, icells ! aice(i,j) > puny
i = indxi(ij)
j = indxj(ij)
if (aicen(i,j,n) > puny .and. &
hicen(ij,n) > hin_max(n)) then
shiftflag = .true.
donor(ij,n) = n
daice(ij,n) = aicen(i,j,n)
dvice(ij,n) = vicen(i,j,n)
endif
enddo ! ij
if (shiftflag) then
!-----------------------------------------------------------------
! shift ice between categories
!-----------------------------------------------------------------
call shift_ice (nx_block, ny_block, &
indxi, indxj, &
icells, &
ntrcr, trcr_depend, &
aicen, trcrn, &
vicen, vsnon, &
eicen, esnon, &
hicen, donor, &
daice, dvice, &
l_stop, &
istop, jstop)
if (l_stop) return
!-----------------------------------------------------------------
! reset shift parameters
!-----------------------------------------------------------------
do ij = 1, icells ! aice(i,j) > puny
donor(ij,n) = 0
daice(ij,n) = c0
dvice(ij,n) = c0
enddo
endif ! shiftflag
enddo ! n
!-----------------------------------------------------------------
! Move thick categories down
!-----------------------------------------------------------------
do n = ncat-1, 1, -1 ! loop over category boundaries
!-----------------------------------------------------------------
! identify thicknesses that are too small
!-----------------------------------------------------------------
shiftflag = .false.
do ij = 1, icells ! aice(i,j) > puny
i = indxi(ij)
j = indxj(ij)
if (aicen(i,j,n+1) > puny .and. &
hicen(ij,n+1) <= hin_max(n)) then
shiftflag = .true.
donor(ij,n) = n+1
daice(ij,n) = aicen(i,j,n+1)
dvice(ij,n) = vicen(i,j,n+1)
endif
enddo ! ij
if (shiftflag) then
!-----------------------------------------------------------------
! shift ice between categories
!-----------------------------------------------------------------
call shift_ice (nx_block, ny_block, &
indxi, indxj, &
icells, &
ntrcr, trcr_depend, &
aicen, trcrn, &
vicen, vsnon, &
eicen, esnon, &
hicen, donor, &
daice, dvice, &
l_stop, &
istop, jstop)
if (l_stop) return
!-----------------------------------------------------------------
! reset shift parameters
!-----------------------------------------------------------------
do ij = 1, icells ! aice(i,j) > puny
donor(ij,n) = 0
daice(ij,n) = c0
dvice(ij,n) = c0
enddo
endif ! shiftflag
enddo ! n
end subroutine rebin
!=======================================================================
!BOP
!
! !IROUTINE: reduce_area - reduce area when ice melts for special case ncat=1
!
! !INTERFACE:
!
subroutine reduce_area (nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
tmask, &
aicen, vicen, &
aicen_init,vicen_init)
!
! !DESCRIPTION:
!
! Reduce area when ice melts for special case of ncat=1
!
! Use CSM 1.0-like method of reducing ice area
! when melting occurs: assume only half the ice volume
! change goes to thickness decrease, the other half
! to reduction in ice fraction
!
! !REVISION HISTORY:
!
! authors: C. M. Bitz, UW
! modified by: Elizabeth C. Hunke, 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
logical (kind=log_kind), dimension (nx_block,ny_block), &
intent(in) :: &
tmask ! land/boundary mask, thickness (T-cell)
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(inout) :: &
aicen , & ! concentration of ice
vicen ! volume per unit area of ice (m)
real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: &
aicen_init, & ! old ice area for category 1 (m)
vicen_init ! old ice volume for category 1 (m)
!
!EOP
!
integer (kind=int_kind) :: &
i, j ! horizontal indices
real (kind=dbl_kind) :: &
hi0 , & ! initial hi
hi1 , & ! current hi
dhi ! hi1 - hi0
do j = jlo, jhi
do i = ilo, ihi
if (tmask(i,j)) then
hi0 = c0
if (aicen_init(i,j) > c0) &
hi0 = vicen_init(i,j) / aicen_init(i,j)
hi1 = c0
if (aicen(i,j) > c0) &
hi1 = vicen(i,j) / aicen(i,j)
! make sure thickness of cat 1 is at least hin_max(0)
if (hi1 <= hin_max(0) .and. hin_max(0) > c0 ) then
aicen(i,j) = vicen(i,j) / hin_max(0)
hi1 = hin_max(0)
endif
if (aicen(i,j) > c0) then
dhi = hi1 - hi0
if (dhi < c0) then
hi1 = vicen(i,j) / aicen(i,j)
aicen(i,j) = c2 * vicen(i,j) / (hi1 + hi0)
endif
endif
endif ! tmask
enddo ! i
enddo ! j
end subroutine reduce_area
!=======================================================================
!BOP
!
! !IROUTINE: shift_ice - shift ice across category boundaries
!
! !INTERFACE:
!
subroutine shift_ice (nx_block, ny_block, & 3,1
indxi, indxj, &
icells, &
ntrcr, trcr_depend, &
aicen, trcrn, &
vicen, vsnon, &
eicen, esnon, &
hicen, donor, &
daice, dvice, &
l_stop, &
istop, jstop)
!
! !DESCRIPTION:
!
! Shift ice across category boundaries, conserving area, volume, and
! energy.
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
icells , & ! number of ocean/ice cells
ntrcr ! number of tracers in use
integer (kind=int_kind), dimension (nx_block*ny_block), &
intent(in) :: &
indxi , & ! compressed indices in i/j directions
indxj
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(inout) :: &
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(inout) :: &
trcrn ! ice 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)
! NOTE: Third index of donor, daice, dvice should be ncat-1,
! except that compilers would have trouble when ncat = 1
integer (kind=int_kind), dimension(icells,ncat), &
intent(in) :: &
donor ! donor category index
real (kind=dbl_kind), dimension(icells,ncat), &
intent(inout) :: &
daice , & ! ice area transferred across boundary
dvice , & ! ice volume transferred across boundary
hicen ! ice thickness for each cat (m)
logical (kind=log_kind), intent(out) :: &
l_stop ! if true, abort on return
integer (kind=int_kind), intent(out) :: &
istop, jstop ! indices of grid cell where model aborts
!
!EOP
!
integer (kind=int_kind) :: &
i, j, m , & ! horizontal indices
n , & ! thickness category index
nr , & ! receiver category
nd , & ! donor category
k , & ! ice layer index
it , & ! tracer index
ilo,ihi,jlo,jhi ! beginning and end of physical domain
real (kind=dbl_kind), dimension(icells,max_ntrcr,ncat) :: &
atrcrn ! aicen*trcrn
real (kind=dbl_kind) :: &
dvsnow , & ! snow volume transferred
desnow , & ! snow energy transferred
deice , & ! ice energy transferred
datrcr ! aicen*train transferred
integer (kind=int_kind), dimension (icells) :: &
indxii , & ! compressed indices for i/j directions
indxjj , &
indxij
integer (kind=int_kind) :: &
ishift , & ! number of cells with ice to transfer
ij ! combined i/j horizontal index
logical (kind=log_kind) :: &
daice_negative , & ! true if daice < -puny
dvice_negative , & ! true if dvice < -puny
daice_greater_aicen, & ! true if daice > aicen
dvice_greater_vicen ! true if dvice > vicen
real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
worka, &
workb
!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------
l_stop = .false.
istop = 0
jstop = 0
worka(:,:) = c0
workb(:,:) = c0
!-----------------------------------------------------------------
! Define variables equal to aicen*trcrn, vicen*trcrn, vsnon*trcrn
!-----------------------------------------------------------------
do n = 1, ncat
do it = 1, ntrcr
if (trcr_depend(it) == 0) then ! ice area tracer
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
atrcrn(ij,it,n) = aicen(i,j,n)*trcrn(i,j,it,n)
enddo
elseif (trcr_depend(it) ==1) then ! ice volume tracer
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
atrcrn(ij,it,n) = vicen(i,j,n)*trcrn(i,j,it,n)
enddo
elseif (trcr_depend(it) ==2) then ! snow volume tracer
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
atrcrn(ij,it,n) = vsnon(i,j,n)*trcrn(i,j,it,n)
enddo
endif
enddo
enddo
!-----------------------------------------------------------------
! Check for daice or dvice out of range, allowing for roundoff error
!-----------------------------------------------------------------
do n = 1, ncat-1
daice_negative = .false.
dvice_negative = .false.
daice_greater_aicen = .false.
dvice_greater_vicen = .false.
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (donor(ij,n) > 0) then
nd = donor(ij,n)
if (daice(ij,n) < c0) then
if (daice(ij,n) > -puny*aicen(i,j,nd)) then
daice(ij,n) = c0 ! shift no ice
dvice(ij,n) = c0
else
daice_negative = .true.
endif
endif
if (dvice(ij,n) < c0) then
if (dvice(ij,n) > -puny*vicen(i,j,nd)) then
daice(ij,n) = c0 ! shift no ice
dvice(ij,n) = c0
else
dvice_negative = .true.
endif
endif
if (daice(ij,n) > aicen(i,j,nd)*(c1-puny)) then
if (daice(ij,n) < aicen(i,j,nd)*(c1+puny)) then
daice(ij,n) = aicen(i,j,nd)
dvice(ij,n) = vicen(i,j,nd)
else
daice_greater_aicen = .true.
endif
endif
if (dvice(ij,n) > vicen(i,j,nd)*(c1-puny)) then
if (dvice(ij,n) < vicen(i,j,nd)*(c1+puny)) then
daice(ij,n) = aicen(i,j,nd)
dvice(ij,n) = vicen(i,j,nd)
else
dvice_greater_vicen = .true.
endif
endif
endif ! donor > 0
enddo ! ij
!-----------------------------------------------------------------
! error messages
!-----------------------------------------------------------------
if (daice_negative) then
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (donor(ij,n) > 0 .and. &
daice(ij,n) <= -puny*aicen(i,j,nd)) then
write(nu_diag,*) ' '
write(nu_diag,*) 'shift_ice: negative daice'
write(nu_diag,*) 'i, j:', i, j
write(nu_diag,*) 'boundary, donor cat:', n, nd
write(nu_diag,*) 'daice =', daice(ij,n)
write(nu_diag,*) 'dvice =', dvice(ij,n)
l_stop = .true.
istop = i
jstop = j
endif
enddo
endif
if (l_stop) return
if (dvice_negative) then
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (donor(ij,n) > 0 .and. &
dvice(ij,n) <= -puny*vicen(i,j,nd)) then
write(nu_diag,*) ' '
write(nu_diag,*) 'shift_ice: negative dvice'
write(nu_diag,*) 'i, j:', i, j
write(nu_diag,*) 'boundary, donor cat:', n, nd
write(nu_diag,*) 'daice =', daice(ij,n)
write(nu_diag,*) 'dvice =', dvice(ij,n)
l_stop = .true.
istop = i
jstop = j
endif
enddo
endif
if (l_stop) return
if (daice_greater_aicen) then
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (donor(ij,n) > 0) then
nd = donor(ij,n)
if (daice(ij,n) >= aicen(i,j,nd)*(c1+puny)) then
write(nu_diag,*) ' '
write(nu_diag,*) 'shift_ice: daice > aicen'
write(nu_diag,*) 'i, j:', i, j
write(nu_diag,*) 'boundary, donor cat:', n, nd
write(nu_diag,*) 'daice =', daice(ij,n)
write(nu_diag,*) 'aicen =', aicen(i,j,nd)
l_stop = .true.
istop = i
jstop = j
endif
endif
enddo
endif
if (l_stop) return
if (dvice_greater_vicen) then
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (donor(ij,n) > 0) then
nd = donor(ij,n)
if (dvice(ij,n) >= vicen(i,j,nd)*(c1+puny)) then
write(nu_diag,*) ' '
write(nu_diag,*) 'shift_ice: dvice > vicen'
write(nu_diag,*) 'i, j:', i, j
write(nu_diag,*) 'boundary, donor cat:', n, nd
write(nu_diag,*) 'dvice =', dvice(ij,n)
write(nu_diag,*) 'vicen =', vicen(i,j,nd)
l_stop = .true.
istop = i
jstop = j
endif
endif
enddo
endif
if (l_stop) return
!-----------------------------------------------------------------
! transfer volume and energy between categories
!-----------------------------------------------------------------
ishift = 0
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (daice(ij,n) > c0) then ! daice(n) can be < puny
ishift = ishift + 1
indxii(ishift) = i
indxjj(ishift) = j
indxij(ishift) = ij
endif ! tmask
enddo
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, ishift
i = indxii(ij)
j = indxjj(ij)
m = indxij(ij)
nd = donor(m,n)
worka(i,j) = dvice(m,n) / vicen(i,j,nd)
if (nd == n) then
nr = nd+1
else ! nd = n+1
nr = n
endif
aicen(i,j,nd) = aicen(i,j,nd) - daice(m,n)
aicen(i,j,nr) = aicen(i,j,nr) + daice(m,n)
vicen(i,j,nd) = vicen(i,j,nd) - dvice(m,n)
vicen(i,j,nr) = vicen(i,j,nr) + dvice(m,n)
dvsnow = vsnon(i,j,nd) * worka(i,j)
vsnon(i,j,nd) = vsnon(i,j,nd) - dvsnow
vsnon(i,j,nr) = vsnon(i,j,nr) + dvsnow
workb(i,j) = dvsnow
enddo ! ij
do it = 1, ntrcr
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, ishift
i = indxii(ij)
j = indxjj(ij)
m = indxij(ij)
nd = donor(m,n)
if (nd == n) then
nr = nd+1
else ! nd = n+1
nr = n
endif
if (trcr_depend(it) == 0) then
datrcr = daice(m,n)*trcrn(i,j,it,nd)
elseif (trcr_depend(it) == 1) then
datrcr = dvice(m,n)*trcrn(i,j,it,nd)
elseif (trcr_depend(it) == 2) then
datrcr = workb(i,j) *trcrn(i,j,it,nd)
endif
atrcrn(m,it,nd) = atrcrn(m,it,nd) - datrcr
atrcrn(m,it,nr) = atrcrn(m,it,nr) + datrcr
enddo ! ij
enddo ! ntrcr
do k = 1, nilyr
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, ishift
i = indxii(ij)
j = indxjj(ij)
m = indxij(ij)
nd = donor(m,n)
if (nd == n) then
nr = nd+1
else ! nd = n+1
nr = n
endif
deice = eicen(i,j,ilyr1(nd)+k-1) * worka(i,j)
eicen(i,j,ilyr1(nd)+k-1) = &
eicen(i,j,ilyr1(nd)+k-1) - deice
eicen(i,j,ilyr1(nr)+k-1) = &
eicen(i,j,ilyr1(nr)+k-1) + deice
enddo ! ij
enddo ! nilyr
do k = 1, nslyr
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, ishift
i = indxii(ij)
j = indxjj(ij)
m = indxij(ij)
nd = donor(m,n)
if (nd == n) then
nr = nd+1
else ! nd = n+1
nr = n
endif
desnow = esnon(i,j,slyr1(nd)+k-1) * worka(i,j)
esnon(i,j,slyr1(nd)+k-1) = &
esnon(i,j,slyr1(nd)+k-1) - desnow
esnon(i,j,slyr1(nr)+k-1) = &
esnon(i,j,slyr1(nr)+k-1) + desnow
enddo ! ij
enddo ! nslyr
enddo ! boundaries, 1 to ncat-1
!-----------------------------------------------------------------
! Update ice thickness and tracers
!-----------------------------------------------------------------
do n = 1, ncat
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (aicen(i,j,n) > puny) then
hicen(ij,n) = vicen (i,j,n) / aicen(i,j,n)
else
hicen(ij,n) = c0
endif
enddo
call compute_tracers (nx_block, ny_block, &
icells, indxi, indxj, &
ntrcr, trcr_depend, &
atrcrn(:,:,n), aicen(:,:, n), &
vicen (:,:, n), vsnon(:,:, n), &
trcrn(:,:,:,n))
enddo ! ncat
end subroutine shift_ice
!=======================================================================
!BOP
!
! !IROUTINE: column_sum - sum field over all ice categories
!
! !INTERFACE:
!
subroutine column_sum (nx_block, ny_block, & 18
icells, indxi, indxj, &
nsum, &
xin, xout)
!
! !DESCRIPTION:
!
! For each grid cell, sum field over all ice categories.
!
! !REVISION HISTORY:
!
! author: William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
nsum , & ! number of categories/layers
icells ! number of ice/ocean grid cells
integer (kind=int_kind), dimension (nx_block*ny_block), &
intent(in) :: &
indxi, indxj ! compressed i/j indices
real (kind=dbl_kind), dimension (nx_block,ny_block,nsum), &
intent(in) :: &
xin ! input field
real (kind=dbl_kind), dimension (icells), intent(out) :: &
xout ! output field
!
!EOP
!
integer (kind=int_kind) :: &
i, j, ij , & ! horizontal indices
n ! category/layer index
do ij = 1, icells
xout(ij) = c0
enddo
do n = 1, nsum
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
xout(ij) = xout(ij) + xin(i,j,n)
enddo ! ij
enddo ! n
end subroutine column_sum
!=======================================================================
!BOP
!
! !IROUTINE: column_conservation_check
!
! !INTERFACE:
!
subroutine column_conservation_check (nx_block, ny_block, & 9
icells, indxi, indxj, &
fieldid, &
x1, x2, &
max_err, l_stop, &
istop, jstop)
!
! !DESCRIPTION:
!
! For each physical grid cell, check that initial and final values
! of a conserved field are equal to within a small value.
!
! !REVISION HISTORY:
!
! author: William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
icells ! number of ice/ocean grid cells
integer (kind=int_kind), dimension (nx_block*ny_block), &
intent(in) :: &
indxi, indxj ! compressed i/j indices
real (kind=dbl_kind), dimension(icells), intent(in) :: &
x1 , & ! initial field
x2 ! final field
real (kind=dbl_kind), intent(in) :: &
max_err ! max allowed error
character (len=char_len), intent(in) :: &
fieldid ! field identifier
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) :: &
ij ! horizontal indices
do ij = 1, icells
if (abs (x2(ij)-x1(ij)) > max_err) then
l_stop = .true.
istop = indxi(ij)
jstop = indxj(ij)
write (nu_diag,*) ' '
write (nu_diag,*) 'Conservation error: ', trim(fieldid)
write (nu_diag,*) 'i, j =', istop, jstop
write (nu_diag,*) 'Initial value =', x1(ij)
write (nu_diag,*) 'Final value =', x2(ij)
write (nu_diag,*) 'Difference =', x2(ij) - x1(ij)
endif
enddo
end subroutine column_conservation_check
!=======================================================================
!BOP
!
! !IROUTINE: compute_tracers - compute tracer fields
!
! !INTERFACE:
!
subroutine compute_tracers (nx_block, ny_block, & 4,1
icells, indxi, indxj, &
ntrcr, trcr_depend, &
atrcrn, aicen, &
vicen, vsnon, &
trcrn)
!
! !DESCRIPTION:
!
! Compute tracer fields.
! Given atrcrn = aicen*trcrn (or vicen*trcrn, vsnon*trcrn), compute trcrn.
!
! !REVISION HISTORY:
!
! author: William H. Lipscomb, LANL
!
! !USES:
!
use ice_state, only: nt_Tsfc
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
icells , & ! number of ice/ocean grid cells
ntrcr ! number of tracers in use
integer (kind=int_kind), dimension (nx_block*ny_block), &
intent(in) :: &
indxi, indxj ! compressed i/j indices
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 (icells,max_ntrcr), &
intent(in) :: &
atrcrn ! aicen*trcrn or vicen*trcrn or vsnon*trcrn
real (kind=dbl_kind), dimension (nx_block,ny_block), &
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), &
intent(out) :: &
trcrn ! ice tracers
!
!EOP
!
integer (kind=int_kind) :: &
i, j, it, ij ! counting indices
trcrn(:,:,:) = c0
!-----------------------------------------------------------------
! Compute new tracers
!-----------------------------------------------------------------
do it = 1, ntrcr
if (it == nt_Tsfc) then ! surface temperature
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (aicen(i,j) > puny) then
trcrn(i,j,it) = atrcrn(ij,it) / aicen(i,j)
else
trcrn(i,j,it) = Tocnfrz
endif
enddo
elseif (trcr_depend(it) == 0) then ! ice area tracers
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (aicen(i,j) > puny) then
trcrn(i,j,it) = atrcrn(ij,it) / aicen(i,j)
else
trcrn(i,j,it) = c0
endif
enddo
elseif (trcr_depend(it) == 1) then ! ice volume tracers
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (vicen(i,j) > puny) then
trcrn(i,j,it) = atrcrn(ij,it) / vicen(i,j)
else
trcrn(i,j,it) = c0
endif
enddo
elseif (trcr_depend(it) == 2) then ! snow volume tracers
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (vsnon(i,j) > puny) then
trcrn(i,j,it) = atrcrn(ij,it) / vsnon(i,j)
else
trcrn(i,j,it) = c0
endif
enddo
endif ! trcr_depend
enddo ! ntrcr
end subroutine compute_tracers
!=======================================================================
!BOP
!
! !IROUTINE: cleanup_itd - rebin if needed, eliminate small ice areas,
! and aggregate over categories
!
! !INTERFACE:
!
subroutine cleanup_itd (nx_block, ny_block, & 2,4
ilo, ihi, jlo, jhi, &
dt, ntrcr, &
aicen, trcrn, &
vicen, vsnon, &
eicen, esnon, &
aice0, aice, &
trcr_depend, fresh, &
fsalt, fhocn, &
fsoot, tr_aero, &
heat_capacity, l_stop, &
istop, jstop, &
limit_aice_in)
!
! !DESCRIPTION:
!
! Cleanup subroutine that rebins thickness categories if necessary,
! eliminates very small ice areas while conserving mass and energy,
! aggregates state variables, and does a boundary call.
! It is a good idea to call this subroutine after the thermodynamics
! (thermo_vertical/thermo_itd) and again after the dynamics
! (evp/transport/ridging).
!
! !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
ntrcr ! number of tracers in use
real (kind=dbl_kind), intent(in) :: &
dt ! time step
real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), &
intent(inout) :: &
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(inout) :: &
trcrn ! ice 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)
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(inout) :: &
aice , & ! total ice concentration
aice0 ! concentration of open water
integer (kind=int_kind), dimension(max_ntrcr), intent(in) :: &
trcr_depend ! tracer dependency information
logical (kind=log_kind), intent(in) :: &
tr_aero, &
heat_capacity ! if false, ice and snow have zero heat capacity
logical (kind=log_kind), intent(out) :: &
l_stop ! if true, abort on return
integer (kind=int_kind), intent(out) :: &
istop, jstop ! indices of grid cell where model aborts
! ice-ocean fluxes (required for strict conservation)
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(inout), optional :: &
fresh , & ! fresh water flux to ocean (kg/m^2/s)
fsalt , & ! salt flux to ocean (kg/m^2/s)
fhocn ! net heat flux to ocean (W/m^2)
real (kind=dbl_kind), dimension (nx_block,ny_block,n_aeromx), &
intent(inout), optional :: &
fsoot ! soot flux to ocean (kg/m^2/s)
logical (kind=log_kind), intent(in), optional :: &
limit_aice_in ! if false, allow aice to be out of bounds
! may want to allow this for unit tests
!
!EOP
!
integer (kind=int_kind) :: &
i, j , & ! horizontal indices
n , & ! category index
icells ! number of grid cells with ice
integer (kind=int_kind), dimension (nx_block*ny_block) :: &
indxi, indxj ! compressed i/j indices
real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
dfresh , & ! zapped fresh water flux (kg/m^2/s)
dfsalt , & ! zapped salt flux (kg/m^2/s)
dfhocn ! zapped energy flux ( W/m^2)
real (kind=dbl_kind), dimension (nx_block,ny_block,n_aeromx) :: &
dfsoot ! zapped soot flux (kg/m^2/s)
logical (kind=log_kind) :: &
limit_aice ! if true, check for aice out of bounds
!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------
if (present(limit_aice_in)) then
limit_aice = limit_aice_in
else
limit_aice = .true.
endif
l_stop = .false.
istop = 0
jstop = 0
!-----------------------------------------------------------------
! Compute total ice area.
!-----------------------------------------------------------------
call aggregate_area (nx_block, ny_block, &
aicen(:,:,:), &
aice, aice0)
if (limit_aice) then ! check for aice out of bounds
do j = jlo,jhi
do i = ilo,ihi
if (aice(i,j) > c1+puny .or. aice(i,j) < -puny) then
l_stop = .true.
istop = i
jstop = j
endif
enddo
enddo
if (l_stop) then ! area out of bounds
i = istop
j = jstop
write(nu_diag,*) ' '
write(nu_diag,*) 'aggregate ice area out of bounds'
write(nu_diag,*) 'my_task, i, j, aice:', &
my_task, i, j, aice(i,j)
do n = 1, ncat
write(nu_diag,*) 'n, aicen:', n, aicen(i,j,n)
enddo
return
endif ! l_stop
endif ! limit_aice
!-----------------------------------------------------------------
! Identify grid cells with ice.
!-----------------------------------------------------------------
icells = 0
do j = jlo,jhi
do i = ilo,ihi
if (aice(i,j) > puny) then
icells = icells + 1
indxi(icells) = i
indxj(icells) = j
endif
enddo
enddo
!-----------------------------------------------------------------
! Make sure ice in each category is within its thickness bounds.
! NOTE: The rebin subroutine is needed only in the rare cases
! when the linear_itd subroutine cannot transfer ice
! correctly (e.g., very fast ice growth).
!-----------------------------------------------------------------
call rebin (nx_block, ny_block, &
icells, indxi, indxj, &
ntrcr, trcr_depend, &
aicen(:,:,:), trcrn(:,:,:,:), &
vicen(:,:,:), vsnon(:,:,:), &
eicen(:,:,:), esnon(:,:,:), &
l_stop, &
istop, jstop)
if (l_stop) return
!-----------------------------------------------------------------
! Zero out ice categories with very small areas.
!-----------------------------------------------------------------
if (limit_aice) then
call zap_small_areas (nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
dt, ntrcr, &
aice, aice0, &
aicen(:,:,:), trcrn(:,:,:,:), &
vicen(:,:,:), vsnon(:,:,:), &
eicen(:,:,:), esnon(:,:,:), &
dfresh, dfsalt, &
dfhocn, dfsoot, &
tr_aero, &
l_stop, &
istop, jstop)
if (l_stop) return
endif ! l_limit_aice
!-------------------------------------------------------------------
! Update ice-ocean fluxes for strict conservation
!-------------------------------------------------------------------
if (present(fresh)) &
fresh (:,:) = fresh(:,:) + dfresh(:,:)
if (present(fsalt)) &
fsalt (:,:) = fsalt(:,:) + dfsalt(:,:)
if (present(fhocn)) &
fhocn (:,:) = fhocn(:,:) + dfhocn(:,:)
if (present(fsoot)) &
fsoot (:,:,:) = fsoot(:,:,:) + dfsoot(:,:,:)
!----------------------------------------------------------------
! If using zero-layer model (no heat capacity), check that the
! energy of snow and ice is correct.
!----------------------------------------------------------------
if (.not. heat_capacity) then
call zerolayer_check(nx_block, ny_block, &
icells, indxi, indxj, &
aicen, &
vicen, vsnon, &
eicen, esnon, &
l_stop, &
istop, jstop)
endif
end subroutine cleanup_itd
!=======================================================================
!BOP
!
! !IROUTINE: zap_small_areas - eliminate very small ice areas
!
! !INTERFACE:
!
subroutine zap_small_areas (nx_block, ny_block, & 1,1
ilo, ihi, jlo, jhi, &
dt, ntrcr, &
aice, aice0, &
aicen, trcrn, &
vicen, vsnon, &
eicen, esnon, &
dfresh, dfsalt, &
dfhocn, dfsoot, &
tr_aero, &
l_stop, &
istop, jstop)
!
! !DESCRIPTION:
!
! For each ice category in each grid cell, remove ice if the fractional
! area is less than puny.
!
! !REVISION HISTORY:
!
! author: William H. Lipscomb, LANL
!
! !USES:
!
use ice_state, only: nt_Tsfc, nt_aero
!
! !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
ntrcr ! number of tracers in use
real (kind=dbl_kind), intent(in) :: &
dt ! time step
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(inout) :: &
aice , & ! total ice concentration
aice0 ! concentration of open water
real (kind=dbl_kind), dimension(nx_block,ny_block,ncat), &
intent(inout) :: &
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,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)
real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrcr,ncat), &
intent(inout) :: &
trcrn ! ice tracers
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(out) :: &
dfresh , & ! zapped fresh water flux (kg/m^2/s)
dfsalt , & ! zapped salt flux (kg/m^2/s)
dfhocn ! zapped energy flux ( W/m^2)
real (kind=dbl_kind), dimension (nx_block,ny_block,n_aeromx), &
intent(out) :: &
dfsoot ! zapped soot flux (kg/m^2/s)
logical (kind=log_kind), intent(in) :: &
tr_aero
logical (kind=log_kind), intent(out) :: &
l_stop ! if true, abort on return
integer (kind=int_kind), intent(out) :: &
istop, jstop ! indices of grid cell where model aborts
!
!EOP
!
integer (kind=int_kind) :: &
i,j, n, k, it , & ! counting indices
icells , & ! number of cells with ice to zap
ij ! combined i/j horizontal index
integer (kind=int_kind), dimension (nx_block*ny_block) :: &
indxi , & ! compressed indices for i/j directions
indxj
real (kind=dbl_kind) :: xtmp ! temporary variable
!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------
l_stop = .false.
istop = 0
jstop = 0
dfresh(:,:) = c0
dfsalt(:,:) = c0
dfhocn(:,:) = c0
dfsoot(:,:,:) = c0
!-----------------------------------------------------------------
! Zap categories with very small areas.
!-----------------------------------------------------------------
do n = 1, ncat
!-----------------------------------------------------------------
! Count categories to be zapped.
!-----------------------------------------------------------------
icells = 0
do j = jlo, jhi
do i = ilo, ihi
if (aicen(i,j,n) < -puny) then
write (nu_diag,*) 'Zap ice: negative ice area'
write (nu_diag,*) 'i, j, n, aicen =', &
i, j, n, aicen(i,j,n)
l_stop = .true.
istop = i
jstop = j
return
elseif ((aicen(i,j,n) >= -puny .and. aicen(i,j,n) < c0) .or. &
(aicen(i,j,n) > c0 .and. aicen(i,j,n) <= puny)) then
icells = icells + 1
indxi(icells) = i
indxj(icells) = j
endif
enddo
enddo
!-----------------------------------------------------------------
! Zap ice energy and use ocean heat to melt ice
!-----------------------------------------------------------------
do k = 1, nilyr
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
xtmp = eicen(i,j,ilyr1(n)+k-1) / dt ! < 0
dfhocn(i,j) = dfhocn(i,j) + xtmp
eicen(i,j,ilyr1(n)+k-1) = c0
enddo ! ij
enddo ! k
!-----------------------------------------------------------------
! Zap snow energy and use ocean heat to melt snow
!-----------------------------------------------------------------
do k = 1, nslyr
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
xtmp = esnon(i,j,slyr1(n)+k-1) / dt ! < 0
dfhocn(i,j) = dfhocn(i,j) + xtmp
esnon(i,j,slyr1(n)+k-1) = c0
enddo ! ij
enddo ! k
!-----------------------------------------------------------------
! Zap ice and snow volume, add water and salt to ocean
!-----------------------------------------------------------------
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt
dfresh(i,j) = dfresh(i,j) + xtmp
xtmp = rhoi*vicen(i,j,n)*ice_ref_salinity*p001 / dt
dfsalt(i,j) = dfsalt(i,j) + xtmp
aice0(i,j) = aice0(i,j) + aicen(i,j,n)
aicen(i,j,n) = c0
vicen(i,j,n) = c0
vsnon(i,j,n) = c0
trcrn(i,j,nt_Tsfc,n) = Tocnfrz
enddo ! ij
if (tr_aero) then
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
do it=1,n_aero
xtmp &
= (vsnon(i,j,n)*(trcrn(i,j,nt_aero +4*(it-1),n) &
+trcrn(i,j,nt_aero+1+4*(it-1),n)) &
+ vicen(i,j,n)*(trcrn(i,j,nt_aero+2+4*(it-1),n) &
+trcrn(i,j,nt_aero+3+4*(it-1),n))) &
/ dt
dfsoot(i,j,it) = dfsoot(i,j,it) + xtmp
enddo ! n
enddo ! ij
endif
!-----------------------------------------------------------------
! Zap tracers
!-----------------------------------------------------------------
if (ntrcr >= 2) then
do it = 1, ntrcr ! this assumes nt_Tsfc = 1
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
trcrn(i,j,it,n) = c0
enddo
enddo
endif
enddo ! n
!-----------------------------------------------------------------
! Count cells with excess ice (aice > c1) due to roundoff errors.
! Zap a little ice in each category so that aice = c1.
!-----------------------------------------------------------------
icells = 0
do j = jlo, jhi
do i = ilo, ihi
if (aice(i,j) > (c1+puny)) then
write (nu_diag,*) 'Zap ice: excess ice area'
write (nu_diag,*) 'i, j, aice =', &
i, j, aice(i,j)
l_stop = .true.
istop = i
jstop = j
return
elseif (aice(i,j) > c1 .and. aice(i,j) < (c1+puny)) then
icells = icells + 1
indxi(icells) = i
indxj(icells) = j
endif
enddo
enddo
do n = 1, ncat
!-----------------------------------------------------------------
! Zap ice energy and use ocean heat to melt ice
!-----------------------------------------------------------------
do k = 1, nilyr
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
xtmp = eicen(i,j,ilyr1(n)+k-1) &
* (aice(i,j)-c1)/aice(i,j) / dt ! < 0
dfhocn(i,j) = dfhocn(i,j) + xtmp
eicen(i,j,ilyr1(n)+k-1) = eicen(i,j,ilyr1(n)+k-1) &
* (c1/aice(i,j))
enddo ! ij
enddo ! k
!-----------------------------------------------------------------
! Zap snow energy and use ocean heat to melt snow
!-----------------------------------------------------------------
do k = 1, nslyr
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
xtmp = esnon(i,j,slyr1(n)+k-1) &
* (aice(i,j)-c1)/aice(i,j) / dt ! < 0
dfhocn(i,j) = dfhocn(i,j) + xtmp
esnon(i,j,slyr1(n)+k-1) = esnon(i,j,slyr1(n)+k-1) &
*(c1/aice(i,j))
enddo ! ij
enddo ! k
!-----------------------------------------------------------------
! Zap ice and snow volume, add water and salt to ocean
!-----------------------------------------------------------------
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) &
* (aice(i,j)-c1)/aice(i,j) / dt
dfresh(i,j) = dfresh(i,j) + xtmp
xtmp = rhoi*vicen(i,j,n)*ice_ref_salinity*p001 &
* (aice(i,j)-c1)/aice(i,j) / dt
dfsalt(i,j) = dfsalt(i,j) + xtmp
aicen(i,j,n) = aicen(i,j,n) * (c1/aice(i,j))
vicen(i,j,n) = vicen(i,j,n) * (c1/aice(i,j))
vsnon(i,j,n) = vsnon(i,j,n) * (c1/aice(i,j))
enddo ! ij
! Note: Tracers are unchanged.
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
if (tr_aero) then
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
do it=1,n_aero
xtmp &
= (vsnon(i,j,n)*(trcrn(i,j,nt_aero +4*(it-1),n) &
+trcrn(i,j,nt_aero+1+4*(it-1),n)) &
+ vicen(i,j,n)*(trcrn(i,j,nt_aero+2+4*(it-1),n) &
+trcrn(i,j,nt_aero+3+4*(it-1),n))) &
* (aice(i,j)-c1)/aice(i,j) / dt
dfsoot(i,j,it) = dfsoot(i,j,it) + xtmp
enddo ! n
enddo ! ij
endif
enddo ! n
!-----------------------------------------------------------------
! Correct aice
!-----------------------------------------------------------------
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
aice(i,j) = c1
aice0(i,j) = c0
enddo
end subroutine zap_small_areas
!=======================================================================
!BOP
!
! !IROUTINE: zerolayer_check - check that snow and ice energy is
! correct when using zero layer thermodynamics
!
! !INTERFACE:
!
subroutine zerolayer_check (nx_block, ny_block, & 1
icells, indxi, indxj, &
aicen, &
vicen, vsnon, &
eicen, esnon, &
l_stop, &
istop, jstop)
!
! !DESCRIPTION:
!
! Checks that the snow and ice energy in the zero layer thermodynamics
! model still agrees with the snow and ice volume.
! If there is an error, the model will abort.
! This subroutine is only called if heat_capacity = .false.
!
! !REVISION HISTORY:
!
! author: Alison McLaren, Met Office
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
icells ! number of grid cells with ice
integer (kind=int_kind), dimension (nx_block*ny_block), &
intent(in) :: &
indxi, indxj ! compressed i/j indices
real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), &
intent(inout) :: &
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,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)
logical (kind=log_kind), intent(out) :: &
l_stop ! if true, abort on return
integer (kind=int_kind), intent(out) :: &
istop, jstop ! indices of grid cell where model aborts
!
!EOP
!
integer (kind=int_kind) :: &
i, j , & ! horizontal indices
n , & ! category index
ij ! combined horizontal index
real (kind=dbl_kind), parameter :: &
max_error = puny*Lfresh*rhos ! max error in zero layer energy check
! (so max volume error = puny)
logical (kind=log_kind) :: &
ice_energy_correct , & ! zero layer ice energy check
snow_energy_correct ! zero layer snow energy check
real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
worka, &
workb
!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------
l_stop = .false.
istop = 0
jstop = 0
worka(:,:) = c0
workb(:,:) = c0
!----------------------------------------------------------------
! Calculate difference between ice and snow energies and the
! energy values derived from the ice and snow volumes
!----------------------------------------------------------------
ice_energy_correct = .true.
snow_energy_correct = .true.
do n=1,ncat
do ij=1,icells
i=indxi(ij)
j=indxj(ij)
worka(i,j) = eicen(i,j,n) + rhoi * Lfresh * vicen(i,j,n)
workb(i,j) = esnon(i,j,n) + rhos * Lfresh * vsnon(i,j,n)
if(abs(worka(i,j)) > max_error) then
ice_energy_correct = .false.
endif
if(abs(workb(i,j)) > max_error) then
snow_energy_correct = .false.
endif
enddo
!----------------------------------------------------------------
! If there is a problem, abort with error message
!----------------------------------------------------------------
if (.not. ice_energy_correct) then
do ij=1,icells
i=indxi(ij)
j=indxj(ij)
if(abs(worka(i,j)) > max_error) then
write(nu_diag,*) ' '
write(nu_diag,*) &
'zerolayer check - wrong ice energy'
write(nu_diag,*) 'i, j, n:', i,j,n
write(nu_diag,*) 'eicen =', eicen(i,j,n)
write(nu_diag,*) 'error=', worka(i,j)
write(nu_diag,*) 'vicen =', vicen(i,j,n)
write(nu_diag,*) 'aicen =', aicen(i,j,n)
l_stop = .true.
istop = i
jstop = j
endif
enddo
endif
if (l_stop) return
if (.not. snow_energy_correct) then
do ij=1,icells
i=indxi(ij)
j=indxj(ij)
if(abs(workb(i,j)) > max_error) then
write(nu_diag,*) ' '
write(nu_diag,*) &
'zerolayer_check - wrong snow energy'
write(nu_diag,*) 'i, j, n:', i,j,n
write(nu_diag,*) 'esnon =', esnon(i,j,n)
write(nu_diag,*) 'error=', workb(i,j)
write(nu_diag,*) 'vsnon =', vsnon(i,j,n)
write(nu_diag,*) 'aicen =', aicen(i,j,n)
l_stop = .true.
istop = i
jstop = j
return
endif
enddo
endif
enddo ! ncat
end subroutine zerolayer_check
!=======================================================================
end module ice_itd
!=======================================================================