!=======================================================================
!BOP
!
! !MODULE: ice_therm_itd - thermo calculations after call to coupler
!
! !DESCRIPTION:
!
! Thermo calculations after call to coupler, related to ITD:
! ice thickness redistribution, lateral growth and melting.
!
! NOTE: The thermodynamic calculation is split in two for load balancing.
! First ice_therm_vertical computes vertical growth rates and coupler
! fluxes. Then ice_therm_itd does thermodynamic calculations not
! needed for coupling.
!
! !REVISION HISTORY:
! SVN:$Id: ice_therm_itd.F90 48 2007-01-09 23:57:33Z eclare $
!
! authors William H. Lipscomb, LANL
! C. M. Bitz, UW
! Elizabeth C. Hunke, LANL
!
! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
! 2004: Block structure added by William Lipscomb.
! 2006: Streamlined for efficiency by Elizabeth Hunke
!
! !INTERFACE:
!
module ice_therm_itd 3,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
real (kind=dbl_kind), parameter, private :: &
hfrazilmin = 0.05_dbl_kind ! min thickness of new frazil ice (m)
!=======================================================================
contains
!=======================================================================
!BOP
!
! !IROUTINE: linear_itd - ITD scheme that shifts ice among categories
!
! !INTERFACE:
!
subroutine linear_itd (nx_block, ny_block, & 1,17
icells, indxi, indxj, &
ntrcr, trcr_depend, &
aicen_init, vicen_init, &
aicen, trcrn, &
vicen, vsnon, &
eicen, esnon, &
aice, aice0, &
l_stop, &
istop, jstop)
! See Lipscomb, W. H. Remapping the thickness distribution in sea
! ice models. 2001, J. Geophys. Res., Vol 106, 13989--14000.
!
! Using the thermodynamic "velocities", interpolate to find the
! velocities in thickness space at the category boundaries, and
! compute the new locations of the boundaries. Then for each
! category, compute the thickness distribution function, g(h),
! between hL and hR, the left and right boundaries of the category.
! Assume g(h) is a linear polynomial that satisfies two conditions:
!
! (1) The ice area implied by g(h) equals aicen(n).
! (2) The ice volume implied by g(h) equals aicen(n)*hicen(n).
!
! Given g(h), at each boundary compute the ice area and volume lying
! between the original and new boundary locations. Transfer area
! and volume across each boundary in the appropriate direction, thus
! restoring the original boundaries.
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb, LANL
! Elizabeth C. Hunke, LANL
!
! !USES:
!
use ice_itd, only: hin_max, hi_min, aggregate_area, shift_ice, &
column_sum, column_conservation_check
!
! !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(in) :: &
aicen_init, & ! initial ice concentration (before vertical thermo)
vicen_init ! initial ice volume (m)
real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), &
intent(inout) :: &
aicen , & ! ice concentration
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 , & ! concentration of ice
aice0 ! concentration of open water
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, nd , & ! category indices
k ! ice layer index
real (kind=dbl_kind) :: &
slope , & ! rate of change of dhice with hice
dh0 , & ! change in ice thickness at h = 0
da0 , & ! area melting from category 1
damax , & ! max allowed reduction in category 1 area
etamin, etamax,& ! left and right limits of integration
x1 , & ! etamax - etamin
x2 , & ! (etamax^2 - etamin^2) / 2
x3 , & ! (etamax^3 - etamin^3) / 3
wk1, wk2 ! temporary variables
real (kind=dbl_kind), dimension(icells,0:ncat) :: &
hbnew ! new boundary locations
real (kind=dbl_kind), dimension(icells) :: &
work ! temporary work array (for hbnew)
integer (kind=int_kind), dimension(icells) :: &
indxii, indxjj,& ! compressed i/j indices
indxij ! compressed 1D i/j indices
real (kind=dbl_kind), dimension(:,:), allocatable :: &
g0 , & ! constant coefficient in g(h)
g1 , & ! linear coefficient in g(h)
hL , & ! left end of range over which g(h) > 0
hR ! right end of range over which g(h) > 0
real (kind=dbl_kind), dimension(icells,ncat) :: &
hicen , & ! ice thickness for each cat (m)
hicen_init , & ! initial ice thickness for each cat (pre-thermo)
dhicen , & ! thickness change for remapping (m)
daice , & ! ice area transferred across boundary
dvice ! ice volume transferred across boundary
real (kind=dbl_kind), dimension(icells) :: &
vice_init, vice_final, & ! ice volume summed over categories
vsno_init, vsno_final, & ! snow volume summed over categories
eice_init, eice_final, & ! ice energy summed over categories
esno_init, esno_final ! snow energy summed over categories
! 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) :: &
donor ! donor category index
logical (kind=log_kind), dimension(icells) :: &
remap_flag ! remap ITD if remap_flag(ij) is true
character (len=char_len) :: &
fieldid ! field identifier
logical (kind=log_kind), parameter :: &
l_conservation_check = .true. ! if true, check conservation
! l_conservation_check = .false. ! if true, check conservation
! (useful for debugging)
integer (kind=int_kind) :: &
iflag , & ! number of grid cells with remap_flag = .true.
ij, m ! combined horizontal indices
!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------
l_stop = .false.
istop = 0
jstop = 0
hin_max(ncat) = 999.9_dbl_kind ! arbitrary big number
!-----------------------------------------------------------------
! Compute volume and energy sums that linear remapping should
! conserve.
!-----------------------------------------------------------------
if (l_conservation_check) then
call column_sum (nx_block, ny_block, &
icells, indxi, indxj, &
ncat, &
vicen, vice_init)
call column_sum (nx_block, ny_block, &
icells, indxi, indxj, &
ncat, &
vsnon, vsno_init)
call column_sum (nx_block, ny_block, &
icells, indxi, indxj, &
ntilyr, &
eicen, eice_init)
call column_sum (nx_block, ny_block, &
icells, indxi, indxj, &
ntslyr, &
esnon, esno_init)
endif
!-----------------------------------------------------------------
! Initialize remapping flag.
! Remapping is done wherever remap_flag = .true.
! In rare cases the category boundaries may shift too far for the
! remapping algorithm to work, and remap_flag is set to .false.
! In these cases the simpler 'rebin' subroutine will shift ice
! between categories if needed.
!-----------------------------------------------------------------
do ij = 1, icells
remap_flag(ij) = .true.
enddo
!-----------------------------------------------------------------
! Compute thickness change in each category.
!-----------------------------------------------------------------
do n = 1, ncat
do ij = 1, icells ! aice(i,j) > puny
i = indxi(ij)
j = indxj(ij)
if (aicen_init(i,j,n) > puny) then
hicen_init (ij,n) = vicen_init(i,j,n) / &
aicen_init(i,j,n)
else
hicen_init(ij,n) = c0
endif ! aicen_init > puny
if (aicen(i,j,n) > puny) then
hicen (ij,n) = vicen(i,j,n) / aicen(i,j,n)
dhicen(ij,n) = hicen(ij,n) - hicen_init(ij,n)
else
hicen (ij,n) = c0
dhicen(ij,n) = c0
endif ! aicen > puny
enddo ! ij
enddo ! n
!-----------------------------------------------------------------
! Compute new category boundaries, hbnew, based on changes in
! ice thickness from vertical thermodynamics.
!-----------------------------------------------------------------
do ij = 1, icells ! aice(i,j) > puny
hbnew(ij,0) = hin_max(0)
enddo
do n = 1, ncat-1
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells ! aice(i,j) > puny
i = indxi(ij)
j = indxj(ij)
if (hicen_init(ij,n) > puny .and. &
hicen_init(ij,n+1) > puny) then
! interpolate between adjacent category growth rates
slope = (dhicen(ij,n+1) - dhicen(ij,n)) / &
(hicen_init(ij,n+1) - hicen_init(ij,n))
hbnew(ij,n) = hin_max(n) + dhicen(ij,n) &
+ slope * (hin_max(n) - hicen_init(ij,n))
elseif (hicen_init(ij,n) > puny) then ! hicen_init(n+1)=0
hbnew(ij,n) = hin_max(n) + dhicen(ij,n)
elseif (hicen_init(ij,n+1) > puny) then ! hicen_init(n)=0
hbnew(ij,n) = hin_max(n) + dhicen(ij,n+1)
else
hbnew(ij,n) = hin_max(n)
endif
!-----------------------------------------------------------------
! Check that each boundary lies between adjacent values of hicen.
! If not, set remap_flag = .false.
! Write diagnosis outputs if remap_flag was changed to false
!-----------------------------------------------------------------
if (aicen(i,j,n) > puny .and. &
hicen(ij,n) >= hbnew(ij,n)) then
remap_flag(ij) = .false.
write(nu_diag,*) my_task,':',i,j, &
'ITD: hicen(n) > hbnew(n)'
write(nu_diag,*) 'cat ',n
write(nu_diag,*) my_task,':',i,j, &
'hicen(n) =', hicen(ij,n)
write(nu_diag,*) my_task,':',i,j, &
'hbnew(n) =', hbnew(ij,n)
elseif (aicen(i,j,n+1) > puny .and. &
hicen(ij,n+1) <= hbnew(ij,n)) then
remap_flag(ij) = .false.
write(nu_diag,*) my_task,':',i,j, &
'ITD: hicen(n+1) < hbnew(n)'
write(nu_diag,*) 'cat ',n
write(nu_diag,*) my_task,':',i,j, &
'hicen(n+1) =', hicen(ij,n+1)
write(nu_diag,*) my_task,':',i,j, &
'hbnew(n) =', hbnew(ij,n)
endif
!-----------------------------------------------------------------
! Check that hbnew(n) lies between hin_max(n-1) and hin_max(n+1).
! If not, set remap_flag = .false.
! (In principle we could allow this, but it would make the code
! more complicated.)
! Write diagnosis outputs if remap_flag was changed to false
!-----------------------------------------------------------------
if (hbnew(ij,n) > hin_max(n+1)) then
remap_flag(ij) = .false.
write(nu_diag,*) my_task,':',i,j, &
'ITD hbnew(n) > hin_max(n+1)'
write(nu_diag,*) 'cat ',n
write(nu_diag,*) my_task,':',i,j, &
'hbnew(n) =', hbnew(ij,n)
write(nu_diag,*) my_task,':',i,j, &
'hin_max(n+1) =', hin_max(n+1)
endif
if (hbnew(ij,n) < hin_max(n-1)) then
remap_flag(ij) = .false.
write(nu_diag,*) my_task,':',i,j, &
'ITD: hbnew(n) < hin_max(n-1)'
write(nu_diag,*) 'cat ',n
write(nu_diag,*) my_task,':',i,j, &
'hbnew(n) =', hbnew(ij,n)
write(nu_diag,*) my_task,':',i,j, &
'hin_max(n-1) =', hin_max(n-1)
endif
enddo ! ij
enddo ! boundaries, 1 to ncat-1
!-----------------------------------------------------------------
! Compute hbnew(ncat)
!-----------------------------------------------------------------
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells ! aice(i,j) > puny
i = indxi(ij)
j = indxj(ij)
if (aicen(i,j,ncat) > puny) then
hbnew(ij,ncat) = c3*hicen(ij,ncat) - c2*hbnew(ij,ncat-1)
else
hbnew(ij,ncat) = hin_max(ncat)
endif
hbnew(ij,ncat) = max(hbnew(ij,ncat),hin_max(ncat-1))
enddo
!-----------------------------------------------------------------
! Identify cells where the ITD is to be remapped
!-----------------------------------------------------------------
iflag = 0
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
if (remap_flag(ij)) then
iflag = iflag + 1
indxii(iflag) = i
indxjj(iflag) = j
indxij(iflag) = ij
endif
enddo
allocate(g0(iflag,ncat))
allocate(g1(iflag,ncat))
allocate(hL(iflag,ncat))
allocate(hR(iflag,ncat))
!-----------------------------------------------------------------
! Compute g(h) for category 1 at start of time step
! (hicen = hicen_init)
!-----------------------------------------------------------------
do ij = 1, icells ! aice(i,j) > puny
work(ij) = hin_max(1)
enddo
call fit_line(nx_block, ny_block, &
iflag, icells, &
indxii, indxjj, indxij, &
aicen(:,:,1), hicen_init(:,1), &
hbnew(:,0), work (:), &
g0 (:,1), g1 (:,1), &
hL (:,1), hR (:,1))
!-----------------------------------------------------------------
! Find area lost due to melting of thin (category 1) ice
!-----------------------------------------------------------------
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, iflag ! remap_flag = .true.
i = indxii(ij)
j = indxjj(ij)
m = indxij(ij)
if (aicen(i,j,1) > puny) then
dh0 = dhicen(m,1)
if (dh0 < c0) then ! remove area from category 1
dh0 = min(-dh0,hin_max(1)) ! dh0 --> |dh0|
!-----------------------------------------------------------------
! Integrate g(1) from 0 to dh0 to estimate area melted
!-----------------------------------------------------------------
! right integration limit (left limit = 0)
etamax = min(dh0,hR(ij,1)) - hL(ij,1)
if (etamax > c0) then
x1 = etamax
x2 = p5 * etamax*etamax
da0 = g1(ij,1)*x2 + g0(ij,1)*x1 ! ice area removed
! constrain new thickness <= hicen_init
damax = aicen(i,j,1) &
* (c1-hicen(m,1)/hicen_init(m,1)) ! damax > 0
da0 = min (da0, damax)
! remove area, conserving volume
hicen(m,1) = hicen(m,1) &
* aicen(i,j,1) / (aicen(i,j,1)-da0)
aicen(i,j,1) = aicen(i,j,1) - da0
endif ! etamax > 0
else ! dh0 >= 0
hbnew(m,0) = min(dh0,hin_max(1)) ! shift hbnew(0) to right
endif
endif ! aicen(i,j,1) > puny
enddo ! ij
!-----------------------------------------------------------------
! Compute g(h) for each ice thickness category.
!-----------------------------------------------------------------
do n = 1, ncat
call fit_line(nx_block, ny_block, &
iflag, icells, &
indxii, indxjj, indxij, &
aicen(:,:,n), hicen(:,n), &
hbnew(:,n-1), hbnew(:,n), &
g0 (:,n), g1 (:,n), &
hL (:,n), hR (:,n))
enddo
!-----------------------------------------------------------------
! Compute area and volume to be shifted across each boundary.
!-----------------------------------------------------------------
do n = 1, ncat
do ij = 1, icells
donor(ij,n) = 0
daice(ij,n) = c0
dvice(ij,n) = c0
enddo
enddo
do n = 1, ncat-1
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, iflag ! remap_flag = .true.
i = indxii(ij)
j = indxjj(ij)
m = indxij(ij)
if (hbnew(m,n) > hin_max(n)) then ! transfer from n to n+1
! left and right integration limits in eta space
etamin = max(hin_max(n), hL(ij,n)) - hL(ij,n)
etamax = min(hbnew(m,n), hR(ij,n)) - hL(ij,n)
donor(m,n) = n
else ! hbnew(n) <= hin_max(n); transfer from n+1 to n
! left and right integration limits in eta space
etamin = c0
etamax = min(hin_max(n), hR(ij,n+1)) - hL(ij,n+1)
donor(m,n) = n+1
endif ! hbnew(n) > hin_max(n)
if (etamax > etamin) then
x1 = etamax - etamin
wk1 = etamin*etamin
wk2 = etamax*etamax
x2 = p5 * (wk2 - wk1)
wk1 = wk1*etamin
wk2 = wk2*etamax
x3 = p333 * (wk2 - wk1)
nd = donor(m,n)
daice(m,n) = g1(ij,nd)*x2 + g0(ij,nd)*x1
dvice(m,n) = g1(ij,nd)*x3 + g0(ij,nd)*x2 &
+ daice(m,n)*hL(ij,nd)
endif ! etamax > etamin
! If daice or dvice is very small, shift no ice.
nd = donor(m,n)
if (daice(m,n) < aicen(i,j,nd)*puny) then
daice(m,n) = c0
dvice(m,n) = c0
donor(m,n) = 0
endif
if (dvice(m,n) < vicen(i,j,nd)*puny) then
daice(m,n) = c0
dvice(m,n) = c0
donor(m,n) = 0
endif
! If daice is close to aicen or dvice is close to vicen,
! shift entire category
if (daice(m,n) > aicen(i,j,nd)*(c1-puny)) then
daice(m,n) = aicen(i,j,nd)
dvice(m,n) = vicen(i,j,nd)
endif
if (dvice(m,n) > vicen(i,j,nd)*(c1-puny)) then
daice(m,n) = aicen(i,j,nd)
dvice(m,n) = vicen(i,j,nd)
endif
enddo ! ij
enddo ! boundaries, 1 to ncat-1
deallocate(g0)
deallocate(g1)
deallocate(hL)
deallocate(hR)
!-----------------------------------------------------------------
! Shift ice between categories as necessary
!-----------------------------------------------------------------
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
!-----------------------------------------------------------------
! Make sure hice(i,j,1) >= minimum ice thickness hi_min.
!-----------------------------------------------------------------
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, iflag ! remap_flag = .true.
i = indxii(ij)
j = indxjj(ij)
m = indxij(ij)
if (hi_min > c0 .and. &
aicen(i,j,1) > puny .and. hicen(m,1) < hi_min) then
aicen(i,j,1) = aicen(i,j,1) * hicen(m,1)/hi_min
hicen(m,1) = hi_min
endif
enddo ! ij
!-----------------------------------------------------------------
! Update fractional ice area in each grid cell.
!-----------------------------------------------------------------
call aggregate_area (nx_block, ny_block, &
aicen, &
aice, aice0)
if (l_stop) return
!-----------------------------------------------------------------
! Check volume and energy conservation.
!-----------------------------------------------------------------
if (l_conservation_check) then
call column_sum (nx_block, ny_block, &
icells, indxi, indxj, &
ncat, &
vicen, vice_final)
fieldid = 'vice, ITD remap'
call column_conservation_check (nx_block, ny_block, &
icells, indxi, indxj, &
fieldid, &
vice_init, vice_final, &
puny, l_stop, &
istop, jstop)
if (l_stop) return
call column_sum (nx_block, ny_block, &
icells, indxi, indxj, &
ncat, &
vsnon, vsno_final)
fieldid = 'vsno, ITD remap'
call column_conservation_check (nx_block, ny_block, &
icells, indxi, indxj, &
fieldid, &
vsno_init, vsno_final, &
puny, l_stop, &
istop, jstop)
if (l_stop) return
call column_sum (nx_block, ny_block, &
icells, indxi, indxj, &
ntilyr, &
eicen, eice_final)
fieldid = 'eice, ITD remap'
call column_conservation_check (nx_block, ny_block, &
icells, indxi, indxj, &
fieldid, &
eice_init, eice_final, &
puny*Lfresh*rhoi, &
l_stop, &
istop, jstop)
if (l_stop) return
call column_sum (nx_block, ny_block, &
icells, indxi, indxj, &
ntslyr, &
esnon, esno_final)
fieldid = 'esno, ITD remap'
call column_conservation_check (nx_block, ny_block, &
icells, indxi, indxj, &
fieldid, &
esno_init, esno_final, &
puny*Lfresh*rhos, &
l_stop, &
istop, jstop)
if (l_stop) return
endif ! conservation check
end subroutine linear_itd
!=======================================================================
!BOP
!
! !IROUTINE: fit_line - fit g(h) with a line using area, volume constraints
!
! !INTERFACE:
!
subroutine fit_line (nx_block, ny_block, & 2
iflag, icells, &
indxii, indxjj, indxij, &
aicen, hice, &
hbL, hbR, &
g0, g1, &
hL, hR)
!
! !DESCRIPTION:
!
! Fit g(h) with a line, satisfying area and volume constraints.
! To reduce roundoff errors caused by large values of g0 and g1,
! we actually compute g(eta), where eta = h - hL, and hL is the
! left boundary.
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb, LANL
! 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
iflag ! number of grid cells with remap_flag = .true.
integer (kind=int_kind), dimension (icells), &
intent(in) :: &
indxii, indxjj, & ! compressed i/j indices (from iflag)
indxij ! compressed i/j indices (from icells)
real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
aicen ! concentration of ice
real (kind=dbl_kind), dimension (icells), intent(in) :: &
hbL, hbR , & ! left and right category boundaries
hice ! ice thickness
real (kind=dbl_kind), dimension (iflag), intent(out):: &
g0, g1 , & ! coefficients in linear equation for g(eta)
hL , & ! min value of range over which g(h) > 0
hR ! max value of range over which g(h) > 0
!
!EOP
!
integer (kind=int_kind) :: &
i,j , & ! horizontal indices
ij, m ! combined horizontal indices
real (kind=dbl_kind) :: &
h13 , & ! hbL + 1/3 * (hbR - hbL)
h23 , & ! hbL + 2/3 * (hbR - hbL)
dhr , & ! 1 / (hR - hL)
wk1, wk2 ! temporary variables
!-----------------------------------------------------------------
! Compute g0, g1, hL, and hR for each category to be remapped.
!-----------------------------------------------------------------
do ij = 1, iflag
i = indxii(ij)
j = indxjj(ij)
m = indxij(ij)
if (aicen(i,j) > puny .and. hbR(m) - hbL(m) > puny) then
! Initialize hL and hR
hL(ij) = hbL(m)
hR(ij) = hbR(m)
! Change hL or hR if hicen(n) falls outside central third of range
h13 = p333 * (c2*hL(ij) + hR(ij))
h23 = p333 * (hL(ij) + c2*hR(ij))
if (hice(m) < h13) then
hR(ij) = c3*hice(m) - c2*hL(ij)
elseif (hice(m) > h23) then
hL(ij) = c3*hice(m) - c2*hR(ij)
endif
! Compute coefficients of g(eta) = g0 + g1*eta
dhr = c1 / (hR(ij) - hL(ij))
wk1 = c6 * aicen(i,j) * dhr
wk2 = (hice(m) - hL(ij)) * dhr
g0(ij) = wk1 * (p666 - wk2)
g1(ij) = c2*dhr * wk1 * (wk2 - p5)
else
g0(ij) = c0
g1(ij) = c0
hL(ij) = c0
hR(ij) = c0
endif ! aicen > puny
enddo ! ij
end subroutine fit_line
!=======================================================================
!BOP
!
! !ROUTINE: add_new_ice - add frazil ice to ice thickness distribution
!
! !DESCRIPTION:
!
! Given the volume of new ice grown in open water, compute its area
! and thickness and add it to the appropriate category or categories.
!
! NOTE: Usually all the new ice is added to category 1. An exception is
! made if there is no open water or if the new ice is too thick
! for category 1, in which case ice is distributed evenly over the
! entire cell. Subroutine rebin should be called in case the ice
! thickness lies outside category bounds after new ice formation.
!
! !REVISION HISTORY:
!
! authors William H. Lipscomb, LANL
! Elizabeth C. Hunke, LANL
!
! !INTERFACE:
!
subroutine add_new_ice (nx_block, ny_block, & 1,6
ntrcr, icells, &
indxi, indxj, &
tmask, dt, &
aicen, trcrn, &
vicen, eicen, &
aice0, aice, &
frzmlt, frazil, &
frz_onset, yday, &
fresh, fsalt, &
Tf, l_stop, &
istop, jstop)
!
! !USES:
!
use ice_itd, only: hin_max, ilyr1, column_sum, &
column_conservation_check
use ice_state, only: nt_Tsfc, nt_iage, nt_FY, nt_aero, &
tr_iage, tr_FY, tr_aero, &
nt_alvl, nt_vlvl, tr_lvl
use ice_flux, only: update_ocn_f
! !INPUT/OUTPUT PARAMETERS:
!
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
ntrcr , & ! number of tracers in use
icells ! number of ice/ocean grid cells
integer (kind=int_kind), dimension (nx_block*ny_block), &
intent(in) :: &
indxi, indxj ! compressed i/j indices
logical (kind=log_kind), dimension (nx_block,ny_block), &
intent(in) :: &
tmask ! land/boundary mask, thickness (T-cell)
real (kind=dbl_kind), intent(in) :: &
dt ! time step (s)
real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
aice , & ! total concentration of ice
frzmlt, & ! freezing/melting potential (W/m^2)
Tf ! freezing temperature (C)
real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), &
intent(inout) :: &
aicen , & ! concentration of ice
vicen ! volume per unit area of ice (m)
real (kind=dbl_kind), dimension (nx_block,ny_block,max_ntrcr,ncat), &
intent(inout) :: &
trcrn ! ice tracers
! 1: surface temperature
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), &
intent(inout) :: &
aice0 , & ! concentration of open water
frazil , & ! frazil ice growth (m/step-->cm/day)
fresh , & ! fresh water flux to ocean (kg/m^2/s)
fsalt ! salt flux to ocean (kg/m^2/s)
real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(inout), optional :: &
frz_onset ! day of year that freezing begins (congel or frazil)
real (kind=dbl_kind), intent(in), optional :: &
yday ! day of year
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 , & ! ice category index
k , & ! ice layer index
it ! aerosol tracer index
real (kind=dbl_kind), dimension (icells) :: &
ai0new , & ! area of new ice added to cat 1
vi0new , & ! volume of new ice added to cat 1
hsurp , & ! thickness of new ice added to each cat
vlyr ! ice layer volume
real (kind=dbl_kind), dimension (icells) :: &
vice_init, vice_final ! ice volume summed over categories
real (kind=dbl_kind) :: &
fnew , & ! heat flx to open water for new ice (W/m^2)
hi0new , & ! thickness of new ice
hi0max , & ! max allowed thickness of new ice
qi0(nilyr) , & ! frazil ice enthalpy
qi0av , & ! mean value of qi0 for new ice (J kg-1)
vsurp , & ! volume of new ice added to each cat
area1 , & ! starting fractional area of existing ice
vice1 , & ! starting volume of existing ice
rnilyr , & ! real(nilyr)
dfresh , & ! change in fresh
dfsalt , & ! change in fsalt
vtmp
integer (kind=int_kind) :: &
jcells, kcells , & ! grid cell counters
ij, m ! combined i/j horizontal indices
integer (kind=int_kind), dimension (icells) :: &
indxij2, indxij3 , & ! compressed i/j indices
indxi2, indxj2 , &
indxi3, indxj3
character (len=char_len) :: &
fieldid ! field identifier
l_stop = .false.
istop = 0
jstop = 0
jcells = 0
kcells = 0
if (ncat > 1) then
hi0max = hin_max(1)*0.9_dbl_kind ! not too close to boundary
else
hi0max = bignum ! big number
endif
! initial ice volume in each grid cell
call column_sum (nx_block, ny_block, &
icells, indxi, indxj, &
ncat, &
vicen, vice_init)
!-----------------------------------------------------------------
! Compute average enthalpy of new ice.
!
! POP assumes new ice is fresh. Otherwise, it would be better
! to do something like this:
! qi0(i,j,k) = -rhoi * (cp_ice*(Tmlt(k)-Tf(i,j))
! + Lfresh*(1.-Tmlt(k)/Tf(i,j)) - cp_ocn*Tmlt(k))
!-----------------------------------------------------------------
rnilyr = real(nilyr,kind=dbl_kind)
qi0av = c0
do k = 1, nilyr
qi0(k) = -rhoi*Lfresh ! note sign convention, qi < 0
qi0av = qi0av + qi0(k)
enddo
qi0av = qi0av/rnilyr
!-----------------------------------------------------------------
! Compute the volume, area, and thickness of new ice.
!-----------------------------------------------------------------
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
fnew = max (frzmlt(i,j), c0) ! fnew > 0 iff frzmlt > 0
vi0new(ij) = -fnew*dt / qi0av ! note sign convention, qi < 0
! increment ice volume
vice_init(ij) = vice_init(ij) + vi0new(ij)
! history diagnostics
frazil(i,j) = vi0new(ij)
if (present(frz_onset) .and. present(yday)) then
if (frazil(i,j) > puny .and. frz_onset(i,j) < puny) &
frz_onset(i,j) = yday
endif
!-----------------------------------------------------------------
! Update fresh water and salt fluxes.
!
! NOTE: POP assumes fresh water and salt flux due to frzmlt > 0
! is NOT included in fluxes fresh and fsalt.
!-----------------------------------------------------------------
if (update_ocn_f) then
dfresh = -rhoi*vi0new(ij)/dt
dfsalt = ice_ref_salinity*p001*dfresh
fresh(i,j) = fresh(i,j) + dfresh
fsalt(i,j) = fsalt(i,j) + dfsalt
endif
!-----------------------------------------------------------------
! Decide how to distribute the new ice.
!-----------------------------------------------------------------
hsurp(ij) = c0
ai0new(ij) = c0
if (vi0new(ij) > c0) then
! new ice area and thickness
! hin_max(0) < new ice thickness < hin_max(1)
if (aice0(i,j) > puny) then
hi0new = max(vi0new(ij)/aice0(i,j), hfrazilmin)
if (hi0new > hi0max .and. aice0(i,j)+puny < c1) then
! distribute excess volume over all categories (below)
hi0new = hi0max
ai0new(ij) = aice0(i,j)
vsurp = vi0new(ij) - ai0new(ij)*hi0new
hsurp(ij) = vsurp / aice(i,j)
vi0new(ij) = ai0new(ij)*hi0new
else
! put ice in a single category, with hsurp = 0
ai0new(ij) = vi0new(ij)/hi0new
endif
else ! aice0 < puny
hsurp(ij) = vi0new(ij)/aice(i,j) ! new thickness in each cat
vi0new(ij) = c0
endif ! aice0 > puny
endif ! vi0new > puny
!-----------------------------------------------------------------
! Identify grid cells receiving new ice.
!-----------------------------------------------------------------
i = indxi(ij)
j = indxj(ij)
if (vi0new(ij) > c0) then ! add ice to category 1
jcells = jcells + 1
indxi2(jcells) = i
indxj2(jcells) = j
indxij2(jcells) = ij
endif
if (hsurp(ij) > c0) then ! add ice to all categories
kcells = kcells + 1
indxi3(kcells) = i
indxj3(kcells) = j
indxij3(kcells) = ij
endif
enddo ! ij
!-----------------------------------------------------------------
! Distribute excess ice volume among ice categories by increasing
! ice thickness, leaving ice area unchanged.
!
! NOTE: If new ice contains globally conserved tracers
! (e.g., isotopes from seawater), code must be added here.
!-----------------------------------------------------------------
do n = 1, ncat
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, kcells
i = indxi3(ij)
j = indxj3(ij)
m = indxij3(ij)
vsurp = hsurp(m) * aicen(i,j,n)
! update ice age due to freezing (new ice age = dt)
vtmp = vicen(i,j,n) + vsurp
if (vtmp > puny) then
if (tr_iage) trcrn(i,j,nt_iage,n) &
= (trcrn(i,j,nt_iage,n)*vicen(i,j,n) + dt*vsurp) &
/ vtmp
if (tr_lvl) &
trcrn(i,j,nt_vlvl,n) = &
(trcrn(i,j,nt_vlvl,n)*vicen(i,j,n) + &
trcrn(i,j,nt_alvl,n)*vsurp) / vtmp
if (tr_aero) then
do it=1,n_aero
trcrn(i,j,nt_aero+2+4*(it-1),n) &
= trcrn(i,j,nt_aero+2+4*(it-1),n)*vicen(i,j,n) &
/ vtmp
trcrn(i,j,nt_aero+3+4*(it-1),n) &
= trcrn(i,j,nt_aero+3+4*(it-1),n)*vicen(i,j,n) &
/ vtmp
enddo
endif
endif
! update category volumes
vicen(i,j,n) = vicen(i,j,n) + vsurp
vlyr(m) = vsurp/rnilyr
enddo ! ij
do k = 1, nilyr
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, kcells
i = indxi3(ij)
j = indxj3(ij)
m = indxij3(ij)
eicen(i,j,ilyr1(n)+k-1) = &
eicen(i,j,ilyr1(n)+k-1) + qi0(k)*vlyr(m)
enddo ! ij
enddo ! k
enddo ! n
!-----------------------------------------------------------------
! Combine new ice grown in open water with category 1 ice.
! Assume that vsnon and esnon are unchanged.
!-----------------------------------------------------------------
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, jcells
i = indxi2(ij)
j = indxj2(ij)
m = indxij2(ij)
area1 = aicen(i,j,1) ! save
vice1 = vicen(i,j,1) ! save
aicen(i,j,1) = aicen(i,j,1) + ai0new(m)
aice0(i,j) = aice0(i,j) - ai0new(m)
vicen(i,j,1) = vicen(i,j,1) + vi0new(m)
trcrn(i,j,nt_Tsfc,1) = (Tf(i,j)*ai0new(m) + trcrn(i,j,nt_Tsfc,1)*area1) &
/ aicen(i,j,1)
trcrn(i,j,nt_Tsfc,1) = min (trcrn(i,j,nt_Tsfc,1), c0)
if (vicen(i,j,1) > puny) then
if (tr_iage) trcrn(i,j,nt_iage,1) = &
(trcrn(i,j,nt_iage,1)*vice1 + dt*vi0new(m))/vicen(i,j,1)
if (tr_aero) then
do it=1,n_aero
trcrn(i,j,nt_aero+2+4*(it-1),1) = &
trcrn(i,j,nt_aero+2+4*(it-1),1)*vice1/vicen(i,j,1)
trcrn(i,j,nt_aero+3+4*(it-1),1) = &
trcrn(i,j,nt_aero+3+4*(it-1),1)*vice1/vicen(i,j,1)
enddo
endif
endif
if (tr_lvl .and. aicen(i,j,1) > puny) then
trcrn(i,j,nt_alvl,1) = &
(trcrn(i,j,nt_alvl,1)*area1 + ai0new(m))/aicen(i,j,1)
trcrn(i,j,nt_vlvl,1) = &
(trcrn(i,j,nt_vlvl,1)*vice1 + vi0new(m))/vicen(i,j,1)
endif
if (tr_FY) &
trcrn(i,j,nt_FY,1) = (ai0new(m) + trcrn(i,j,nt_FY,1)*area1) &
/ aicen(i,j,1)
vlyr(m) = vi0new(m) / rnilyr
enddo ! ij
do k = 1, nilyr
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, jcells
i = indxi2(ij)
j = indxj2(ij)
m = indxij2(ij)
eicen(i,j,k) = eicen(i,j,k) + qi0(k)*vlyr(m)
enddo
enddo
call column_sum (nx_block, ny_block, &
icells, indxi, indxj, &
ncat, &
vicen, vice_final)
fieldid = 'vice, add_new_ice'
call column_conservation_check (nx_block, ny_block, &
icells, indxi, indxj, &
fieldid, &
vice_init, vice_final, &
puny, l_stop, &
istop, jstop)
if (l_stop) return
end subroutine add_new_ice
!=======================================================================
!BOP
!
! !ROUTINE: lateral_melt - melt ice laterally
!
! !DESCRIPTION:
!
! Given the fraction of ice melting laterally in each grid cell
! (computed in subroutine frzmlt_bottom_lateral), melt ice.
!
! !REVISION HISTORY:
!
! author: C. M. Bitz, UW
! 2003: Modified by William H. Lipscomb and Elizabeth C. Hunke, LANL
!
! !INTERFACE:
!
subroutine lateral_melt (nx_block, ny_block, & 1,2
ilo, ihi, jlo, jhi, &
dt, &
fresh, fsalt, &
fhocn, fsoot, &
rside, meltl, &
aicen, vicen, &
vsnon, eicen, &
esnon, trcrn)
!
! !USES:
!
use ice_itd, only: ilyr1, slyr1
use ice_state, only: nt_aero, tr_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
real (kind=dbl_kind), intent(in) :: &
dt ! time step (s)
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(in) :: &
trcrn ! energy of melting for each snow layer (J/m^2)
real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: &
rside ! fraction of ice that melts laterally
real (kind=dbl_kind), dimension(nx_block,ny_block), &
intent(inout) :: &
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)
meltl ! lateral ice melt (m/step-->cm/day)
real (kind=dbl_kind), dimension(nx_block,ny_block,n_aeromx), &
intent(inout) :: &
fsoot !
!
!EOP
!
integer (kind=int_kind) :: &
i, j , & ! horizontal indices
n , & ! thickness category index
k , & ! layer index
ij , & ! horizontal index, combines i and j loops
icells , & ! number of cells with aice > puny
it ! tracer index for aerosols
integer (kind=int_kind), dimension(nx_block*ny_block) :: &
indxi, indxj ! compressed indices for cells with aice > puny
real (kind=dbl_kind) :: &
dfhocn , & ! change in fhocn
dfresh , & ! change in fresh
dfsalt ! change in fsalt
do n = 1, ncat
!-----------------------------------------------------------------
! Identify grid cells with lateral melting.
!-----------------------------------------------------------------
icells = 0
do j = jlo, jhi
do i = ilo, ihi
if (rside(i,j) > c0) then
icells = icells + 1
indxi(icells) = i
indxj(icells) = j
endif
enddo ! i
enddo ! j
!-----------------------------------------------------------------
! Melt the ice and increment fluxes.
!-----------------------------------------------------------------
if (tr_aero) then
do k=1,n_aero
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
fsoot(i,j,k) = fsoot(i,j,k) &
+ (vsnon(i,j,n) &
*(trcrn(i,j,nt_aero +4*(k-1),n) &
+trcrn(i,j,nt_aero+1+4*(k-1),n)) &
+ vicen(i,j,n) &
*(trcrn(i,j,nt_aero+2+4*(k-1),n) &
+trcrn(i,j,nt_aero+3+4*(k-1),n))) &
* rside(i,j) / dt
enddo
enddo
endif
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
! fluxes to coupler
! dfresh > 0, dfsalt > 0
dfresh = (rhos*vsnon(i,j,n) + rhoi*vicen(i,j,n)) &
* rside(i,j) / dt
dfsalt = rhoi*vicen(i,j,n)*ice_ref_salinity*p001 &
* rside(i,j) / dt
fresh(i,j) = fresh(i,j) + dfresh
fsalt(i,j) = fsalt(i,j) + dfsalt
! history diagnostics
meltl(i,j) = meltl(i,j) + vicen(i,j,n)*rside(i,j)
! state variables
aicen(i,j,n) = aicen(i,j,n) * (c1 - rside(i,j))
vicen(i,j,n) = vicen(i,j,n) * (c1 - rside(i,j))
vsnon(i,j,n) = vsnon(i,j,n) * (c1 - rside(i,j))
enddo ! ij
do k = 1, nilyr
!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
do ij = 1, icells
i = indxi(ij)
j = indxj(ij)
! heat flux to coupler for ice melt (dfhocn < 0)
dfhocn = eicen(i,j,ilyr1(n)+k-1)*rside(i,j) / dt
fhocn(i,j) = fhocn(i,j) + dfhocn
! ice energy
eicen(i,j,ilyr1(n)+k-1) = eicen(i,j,ilyr1(n)+k-1) &
* (c1 - rside(i,j))
enddo ! ij
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)
! heat flux to coupler for snow melt (dfhocn < 0)
dfhocn = esnon(i,j,slyr1(n)+k-1)*rside(i,j) / dt
fhocn(i,j) = fhocn(i,j) + dfhocn
! snow energy
esnon(i,j,slyr1(n)+k-1) = esnon(i,j,slyr1(n)+k-1) &
* (c1 - rside(i,j))
enddo ! ij
enddo ! nslyr
enddo ! n
end subroutine lateral_melt
!=======================================================================
end module ice_therm_itd
!=======================================================================