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
00030
00031
00032
00033
00034
00035 module ice_itd
00036
00037
00038
00039 use ice_kinds_mod
00040 use ice_communicate, only: my_task, master_task
00041 use ice_domain_size
00042 use ice_constants
00043 use ice_fileunits
00044 use ice_exit
00045
00046
00047
00048 implicit none
00049 save
00050
00051 integer (kind=int_kind) ::
00052 kitd ,
00053
00054
00055 kcatbound , &
00056
00057
00058 ilyr1 (ncat), &
00059 ilyrn (ncat), &
00060 slyr1 (ncat), &
00061 slyrn (ncat)
00062
00063 real (kind=dbl_kind), parameter ::
00064 hi_min = p01
00065
00066 real (kind=dbl_kind) ::
00067 hin_max(0:ncat)
00068
00069 character (len=35) :: c_hi_range(ncat)
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088 contains
00089
00090
00091
00092
00093
00094
00095
00096
00097 subroutine init_itd
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114 integer (kind=int_kind) ::
00115 n
00116
00117 real (kind=dbl_kind) ::
00118 cc1, cc2, cc3,
00119 x1 ,
00120 rn ,
00121 rncat ,
00122 d1 ,
00123 d2
00124
00125 real (kind=dbl_kind), dimension(5) :: wmo5
00126 real (kind=dbl_kind), dimension(6) :: wmo6
00127 real (kind=dbl_kind), dimension(7) :: wmo7
00128
00129 character(len=8) :: c_hinmax1,c_hinmax2
00130 character(len=2) :: c_nc
00131
00132 rncat = real(ncat, kind=dbl_kind)
00133 d1 = 3.0_dbl_kind / rncat
00134 d2 = 0.5_dbl_kind / rncat
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162 if (kcatbound == 0) then
00163
00164 if (kitd == 1) then
00165
00166 cc1 = c3/rncat
00167 cc2 = c15*cc1
00168 cc3 = c3
00169
00170 hin_max(0) = c0
00171 else
00172
00173 cc1 = max(1.1_dbl_kind/rncat,c1*hi_min)
00174 cc2 = c25*cc1
00175 cc3 = 2.25_dbl_kind
00176
00177
00178
00179 hin_max(0) = hi_min
00180 endif
00181
00182 do n = 1, ncat
00183 x1 = real(n-1,kind=dbl_kind) / rncat
00184 hin_max(n) = hin_max(n-1) &
00185 + cc1 + cc2*(c1 + tanh(cc3*(x1-c1)))
00186 enddo
00187
00188 elseif (kcatbound == 1) then
00189
00190 hin_max(0) = c0
00191 do n = 1, ncat
00192 rn = real(n, kind=dbl_kind)
00193 hin_max(n) = rn * (d1 + (rn-c1)*d2)
00194 enddo
00195
00196 elseif (kcatbound == 2) then
00197
00198 if (ncat == 5) then
00199
00200 data wmo5 / 0.30_dbl_kind, 0.70_dbl_kind, &
00201 1.20_dbl_kind, 2.00_dbl_kind, &
00202 999._dbl_kind /
00203 hin_max(0) = c0
00204 do n = 1, ncat
00205 hin_max(n) = wmo5(n)
00206 enddo
00207 elseif (ncat == 6) then
00208
00209 data wmo6 / 0.15_dbl_kind, &
00210 0.30_dbl_kind, 0.70_dbl_kind, &
00211 1.20_dbl_kind, 2.00_dbl_kind, &
00212 999._dbl_kind /
00213 hin_max(0) = c0
00214 do n = 1, ncat
00215 hin_max(n) = wmo6(n)
00216 enddo
00217 elseif (ncat == 7) then
00218
00219 data wmo7 / 0.10_dbl_kind, 0.15_dbl_kind, &
00220 0.30_dbl_kind, 0.70_dbl_kind, &
00221 1.20_dbl_kind, 2.00_dbl_kind, &
00222 999._dbl_kind /
00223 hin_max(0) = c0
00224 do n = 1, ncat
00225 hin_max(n) = wmo7(n)
00226 enddo
00227 else
00228 write (nu_diag,*) 'kcatbound=3 (WMO) must have ncat=5, 6 or 7'
00229 stop
00230 endif
00231
00232 endif
00233
00234 if (my_task == master_task) then
00235 write (nu_diag,*) ' '
00236 write (nu_diag,*) 'hin_max(n-1) < Cat n < hin_max(n)'
00237 do n = 1, ncat
00238 write (nu_diag,*) hin_max(n-1),' < Cat ',n, ' < ',hin_max(n)
00239
00240 write (c_nc, '(i2)') n
00241
00242
00243 write (c_hinmax1, '(f5.3)') hin_max(n-1)
00244 write (c_hinmax2, '(f5.3)') hin_max(n)
00245
00246
00247 c_hi_range(n)=c_hinmax1//'m < hi Cat '//c_nc//' < '// &
00248 c_hinmax2//'m'
00249 enddo
00250 write (nu_diag,*) ' '
00251 endif
00252
00253
00254
00255
00256 ilyr1(1) = 1
00257 ilyrn(1) = nilyr
00258 do n = 2, ncat
00259 ilyr1(n) = ilyrn(n-1) + 1
00260 ilyrn(n) = ilyrn(n-1) + nilyr
00261 enddo
00262
00263 slyr1(1) = 1
00264 slyrn(1) = nslyr
00265 do n = 2, ncat
00266 slyr1(n) = slyrn(n-1) + 1
00267 slyrn(n) = slyrn(n-1) + nslyr
00268 enddo
00269
00270 end subroutine init_itd
00271
00272
00273
00274
00275
00276
00277
00278
00279 subroutine aggregate (nx_block, ny_block, &
00280 aicen, trcrn, &
00281 vicen, vsnon, &
00282 eicen, esnon, &
00283 aice, trcr, &
00284 vice, vsno, &
00285 eice, esno, &
00286 aice0, &
00287 tmask, trcr_depend)
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302 integer (kind=int_kind), intent(in) ::
00303 nx_block, ny_block
00304
00305 real (kind=dbl_kind), dimension (nx_block,ny_block,ncat),
00306 intent(in) ::
00307 aicen ,
00308 vicen ,
00309 vsnon
00310
00311 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat),
00312 intent(in) ::
00313 trcrn
00314
00315 real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr),
00316 intent(in) ::
00317 eicen
00318
00319 real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr),
00320 intent(in) ::
00321 esnon
00322
00323 logical (kind=log_kind), dimension (nx_block,ny_block),
00324 intent(in) ::
00325 tmask
00326
00327 integer (kind=int_kind), dimension (ntrcr), intent(in) ::
00328 trcr_depend
00329
00330 real (kind=dbl_kind), dimension (nx_block,ny_block),
00331 intent(out) ::
00332 aice ,
00333 vice ,
00334 vsno ,
00335 eice ,
00336 esno ,
00337 aice0
00338
00339 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr),
00340 intent(out) ::
00341 trcr
00342
00343
00344
00345 integer (kind=int_kind) ::
00346 icells
00347
00348 integer (kind=int_kind), dimension (nx_block*ny_block) ::
00349 indxi,
00350 indxj
00351
00352 integer (kind=int_kind) ::
00353 i, j, k, n, it,
00354 ij
00355
00356 real (kind=dbl_kind), dimension (:,:), allocatable ::
00357 atrcr
00358
00359
00360
00361
00362
00363 icells = 0
00364 do j = 1, ny_block
00365 do i = 1, nx_block
00366 if (tmask(i,j)) then
00367 icells = icells + 1
00368 indxi(icells) = i
00369 indxj(icells) = j
00370 endif
00371
00372 aice0(i,j) = c1
00373 aice (i,j) = c0
00374 vice (i,j) = c0
00375 vsno (i,j) = c0
00376 eice (i,j) = c0
00377 esno (i,j) = c0
00378 enddo
00379 enddo
00380
00381
00382 allocate (atrcr(icells,ntrcr))
00383
00384
00385
00386
00387
00388 atrcr(:,:) = c0
00389
00390 do n = 1, ncat
00391
00392
00393
00394
00395 do ij = 1, icells
00396 i = indxi(ij)
00397 j = indxj(ij)
00398 aice(i,j) = aice(i,j) + aicen(i,j,n)
00399 vice(i,j) = vice(i,j) + vicen(i,j,n)
00400 vsno(i,j) = vsno(i,j) + vsnon(i,j,n)
00401 enddo
00402
00403 do it = 1, ntrcr
00404 if (trcr_depend(it) == 0) then
00405
00406
00407
00408 do ij = 1, icells
00409 i = indxi(ij)
00410 j = indxj(ij)
00411 atrcr(ij,it) = atrcr(ij,it) &
00412 + trcrn(i,j,it,n)*aicen(i,j,n)
00413 enddo
00414
00415 elseif (trcr_depend(it) == 1) then
00416
00417
00418
00419 do ij = 1, icells
00420 i = indxi(ij)
00421 j = indxj(ij)
00422 atrcr(ij,it) = atrcr(ij,it) &
00423 + trcrn(i,j,it,n)*vicen(i,j,n)
00424 enddo
00425
00426 elseif (trcr_depend(it) ==2) then
00427
00428
00429
00430
00431 do ij = 1, icells
00432 i = indxi(ij)
00433 j = indxj(ij)
00434 atrcr(ij,it) = atrcr(ij,it) &
00435 + trcrn(i,j,it,n)*vsnon(i,j,n)
00436 enddo
00437
00438 endif
00439 enddo
00440
00441 do k = 1, nilyr
00442
00443
00444
00445 do ij = 1, icells
00446 i = indxi(ij)
00447 j = indxj(ij)
00448 eice(i,j) = eice(i,j) + eicen(i,j,ilyr1(n)+k-1)
00449 enddo
00450 enddo
00451
00452 do k = 1, nslyr
00453
00454
00455
00456 do ij = 1, icells
00457 i = indxi(ij)
00458 j = indxj(ij)
00459 esno(i,j) = esno(i,j) + esnon(i,j,slyr1(n)+k-1)
00460 enddo
00461 enddo
00462
00463 enddo
00464
00465
00466
00467 do ij = 1, icells
00468 i = indxi(ij)
00469 j = indxj(ij)
00470 aice0(i,j) = max (c1 - aice(i,j), c0)
00471 enddo
00472
00473
00474
00475 call compute_tracers (nx_block, ny_block, &
00476 icells, indxi, indxj, &
00477 trcr_depend, &
00478 atrcr(:,:), aice(:,:), &
00479 vice (:,:), vsno(:,:), &
00480 trcr(:,:,:))
00481
00482 deallocate (atrcr)
00483
00484 end subroutine aggregate
00485
00486
00487
00488
00489
00490
00491
00492
00493 subroutine aggregate_area (nx_block, ny_block, &
00494 aicen, aice, aice0)
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510 integer (kind=int_kind), intent(in) ::
00511 nx_block, ny_block
00512
00513 real (kind=dbl_kind), dimension (:,:,:), intent(in) ::
00514 aicen
00515
00516 real (kind=dbl_kind), dimension (:,:), intent(inout) ::
00517 aice,
00518 aice0
00519
00520
00521
00522 integer (kind=int_kind) :: i, j, n
00523
00524
00525
00526
00527
00528 aice(:,:) = c0
00529
00530 do n = 1, ncat
00531 do j = 1, ny_block
00532 do i = 1, nx_block
00533 aice(i,j) = aice(i,j) + aicen(i,j,n)
00534 enddo
00535 enddo
00536 enddo
00537
00538 do j = 1, ny_block
00539 do i = 1, nx_block
00540
00541
00542 aice0(i,j) = max (c1 - aice(i,j), c0)
00543
00544 enddo
00545 enddo
00546
00547 end subroutine aggregate_area
00548
00549
00550
00551
00552
00553
00554
00555
00556 subroutine rebin (nx_block, ny_block, &
00557 icells, indxi, indxj, &
00558 trcr_depend, &
00559 aicen, trcrn, &
00560 vicen, vsnon, &
00561 eicen, esnon, &
00562 l_stop, &
00563 istop, jstop)
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577 integer (kind=int_kind), intent(in) ::
00578 nx_block, ny_block,
00579 icells
00580
00581 integer (kind=int_kind), dimension (nx_block*ny_block),
00582 intent(in) ::
00583 indxi, indxj
00584
00585 integer (kind=int_kind), dimension (ntrcr), intent(in) ::
00586 trcr_depend
00587
00588 real (kind=dbl_kind), dimension (nx_block,ny_block,ncat),
00589 intent(inout) ::
00590 aicen ,
00591 vicen ,
00592 vsnon
00593
00594 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat),
00595 intent(inout) ::
00596 trcrn
00597
00598 real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr),
00599 intent(inout) ::
00600 eicen
00601
00602 real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr),
00603 intent(inout) ::
00604 esnon
00605
00606 logical (kind=log_kind), intent(out) ::
00607 l_stop
00608
00609 integer (kind=int_kind), intent(out) ::
00610 istop, jstop
00611
00612
00613
00614 integer (kind=int_kind) ::
00615 i,j ,
00616 n ,
00617 ij
00618
00619 logical (kind=log_kind) ::
00620 shiftflag
00621
00622 integer (kind=int_kind), dimension (icells,ncat) ::
00623 donor
00624
00625 real (kind=dbl_kind), dimension (icells,ncat) ::
00626 daice ,
00627 dvice ,
00628 hicen
00629
00630
00631
00632
00633
00634 l_stop = .false.
00635 istop = 0
00636 jstop = 0
00637
00638 do n = 1, ncat
00639 do ij = 1, icells
00640 i = indxi(ij)
00641 j = indxj(ij)
00642
00643 donor(ij,n) = 0
00644 daice(ij,n) = c0
00645 dvice(ij,n) = c0
00646
00647
00648
00649
00650 if (aicen(i,j,n) > puny) then
00651 hicen(ij,n) = vicen(i,j,n) / aicen(i,j,n)
00652 else
00653 hicen(ij,n) = c0
00654 endif
00655 enddo
00656 enddo
00657
00658
00659
00660
00661 do ij = 1, icells
00662 i = indxi(ij)
00663 j = indxj(ij)
00664
00665 if (aicen(i,j,1) > puny) then
00666 if (hicen(ij,1) <= hin_max(0) .and. hin_max(0) > c0 ) then
00667 aicen(i,j,1) = vicen(i,j,1) / hin_max(0)
00668 hicen(ij,1) = hin_max(0)
00669 endif
00670 endif
00671 enddo
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682 do n = 1, ncat-1
00683
00684
00685
00686
00687 shiftflag = .false.
00688 do ij = 1, icells
00689 i = indxi(ij)
00690 j = indxj(ij)
00691
00692 if (aicen(i,j,n) > puny .and. &
00693 hicen(ij,n) > hin_max(n)) then
00694 shiftflag = .true.
00695 donor(ij,n) = n
00696 daice(ij,n) = aicen(i,j,n)
00697 dvice(ij,n) = vicen(i,j,n)
00698 endif
00699 enddo
00700
00701 if (shiftflag) then
00702
00703
00704
00705
00706 call shift_ice (nx_block, ny_block, &
00707 indxi, indxj, &
00708 icells, trcr_depend, &
00709 aicen, trcrn, &
00710 vicen, vsnon, &
00711 eicen, esnon, &
00712 hicen, donor, &
00713 daice, dvice, &
00714 l_stop, &
00715 istop, jstop)
00716
00717 if (l_stop) return
00718
00719
00720
00721
00722
00723 do ij = 1, icells
00724 donor(ij,n) = 0
00725 daice(ij,n) = c0
00726 dvice(ij,n) = c0
00727 enddo
00728
00729 endif
00730
00731 enddo
00732
00733
00734
00735
00736
00737 do n = ncat-1, 1, -1
00738
00739
00740
00741
00742 shiftflag = .false.
00743 do ij = 1, icells
00744 i = indxi(ij)
00745 j = indxj(ij)
00746
00747 if (aicen(i,j,n+1) > puny .and. &
00748 hicen(ij,n+1) <= hin_max(n)) then
00749 shiftflag = .true.
00750 donor(ij,n) = n+1
00751 daice(ij,n) = aicen(i,j,n+1)
00752 dvice(ij,n) = vicen(i,j,n+1)
00753 endif
00754 enddo
00755
00756 if (shiftflag) then
00757
00758
00759
00760
00761 call shift_ice (nx_block, ny_block, &
00762 indxi, indxj, &
00763 icells, trcr_depend, &
00764 aicen, trcrn, &
00765 vicen, vsnon, &
00766 eicen, esnon, &
00767 hicen, donor, &
00768 daice, dvice, &
00769 l_stop, &
00770 istop, jstop)
00771
00772 if (l_stop) return
00773
00774
00775
00776
00777
00778 do ij = 1, icells
00779 donor(ij,n) = 0
00780 daice(ij,n) = c0
00781 dvice(ij,n) = c0
00782 enddo
00783
00784 endif
00785
00786 enddo
00787
00788
00789 end subroutine rebin
00790
00791
00792
00793
00794
00795
00796
00797
00798 subroutine reduce_area (nx_block, ny_block, &
00799 ilo, ihi, jlo, jhi, &
00800 tmask, &
00801 aicen, vicen, &
00802 aicen_init,vicen_init)
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822 integer (kind=int_kind), intent(in) ::
00823 nx_block, ny_block,
00824 ilo,ihi,jlo,jhi
00825
00826 logical (kind=log_kind), dimension (nx_block,ny_block),
00827 intent(in) ::
00828 tmask
00829
00830 real (kind=dbl_kind), dimension (nx_block,ny_block),
00831 intent(inout) ::
00832 aicen ,
00833 vicen
00834
00835 real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) ::
00836 aicen_init,
00837 vicen_init
00838
00839
00840
00841 integer (kind=int_kind) ::
00842 i, j
00843
00844 real (kind=dbl_kind) ::
00845 hi0 ,
00846 hi1 ,
00847 dhi
00848
00849 do j = jlo, jhi
00850 do i = ilo, ihi
00851 if (tmask(i,j)) then
00852
00853 hi0 = c0
00854 if (aicen_init(i,j) > c0) &
00855 hi0 = vicen_init(i,j) / aicen_init(i,j)
00856
00857 hi1 = c0
00858 if (aicen(i,j) > c0) &
00859 hi1 = vicen(i,j) / aicen(i,j)
00860
00861
00862 if (hi1 <= hin_max(0) .and. hin_max(0) > c0 ) then
00863 aicen(i,j) = vicen(i,j) / hin_max(0)
00864 hi1 = hin_max(0)
00865 endif
00866
00867 if (aicen(i,j) > c0) then
00868 dhi = hi1 - hi0
00869 if (dhi < c0) then
00870 hi1 = vicen(i,j) / aicen(i,j)
00871 aicen(i,j) = c2 * vicen(i,j) / (hi1 + hi0)
00872 endif
00873 endif
00874
00875 endif
00876 enddo
00877 enddo
00878
00879 end subroutine reduce_area
00880
00881
00882
00883
00884
00885
00886
00887
00888 subroutine shift_ice (nx_block, ny_block, &
00889 indxi, indxj, &
00890 icells, trcr_depend, &
00891 aicen, trcrn, &
00892 vicen, vsnon, &
00893 eicen, esnon, &
00894 hicen, donor, &
00895 daice, dvice, &
00896 l_stop, &
00897 istop, jstop)
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913 integer (kind=int_kind), intent(in) ::
00914 nx_block, ny_block,
00915 icells
00916
00917 integer (kind=int_kind), dimension (nx_block*ny_block),
00918 intent(in) ::
00919 indxi ,
00920 indxj
00921
00922 integer (kind=int_kind), dimension (ntrcr), intent(in) ::
00923 trcr_depend
00924
00925 real (kind=dbl_kind), dimension (nx_block,ny_block,ncat),
00926 intent(inout) ::
00927 aicen ,
00928 vicen ,
00929 vsnon
00930
00931 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat),
00932 intent(inout) ::
00933 trcrn
00934
00935 real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr),
00936 intent(inout) ::
00937 eicen
00938
00939 real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr),
00940 intent(inout) ::
00941 esnon
00942
00943
00944
00945 integer (kind=int_kind), dimension(icells,ncat),
00946 intent(in) ::
00947 donor
00948
00949 real (kind=dbl_kind), dimension(icells,ncat),
00950 intent(inout) ::
00951 daice ,
00952 dvice ,
00953 hicen
00954
00955 logical (kind=log_kind), intent(out) ::
00956 l_stop
00957
00958 integer (kind=int_kind), intent(out) ::
00959 istop, jstop
00960
00961
00962
00963 integer (kind=int_kind) ::
00964 i, j, m ,
00965 n ,
00966 nr ,
00967 nd ,
00968 k ,
00969 it ,
00970 ilo,ihi,jlo,jhi
00971
00972 real (kind=dbl_kind), dimension(icells,ntrcr,ncat) ::
00973 atrcrn
00974
00975 real (kind=dbl_kind) ::
00976 dvsnow ,
00977 desnow ,
00978 deice ,
00979 datrcr
00980
00981 integer (kind=int_kind), dimension (icells) ::
00982 indxii ,
00983 indxjj ,
00984 indxij
00985
00986 integer (kind=int_kind) ::
00987 ishift ,
00988 ij
00989
00990 logical (kind=log_kind) ::
00991 daice_negative ,
00992 dvice_negative ,
00993 daice_greater_aicen,
00994 dvice_greater_vicen
00995
00996 real (kind=dbl_kind), dimension (nx_block,ny_block) ::
00997 worka,
00998 workb
00999
01000
01001
01002
01003
01004 l_stop = .false.
01005 istop = 0
01006 jstop = 0
01007
01008 worka(:,:) = c0
01009 workb(:,:) = c0
01010
01011
01012
01013
01014
01015 do n = 1, ncat
01016 do it = 1, ntrcr
01017 if (trcr_depend(it) == 0) then
01018 do ij = 1, icells
01019 i = indxi(ij)
01020 j = indxj(ij)
01021 atrcrn(ij,it,n) = aicen(i,j,n)*trcrn(i,j,it,n)
01022 enddo
01023 elseif (trcr_depend(it) ==1) then
01024 do ij = 1, icells
01025 i = indxi(ij)
01026 j = indxj(ij)
01027 atrcrn(ij,it,n) = vicen(i,j,n)*trcrn(i,j,it,n)
01028 enddo
01029 elseif (trcr_depend(it) ==2) then
01030 do ij = 1, icells
01031 i = indxi(ij)
01032 j = indxj(ij)
01033 atrcrn(ij,it,n) = vsnon(i,j,n)*trcrn(i,j,it,n)
01034 enddo
01035 endif
01036 enddo
01037 enddo
01038
01039
01040
01041
01042
01043 do n = 1, ncat-1
01044
01045 daice_negative = .false.
01046 dvice_negative = .false.
01047 daice_greater_aicen = .false.
01048 dvice_greater_vicen = .false.
01049
01050
01051 do ij = 1, icells
01052 i = indxi(ij)
01053 j = indxj(ij)
01054
01055 if (donor(ij,n) > 0) then
01056 nd = donor(ij,n)
01057
01058 if (daice(ij,n) < c0) then
01059 if (daice(ij,n) > -puny*aicen(i,j,nd)) then
01060 daice(ij,n) = c0
01061 dvice(ij,n) = c0
01062 else
01063 daice_negative = .true.
01064 endif
01065 endif
01066
01067 if (dvice(ij,n) < c0) then
01068 if (dvice(ij,n) > -puny*vicen(i,j,nd)) then
01069 daice(ij,n) = c0
01070 dvice(ij,n) = c0
01071 else
01072 dvice_negative = .true.
01073 endif
01074 endif
01075
01076 if (daice(ij,n) > aicen(i,j,nd)*(c1-puny)) then
01077 if (daice(ij,n) < aicen(i,j,nd)*(c1+puny)) then
01078 daice(ij,n) = aicen(i,j,nd)
01079 dvice(ij,n) = vicen(i,j,nd)
01080 else
01081 daice_greater_aicen = .true.
01082 endif
01083 endif
01084
01085 if (dvice(ij,n) > vicen(i,j,nd)*(c1-puny)) then
01086 if (dvice(ij,n) < vicen(i,j,nd)*(c1+puny)) then
01087 daice(ij,n) = aicen(i,j,nd)
01088 dvice(ij,n) = vicen(i,j,nd)
01089 else
01090 dvice_greater_vicen = .true.
01091 endif
01092 endif
01093
01094 endif
01095 enddo
01096
01097
01098
01099
01100
01101 if (daice_negative) then
01102 do ij = 1, icells
01103 i = indxi(ij)
01104 j = indxj(ij)
01105
01106 if (donor(ij,n) > 0 .and. &
01107 daice(ij,n) <= -puny*aicen(i,j,nd)) then
01108 write(nu_diag,*) ' '
01109 write(nu_diag,*) 'shift_ice: negative daice'
01110 write(nu_diag,*) 'i, j:', i, j
01111 write(nu_diag,*) 'boundary, donor cat:', n, nd
01112 write(nu_diag,*) 'daice =', daice(ij,n)
01113 write(nu_diag,*) 'dvice =', dvice(ij,n)
01114 l_stop = .true.
01115 istop = i
01116 jstop = j
01117 endif
01118 enddo
01119 endif
01120 if (l_stop) return
01121
01122 if (dvice_negative) then
01123 do ij = 1, icells
01124 i = indxi(ij)
01125 j = indxj(ij)
01126
01127 if (donor(ij,n) > 0 .and. &
01128 dvice(ij,n) <= -puny*vicen(i,j,nd)) then
01129 write(nu_diag,*) ' '
01130 write(nu_diag,*) 'shift_ice: negative dvice'
01131 write(nu_diag,*) 'i, j:', i, j
01132 write(nu_diag,*) 'boundary, donor cat:', n, nd
01133 write(nu_diag,*) 'daice =', daice(ij,n)
01134 write(nu_diag,*) 'dvice =', dvice(ij,n)
01135 l_stop = .true.
01136 istop = i
01137 jstop = j
01138 endif
01139 enddo
01140 endif
01141 if (l_stop) return
01142
01143 if (daice_greater_aicen) then
01144 do ij = 1, icells
01145 i = indxi(ij)
01146 j = indxj(ij)
01147
01148 if (donor(ij,n) > 0) then
01149 nd = donor(ij,n)
01150 if (daice(ij,n) >= aicen(i,j,nd)*(c1+puny)) then
01151 write(nu_diag,*) ' '
01152 write(nu_diag,*) 'shift_ice: daice > aicen'
01153 write(nu_diag,*) 'i, j:', i, j
01154 write(nu_diag,*) 'boundary, donor cat:', n, nd
01155 write(nu_diag,*) 'daice =', daice(ij,n)
01156 write(nu_diag,*) 'aicen =', aicen(i,j,nd)
01157 l_stop = .true.
01158 istop = i
01159 jstop = j
01160 endif
01161 endif
01162 enddo
01163 endif
01164 if (l_stop) return
01165
01166 if (dvice_greater_vicen) then
01167 do ij = 1, icells
01168 i = indxi(ij)
01169 j = indxj(ij)
01170
01171 if (donor(ij,n) > 0) then
01172 nd = donor(ij,n)
01173 if (dvice(ij,n) >= vicen(i,j,nd)*(c1+puny)) then
01174 write(nu_diag,*) ' '
01175 write(nu_diag,*) 'shift_ice: dvice > vicen'
01176 write(nu_diag,*) 'i, j:', i, j
01177 write(nu_diag,*) 'boundary, donor cat:', n, nd
01178 write(nu_diag,*) 'dvice =', dvice(ij,n)
01179 write(nu_diag,*) 'vicen =', vicen(i,j,nd)
01180 l_stop = .true.
01181 istop = i
01182 jstop = j
01183 endif
01184 endif
01185 enddo
01186 endif
01187 if (l_stop) return
01188
01189
01190
01191
01192
01193 ishift = 0
01194 do ij = 1, icells
01195 i = indxi(ij)
01196 j = indxj(ij)
01197
01198 if (daice(ij,n) > c0) then
01199 ishift = ishift + 1
01200 indxii(ishift) = i
01201 indxjj(ishift) = j
01202 indxij(ishift) = ij
01203 endif
01204 enddo
01205
01206
01207
01208
01209 do ij = 1, ishift
01210 i = indxii(ij)
01211 j = indxjj(ij)
01212 m = indxij(ij)
01213
01214 nd = donor(m,n)
01215 worka(i,j) = dvice(m,n) / vicen(i,j,nd)
01216 if (nd == n) then
01217 nr = nd+1
01218 else
01219 nr = n
01220 endif
01221
01222 aicen(i,j,nd) = aicen(i,j,nd) - daice(m,n)
01223 aicen(i,j,nr) = aicen(i,j,nr) + daice(m,n)
01224
01225 vicen(i,j,nd) = vicen(i,j,nd) - dvice(m,n)
01226 vicen(i,j,nr) = vicen(i,j,nr) + dvice(m,n)
01227
01228 dvsnow = vsnon(i,j,nd) * worka(i,j)
01229 vsnon(i,j,nd) = vsnon(i,j,nd) - dvsnow
01230 vsnon(i,j,nr) = vsnon(i,j,nr) + dvsnow
01231 workb(i,j) = dvsnow
01232
01233 enddo
01234
01235 do it = 1, ntrcr
01236
01237
01238
01239 do ij = 1, ishift
01240 i = indxii(ij)
01241 j = indxjj(ij)
01242 m = indxij(ij)
01243
01244 nd = donor(m,n)
01245 if (nd == n) then
01246 nr = nd+1
01247 else
01248 nr = n
01249 endif
01250
01251 if (trcr_depend(it) == 0) then
01252 datrcr = daice(m,n)*trcrn(i,j,it,nd)
01253 elseif (trcr_depend(it) == 1) then
01254 datrcr = dvice(m,n)*trcrn(i,j,it,nd)
01255 elseif (trcr_depend(it) == 2) then
01256 datrcr = workb(i,j) *trcrn(i,j,it,nd)
01257 endif
01258
01259 atrcrn(m,it,nd) = atrcrn(m,it,nd) - datrcr
01260 atrcrn(m,it,nr) = atrcrn(m,it,nr) + datrcr
01261 enddo
01262 enddo
01263
01264 do k = 1, nilyr
01265
01266
01267
01268 do ij = 1, ishift
01269 i = indxii(ij)
01270 j = indxjj(ij)
01271 m = indxij(ij)
01272
01273 nd = donor(m,n)
01274 if (nd == n) then
01275 nr = nd+1
01276 else
01277 nr = n
01278 endif
01279
01280 deice = eicen(i,j,ilyr1(nd)+k-1) * worka(i,j)
01281 eicen(i,j,ilyr1(nd)+k-1) = &
01282 eicen(i,j,ilyr1(nd)+k-1) - deice
01283 eicen(i,j,ilyr1(nr)+k-1) = &
01284 eicen(i,j,ilyr1(nr)+k-1) + deice
01285 enddo
01286 enddo
01287
01288 do k = 1, nslyr
01289
01290
01291
01292 do ij = 1, ishift
01293 i = indxii(ij)
01294 j = indxjj(ij)
01295 m = indxij(ij)
01296
01297 nd = donor(m,n)
01298 if (nd == n) then
01299 nr = nd+1
01300 else
01301 nr = n
01302 endif
01303
01304 desnow = esnon(i,j,slyr1(nd)+k-1) * worka(i,j)
01305 esnon(i,j,slyr1(nd)+k-1) = &
01306 esnon(i,j,slyr1(nd)+k-1) - desnow
01307 esnon(i,j,slyr1(nr)+k-1) = &
01308 esnon(i,j,slyr1(nr)+k-1) + desnow
01309 enddo
01310 enddo
01311
01312 enddo
01313
01314
01315
01316
01317
01318 do n = 1, ncat
01319
01320 do ij = 1, icells
01321 i = indxi(ij)
01322 j = indxj(ij)
01323
01324 if (aicen(i,j,n) > puny) then
01325 hicen(ij,n) = vicen (i,j,n) / aicen(i,j,n)
01326 else
01327 hicen(ij,n) = c0
01328 endif
01329 enddo
01330
01331 call compute_tracers (nx_block, ny_block, &
01332 icells, indxi, indxj, &
01333 trcr_depend, &
01334 atrcrn(:,:,n), aicen(:,:, n), &
01335 vicen (:,:, n), vsnon(:,:, n), &
01336 trcrn(:,:,:,n))
01337
01338 enddo
01339
01340 end subroutine shift_ice
01341
01342
01343
01344
01345
01346
01347
01348
01349 subroutine column_sum (nx_block, ny_block, &
01350 icells, indxi, indxj, &
01351 nsum, &
01352 xin, xout)
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366 integer (kind=int_kind), intent(in) ::
01367 nx_block, ny_block,
01368 nsum ,
01369 icells
01370
01371 integer (kind=int_kind), dimension (nx_block*ny_block),
01372 intent(in) ::
01373 indxi, indxj
01374
01375 real (kind=dbl_kind), dimension (nx_block,ny_block,nsum),
01376 intent(in) ::
01377 xin
01378
01379 real (kind=dbl_kind), dimension (icells), intent(out) ::
01380 xout
01381
01382
01383
01384 integer (kind=int_kind) ::
01385 i, j, ij ,
01386 n
01387
01388 do ij = 1, icells
01389 xout(ij) = c0
01390 enddo
01391
01392 do n = 1, nsum
01393 do ij = 1, icells
01394 i = indxi(ij)
01395 j = indxj(ij)
01396 xout(ij) = xout(ij) + xin(i,j,n)
01397 enddo
01398 enddo
01399
01400 end subroutine column_sum
01401
01402
01403
01404
01405
01406
01407
01408
01409 subroutine column_conservation_check (nx_block, ny_block, &
01410 icells, indxi, indxj, &
01411 fieldid, &
01412 x1, x2, &
01413 max_err, l_stop, &
01414 istop, jstop)
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429 integer (kind=int_kind), intent(in) ::
01430 nx_block, ny_block,
01431 icells
01432
01433 integer (kind=int_kind), dimension (nx_block*ny_block),
01434 intent(in) ::
01435 indxi, indxj
01436
01437 real (kind=dbl_kind), dimension(icells), intent(in) ::
01438 x1 ,
01439 x2
01440
01441 real (kind=dbl_kind), intent(in) ::
01442 max_err
01443
01444 character (len=char_len), intent(in) ::
01445 fieldid
01446
01447 logical (kind=log_kind), intent(inout) ::
01448 l_stop
01449
01450 integer (kind=int_kind), intent(inout) ::
01451 istop, jstop
01452
01453
01454
01455 integer (kind=int_kind) ::
01456 ij
01457
01458 do ij = 1, icells
01459 if (abs (x2(ij)-x1(ij)) > max_err) then
01460 l_stop = .true.
01461 istop = indxi(ij)
01462 jstop = indxj(ij)
01463
01464 write (nu_diag,*) ' '
01465 write (nu_diag,*) 'Conservation error: ', trim(fieldid)
01466 write (nu_diag,*) 'i, j =', istop, jstop
01467 write (nu_diag,*) 'Initial value =', x1(ij)
01468 write (nu_diag,*) 'Final value =', x2(ij)
01469 write (nu_diag,*) 'Difference =', x2(ij) - x1(ij)
01470 endif
01471 enddo
01472
01473 end subroutine column_conservation_check
01474
01475
01476
01477
01478
01479
01480
01481
01482 subroutine compute_tracers (nx_block, ny_block, &
01483 icells, indxi, indxj, &
01484 trcr_depend, &
01485 atrcrn, aicen, &
01486 vicen, vsnon, &
01487 trcrn)
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500 use ice_state, only: nt_Tsfc
01501
01502
01503
01504 integer (kind=int_kind), intent(in) ::
01505 nx_block, ny_block,
01506 icells
01507
01508 integer (kind=int_kind), dimension (nx_block*ny_block),
01509 intent(in) ::
01510 indxi, indxj
01511
01512 integer (kind=int_kind), dimension (ntrcr), intent(in) ::
01513 trcr_depend
01514
01515 real (kind=dbl_kind), dimension (icells,ntrcr),
01516 intent(in) ::
01517 atrcrn
01518
01519 real (kind=dbl_kind), dimension (nx_block,ny_block),
01520 intent(in) ::
01521 aicen ,
01522 vicen ,
01523 vsnon
01524
01525 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr),
01526 intent(out) ::
01527 trcrn
01528
01529
01530
01531 integer (kind=int_kind) ::
01532 i, j, it, ij
01533
01534
01535 trcrn(:,:,:) = c0
01536
01537
01538
01539
01540
01541 do it = 1, ntrcr
01542 if (it == nt_Tsfc) then
01543 do ij = 1, icells
01544 i = indxi(ij)
01545 j = indxj(ij)
01546 if (aicen(i,j) > puny) then
01547 trcrn(i,j,it) = atrcrn(ij,it) / aicen(i,j)
01548 else
01549 trcrn(i,j,it) = Tocnfrz
01550 endif
01551 enddo
01552
01553 elseif (trcr_depend(it) == 0) then
01554 do ij = 1, icells
01555 i = indxi(ij)
01556 j = indxj(ij)
01557 if (aicen(i,j) > puny) then
01558 trcrn(i,j,it) = atrcrn(ij,it) / aicen(i,j)
01559 else
01560 trcrn(i,j,it) = c0
01561 endif
01562 enddo
01563
01564 elseif (trcr_depend(it) == 1) then
01565 do ij = 1, icells
01566 i = indxi(ij)
01567 j = indxj(ij)
01568 if (vicen(i,j) > puny) then
01569 trcrn(i,j,it) = atrcrn(ij,it) / vicen(i,j)
01570 else
01571 trcrn(i,j,it) = c0
01572 endif
01573 enddo
01574
01575 elseif (trcr_depend(it) == 2) then
01576 do ij = 1, icells
01577 i = indxi(ij)
01578 j = indxj(ij)
01579 if (vsnon(i,j) > puny) then
01580 trcrn(i,j,it) = atrcrn(ij,it) / vsnon(i,j)
01581 else
01582 trcrn(i,j,it) = c0
01583 endif
01584 enddo
01585
01586 endif
01587 enddo
01588
01589 end subroutine compute_tracers
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599 subroutine cleanup_itd (nx_block, ny_block, &
01600 ilo, ihi, jlo, jhi, &
01601 dt, &
01602 aicen, trcrn, &
01603 vicen, vsnon, &
01604 eicen, esnon, &
01605 aice0, aice, &
01606 trcr_depend, fresh, &
01607 fsalt, fhocn, &
01608 fsoot, tr_aero, &
01609 heat_capacity, l_stop, &
01610 istop, jstop, &
01611 limit_aice_in)
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630 integer (kind=int_kind), intent(in) ::
01631 nx_block, ny_block,
01632 ilo,ihi,jlo,jhi
01633
01634 real (kind=dbl_kind), intent(in) ::
01635 dt
01636
01637 real (kind=dbl_kind), dimension (nx_block,ny_block,ncat),
01638 intent(inout) ::
01639 aicen ,
01640 vicen ,
01641 vsnon
01642
01643 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat),
01644 intent(inout) ::
01645 trcrn
01646
01647 real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr),
01648 intent(inout) ::
01649 eicen
01650
01651 real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr),
01652 intent(inout) ::
01653 esnon
01654
01655 real (kind=dbl_kind), dimension (nx_block,ny_block),
01656 intent(inout) ::
01657 aice ,
01658 aice0
01659
01660 integer (kind=int_kind), dimension(ntrcr), intent(in) ::
01661 trcr_depend
01662
01663 logical (kind=log_kind), intent(in) ::
01664 tr_aero,
01665 heat_capacity
01666
01667 logical (kind=log_kind), intent(out) ::
01668 l_stop
01669
01670 integer (kind=int_kind), intent(out) ::
01671 istop, jstop
01672
01673
01674 real (kind=dbl_kind), dimension (nx_block,ny_block),
01675 intent(inout), optional ::
01676 fresh ,
01677 fsalt ,
01678 fhocn
01679
01680 real (kind=dbl_kind), dimension (nx_block,ny_block,n_aeromx),
01681 intent(inout), optional ::
01682 fsoot
01683
01684 logical (kind=log_kind), intent(in), optional ::
01685 limit_aice_in
01686
01687
01688
01689
01690 integer (kind=int_kind) ::
01691 i, j ,
01692 n ,
01693 icells
01694
01695 integer (kind=int_kind), dimension (nx_block*ny_block) ::
01696 indxi, indxj
01697
01698 real (kind=dbl_kind), dimension (nx_block,ny_block) ::
01699 dfresh ,
01700 dfsalt ,
01701 dfhocn
01702
01703 real (kind=dbl_kind), dimension (nx_block,ny_block,n_aeromx) ::
01704 dfsoot
01705
01706 logical (kind=log_kind) ::
01707 limit_aice
01708
01709
01710
01711
01712
01713 if (present(limit_aice_in)) then
01714 limit_aice = limit_aice_in
01715 else
01716 limit_aice = .true.
01717 endif
01718
01719 l_stop = .false.
01720 istop = 0
01721 jstop = 0
01722
01723
01724
01725
01726
01727 call aggregate_area (nx_block, ny_block, &
01728 aicen(:,:,:), &
01729 aice, aice0)
01730
01731
01732 if (limit_aice) then
01733
01734 do j = jlo,jhi
01735 do i = ilo,ihi
01736 if (aice(i,j) > c1+puny .or. aice(i,j) < -puny) then
01737 l_stop = .true.
01738 istop = i
01739 jstop = j
01740 endif
01741 enddo
01742 enddo
01743
01744 if (l_stop) then
01745 i = istop
01746 j = jstop
01747 write(nu_diag,*) ' '
01748 write(nu_diag,*) 'aggregate ice area out of bounds'
01749 write(nu_diag,*) 'my_task, i, j, aice:', &
01750 my_task, i, j, aice(i,j)
01751 do n = 1, ncat
01752 write(nu_diag,*) 'n, aicen:', n, aicen(i,j,n)
01753 enddo
01754 return
01755 endif
01756 endif
01757
01758
01759
01760
01761
01762 icells = 0
01763 do j = jlo,jhi
01764 do i = ilo,ihi
01765 if (aice(i,j) > puny) then
01766 icells = icells + 1
01767 indxi(icells) = i
01768 indxj(icells) = j
01769 endif
01770 enddo
01771 enddo
01772
01773
01774
01775
01776
01777
01778
01779
01780 call rebin (nx_block, ny_block, &
01781 icells, indxi, indxj, &
01782 trcr_depend, &
01783 aicen(:,:,:), trcrn(:,:,:,:), &
01784 vicen(:,:,:), vsnon(:,:,:), &
01785 eicen(:,:,:), esnon(:,:,:), &
01786 l_stop, &
01787 istop, jstop)
01788
01789 if (l_stop) return
01790
01791
01792
01793
01794
01795 if (limit_aice) then
01796 call zap_small_areas (nx_block, ny_block, &
01797 ilo, ihi, jlo, jhi, &
01798 dt, &
01799 aice, aice0, &
01800 aicen(:,:,:), trcrn(:,:,:,:), &
01801 vicen(:,:,:), vsnon(:,:,:), &
01802 eicen(:,:,:), esnon(:,:,:), &
01803 dfresh, dfsalt, &
01804 dfhocn, dfsoot, &
01805 tr_aero, &
01806 l_stop, &
01807 istop, jstop)
01808 if (l_stop) return
01809 endif
01810
01811
01812
01813
01814
01815 if (present(fresh)) &
01816 fresh (:,:) = fresh(:,:) + dfresh(:,:)
01817 if (present(fsalt)) &
01818 fsalt (:,:) = fsalt(:,:) + dfsalt(:,:)
01819 if (present(fhocn)) &
01820 fhocn (:,:) = fhocn(:,:) + dfhocn(:,:)
01821 if (present(fsoot)) &
01822 fsoot (:,:,:) = fsoot(:,:,:) + dfsoot(:,:,:)
01823
01824
01825
01826
01827
01828
01829 if (.not. heat_capacity) then
01830
01831 call zerolayer_check(nx_block, ny_block, &
01832 icells, indxi, indxj, &
01833 aicen, &
01834 vicen, vsnon, &
01835 eicen, esnon, &
01836 l_stop, &
01837 istop, jstop)
01838
01839 endif
01840
01841 end subroutine cleanup_itd
01842
01843
01844
01845
01846
01847
01848
01849
01850 subroutine zap_small_areas (nx_block, ny_block, &
01851 ilo, ihi, jlo, jhi, &
01852 dt, &
01853 aice, aice0, &
01854 aicen, trcrn, &
01855 vicen, vsnon, &
01856 eicen, esnon, &
01857 dfresh, dfsalt, &
01858 dfhocn, dfsoot, &
01859 tr_aero, &
01860 l_stop, &
01861 istop, jstop)
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874 use ice_state, only: nt_Tsfc, nt_aero
01875
01876
01877
01878 integer (kind=int_kind), intent(in) ::
01879 nx_block, ny_block,
01880 ilo,ihi,jlo,jhi
01881
01882 real (kind=dbl_kind), intent(in) ::
01883 dt
01884
01885 real (kind=dbl_kind), dimension (nx_block,ny_block),
01886 intent(inout) ::
01887 aice ,
01888 aice0
01889
01890 real (kind=dbl_kind), dimension(nx_block,ny_block,ncat),
01891 intent(inout) ::
01892 aicen ,
01893 vicen ,
01894 vsnon
01895
01896 real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr),
01897 intent(inout) ::
01898 eicen
01899
01900 real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr),
01901 intent(inout) ::
01902 esnon
01903
01904 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr,ncat),
01905 intent(inout) ::
01906 trcrn
01907
01908 real (kind=dbl_kind), dimension (nx_block,ny_block),
01909 intent(out) ::
01910 dfresh ,
01911 dfsalt ,
01912 dfhocn
01913
01914 real (kind=dbl_kind), dimension (nx_block,ny_block,n_aeromx),
01915 intent(out) ::
01916 dfsoot
01917
01918 logical (kind=log_kind), intent(in) ::
01919 tr_aero
01920
01921 logical (kind=log_kind), intent(out) ::
01922 l_stop
01923
01924 integer (kind=int_kind), intent(out) ::
01925 istop, jstop
01926
01927
01928
01929 integer (kind=int_kind) ::
01930 i,j, n, k, it ,
01931 icells ,
01932 ij
01933
01934 integer (kind=int_kind), dimension (nx_block*ny_block) ::
01935 indxi ,
01936 indxj
01937
01938 real (kind=dbl_kind) :: xtmp
01939
01940
01941
01942
01943
01944 l_stop = .false.
01945 istop = 0
01946 jstop = 0
01947
01948 dfresh(:,:) = c0
01949 dfsalt(:,:) = c0
01950 dfhocn(:,:) = c0
01951 dfsoot(:,:,:) = c0
01952
01953
01954
01955
01956
01957 do n = 1, ncat
01958
01959
01960
01961
01962
01963 icells = 0
01964 do j = jlo, jhi
01965 do i = ilo, ihi
01966 if (aicen(i,j,n) < -puny) then
01967 write (nu_diag,*) 'Zap ice: negative ice area'
01968 write (nu_diag,*) 'i, j, n, aicen =', &
01969 i, j, n, aicen(i,j,n)
01970 l_stop = .true.
01971 istop = i
01972 jstop = j
01973 return
01974 elseif ((aicen(i,j,n) >= -puny .and. aicen(i,j,n) < c0) .or. &
01975 (aicen(i,j,n) > c0 .and. aicen(i,j,n) <= puny)) then
01976 icells = icells + 1
01977 indxi(icells) = i
01978 indxj(icells) = j
01979 endif
01980 enddo
01981 enddo
01982
01983
01984
01985
01986
01987 do k = 1, nilyr
01988
01989
01990
01991 do ij = 1, icells
01992 i = indxi(ij)
01993 j = indxj(ij)
01994
01995 xtmp = eicen(i,j,ilyr1(n)+k-1) / dt
01996 dfhocn(i,j) = dfhocn(i,j) + xtmp
01997 eicen(i,j,ilyr1(n)+k-1) = c0
01998
01999 enddo
02000 enddo
02001
02002
02003
02004
02005
02006 do k = 1, nslyr
02007
02008
02009
02010 do ij = 1, icells
02011 i = indxi(ij)
02012 j = indxj(ij)
02013
02014 xtmp = esnon(i,j,slyr1(n)+k-1) / dt
02015 dfhocn(i,j) = dfhocn(i,j) + xtmp
02016 esnon(i,j,slyr1(n)+k-1) = c0
02017
02018 enddo
02019 enddo
02020
02021
02022
02023
02024
02025
02026
02027
02028 do ij = 1, icells
02029 i = indxi(ij)
02030 j = indxj(ij)
02031
02032 xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt
02033 dfresh(i,j) = dfresh(i,j) + xtmp
02034
02035 xtmp = rhoi*vicen(i,j,n)*ice_ref_salinity*p001 / dt
02036 dfsalt(i,j) = dfsalt(i,j) + xtmp
02037
02038 aice0(i,j) = aice0(i,j) + aicen(i,j,n)
02039 aicen(i,j,n) = c0
02040 vicen(i,j,n) = c0
02041 vsnon(i,j,n) = c0
02042 trcrn(i,j,nt_Tsfc,n) = Tocnfrz
02043
02044 enddo
02045
02046 if (tr_aero) then
02047
02048
02049
02050 do ij = 1, icells
02051 i = indxi(ij)
02052 j = indxj(ij)
02053 do it=1,n_aero
02054 xtmp &
02055 = (vsnon(i,j,n)*(trcrn(i,j,nt_aero +4*(it-1),n) &
02056 +trcrn(i,j,nt_aero+1+4*(it-1),n)) &
02057 + vicen(i,j,n)*(trcrn(i,j,nt_aero+2+4*(it-1),n) &
02058 +trcrn(i,j,nt_aero+3+4*(it-1),n))) &
02059 / dt
02060 dfsoot(i,j,it) = dfsoot(i,j,it) + xtmp
02061 enddo
02062 enddo
02063 endif
02064
02065
02066
02067
02068
02069 if (ntrcr >= 2) then
02070 do it = 1, ntrcr
02071 do ij = 1, icells
02072 i = indxi(ij)
02073 j = indxj(ij)
02074 trcrn(i,j,it,n) = c0
02075 enddo
02076 enddo
02077 endif
02078
02079 enddo
02080
02081
02082
02083
02084
02085
02086 icells = 0
02087 do j = jlo, jhi
02088 do i = ilo, ihi
02089 if (aice(i,j) > (c1+puny)) then
02090 write (nu_diag,*) 'Zap ice: excess ice area'
02091 write (nu_diag,*) 'i, j, aice =', &
02092 i, j, aice(i,j)
02093 l_stop = .true.
02094 istop = i
02095 jstop = j
02096 return
02097 elseif (aice(i,j) > c1 .and. aice(i,j) < (c1+puny)) then
02098 icells = icells + 1
02099 indxi(icells) = i
02100 indxj(icells) = j
02101 endif
02102 enddo
02103 enddo
02104
02105 do n = 1, ncat
02106
02107
02108
02109
02110
02111 do k = 1, nilyr
02112
02113
02114
02115 do ij = 1, icells
02116 i = indxi(ij)
02117 j = indxj(ij)
02118
02119 xtmp = eicen(i,j,ilyr1(n)+k-1) &
02120 * (aice(i,j)-c1)/aice(i,j) / dt
02121 dfhocn(i,j) = dfhocn(i,j) + xtmp
02122 eicen(i,j,ilyr1(n)+k-1) = eicen(i,j,ilyr1(n)+k-1) &
02123 * (c1/aice(i,j))
02124
02125 enddo
02126 enddo
02127
02128
02129
02130
02131
02132 do k = 1, nslyr
02133
02134
02135
02136 do ij = 1, icells
02137 i = indxi(ij)
02138 j = indxj(ij)
02139
02140 xtmp = esnon(i,j,slyr1(n)+k-1) &
02141 * (aice(i,j)-c1)/aice(i,j) / dt
02142 dfhocn(i,j) = dfhocn(i,j) + xtmp
02143 esnon(i,j,slyr1(n)+k-1) = esnon(i,j,slyr1(n)+k-1) &
02144 *(c1/aice(i,j))
02145
02146 enddo
02147 enddo
02148
02149
02150
02151
02152
02153
02154
02155
02156 do ij = 1, icells
02157 i = indxi(ij)
02158 j = indxj(ij)
02159
02160 xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) &
02161 * (aice(i,j)-c1)/aice(i,j) / dt
02162 dfresh(i,j) = dfresh(i,j) + xtmp
02163
02164 xtmp = rhoi*vicen(i,j,n)*ice_ref_salinity*p001 &
02165 * (aice(i,j)-c1)/aice(i,j) / dt
02166 dfsalt(i,j) = dfsalt(i,j) + xtmp
02167
02168 aicen(i,j,n) = aicen(i,j,n) * (c1/aice(i,j))
02169 vicen(i,j,n) = vicen(i,j,n) * (c1/aice(i,j))
02170 vsnon(i,j,n) = vsnon(i,j,n) * (c1/aice(i,j))
02171
02172 enddo
02173
02174
02175
02176
02177
02178
02179 if (tr_aero) then
02180 do ij = 1, icells
02181 i = indxi(ij)
02182 j = indxj(ij)
02183 do it=1,n_aero
02184 xtmp &
02185 = (vsnon(i,j,n)*(trcrn(i,j,nt_aero +4*(it-1),n) &
02186 +trcrn(i,j,nt_aero+1+4*(it-1),n)) &
02187 + vicen(i,j,n)*(trcrn(i,j,nt_aero+2+4*(it-1),n) &
02188 +trcrn(i,j,nt_aero+3+4*(it-1),n))) &
02189 * (aice(i,j)-c1)/aice(i,j) / dt
02190 dfsoot(i,j,it) = dfsoot(i,j,it) + xtmp
02191 enddo
02192 enddo
02193 endif
02194
02195 enddo
02196
02197
02198
02199
02200
02201
02202
02203
02204 do ij = 1, icells
02205 i = indxi(ij)
02206 j = indxj(ij)
02207 aice(i,j) = c1
02208 aice0(i,j) = c0
02209 enddo
02210
02211 end subroutine zap_small_areas
02212
02213
02214
02215
02216
02217
02218
02219
02220
02221 subroutine zerolayer_check (nx_block, ny_block, &
02222 icells, indxi, indxj, &
02223 aicen, &
02224 vicen, vsnon, &
02225 eicen, esnon, &
02226 l_stop, &
02227 istop, jstop)
02228
02229
02230
02231
02232
02233
02234
02235
02236
02237
02238
02239
02240
02241
02242
02243
02244
02245 integer (kind=int_kind), intent(in) ::
02246 nx_block, ny_block,
02247 icells
02248
02249 integer (kind=int_kind), dimension (nx_block*ny_block),
02250 intent(in) ::
02251 indxi, indxj
02252
02253 real (kind=dbl_kind), dimension (nx_block,ny_block,ncat),
02254 intent(inout) ::
02255 aicen ,
02256 vicen ,
02257 vsnon
02258
02259 real (kind=dbl_kind), dimension (nx_block,ny_block,ntilyr),
02260 intent(in) ::
02261 eicen
02262
02263 real (kind=dbl_kind), dimension (nx_block,ny_block,ntslyr),
02264 intent(in) ::
02265 esnon
02266
02267 logical (kind=log_kind), intent(out) ::
02268 l_stop
02269
02270 integer (kind=int_kind), intent(out) ::
02271 istop, jstop
02272
02273
02274
02275 integer (kind=int_kind) ::
02276 i, j ,
02277 n ,
02278 ij
02279
02280 real (kind=dbl_kind), parameter ::
02281 max_error = puny*Lfresh*rhos
02282
02283
02284 logical (kind=log_kind) ::
02285 ice_energy_correct ,
02286 snow_energy_correct
02287
02288 real (kind=dbl_kind), dimension (nx_block,ny_block) ::
02289 worka,
02290 workb
02291
02292
02293
02294
02295
02296 l_stop = .false.
02297 istop = 0
02298 jstop = 0
02299
02300 worka(:,:) = c0
02301 workb(:,:) = c0
02302
02303
02304
02305
02306
02307
02308 ice_energy_correct = .true.
02309 snow_energy_correct = .true.
02310
02311 do n=1,ncat
02312
02313 do ij=1,icells
02314 i=indxi(ij)
02315 j=indxj(ij)
02316
02317 worka(i,j) = eicen(i,j,n) + rhoi * Lfresh * vicen(i,j,n)
02318 workb(i,j) = esnon(i,j,n) + rhos * Lfresh * vsnon(i,j,n)
02319
02320 if(abs(worka(i,j)) > max_error) then
02321 ice_energy_correct = .false.
02322 endif
02323
02324 if(abs(workb(i,j)) > max_error) then
02325 snow_energy_correct = .false.
02326 endif
02327 enddo
02328
02329
02330
02331
02332
02333 if (.not. ice_energy_correct) then
02334
02335 do ij=1,icells
02336 i=indxi(ij)
02337 j=indxj(ij)
02338
02339 if(abs(worka(i,j)) > max_error) then
02340 write(nu_diag,*) ' '
02341 write(nu_diag,*) &
02342 'zerolayer check - wrong ice energy'
02343 write(nu_diag,*) 'i, j, n:', i,j,n
02344 write(nu_diag,*) 'eicen =', eicen(i,j,n)
02345 write(nu_diag,*) 'error=', worka(i,j)
02346 write(nu_diag,*) 'vicen =', vicen(i,j,n)
02347 write(nu_diag,*) 'aicen =', aicen(i,j,n)
02348 l_stop = .true.
02349 istop = i
02350 jstop = j
02351 endif
02352 enddo
02353
02354 endif
02355 if (l_stop) return
02356
02357 if (.not. snow_energy_correct) then
02358
02359 do ij=1,icells
02360 i=indxi(ij)
02361 j=indxj(ij)
02362
02363 if(abs(workb(i,j)) > max_error) then
02364 write(nu_diag,*) ' '
02365 write(nu_diag,*) &
02366 'zerolayer_check - wrong snow energy'
02367 write(nu_diag,*) 'i, j, n:', i,j,n
02368 write(nu_diag,*) 'esnon =', esnon(i,j,n)
02369 write(nu_diag,*) 'error=', workb(i,j)
02370 write(nu_diag,*) 'vsnon =', vsnon(i,j,n)
02371 write(nu_diag,*) 'aicen =', aicen(i,j,n)
02372 l_stop = .true.
02373 istop = i
02374 jstop = j
02375 return
02376 endif
02377 enddo
02378
02379 endif
02380
02381 enddo
02382
02383 end subroutine zerolayer_check
02384
02385
02386
02387 end module ice_itd
02388
02389
02390
02391
02392
02393
02394
02395
02396
02397
02398