My Project
Functions/Subroutines
ice_itd_linear Module Reference

Functions/Subroutines

subroutine linear_itd (hicen_old, hicen)
 
subroutine fit_line (ni, HbL, HbR, hice, g0, g1, hL, hR, remap_flag)
 

Function/Subroutine Documentation

◆ fit_line()

subroutine ice_itd_linear::fit_line ( integer (kind=int_kind), intent(in)  ni,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  HbL,
real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  HbR,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(in)  hice,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  g0,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  g1,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  hL,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out)  hR,
logical (kind=log_kind), dimension (ilo:ihi,jlo:jhi), intent(in)  remap_flag 
)

Definition at line 637 of file ice_itd_linear.f90.

637 !
638 ! !DESCRIPTION:
639 !
640 ! Fit g(h) with a line, satisfying area and volume constraints.
641 ! To reduce roundoff errors caused by large values of g0 and g1,
642 ! we actually compute g(eta), where eta = h - hL, and hL is the
643 ! left boundary.
644 !
645 ! !REVISION HISTORY:
646 !
647 ! authors: William H. Lipscomb, LANL
648 ! Elizabeth C. Hunke, LANL
649 !
650 ! !USES:
651 !
652 ! !INPUT/OUTPUT PARAMETERS:
653 !
654  integer (kind=int_kind), intent(in) :: ni ! category index
655 
656  real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), &
657  intent(in) :: &
658  hbl, hbr ! left and right category boundaries
659 
660  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(in) :: &
661  hice ! ice thickness
662 
663  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out) :: &
664  g0, g1 & ! coefficients in linear equation for g(eta)
665  , hl & ! min value of range over which g(h) > 0
666  , hr ! max value of range over which g(h) > 0
667 
668  logical (kind=log_kind), dimension (ilo:ihi,jlo:jhi), &
669  intent(in) :: &
670  remap_flag
671 !
672 !EOP
673 !
674  integer (kind=int_kind) :: &
675  i,j ! horizontal indices
676 
677  real (kind=dbl_kind) :: &
678  h13 & ! HbL + 1/3 * (HbR - HbL)
679  , h23 & ! HbL + 2/3 * (HbR - HbL)
680  , dhr & ! 1 / (hR - hL)
681  , wk1, wk2 ! temporary variables
682 
683 !! ggao 60162008
684  real (kind=dbl_kind) ::eps
685 !! change end
686 
687 
688  do j = jlo,jhi
689  do i = ilo,ihi
690 
691  if (remap_flag(i,j) .and. aicen(i,j,ni) > puny &
692  .and. hbr(i,j) - hbl(i,j) > puny) then
693 
694  ! Initialize hL and hR
695 
696  hl(i,j) = hbl(i,j)
697  hr(i,j) = hbr(i,j)
698 
699  ! Change hL or hR if hicen(n) falls outside central third of range
700 
701  h13 = p333 * (c2i*hl(i,j) + hr(i,j))
702  h23 = p333 * (hl(i,j) + c2i*hr(i,j))
703  if (hice(i,j) < h13) then
704  hr(i,j) = c3i*hice(i,j) - c2i*hl(i,j)
705  elseif (hice(i,j) > h23) then
706  hl(i,j) = c3i*hice(i,j) - c2i*hr(i,j)
707  endif
708 
709  ! Compute coefficients of g(eta) = g0 + g1*eta
710 
711 ! dhr = c1i / (hR(i,j) - hL(i,j))
712 !! ggao change 0616-2008
713  dhr = c1i / (hr(i,j) - hl(i,j)+epsilon(eps))
714 !! change end
715  wk1 = c6i * aicen(i,j,ni) * dhr
716  wk2 = (hice(i,j) - hl(i,j)) * dhr
717  g0(i,j) = wk1 * (p666 - wk2)
718  g1(i,j) = c2i*dhr * wk1 * (wk2 - p5)
719 
720  else ! remap_flag = .false. or aicen < puny
721  ! or hbR <= hbL
722  hl(i,j) = c0i
723  hr(i,j) = c0i
724  g0(i,j) = c0i
725  g1(i,j) = c0i
726 
727  endif ! aicen > puny
728 
729  enddo ! i
730  enddo ! j
731 
ice_domain::ihi
integer(kind=int_kind) ihi
Definition: ice_domain.f90:101
ice_constants::c0i
real(kind=dbl_kind), parameter c0i
Definition: ice_constants.f90:150
ice_constants::c3i
real(kind=dbl_kind), parameter c3i
Definition: ice_constants.f90:150
ice_constants::c6i
real(kind=dbl_kind), parameter c6i
Definition: ice_constants.f90:150
ice_domain::jlo
integer(kind=int_kind) jlo
Definition: ice_domain.f90:101
ice_domain::ilo
integer(kind=int_kind) ilo
Definition: ice_domain.f90:101
ice_constants::p5
real(kind=dbl_kind), parameter p5
Definition: ice_constants.f90:150
ice_constants::puny
real(kind=dbl_kind), parameter puny
Definition: ice_constants.f90:150
ice_domain::jhi
integer(kind=int_kind) jhi
Definition: ice_domain.f90:101
ice_constants::c2i
real(kind=dbl_kind), parameter c2i
Definition: ice_constants.f90:150
ice_state::aicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save aicen
Definition: ice_state.f90:97
ice_constants::p666
real(kind=dbl_kind), parameter p666
Definition: ice_constants.f90:150
ice_constants::c1i
real(kind=dbl_kind), parameter c1i
Definition: ice_constants.f90:150
ice_constants::p333
real(kind=dbl_kind), parameter p333
Definition: ice_constants.f90:150
Here is the caller graph for this function:

◆ linear_itd()

subroutine ice_itd_linear::linear_itd ( real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,ncat), intent(in)  hicen_old,
real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,ncat), intent(inout)  hicen 
)

Definition at line 69 of file ice_itd_linear.f90.

69 !
70 ! !DESCRIPTION:
71 !
72 ! Ice thickness distribution scheme that shifts ice among categories. \!
73 ! The default scheme is linear remapping, which works as follows. See
74 ! Lipscomb (2001) for more details. \!
75 ! Using the thermodynamic "velocities", interpolate to find the
76 ! velocities in thickness space at the category boundaries, and
77 ! compute the new locations of the boundaries. Then for each
78 ! category, compute the thickness distribution function, g(h),
79 ! between hL and hR, the left and right boundaries of the category.
80 ! Assume g(h) is a linear polynomial that satisfies two conditions: \!
81 ! (1) The ice area implied by g(h) equals aicen(n).
82 ! (2) The ice volume implied by g(h) equals aicen(n)*hicen(n).
83 !
84 ! Given g(h), at each boundary compute the ice area and volume lying
85 ! between the original and new boundary locations. Transfer area
86 ! and volume across each boundary in the appropriate direction, thus
87 ! restoring the original boundaries. See Lipscomb (2001) for details.
88 !
89 ! !REVISION HISTORY:
90 !
91 ! authors: William H. Lipscomb, LANL
92 ! Elizabeth C. Hunke, LANL
93 !
94 ! !USES:
95 !
96 ! !INPUT/OUTPUT PARAMETERS:
97 !
98  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,ncat), &
99  intent(in) :: &
100  hicen_old ! starting value of hicen (m)
101 
102  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,ncat), &
103  intent(inout) :: &
104  hicen ! ice thickness for each cat (m)
105 !
106 !EOP
107 !
108  integer (kind=int_kind) :: &
109  i, j & ! horizontal indices
110  , ni, nd & ! category indices
111  , k ! ice layer index
112 
113  real (kind=dbl_kind) :: &
114  slope &! rate of change of dhice with hice
115  , dh0 &! change in ice thickness at h = 0
116  , da0 &! area melting from category 1
117  , damax &! max allowed reduction in category 1 area
118  , etamin, etamax &! left and right limits of integration
119  , x1 &! etamax - etamin
120  , x2 &! (etamax^2 - etamin^2) / 2
121  , x3 &! (etamax^3 - etamin^3) / 3
122  , wk1, wk2 ! temporary variables
123 
124  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,0:ncat) :: &
125  hbnew ! new boundary locations
126 
127  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi) :: &
128  hb0 & ! hin_max(0)
129  , hb1 ! hin_max(1)
130 
131  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,ncat) :: &
132  dhicen & ! thickness change for remapping (m)
133  , g0 & ! constant coefficient in g(h)
134  , g1 & ! linear coefficient in g(h)
135  , hl & ! left end of range over which g(h) > 0
136  , hr ! right end of range over which g(h) > 0
137 
138  real (kind=dbl_kind), dimension(imt_local, jmt_local) :: &
139  vice_init, vice_final &! ice volume summed over categories
140  , vsno_init, vsno_final &! snow volume summed over categories
141  , eice_init, eice_final &! ice energy summed over categories
142  , esno_init, esno_final ! snow energy summed over categories
143 
144  ! NOTE: Third index of donor, daice, dvice should be ncat-1,
145  ! except that compilers would have trouble when ncat = 1
146  integer (kind=int_kind), dimension(ilo:ihi,jlo:jhi,ncat) :: &
147  donor ! donor category index
148 
149  real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,ncat) :: &
150  daice & ! ice area transferred across boundary
151  , dvice ! ice volume transferred across boundary
152 
153  logical (kind=log_kind), dimension(ilo:ihi,jlo:jhi) :: &
154  remap_flag ! remap ITD if remap_flag(i,j) is true
155 
156  character (len=char_len) :: &
157  fieldid ! field identifier
158 
159  logical (kind=log_kind), parameter :: &
160  l_conservation_check = .true. ! if true, check conservation
161  ! (useful for debugging)
162 
163  integer (kind=int_kind) :: &
164  icells & ! number of grid cells with ice
165  , ij ! combined horizontal index
166 
167  integer (kind=int_kind), dimension (1:(ihi-ilo+1)*(jhi-jlo+1)) :: &
168  indxi & ! compressed i/j indices
169  , indxj
170 
171  logical (kind=log_kind) :: &
172  flag_changed
173 
174 !! ggao 60162008
175  real (kind=dbl_kind) ::eps
176 !! change end
177 
178 
179  !-----------------------------------------------------------------
180  ! Compute volume and energy sums that linear remapping should
181  ! conserve.
182  !-----------------------------------------------------------------
183 
184  if (l_conservation_check) then
185  call column_sum (ncat, vicen, vice_init)
186  call column_sum (ncat, vsnon, vsno_init)
187  call column_sum (ntilay, eicen, eice_init)
188  call column_sum (ncat, esnon, esno_init)
189  endif
190 
191  !-----------------------------------------------------------------
192  ! Compute thickness change in each category.
193  !-----------------------------------------------------------------
194 
195  dhicen = c0i
196  do ni = 1, ncat
197  do j = jlo,jhi
198  do i = ilo,ihi
199  if (aicen(i,j,ni) > puny) then
200  dhicen(i,j,ni) = hicen(i,j,ni) - hicen_old(i,j,ni)
201  endif ! aicen > puny
202  enddo ! i
203  enddo ! j
204  enddo ! n
205 
206  !-----------------------------------------------------------------
207  ! Compute fractional ice area in each grid cell.
208  !-----------------------------------------------------------------
209  call aggregate_area
210 
211  !-----------------------------------------------------------------
212  ! Identify grid cells with ice and initialize remapping flag.
213  ! Remapping is done wherever remap_flag = .true.
214  ! In rare cases the category boundaries may shift too far for the
215  ! remapping algorithm to work, and remap_flag is set to .false.
216  ! In these cases the simpler 'rebin' subroutine will shift ice
217  ! between categories if needed.
218  !-----------------------------------------------------------------
219 
220  icells = 0
221  do j = jlo,jhi
222  do i = ilo,ihi
223  if (aice(i,j) > puny) then
224  remap_flag(i,j) = .true.
225  icells = icells + 1
226  indxi(icells) = i
227  indxj(icells) = j
228  else
229  remap_flag(i,j) = .false.
230  endif
231  enddo
232  enddo
233 
234  !-----------------------------------------------------------------
235  ! Compute new category boundaries, Hbnew, based on changes in
236  ! ice thickness from vertical thermodynamics.
237  !-----------------------------------------------------------------
238 
239  hin_max(ncat) = 999.9_dbl_kind ! arbitrary big number
240  hbnew = c0i
241 
242  do ni = 1, ncat-1
243 
244 !DIR$ CONCURRENT !Cray
245 !cdir nodep !NEC
246 !ocl novrec !Fujitsu
247  do ij = 1, icells ! aice(i,j) > puny
248  i = indxi(ij)
249  j = indxj(ij)
250 
251  if (hicen_old(i,j,ni) > puny .and. &
252  hicen_old(i,j,ni+1) > puny) then
253  ! interpolate between adjacent category growth rates
254  slope = (dhicen(i,j,ni+1)-dhicen(i,j,ni)) / &
255 ! (hicen_old(i,j,ni+1)-hicen_old(i,j,ni))
256 !! ggao change 0616-2008
257  (hicen_old(i,j,ni+1)-hicen_old(i,j,ni)+epsilon(eps))
258  hbnew(i,j,ni) = hin_max(ni) + dhicen(i,j,ni) &
259  + slope * (hin_max(ni) - hicen_old(i,j,ni))
260  elseif (hicen_old(i,j,ni) > puny) then ! hicen_old(ni+1)=0
261  hbnew(i,j,ni) = hin_max(ni) + dhicen(i,j,ni)
262  elseif (hicen_old(i,j,ni+1) > puny) then ! hicen_old(ni)=0
263  hbnew(i,j,ni) = hin_max(ni) + dhicen(i,j,ni+1)
264  else
265  hbnew(i,j,ni) = hin_max(ni)
266  endif
267  enddo ! ij
268 
269  !-----------------------------------------------------------------
270  ! Check that each boundary lies between adjacent values of hicen.
271  ! If not, set remap_flag = .false.
272  !-----------------------------------------------------------------
273 
274  flag_changed = .false.
275  do ij = 1, icells ! aice(i,j) > puny
276  i = indxi(ij)
277  j = indxj(ij)
278 
279  if (aicen(i,j,ni) > puny .and. &
280  hicen(i,j,ni) >= hbnew(i,j,ni)) then
281  remap_flag(i,j) = .false.
282  flag_changed = .true.
283  elseif (aicen(i,j,ni+1) > puny .and. &
284  hicen(i,j,ni+1) <= hbnew(i,j,ni)) then
285  remap_flag(i,j) = .false.
286  flag_changed = .true.
287  endif
288 
289  !-----------------------------------------------------------------
290  ! Check that Hbnew(n) lies between hin_max(n-1) and hin_max(n+1).
291  ! If not, set remap_flag = .false.
292  ! (In principle we could allow this, but it would make the code
293  ! more complicated.)
294  !-----------------------------------------------------------------
295 
296  if (hbnew(i,j,ni) > hin_max(ni+1)) then
297  remap_flag(i,j) = .false.
298  flag_changed = .true.
299  endif
300 
301  if (hbnew(i,j,ni) < hin_max(ni-1)) then
302  remap_flag(i,j) = .false.
303  flag_changed = .true.
304  endif
305 
306  enddo ! ij
307 
308  !-----------------------------------------------------------------
309  ! Write diagnosis outputs if remap_flag was changed to false
310  !-----------------------------------------------------------------
311 
312  if (flag_changed) then
313  do ij = 1, icells ! aice(i,j) > puny
314 
315  i = indxi(ij)
316  j = indxj(ij)
317 
318  if (aicen(i,j,ni) > puny .and. &
319  hicen(i,j,ni) >= hbnew(i,j,ni)) then
320  write(nu_diag,*) 'istep1 = ',istep1
321  write(nu_diag,*) my_task,':',i,j, &
322  'ITD: hicen(ni) > Hbnew(ni)'
323  write(nu_diag,*) 'cat ',ni
324  write(nu_diag,*) my_task,':',i,j, &
325  'hicen(ni) =', hicen(i,j,ni)
326  write(nu_diag,*) my_task,':',i,j, &
327  'Hbnew(ni) =', hbnew(i,j,ni)
328  elseif (aicen(i,j,ni+1) > puny .and. &
329  hicen(i,j,ni+1) <= hbnew(i,j,ni)) then
330  write(nu_diag,*) 'istep1 = ',istep1
331  write(nu_diag,*) my_task,':',i,j, &
332  'ITD: hicen(ni+1) < Hbnew(ni)'
333  write(nu_diag,*) 'cat ',ni
334  write(nu_diag,*) my_task,':',i,j, &
335  'hicen(ni+1) =', hicen(i,j,ni+1)
336  write(nu_diag,*) my_task,':',i,j, &
337  'Hbnew(ni) =', hbnew(i,j,ni)
338  endif
339 
340  if (hbnew(i,j,ni) > hin_max(ni+1)) then
341  write(nu_diag,*) 'istep1 = ',istep1
342  write(nu_diag,*) my_task,':',i,j, &
343  'ITD Hbnew(ni) > hin_max(ni+1)'
344  write(nu_diag,*) 'cat ',ni
345  write(nu_diag,*) my_task,':',i,j, &
346  'Hbnew(ni) =', hbnew(i,j,ni)
347  write(nu_diag,*) my_task,':',i,j, &
348  'hin_max(ni+1) =', hin_max(ni+1)
349  endif
350 
351  if (hbnew(i,j,ni) < hin_max(ni-1)) then
352  write(nu_diag,*) 'istep1 = ',istep1
353  write(nu_diag,*) my_task,':',i,j, &
354  'ITD: Hbnew(ni) < hin_max(ni-1)'
355  write(nu_diag,*) 'cat ',ni
356  write(nu_diag,*) my_task,':',i,j, &
357  'Hbnew(ni) =', hbnew(i,j,ni)
358  write(nu_diag,*) my_task,':',i,j, &
359  'hin_max(ni-1) =', hin_max(ni-1)
360  endif
361 
362  enddo ! ij
363  endif ! flag_changed
364 
365  enddo ! boundaries, 1 to ncat-1
366 
367  !-----------------------------------------------------------------
368  ! Identify cells where the ITD is to be remapped
369  !-----------------------------------------------------------------
370 
371  icells = 0
372  do j = jlo,jhi
373  do i = ilo,ihi
374  if (remap_flag(i,j)) then
375  icells = icells + 1
376  indxi(icells) = i
377  indxj(icells) = j
378  endif
379  enddo
380  enddo
381 
382  !-----------------------------------------------------------------
383  ! Fill arrays with initial boundaries of category 1
384  ! Prescribe Hbnew(0) and Hbnew(ncat)
385  !-----------------------------------------------------------------
386 
387  do j = jlo,jhi
388  do i = ilo,ihi
389  hb0(i,j) = hin_max(0)
390  hb1(i,j) = hin_max(1)
391 
392  hbnew(i,j,0) = c0i
393 
394  if (aicen(i,j,ncat) > puny) then
395  hbnew(i,j,ncat) = c3i*hicen(i,j,ncat) - c2i*hbnew(i,j,ncat-1)
396  else
397  hbnew(i,j,ncat) = hin_max(ncat)
398  endif
399 
400  if (hbnew(i,j,ncat) < hin_max(ncat-1)) &
401  hbnew(i,j,ncat) = hin_max(ncat-1)
402  enddo
403  enddo
404 
405  !-----------------------------------------------------------------
406  ! Compute g(h) for category 1 at start of time step
407  ! (hicen = hicen_old)
408  !-----------------------------------------------------------------
409 
410  call fit_line(1, hb0, hb1, hicen_old(:,:,1), &
411  g0(:,:,1), g1(:,:,1), hl(:,:,1), hr(:,:,1), &
412  remap_flag)
413 
414  !-----------------------------------------------------------------
415  ! Find area lost due to melting of thin (category 1) ice
416  !-----------------------------------------------------------------
417 
418 !DIR$ CONCURRENT !Cray
419 !cdir nodep !NEC
420 !ocl novrec !Fujitsu
421  do ij = 1, icells ! remap_flag = .true.
422  i = indxi(ij)
423  j = indxj(ij)
424 
425  if (aicen(i,j,1) > puny) then
426 
427  dh0 = dhicen(i,j,1)
428 
429  if (dh0 < c0i) then ! remove area from category 1
430  dh0 = min(-dh0,hin_max(1)) ! dh0 --> |dh0|
431 
432  !-----------------------------------------------------------------
433  ! Integrate g(1) from 0 to dh0 to estimate area melted
434  !-----------------------------------------------------------------
435 
436  ! right integration limit (left limit = 0)
437  etamax = min(dh0,hr(i,j,1)) - hl(i,j,1)
438 
439  if (etamax > c0i) then
440  x1 = etamax
441  x2 = p5 * etamax*etamax
442  da0 = g1(i,j,1)*x2 + g0(i,j,1)*x1 ! ice area removed
443 
444  ! constrain new thickness <= hicen_old
445  damax = aicen(i,j,1) &
446  * (c1i-hicen(i,j,1)/hicen_old(i,j,1)) ! damax > 0
447  da0 = min(da0, damax)
448 
449  ! remove area, conserving volume
450  hicen(i,j,1) = hicen(i,j,1) &
451 ! * aicen(i,j,1) / (aicen(i,j,1)-da0)
452 !! gao change 0616-2008
453  * aicen(i,j,1) / (aicen(i,j,1)-da0+epsilon(eps))
454  aicen(i,j,1) = aicen(i,j,1) - da0
455  endif ! etamax > 0
456 
457  else ! dh0 >= 0
458  hbnew(i,j,0) = min(dh0,hin_max(1)) ! shift Hbnew(0) to right
459  endif
460 
461  endif ! aicen(i,j,1) > puny
462  enddo ! ij
463 
464  !-----------------------------------------------------------------
465  ! Compute g(h) for each ice thickness category.
466  !-----------------------------------------------------------------
467 
468  do ni = 1, ncat
469  call fit_line(ni, hbnew(:,:,ni-1), hbnew(:,:,ni), hicen(:,:,ni),&
470  g0(:,:,ni), g1(:,:,ni), hl(:,:,ni), hr(:,:,ni),&
471  remap_flag)
472  enddo
473 
474  !-----------------------------------------------------------------
475  ! Compute area and volume to be shifted across each boundary.
476  !-----------------------------------------------------------------
477 
478  donor(:,:,:) = 0
479  daice(:,:,:) = c0i
480  dvice(:,:,:) = c0i
481 
482  do ni = 1, ncat-1
483 !DIR$ CONCURRENT !Cray
484 !cdir nodep !NEC
485 !ocl novrec !Fujitsu
486  do ij = 1, icells ! remap_flag = .true.
487  i = indxi(ij)
488  j = indxj(ij)
489 
490  if (hbnew(i,j,ni) > hin_max(ni)) then ! transfer from n to n+1
491 
492  ! left and right integration limits in eta space
493  etamin = max(hin_max(ni), hl(i,j,ni)) - hl(i,j,ni)
494  etamax = min(hbnew(i,j,ni), hr(i,j,ni)) - hl(i,j,ni)
495  donor(i,j,ni) = ni
496 
497  else ! Hbnew(n) <= hin_max(n); transfer from n+1 to n
498 
499  ! left and right integration limits in eta space
500  etamin = c0i
501  etamax = min(hin_max(ni), hr(i,j,ni+1)) - hl(i,j,ni+1)
502  donor(i,j,ni) = ni+1
503 
504  endif ! Hbnew(n) > hin_max(n)
505 
506  if (etamax > etamin) then
507  x1 = etamax - etamin
508  wk1 = etamin*etamin
509  wk2 = etamax*etamax
510  x2 = p5 * (wk2 - wk1)
511  wk1 = wk1*etamin
512  wk2 = wk2*etamax
513  x3 = p333 * (wk2 - wk1)
514  nd = donor(i,j,ni)
515  daice(i,j,ni) = g1(i,j,nd)*x2 + g0(i,j,nd)*x1
516  if (daice(i,j,ni) > c0i) then
517  dvice(i,j,ni) = g1(i,j,nd)*x3 + g0(i,j,nd)*x2 &
518  + daice(i,j,ni)*hl(i,j,nd)
519  else
520  daice(i,j,ni) = c0i
521  donor(i,j,ni) = 0
522  endif
523  endif
524 
525 
526  ! If daice or dvice is very small, shift no ice.
527 
528  nd = donor(i,j,ni)
529 
530 !!!===============================
531  IF(nd>0) THEN
532 !!!===============================
533 
534  if (daice(i,j,ni) < aicen(i,j,nd)*puny) then
535  daice(i,j,ni) = c0i
536  dvice(i,j,ni) = c0i
537  donor(i,j,ni) = 0
538  endif
539 
540  if (dvice(i,j,ni) < vicen(i,j,nd)*puny) then
541  daice(i,j,ni) = c0i
542  dvice(i,j,ni) = c0i
543  donor(i,j,ni) = 0
544  endif
545 
546  ! If daice is close to aicen or dvice is close to vicen,
547  ! shift entire category
548 
549  if (daice(i,j,ni) > aicen(i,j,nd)*(c1i-puny)) then
550  daice(i,j,ni) = aicen(i,j,nd)
551  dvice(i,j,ni) = vicen(i,j,nd)
552  endif
553 
554  if (dvice(i,j,ni) > vicen(i,j,nd)*(c1i-puny)) then
555  daice(i,j,ni) = aicen(i,j,nd)
556  dvice(i,j,ni) = vicen(i,j,nd)
557  endif
558 !!! ggao fix the bug 02102008
559  ELSE
560 
561  daice(i,j,ni) = c0i
562  dvice(i,j,ni) = c0i
563  donor(i,j,ni) = 0
564  ENDIF
565 !!! change end
566 
567  enddo ! ij
568  enddo ! boundaries, 1 to ncat-1
569 
570  !-----------------------------------------------------------------
571  ! Shift ice between categories as necessary
572  !-----------------------------------------------------------------
573 
574  call shift_ice (donor, daice, dvice, hicen)
575 
576  !-----------------------------------------------------------------
577  ! Make sure hice(i,j,1) >= minimum ice thickness hi_min.
578  !-----------------------------------------------------------------
579 
580 !DIR$ CONCURRENT !Cray
581 !cdir nodep !NEC
582 !ocl novrec !Fujitsu
583  do ij = 1, icells ! remap_flag = .true.
584  i = indxi(ij)
585  j = indxj(ij)
586  if (hi_min > c0i .and. &
587  aicen(i,j,1) > puny .and. hicen(i,j,1) < hi_min) then
588  aicen(i,j,1) = aicen(i,j,1) * hicen(i,j,1)/hi_min
589  hicen(i,j,1) = hi_min
590  endif
591  enddo ! ij
592 
593  !-----------------------------------------------------------------
594  ! Update ice and open water area.
595  !-----------------------------------------------------------------
596  call aggregate_area
597 
598  !-----------------------------------------------------------------
599  ! Check volume and energy conservation.
600  !-----------------------------------------------------------------
601 
602  if (l_conservation_check) then
603 
604  call column_sum (ncat, vicen, vice_final)
605  fieldid = 'vice, ITD remap'
606  call column_conservation_check (vice_init, vice_final, &
607  puny, fieldid)
608 
609  call column_sum (ncat, vsnon, vsno_final)
610  fieldid = 'vsno, ITD remap'
611  call column_conservation_check (vsno_init, vsno_final, &
612  puny, fieldid)
613 
614  call column_sum (ntilay, eicen, eice_final)
615  fieldid = 'eice, ITD remap'
616  call column_conservation_check (eice_init, eice_final, &
617  puny*lfresh*rhoi, fieldid)
618 
619  call column_sum (ncat, esnon, esno_final)
620  fieldid = 'esno, ITD remap'
621  call column_conservation_check (esno_init, esno_final, &
622  puny*lfresh*rhos, fieldid)
623 
624  endif ! conservation check
625 
ice_state::esnon
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save esnon
Definition: ice_state.f90:97
ice_therm_vertical::hicen_old
real(kind=dbl_kind), dimension(:,:,:), allocatable, save hicen_old
Definition: ice_therm_vertical.f90:84
ice_domain::ihi
integer(kind=int_kind) ihi
Definition: ice_domain.f90:101
ice_constants::c0i
real(kind=dbl_kind), parameter c0i
Definition: ice_constants.f90:150
ice_constants::c3i
real(kind=dbl_kind), parameter c3i
Definition: ice_constants.f90:150
ice_constants::rhos
real(kind=dbl_kind), parameter rhos
Definition: ice_constants.f90:118
ice_constants::lfresh
real(kind=dbl_kind), parameter lfresh
Definition: ice_constants.f90:85
ice_domain::jlo
integer(kind=int_kind) jlo
Definition: ice_domain.f90:101
ice_domain::ilo
integer(kind=int_kind) ilo
Definition: ice_domain.f90:101
ice_constants::p5
real(kind=dbl_kind), parameter p5
Definition: ice_constants.f90:150
ice_constants::puny
real(kind=dbl_kind), parameter puny
Definition: ice_constants.f90:150
ice_domain::jhi
integer(kind=int_kind) jhi
Definition: ice_domain.f90:101
ice_state::vicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save vicen
Definition: ice_state.f90:97
ice_state::aice
real(kind=dbl_kind), dimension(:,:), allocatable, target, save aice
Definition: ice_state.f90:82
ice_state::eicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save eicen
Definition: ice_state.f90:108
ice_model_size::ncat
integer(kind=int_kind), parameter ncat
Definition: ice_model_size.f90:48
ice_constants::c2i
real(kind=dbl_kind), parameter c2i
Definition: ice_constants.f90:150
ice_state::aicen
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save aicen
Definition: ice_state.f90:97
ice_domain::my_task
integer(kind=int_kind), save my_task
Definition: ice_domain.f90:95
ice_constants::c1i
real(kind=dbl_kind), parameter c1i
Definition: ice_constants.f90:150
ice_constants::rhoi
real(kind=dbl_kind), parameter rhoi
Definition: ice_constants.f90:85
ice_model_size::ntilay
integer(kind=int_kind), parameter ntilay
Definition: ice_model_size.f90:48
ice_constants::p333
real(kind=dbl_kind), parameter p333
Definition: ice_constants.f90:150
ice_state::vsnon
real(kind=dbl_kind), dimension(:,:,:), allocatable, target, save vsnon
Definition: ice_state.f90:97
Here is the call graph for this function:
Here is the caller graph for this function:
Generated by   doxygen 1.8.15