00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 module ice_therm_itd
00030
00031
00032
00033 use ice_kinds_mod
00034 use ice_communicate, only: my_task, master_task
00035 use ice_domain_size
00036 use ice_constants
00037 use ice_fileunits
00038
00039
00040
00041 implicit none
00042 save
00043
00044 real (kind=dbl_kind), parameter, private ::
00045 hfrazilmin = 0.05_dbl_kind
00046
00047
00048
00049 contains
00050
00051
00052
00053
00054
00055
00056
00057
00058 subroutine linear_itd (nx_block, ny_block, &
00059 icells, indxi, indxj, &
00060 trcr_depend, &
00061 aicen_init, vicen_init, &
00062 aicen, trcrn, &
00063 vicen, vsnon, &
00064 eicen, esnon, &
00065 aice, aice0, &
00066 l_stop, &
00067 istop, jstop)
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094 use ice_itd, only: hin_max, hi_min, aggregate_area, shift_ice, &
00095 column_sum, column_conservation_check
00096
00097
00098
00099 integer (kind=int_kind), intent(in) ::
00100 nx_block, ny_block,
00101 icells
00102
00103 integer (kind=int_kind), dimension (nx_block*ny_block),
00104 intent(in) ::
00105 indxi, indxj
00106
00107 integer (kind=int_kind), dimension (ntrcr), intent(in) ::
00108 trcr_depend
00109
00110 real (kind=dbl_kind), dimension(nx_block,ny_block,ncat),
00111 intent(in) ::
00112 aicen_init,
00113 vicen_init
00114
00115 real (kind=dbl_kind), dimension (nx_block,ny_block,ncat),
00116 intent(inout) ::
00117 aicen ,
00118 vicen ,
00119 vsnon
00120
00121 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat),
00122 intent(inout) ::
00123 trcrn
00124
00125 real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr),
00126 intent(inout) ::
00127 eicen
00128
00129 real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr),
00130 intent(inout) ::
00131 esnon
00132
00133 real (kind=dbl_kind), dimension (nx_block,ny_block),
00134 intent(inout) ::
00135 aice ,
00136 aice0
00137
00138 logical (kind=log_kind), intent(out) ::
00139 l_stop
00140
00141 integer (kind=int_kind), intent(out) ::
00142 istop, jstop
00143
00144
00145
00146 integer (kind=int_kind) ::
00147 i, j ,
00148 n, nd ,
00149 k
00150
00151 real (kind=dbl_kind) ::
00152 slope ,
00153 dh0 ,
00154 da0 ,
00155 damax ,
00156 etamin, etamax,
00157 x1 ,
00158 x2 ,
00159 x3 ,
00160 wk1, wk2
00161
00162 real (kind=dbl_kind), dimension(icells,0:ncat) ::
00163 hbnew
00164
00165 real (kind=dbl_kind), dimension(icells) ::
00166 work
00167
00168 integer (kind=int_kind), dimension(icells) ::
00169 indxii, indxjj,
00170 indxij
00171
00172 real (kind=dbl_kind), dimension(:,:), allocatable ::
00173 g0 ,
00174 g1 ,
00175 hL ,
00176 hR
00177
00178 real (kind=dbl_kind), dimension(icells,ncat) ::
00179 hicen ,
00180 hicen_init ,
00181 dhicen ,
00182 daice ,
00183 dvice
00184
00185 real (kind=dbl_kind), dimension(icells) ::
00186 vice_init, vice_final,
00187 vsno_init, vsno_final,
00188 eice_init, eice_final,
00189 esno_init, esno_final
00190
00191
00192
00193 integer (kind=int_kind), dimension(icells,ncat) ::
00194 donor
00195
00196 logical (kind=log_kind), dimension(icells) ::
00197 remap_flag
00198
00199 character (len=char_len) ::
00200 fieldid
00201
00202 logical (kind=log_kind), parameter ::
00203 l_conservation_check = .true.
00204
00205
00206
00207 integer (kind=int_kind) ::
00208 iflag ,
00209 ij, m
00210
00211
00212
00213
00214
00215 l_stop = .false.
00216 istop = 0
00217 jstop = 0
00218
00219 hin_max(ncat) = 999.9_dbl_kind
00220
00221
00222
00223
00224
00225
00226 if (l_conservation_check) then
00227 call column_sum (nx_block, ny_block, &
00228 icells, indxi, indxj, &
00229 ncat, &
00230 vicen, vice_init)
00231 call column_sum (nx_block, ny_block, &
00232 icells, indxi, indxj, &
00233 ncat, &
00234 vsnon, vsno_init)
00235 call column_sum (nx_block, ny_block, &
00236 icells, indxi, indxj, &
00237 ntilyr, &
00238 eicen, eice_init)
00239 call column_sum (nx_block, ny_block, &
00240 icells, indxi, indxj, &
00241 ntslyr, &
00242 esnon, esno_init)
00243 endif
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254 do ij = 1, icells
00255 remap_flag(ij) = .true.
00256 enddo
00257
00258
00259
00260
00261
00262 do n = 1, ncat
00263 do ij = 1, icells
00264 i = indxi(ij)
00265 j = indxj(ij)
00266
00267 if (aicen_init(i,j,n) > puny) then
00268 hicen_init (ij,n) = vicen_init(i,j,n) / &
00269 aicen_init(i,j,n)
00270 else
00271 hicen_init(ij,n) = c0
00272 endif
00273
00274 if (aicen(i,j,n) > puny) then
00275 hicen (ij,n) = vicen(i,j,n) / aicen(i,j,n)
00276 dhicen(ij,n) = hicen(ij,n) - hicen_init(ij,n)
00277 else
00278 hicen (ij,n) = c0
00279 dhicen(ij,n) = c0
00280 endif
00281
00282 enddo
00283 enddo
00284
00285
00286
00287
00288
00289
00290 do ij = 1, icells
00291 hbnew(ij,0) = hin_max(0)
00292 enddo
00293
00294 do n = 1, ncat-1
00295
00296
00297
00298
00299 do ij = 1, icells
00300 i = indxi(ij)
00301 j = indxj(ij)
00302
00303 if (hicen_init(ij,n) > puny .and. &
00304 hicen_init(ij,n+1) > puny) then
00305
00306 slope = (dhicen(ij,n+1) - dhicen(ij,n)) / &
00307 (hicen_init(ij,n+1) - hicen_init(ij,n))
00308 hbnew(ij,n) = hin_max(n) + dhicen(ij,n) &
00309 + slope * (hin_max(n) - hicen_init(ij,n))
00310 elseif (hicen_init(ij,n) > puny) then
00311 hbnew(ij,n) = hin_max(n) + dhicen(ij,n)
00312 elseif (hicen_init(ij,n+1) > puny) then
00313 hbnew(ij,n) = hin_max(n) + dhicen(ij,n+1)
00314 else
00315 hbnew(ij,n) = hin_max(n)
00316 endif
00317
00318
00319
00320
00321
00322
00323
00324 if (aicen(i,j,n) > puny .and. &
00325 hicen(ij,n) >= hbnew(ij,n)) then
00326 remap_flag(ij) = .false.
00327
00328 write(nu_diag,*) my_task,':',i,j, &
00329 'ITD: hicen(n) > hbnew(n)'
00330 write(nu_diag,*) 'cat ',n
00331 write(nu_diag,*) my_task,':',i,j, &
00332 'hicen(n) =', hicen(ij,n)
00333 write(nu_diag,*) my_task,':',i,j, &
00334 'hbnew(n) =', hbnew(ij,n)
00335
00336 elseif (aicen(i,j,n+1) > puny .and. &
00337 hicen(ij,n+1) <= hbnew(ij,n)) then
00338 remap_flag(ij) = .false.
00339
00340 write(nu_diag,*) my_task,':',i,j, &
00341 'ITD: hicen(n+1) < hbnew(n)'
00342 write(nu_diag,*) 'cat ',n
00343 write(nu_diag,*) my_task,':',i,j, &
00344 'hicen(n+1) =', hicen(ij,n+1)
00345 write(nu_diag,*) my_task,':',i,j, &
00346 'hbnew(n) =', hbnew(ij,n)
00347 endif
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357 if (hbnew(ij,n) > hin_max(n+1)) then
00358 remap_flag(ij) = .false.
00359
00360 write(nu_diag,*) my_task,':',i,j, &
00361 'ITD hbnew(n) > hin_max(n+1)'
00362 write(nu_diag,*) 'cat ',n
00363 write(nu_diag,*) my_task,':',i,j, &
00364 'hbnew(n) =', hbnew(ij,n)
00365 write(nu_diag,*) my_task,':',i,j, &
00366 'hin_max(n+1) =', hin_max(n+1)
00367 endif
00368
00369 if (hbnew(ij,n) < hin_max(n-1)) then
00370 remap_flag(ij) = .false.
00371
00372 write(nu_diag,*) my_task,':',i,j, &
00373 'ITD: hbnew(n) < hin_max(n-1)'
00374 write(nu_diag,*) 'cat ',n
00375 write(nu_diag,*) my_task,':',i,j, &
00376 'hbnew(n) =', hbnew(ij,n)
00377 write(nu_diag,*) my_task,':',i,j, &
00378 'hin_max(n-1) =', hin_max(n-1)
00379 endif
00380
00381 enddo
00382
00383 enddo
00384
00385
00386
00387
00388
00389
00390
00391
00392 do ij = 1, icells
00393 i = indxi(ij)
00394 j = indxj(ij)
00395 if (aicen(i,j,ncat) > puny) then
00396 hbnew(ij,ncat) = c3*hicen(ij,ncat) - c2*hbnew(ij,ncat-1)
00397 else
00398 hbnew(ij,ncat) = hin_max(ncat)
00399 endif
00400 hbnew(ij,ncat) = max(hbnew(ij,ncat),hin_max(ncat-1))
00401 enddo
00402
00403
00404
00405
00406
00407 iflag = 0
00408 do ij = 1, icells
00409 i = indxi(ij)
00410 j = indxj(ij)
00411 if (remap_flag(ij)) then
00412 iflag = iflag + 1
00413 indxii(iflag) = i
00414 indxjj(iflag) = j
00415 indxij(iflag) = ij
00416 endif
00417 enddo
00418
00419 allocate(g0(iflag,ncat))
00420 allocate(g1(iflag,ncat))
00421 allocate(hL(iflag,ncat))
00422 allocate(hR(iflag,ncat))
00423
00424
00425
00426
00427
00428
00429 do ij = 1, icells
00430 work(ij) = hin_max(1)
00431 enddo
00432
00433 call fit_line(nx_block, ny_block, &
00434 iflag, icells, &
00435 indxii, indxjj, indxij, &
00436 aicen(:,:,1), hicen_init(:,1), &
00437 hbnew(:,0), work (:), &
00438 g0 (:,1), g1 (:,1), &
00439 hL (:,1), hR (:,1))
00440
00441
00442
00443
00444
00445
00446
00447
00448 do ij = 1, iflag
00449 i = indxii(ij)
00450 j = indxjj(ij)
00451 m = indxij(ij)
00452
00453 if (aicen(i,j,1) > puny) then
00454
00455 dh0 = dhicen(m,1)
00456
00457 if (dh0 < c0) then
00458 dh0 = min(-dh0,hin_max(1))
00459
00460
00461
00462
00463
00464
00465 etamax = min(dh0,hR(ij,1)) - hL(ij,1)
00466
00467 if (etamax > c0) then
00468 x1 = etamax
00469 x2 = p5 * etamax*etamax
00470 da0 = g1(ij,1)*x2 + g0(ij,1)*x1
00471
00472
00473 damax = aicen(i,j,1) &
00474 * (c1-hicen(m,1)/hicen_init(m,1))
00475 da0 = min (da0, damax)
00476
00477
00478 hicen(m,1) = hicen(m,1) &
00479 * aicen(i,j,1) / (aicen(i,j,1)-da0)
00480 aicen(i,j,1) = aicen(i,j,1) - da0
00481 endif
00482
00483 else
00484 hbnew(m,0) = min(dh0,hin_max(1))
00485 endif
00486
00487 endif
00488 enddo
00489
00490
00491
00492
00493
00494 do n = 1, ncat
00495 call fit_line(nx_block, ny_block, &
00496 iflag, icells, &
00497 indxii, indxjj, indxij, &
00498 aicen(:,:,n), hicen(:,n), &
00499 hbnew(:,n-1), hbnew(:,n), &
00500 g0 (:,n), g1 (:,n), &
00501 hL (:,n), hR (:,n))
00502
00503 enddo
00504
00505
00506
00507
00508
00509 do n = 1, ncat
00510 do ij = 1, icells
00511 donor(ij,n) = 0
00512 daice(ij,n) = c0
00513 dvice(ij,n) = c0
00514 enddo
00515 enddo
00516
00517 do n = 1, ncat-1
00518
00519
00520
00521
00522 do ij = 1, iflag
00523 i = indxii(ij)
00524 j = indxjj(ij)
00525 m = indxij(ij)
00526
00527 if (hbnew(m,n) > hin_max(n)) then
00528
00529
00530 etamin = max(hin_max(n), hL(ij,n)) - hL(ij,n)
00531 etamax = min(hbnew(m,n), hR(ij,n)) - hL(ij,n)
00532 donor(m,n) = n
00533
00534 else
00535
00536
00537 etamin = c0
00538 etamax = min(hin_max(n), hR(ij,n+1)) - hL(ij,n+1)
00539 donor(m,n) = n+1
00540
00541 endif
00542
00543 if (etamax > etamin) then
00544 x1 = etamax - etamin
00545 wk1 = etamin*etamin
00546 wk2 = etamax*etamax
00547 x2 = p5 * (wk2 - wk1)
00548 wk1 = wk1*etamin
00549 wk2 = wk2*etamax
00550 x3 = p333 * (wk2 - wk1)
00551 nd = donor(m,n)
00552 daice(m,n) = g1(ij,nd)*x2 + g0(ij,nd)*x1
00553 dvice(m,n) = g1(ij,nd)*x3 + g0(ij,nd)*x2 &
00554 + daice(m,n)*hL(ij,nd)
00555 endif
00556
00557
00558
00559 nd = donor(m,n)
00560
00561 if (daice(m,n) < aicen(i,j,nd)*puny) then
00562 daice(m,n) = c0
00563 dvice(m,n) = c0
00564 donor(m,n) = 0
00565 endif
00566
00567 if (dvice(m,n) < vicen(i,j,nd)*puny) then
00568 daice(m,n) = c0
00569 dvice(m,n) = c0
00570 donor(m,n) = 0
00571 endif
00572
00573
00574
00575
00576 if (daice(m,n) > aicen(i,j,nd)*(c1-puny)) then
00577 daice(m,n) = aicen(i,j,nd)
00578 dvice(m,n) = vicen(i,j,nd)
00579 endif
00580
00581 if (dvice(m,n) > vicen(i,j,nd)*(c1-puny)) then
00582 daice(m,n) = aicen(i,j,nd)
00583 dvice(m,n) = vicen(i,j,nd)
00584 endif
00585
00586 enddo
00587 enddo
00588
00589 deallocate(g0)
00590 deallocate(g1)
00591 deallocate(hL)
00592 deallocate(hR)
00593
00594
00595
00596
00597
00598 call shift_ice (nx_block, ny_block, &
00599 indxi, indxj, &
00600 icells, trcr_depend, &
00601 aicen, trcrn, &
00602 vicen, vsnon, &
00603 eicen, esnon, &
00604 hicen, donor, &
00605 daice, dvice, &
00606 l_stop, &
00607 istop, jstop)
00608
00609 if (l_stop) return
00610
00611
00612
00613
00614
00615
00616
00617
00618 do ij = 1, iflag
00619 i = indxii(ij)
00620 j = indxjj(ij)
00621 m = indxij(ij)
00622 if (hi_min > c0 .and. &
00623 aicen(i,j,1) > puny .and. hicen(m,1) < hi_min) then
00624 aicen(i,j,1) = aicen(i,j,1) * hicen(m,1)/hi_min
00625 hicen(m,1) = hi_min
00626 endif
00627 enddo
00628
00629
00630
00631
00632 call aggregate_area (nx_block, ny_block, &
00633 aicen, &
00634 aice, aice0)
00635
00636 if (l_stop) return
00637
00638
00639
00640
00641
00642 if (l_conservation_check) then
00643
00644 call column_sum (nx_block, ny_block, &
00645 icells, indxi, indxj, &
00646 ncat, &
00647 vicen, vice_final)
00648 fieldid = 'vice, ITD remap'
00649 call column_conservation_check (nx_block, ny_block, &
00650 icells, indxi, indxj, &
00651 fieldid, &
00652 vice_init, vice_final, &
00653 puny, l_stop, &
00654 istop, jstop)
00655 if (l_stop) return
00656
00657 call column_sum (nx_block, ny_block, &
00658 icells, indxi, indxj, &
00659 ncat, &
00660 vsnon, vsno_final)
00661 fieldid = 'vsno, ITD remap'
00662 call column_conservation_check (nx_block, ny_block, &
00663 icells, indxi, indxj, &
00664 fieldid, &
00665 vsno_init, vsno_final, &
00666 puny, l_stop, &
00667 istop, jstop)
00668 if (l_stop) return
00669
00670 call column_sum (nx_block, ny_block, &
00671 icells, indxi, indxj, &
00672 ntilyr, &
00673 eicen, eice_final)
00674 fieldid = 'eice, ITD remap'
00675 call column_conservation_check (nx_block, ny_block, &
00676 icells, indxi, indxj, &
00677 fieldid, &
00678 eice_init, eice_final, &
00679 puny*Lfresh*rhoi, &
00680 l_stop, &
00681 istop, jstop)
00682 if (l_stop) return
00683
00684 call column_sum (nx_block, ny_block, &
00685 icells, indxi, indxj, &
00686 ntslyr, &
00687 esnon, esno_final)
00688 fieldid = 'esno, ITD remap'
00689 call column_conservation_check (nx_block, ny_block, &
00690 icells, indxi, indxj, &
00691 fieldid, &
00692 esno_init, esno_final, &
00693 puny*Lfresh*rhos, &
00694 l_stop, &
00695 istop, jstop)
00696 if (l_stop) return
00697
00698 endif
00699
00700 end subroutine linear_itd
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710 subroutine fit_line (nx_block, ny_block, &
00711 iflag, icells, &
00712 indxii, indxjj, indxij, &
00713 aicen, hice, &
00714 hbL, hbR, &
00715 g0, g1, &
00716 hL, hR)
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734 integer (kind=int_kind), intent(in) ::
00735 nx_block, ny_block,
00736 icells ,
00737 iflag
00738
00739 integer (kind=int_kind), dimension (icells),
00740 intent(in) ::
00741 indxii, indxjj,
00742 indxij
00743
00744 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
00745 aicen
00746
00747 real (kind=dbl_kind), dimension (icells), intent(in) ::
00748 hbL, hbR ,
00749 hice
00750
00751 real (kind=dbl_kind), dimension (iflag), intent(out)::
00752 g0, g1 ,
00753 hL ,
00754 hR
00755
00756
00757
00758 integer (kind=int_kind) ::
00759 i,j ,
00760 ij, m
00761
00762 real (kind=dbl_kind) ::
00763 h13 ,
00764 h23 ,
00765 dhr ,
00766 wk1, wk2
00767
00768
00769
00770
00771
00772 do ij = 1, iflag
00773 i = indxii(ij)
00774 j = indxjj(ij)
00775 m = indxij(ij)
00776
00777 if (aicen(i,j) > puny .and. hbR(m) - hbL(m) > puny) then
00778
00779
00780
00781 hL(ij) = hbL(m)
00782 hR(ij) = hbR(m)
00783
00784
00785
00786 h13 = p333 * (c2*hL(ij) + hR(ij))
00787 h23 = p333 * (hL(ij) + c2*hR(ij))
00788 if (hice(m) < h13) then
00789 hR(ij) = c3*hice(m) - c2*hL(ij)
00790 elseif (hice(m) > h23) then
00791 hL(ij) = c3*hice(m) - c2*hR(ij)
00792 endif
00793
00794
00795
00796 dhr = c1 / (hR(ij) - hL(ij))
00797 wk1 = c6 * aicen(i,j) * dhr
00798 wk2 = (hice(m) - hL(ij)) * dhr
00799 g0(ij) = wk1 * (p666 - wk2)
00800 g1(ij) = c2*dhr * wk1 * (wk2 - p5)
00801
00802 else
00803
00804 g0(ij) = c0
00805 g1(ij) = c0
00806 hL(ij) = c0
00807 hR(ij) = c0
00808
00809 endif
00810
00811 enddo
00812
00813 end subroutine fit_line
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838 subroutine add_new_ice (nx_block, ny_block, &
00839 icells, &
00840 indxi, indxj, &
00841 tmask, dt, &
00842 aicen, trcrn, &
00843 vicen, eicen, &
00844 aice0, aice, &
00845 frzmlt, frazil, &
00846 frz_onset, yday, &
00847 fresh, fsalt, &
00848 Tf, l_stop, &
00849 istop, jstop)
00850
00851
00852
00853 use ice_itd, only: hin_max, ilyr1, column_sum, &
00854 column_conservation_check
00855 use ice_state, only: nt_Tsfc, nt_iage, nt_FY, nt_aero, &
00856 tr_iage, tr_FY, tr_aero
00857 use ice_flux, only: update_ocn_f
00858
00859
00860
00861 integer (kind=int_kind), intent(in) ::
00862 nx_block, ny_block,
00863 icells
00864
00865 integer (kind=int_kind), dimension (nx_block*ny_block),
00866 intent(in) ::
00867 indxi, indxj
00868
00869 logical (kind=log_kind), dimension (nx_block,ny_block),
00870 intent(in) ::
00871 tmask
00872
00873 real (kind=dbl_kind), intent(in) ::
00874 dt
00875
00876 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
00877 aice ,
00878 frzmlt,
00879 Tf
00880
00881 real (kind=dbl_kind), dimension (nx_block,ny_block,ncat),
00882 intent(inout) ::
00883 aicen ,
00884 vicen
00885
00886 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat),
00887 intent(inout) ::
00888 trcrn
00889
00890
00891 real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr),
00892 intent(inout) ::
00893 eicen
00894
00895 real (kind=dbl_kind), dimension (nx_block,ny_block),
00896 intent(inout) ::
00897 aice0 ,
00898 frazil ,
00899 fresh ,
00900 fsalt
00901
00902 real (kind=dbl_kind), dimension (nx_block,ny_block),
00903 intent(inout), optional ::
00904 frz_onset
00905
00906 real (kind=dbl_kind), intent(in), optional ::
00907 yday
00908
00909 logical (kind=log_kind), intent(out) ::
00910 l_stop
00911
00912 integer (kind=int_kind), intent(out) ::
00913 istop, jstop
00914
00915
00916
00917 integer (kind=int_kind) ::
00918 i, j ,
00919 n ,
00920 k ,
00921 it
00922
00923 real (kind=dbl_kind), dimension (icells) ::
00924 ai0new ,
00925 vi0new ,
00926 hsurp ,
00927 vlyr
00928
00929 real (kind=dbl_kind), dimension (icells) ::
00930 vice_init, vice_final
00931
00932 real (kind=dbl_kind) ::
00933 fnew ,
00934 hi0new ,
00935 hi0max ,
00936 qi0(nilyr) ,
00937 qi0av ,
00938 vsurp ,
00939 area1 ,
00940 vice1 ,
00941 rnilyr ,
00942 dfresh ,
00943 dfsalt ,
00944 vtmp
00945
00946 integer (kind=int_kind) ::
00947 jcells, kcells ,
00948 ij, m
00949
00950 integer (kind=int_kind), dimension (icells) ::
00951 indxij2, indxij3 ,
00952 indxi2, indxj2 ,
00953 indxi3, indxj3
00954
00955 character (len=char_len) ::
00956 fieldid
00957
00958 l_stop = .false.
00959 istop = 0
00960 jstop = 0
00961
00962 jcells = 0
00963 kcells = 0
00964
00965 if (ncat > 1) then
00966 hi0max = hin_max(1)*0.9_dbl_kind
00967 else
00968 hi0max = bignum
00969 endif
00970
00971
00972 call column_sum (nx_block, ny_block, &
00973 icells, indxi, indxj, &
00974 ncat, &
00975 vicen, vice_init)
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986 rnilyr = real(nilyr,kind=dbl_kind)
00987 qi0av = c0
00988 do k = 1, nilyr
00989 qi0(k) = -rhoi*Lfresh
00990 qi0av = qi0av + qi0(k)
00991 enddo
00992 qi0av = qi0av/rnilyr
00993
00994
00995
00996
00997
00998
00999
01000
01001 do ij = 1, icells
01002 i = indxi(ij)
01003 j = indxj(ij)
01004
01005 fnew = max (frzmlt(i,j), c0)
01006 vi0new(ij) = -fnew*dt / qi0av
01007
01008
01009 vice_init(ij) = vice_init(ij) + vi0new(ij)
01010
01011
01012 frazil(i,j) = vi0new(ij)
01013
01014 if (present(frz_onset) .and. present(yday)) then
01015 if (frazil(i,j) > puny .and. frz_onset(i,j) < puny) &
01016 frz_onset(i,j) = yday
01017 endif
01018
01019
01020
01021
01022
01023
01024
01025
01026 if (update_ocn_f) then
01027 dfresh = -rhoi*vi0new(ij)/dt
01028 dfsalt = ice_ref_salinity*p001*dfresh
01029 fresh(i,j) = fresh(i,j) + dfresh
01030 fsalt(i,j) = fsalt(i,j) + dfsalt
01031 endif
01032
01033
01034
01035
01036
01037 hsurp(ij) = c0
01038 ai0new(ij) = c0
01039
01040 if (vi0new(ij) > c0) then
01041
01042
01043
01044 if (aice0(i,j) > puny) then
01045 hi0new = max(vi0new(ij)/aice0(i,j), hfrazilmin)
01046 if (hi0new > hi0max .and. aice0(i,j)+puny < c1) then
01047
01048 hi0new = hi0max
01049 ai0new(ij) = aice0(i,j)
01050 vsurp = vi0new(ij) - ai0new(ij)*hi0new
01051 hsurp(ij) = vsurp / aice(i,j)
01052 vi0new(ij) = ai0new(ij)*hi0new
01053 else
01054
01055 ai0new(ij) = vi0new(ij)/hi0new
01056 endif
01057 else
01058 hsurp(ij) = vi0new(ij)/aice(i,j)
01059 vi0new(ij) = c0
01060 endif
01061 endif
01062
01063
01064
01065
01066
01067 i = indxi(ij)
01068 j = indxj(ij)
01069
01070 if (vi0new(ij) > c0) then
01071 jcells = jcells + 1
01072 indxi2(jcells) = i
01073 indxj2(jcells) = j
01074 indxij2(jcells) = ij
01075 endif
01076
01077 if (hsurp(ij) > c0) then
01078 kcells = kcells + 1
01079 indxi3(kcells) = i
01080 indxj3(kcells) = j
01081 indxij3(kcells) = ij
01082 endif
01083
01084 enddo
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094 do n = 1, ncat
01095
01096
01097
01098
01099 do ij = 1, kcells
01100 i = indxi3(ij)
01101 j = indxj3(ij)
01102 m = indxij3(ij)
01103
01104 vsurp = hsurp(m) * aicen(i,j,n)
01105
01106
01107 vtmp = vicen(i,j,n) + vsurp
01108 if (vtmp > puny) then
01109
01110 if (tr_iage) trcrn(i,j,nt_iage,n) &
01111 = (trcrn(i,j,nt_iage,n)*vicen(i,j,n) + dt*vsurp) &
01112 / vtmp
01113
01114 if (tr_aero) then
01115 do it=1,n_aero
01116 trcrn(i,j,nt_aero+2+4*(it-1),n) &
01117 = trcrn(i,j,nt_aero+2+4*(it-1),n)*vicen(i,j,n) &
01118 / vtmp
01119 trcrn(i,j,nt_aero+3+4*(it-1),n) &
01120 = trcrn(i,j,nt_aero+3+4*(it-1),n)*vicen(i,j,n) &
01121 / vtmp
01122 enddo
01123 endif
01124
01125 endif
01126
01127
01128 vicen(i,j,n) = vicen(i,j,n) + vsurp
01129 vlyr(m) = vsurp/rnilyr
01130
01131 enddo
01132
01133 do k = 1, nilyr
01134
01135
01136
01137 do ij = 1, kcells
01138 i = indxi3(ij)
01139 j = indxj3(ij)
01140 m = indxij3(ij)
01141
01142 eicen(i,j,ilyr1(n)+k-1) = &
01143 eicen(i,j,ilyr1(n)+k-1) + qi0(k)*vlyr(m)
01144 enddo
01145 enddo
01146
01147 enddo
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157 do ij = 1, jcells
01158 i = indxi2(ij)
01159 j = indxj2(ij)
01160 m = indxij2(ij)
01161
01162 area1 = aicen(i,j,1)
01163 vice1 = vicen(i,j,1)
01164 aicen(i,j,1) = aicen(i,j,1) + ai0new(m)
01165 aice0(i,j) = aice0(i,j) - ai0new(m)
01166 vicen(i,j,1) = vicen(i,j,1) + vi0new(m)
01167 trcrn(i,j,nt_Tsfc,1) = (Tf(i,j)*ai0new(m) + trcrn(i,j,nt_Tsfc,1)*area1) &
01168 / aicen(i,j,1)
01169 trcrn(i,j,nt_Tsfc,1) = min (trcrn(i,j,nt_Tsfc,1), c0)
01170
01171 if (vicen(i,j,1) > puny) then
01172
01173 if (tr_iage) trcrn(i,j,nt_iage,1) = &
01174 (trcrn(i,j,nt_iage,1)*vice1 + dt*vi0new(m))/vicen(i,j,1)
01175
01176 if (tr_aero) then
01177 do it=1,n_aero
01178 trcrn(i,j,nt_aero+2+4*(it-1),1) = &
01179 trcrn(i,j,nt_aero+2+4*(it-1),1)*vice1/vicen(i,j,1)
01180 trcrn(i,j,nt_aero+3+4*(it-1),1) = &
01181 trcrn(i,j,nt_aero+3+4*(it-1),1)*vice1/vicen(i,j,1)
01182 enddo
01183 endif
01184
01185 endif
01186
01187 if (tr_FY) &
01188 trcrn(i,j,nt_FY,1) = (ai0new(m) + trcrn(i,j,nt_FY,1)*area1) &
01189 / aicen(i,j,1)
01190
01191 vlyr(m) = vi0new(m) / rnilyr
01192 enddo
01193
01194 do k = 1, nilyr
01195
01196
01197
01198 do ij = 1, jcells
01199 i = indxi2(ij)
01200 j = indxj2(ij)
01201 m = indxij2(ij)
01202 eicen(i,j,k) = eicen(i,j,k) + qi0(k)*vlyr(m)
01203 enddo
01204 enddo
01205
01206 call column_sum (nx_block, ny_block, &
01207 icells, indxi, indxj, &
01208 ncat, &
01209 vicen, vice_final)
01210
01211 fieldid = 'vice, add_new_ice'
01212 call column_conservation_check (nx_block, ny_block, &
01213 icells, indxi, indxj, &
01214 fieldid, &
01215 vice_init, vice_final, &
01216 puny, l_stop, &
01217 istop, jstop)
01218 if (l_stop) return
01219
01220 end subroutine add_new_ice
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239 subroutine lateral_melt (nx_block, ny_block, &
01240 ilo, ihi, jlo, jhi, &
01241 dt, &
01242 fresh, fsalt, &
01243 fhocn, fsoot, &
01244 rside, meltl, &
01245 aicen, vicen, &
01246 vsnon, eicen, &
01247 esnon, trcrn)
01248
01249
01250
01251 use ice_itd, only: ilyr1, slyr1
01252 use ice_state, only: nt_aero, tr_aero
01253
01254
01255
01256 integer (kind=int_kind), intent(in) ::
01257 nx_block, ny_block,
01258 ilo,ihi,jlo,jhi
01259
01260 real (kind=dbl_kind), intent(in) ::
01261 dt
01262
01263 real (kind=dbl_kind), dimension (nx_block,ny_block,ncat),
01264 intent(inout) ::
01265 aicen ,
01266 vicen ,
01267 vsnon
01268
01269 real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr),
01270 intent(inout) ::
01271 eicen
01272
01273 real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr),
01274 intent(inout) ::
01275 esnon
01276
01277 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat),
01278 intent(in) ::
01279 trcrn
01280
01281 real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) ::
01282 rside
01283
01284 real (kind=dbl_kind), dimension(nx_block,ny_block),
01285 intent(inout) ::
01286 fresh ,
01287 fsalt ,
01288 fhocn ,
01289 meltl
01290
01291 real (kind=dbl_kind), dimension(nx_block,ny_block,n_aeromx),
01292 intent(inout) ::
01293 fsoot
01294
01295
01296
01297 integer (kind=int_kind) ::
01298 i, j ,
01299 n ,
01300 k ,
01301 ij ,
01302 icells ,
01303 it
01304
01305 integer (kind=int_kind), dimension(nx_block*ny_block) ::
01306 indxi, indxj
01307
01308 real (kind=dbl_kind) ::
01309 dfhocn ,
01310 dfresh ,
01311 dfsalt
01312
01313 do n = 1, ncat
01314
01315
01316
01317
01318
01319 icells = 0
01320 do j = jlo, jhi
01321 do i = ilo, ihi
01322 if (rside(i,j) > c0) then
01323 icells = icells + 1
01324 indxi(icells) = i
01325 indxj(icells) = j
01326 endif
01327 enddo
01328 enddo
01329
01330
01331
01332
01333
01334 if (tr_aero) then
01335 do k=1,n_aero
01336 do ij = 1, icells
01337 i = indxi(ij)
01338 j = indxj(ij)
01339 fsoot(i,j,k) = fsoot(i,j,k) &
01340 + (vsnon(i,j,n) &
01341 *(trcrn(i,j,nt_aero +4*(k-1),n) &
01342 +trcrn(i,j,nt_aero+1+4*(k-1),n)) &
01343 + vicen(i,j,n) &
01344 *(trcrn(i,j,nt_aero+2+4*(k-1),n) &
01345 +trcrn(i,j,nt_aero+3+4*(k-1),n))) &
01346 * rside(i,j) / dt
01347 enddo
01348 enddo
01349 endif
01350
01351
01352
01353
01354 do ij = 1, icells
01355 i = indxi(ij)
01356 j = indxj(ij)
01357
01358
01359
01360
01361 dfresh = (rhos*vsnon(i,j,n) + rhoi*vicen(i,j,n)) &
01362 * rside(i,j) / dt
01363 dfsalt = rhoi*vicen(i,j,n)*ice_ref_salinity*p001 &
01364 * rside(i,j) / dt
01365
01366 fresh(i,j) = fresh(i,j) + dfresh
01367 fsalt(i,j) = fsalt(i,j) + dfsalt
01368
01369
01370 meltl(i,j) = meltl(i,j) + vicen(i,j,n)*rside(i,j)
01371
01372
01373 aicen(i,j,n) = aicen(i,j,n) * (c1 - rside(i,j))
01374 vicen(i,j,n) = vicen(i,j,n) * (c1 - rside(i,j))
01375 vsnon(i,j,n) = vsnon(i,j,n) * (c1 - rside(i,j))
01376
01377 enddo
01378
01379 do k = 1, nilyr
01380
01381
01382
01383 do ij = 1, icells
01384 i = indxi(ij)
01385 j = indxj(ij)
01386
01387
01388
01389 dfhocn = eicen(i,j,ilyr1(n)+k-1)*rside(i,j) / dt
01390 fhocn(i,j) = fhocn(i,j) + dfhocn
01391
01392
01393 eicen(i,j,ilyr1(n)+k-1) = eicen(i,j,ilyr1(n)+k-1) &
01394 * (c1 - rside(i,j))
01395 enddo
01396 enddo
01397
01398 do k = 1, nslyr
01399
01400
01401
01402 do ij = 1, icells
01403 i = indxi(ij)
01404 j = indxj(ij)
01405
01406
01407
01408 dfhocn = esnon(i,j,slyr1(n)+k-1)*rside(i,j) / dt
01409 fhocn(i,j) = fhocn(i,j) + dfhocn
01410
01411
01412 esnon(i,j,slyr1(n)+k-1) = esnon(i,j,slyr1(n)+k-1) &
01413 * (c1 - rside(i,j))
01414 enddo
01415 enddo
01416
01417 enddo
01418
01419 end subroutine lateral_melt
01420
01421
01422
01423 end module ice_therm_itd
01424
01425