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 module ice_therm_vertical
00031
00032
00033
00034 use ice_kinds_mod
00035 use ice_domain_size, only: ncat, nilyr, nslyr, ntilyr, ntslyr, ntrcr
00036 use ice_constants
00037 use ice_fileunits, only: nu_diag
00038 use ice_prescribed_mod, only: prescribed_ice
00039
00040
00041
00042 implicit none
00043 save
00044
00045 real (kind=dbl_kind), parameter ::
00046 saltmax = 3.2_dbl_kind,
00047 hs_min = 1.e-4_dbl_kind,
00048 betak = 0.13_dbl_kind,
00049 kimin = 0.10_dbl_kind
00050
00051 real (kind=dbl_kind), dimension(nilyr+1) ::
00052 salin ,
00053 Tmlt
00054
00055
00056 real (kind=dbl_kind) ::
00057 ustar_scale
00058
00059 real (kind=dbl_kind), parameter, private ::
00060 ferrmax = 1.0e-3_dbl_kind
00061
00062
00063 character (char_len) :: stoplabel
00064
00065 logical (kind=log_kind) ::
00066 l_brine
00067
00068 logical (kind=log_kind) ::
00069 heat_capacity,
00070
00071 calc_Tsfc
00072
00073
00074
00075
00076
00077 contains
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098 subroutine thermo_vertical (nx_block, ny_block, &
00099 dt, icells, &
00100 indxi, indxj, &
00101 aicen, trcrn, &
00102 vicen, vsnon, &
00103 eicen, esnon, &
00104 flw, potT, &
00105 Qa, rhoa, &
00106 fsnow, &
00107 fbot, Tbot, &
00108 lhcoef, shcoef, &
00109 fswsfc, fswint, &
00110 fswthrun, &
00111 Sswabs, Iswabs, &
00112 fsurfn, fcondtopn, &
00113 fsensn, flatn, &
00114 fswabsn, flwoutn, &
00115 evapn, freshn, &
00116 fsaltn, fhocnn, &
00117 meltt, melts, &
00118 meltb, &
00119 congel, snoice, &
00120 mlt_onset, frz_onset, &
00121 yday, l_stop, &
00122 istop, jstop)
00123
00124
00125
00126
00127 use ice_communicate, only: my_task, master_task
00128 use ice_calendar, only: istep1
00129 use ice_exit
00130 use ice_ocean
00131 use ice_itd, only: ilyr1, slyr1, ilyrn, slyrn
00132 use ice_state, only: nt_Tsfc, nt_iage, tr_iage
00133
00134
00135
00136 integer (kind=int_kind), intent(in) ::
00137 nx_block, ny_block,
00138 icells
00139
00140 integer (kind=int_kind), dimension (nx_block*ny_block),
00141 intent(in) ::
00142 indxi, indxj
00143
00144 real (kind=dbl_kind), intent(in) ::
00145 dt
00146
00147
00148 real (kind=dbl_kind), dimension (nx_block,ny_block),
00149 intent(inout) ::
00150 aicen ,
00151 vicen ,
00152 vsnon
00153
00154 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr),
00155 intent(inout) ::
00156 trcrn
00157
00158 real (kind=dbl_kind), dimension(nx_block,ny_block,nilyr),
00159 intent(inout) ::
00160 eicen
00161
00162 real (kind=dbl_kind), dimension(nx_block,ny_block,nslyr),
00163 intent(inout) ::
00164 esnon
00165
00166
00167 real (kind=dbl_kind), dimension (nx_block,ny_block),
00168 intent(in) ::
00169 flw ,
00170 potT ,
00171 Qa ,
00172 rhoa ,
00173 fsnow ,
00174 shcoef ,
00175 lhcoef
00176
00177 real (kind=dbl_kind), dimension (nx_block,ny_block),
00178 intent(inout) ::
00179 fswsfc ,
00180 fswint ,
00181 fswthrun
00182
00183 real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr),
00184 intent(inout) ::
00185 Sswabs
00186
00187 real (kind=dbl_kind), dimension (nx_block,ny_block,nilyr),
00188 intent(inout) ::
00189 Iswabs
00190
00191
00192 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
00193 fbot ,
00194 Tbot
00195
00196
00197 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out)::
00198 fsensn ,
00199 fswabsn ,
00200 flwoutn ,
00201 evapn
00202
00203
00204 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout)::
00205 flatn ,
00206 fsurfn ,
00207 fcondtopn
00208
00209
00210 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out)::
00211 freshn ,
00212 fsaltn ,
00213 fhocnn
00214
00215
00216 real (kind=dbl_kind), dimension(nx_block,ny_block),
00217 intent(inout)::
00218 meltt ,
00219 melts ,
00220 meltb ,
00221 congel ,
00222 snoice ,
00223 mlt_onset,
00224 frz_onset
00225
00226 real (kind=dbl_kind), intent(in) ::
00227 yday
00228
00229 logical (kind=log_kind), intent(out) ::
00230 l_stop
00231
00232 integer (kind=int_kind), intent(out) ::
00233 istop, jstop
00234
00235
00236
00237 integer (kind=int_kind) ::
00238 i, j ,
00239 ij ,
00240 ilo,ihi,jlo,jhi,
00241 k ,
00242 il1, il2 ,
00243 sl1, sl2
00244
00245 real (kind=dbl_kind) ::
00246 dhi ,
00247 dhs
00248
00249
00250
00251 real (kind=dbl_kind), dimension (icells) ::
00252 hilyr ,
00253 hslyr ,
00254 Tsf ,
00255 hin ,
00256 hsn ,
00257 hsn_new ,
00258 worki ,
00259 works
00260
00261 real (kind=dbl_kind), dimension (icells,nilyr) ::
00262 qin ,
00263 Tin
00264
00265 real (kind=dbl_kind), dimension (icells,nslyr) ::
00266 qsn ,
00267 Tsn
00268
00269
00270
00271 real (kind=dbl_kind), dimension (icells) ::
00272 fcondbot ,
00273 einit ,
00274 efinal
00275
00276
00277 real (kind=dbl_kind), dimension (nx_block,ny_block) ::
00278 Tsfcn,
00279 iage
00280
00281
00282
00283
00284
00285 l_stop = .false.
00286 istop = 0
00287 jstop = 0
00288
00289 do j=1, ny_block
00290 do i=1, nx_block
00291 fsensn (i,j) = c0
00292 fswabsn(i,j) = c0
00293 flwoutn(i,j) = c0
00294 evapn (i,j) = c0
00295
00296 freshn (i,j) = c0
00297 fsaltn (i,j) = c0
00298 fhocnn (i,j) = c0
00299
00300 meltt (i,j) = c0
00301 meltb (i,j) = c0
00302 melts (i,j) = c0
00303 congel (i,j) = c0
00304 snoice (i,j) = c0
00305
00306 Tsfcn(i,j) = trcrn(i,j,nt_Tsfc)
00307 if (tr_iage) iage(i,j) = trcrn(i,j,nt_iage)
00308 enddo
00309 enddo
00310
00311 if (calc_Tsfc) then
00312 do j=1, ny_block
00313 do i=1, nx_block
00314 flatn (i,j) = c0
00315 fsurfn (i,j) = c0
00316 fcondtopn(i,j) = c0
00317 enddo
00318 enddo
00319 endif
00320
00321
00322
00323
00324
00325 call init_vertical_profile (nx_block, ny_block, &
00326 my_task, istep1, &
00327 icells, &
00328 indxi, indxj, &
00329 aicen(:,:), &
00330 vicen(:,:), vsnon(:,:), &
00331 Tsfcn(:,:), &
00332 eicen(:,:,:), esnon(:,:,:), &
00333 hin, hilyr, &
00334 hsn, hslyr, &
00335 qin, Tin, &
00336 qsn, Tsn, &
00337 Tsf, einit, &
00338 l_stop, &
00339 istop, jstop)
00340
00341 if (l_stop) return
00342
00343 do ij = 1, icells
00344
00345 worki(ij) = hin(ij)
00346 works(ij) = hsn(ij)
00347 enddo
00348
00349
00350
00351
00352
00353
00354 if (heat_capacity) then
00355
00356 call temperature_changes(nx_block, ny_block, &
00357 my_task, istep1, &
00358 dt, icells, &
00359 indxi, indxj, &
00360 rhoa, flw, &
00361 potT, Qa, &
00362 shcoef, lhcoef, &
00363 fswsfc, fswint, &
00364 fswthrun, Sswabs(:,:,:),&
00365 Iswabs, &
00366 hilyr, hslyr, &
00367 qin, Tin, &
00368 qsn, Tsn, &
00369 Tsf, Tbot, &
00370 fsensn, flatn, &
00371 fswabsn, flwoutn, &
00372 fsurfn, &
00373 fcondtopn, fcondbot, &
00374 einit, l_stop, &
00375 istop, jstop)
00376
00377 else
00378
00379 if (calc_Tsfc) then
00380
00381 call zerolayer_temperature(nx_block, ny_block, &
00382 my_task, istep1, &
00383 dt, icells, &
00384 indxi, indxj, &
00385 rhoa, flw, &
00386 potT, Qa, &
00387 shcoef, lhcoef, &
00388 fswsfc, fswthrun, &
00389 hilyr, hslyr, &
00390 Tsf, Tbot, &
00391 fsensn, flatn, &
00392 fswabsn, flwoutn, &
00393 fsurfn, &
00394 fcondtopn, fcondbot, &
00395 l_stop, &
00396 istop, jstop)
00397
00398 else
00399
00400
00401
00402
00403
00404
00405 do ij = 1, icells
00406 i = indxi(ij)
00407 j = indxj(ij)
00408 fcondbot(ij) = fcondtopn(i,j)
00409 enddo
00410
00411 endif
00412
00413 endif
00414
00415 if (l_stop) return
00416
00417
00418
00419
00420
00421
00422
00423 call thickness_changes(nx_block, ny_block, &
00424 dt, &
00425 yday, icells, &
00426 indxi, indxj, &
00427 efinal, &
00428 hin, hilyr, &
00429 hsn, hslyr, &
00430 qin, qsn, &
00431 fbot, Tbot, &
00432 flatn, fsurfn, &
00433 fcondtopn, fcondbot, &
00434 fsnow, hsn_new, &
00435 fhocnn, evapn, &
00436 meltt, melts, &
00437 meltb, iage, &
00438 congel, snoice, &
00439 mlt_onset, frz_onset)
00440
00441
00442
00443
00444
00445
00446 call conservation_check_vthermo(nx_block, ny_block, &
00447 my_task, istep1, &
00448 dt, icells, &
00449 indxi, indxj, &
00450 fsurfn, flatn, &
00451 fhocnn, fswint, &
00452 fsnow, &
00453 einit, efinal, &
00454 l_stop, &
00455 istop, jstop)
00456
00457 if (l_stop) return
00458
00459
00460
00461
00462
00463 if (prescribed_ice) then
00464 do ij = 1, icells
00465 i = indxi(ij)
00466 j = indxj(ij)
00467 hin(ij) = worki(ij)
00468 fhocnn(i,j) = c0
00469 enddo
00470 endif
00471
00472
00473
00474
00475
00476
00477 do ij = 1, icells
00478 i = indxi(ij)
00479 j = indxj(ij)
00480
00481 dhi = hin(ij) - worki(ij)
00482 dhs = hsn(ij) - works(ij)
00483
00484 freshn(i,j) = evapn(i,j) - &
00485 (rhoi*dhi + rhos*(dhs-hsn_new(ij))) / dt
00486 fsaltn(i,j) = -rhoi*dhi*ice_ref_salinity*p001/dt
00487
00488 enddo
00489
00490
00491
00492
00493
00494
00495
00496 call update_state_vthermo(nx_block, ny_block, &
00497 icells, &
00498 indxi, indxj, &
00499 Tbot, Tsf, &
00500 hin, hsn, &
00501 qin, qsn, &
00502 aicen(:,:), &
00503 vicen(:,:), vsnon(:,:), &
00504 Tsfcn(:,:), &
00505 eicen(:,:,:), esnon(:,:,:))
00506
00507
00508
00509
00510
00511 do j = 1, ny_block
00512 do i = 1, nx_block
00513 trcrn(i,j,nt_Tsfc) = Tsfcn(i,j)
00514 if (tr_iage) trcrn(i,j,nt_iage) = iage(i,j)
00515 enddo
00516 enddo
00517
00518 end subroutine thermo_vertical
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536 subroutine init_thermo_vertical
00537
00538
00539
00540
00541
00542
00543
00544 real (kind=dbl_kind), parameter ::
00545 nsal = 0.407_dbl_kind,
00546 msal = 0.573_dbl_kind,
00547 min_salin = 0.1_dbl_kind
00548
00549 integer (kind=int_kind) :: k
00550 real (kind=dbl_kind) :: zn
00551
00552
00553
00554
00555
00556
00557
00558
00559 if (saltmax > min_salin .and. heat_capacity) then
00560 l_brine = .true.
00561 else
00562 l_brine = .false.
00563 endif
00564
00565
00566
00567
00568
00569 if (l_brine) then
00570 do k = 1, nilyr
00571 zn = (real(k,kind=dbl_kind)-p5) /
00572 real(nilyr,kind=dbl_kind)
00573 salin(k)=(saltmax/c2)*(c1-cos(pi*zn**(nsal/(msal+zn))))
00574
00575 Tmlt(k) = -salin(k)*depressT
00576 enddo
00577 salin(nilyr+1) = saltmax
00578 Tmlt(nilyr+1) = -salin(nilyr+1)*depressT
00579 else
00580 do k = 1, nilyr+1
00581 salin(k) = c0
00582 Tmlt(k) = c0
00583 enddo
00584 endif
00585
00586 ustar_scale = c1
00587
00588 end subroutine init_thermo_vertical
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609 subroutine frzmlt_bottom_lateral (nx_block, ny_block, &
00610 ilo, ihi, jlo, jhi, &
00611 dt, &
00612 aice, frzmlt, &
00613 eicen, esnon, &
00614 sst, Tf, &
00615 strocnxT, strocnyT, &
00616 Tbot, fbot, &
00617 rside)
00618
00619
00620
00621 use ice_itd, only: ilyr1, slyr1
00622
00623
00624
00625 integer (kind=int_kind), intent(in) ::
00626 nx_block, ny_block,
00627 ilo,ihi,jlo,jhi
00628
00629 real (kind=dbl_kind), intent(in) ::
00630 dt
00631
00632 real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) ::
00633 aice ,
00634 frzmlt ,
00635 sst ,
00636 Tf ,
00637 strocnxT,
00638 strocnyT
00639
00640 real (kind=dbl_kind), dimension(nx_block,ny_block,ntilyr),
00641 intent(in) ::
00642 eicen
00643
00644 real (kind=dbl_kind), dimension(nx_block,ny_block,ntslyr),
00645 intent(in) ::
00646 esnon
00647
00648 real (kind=dbl_kind), dimension(nx_block,ny_block),
00649 intent(out) ::
00650 Tbot ,
00651 fbot ,
00652 rside
00653
00654
00655
00656 integer (kind=int_kind) ::
00657 i, j ,
00658 n ,
00659 k ,
00660 ij ,
00661 imelt
00662
00663 integer (kind=int_kind), dimension (nx_block*ny_block) ::
00664 indxi, indxj
00665
00666 real (kind=dbl_kind), dimension (:), allocatable ::
00667 etot ,
00668 fside
00669
00670 real (kind=dbl_kind) ::
00671 deltaT ,
00672 ustar ,
00673 wlat ,
00674 xtmp
00675
00676
00677
00678
00679
00680 real (kind=dbl_kind), parameter ::
00681 cpchr = -cp_ocn*rhow*0.006_dbl_kind,
00682 ustar_min = 5.e-3_dbl_kind
00683
00684
00685
00686
00687 real (kind=dbl_kind), parameter ::
00688 floediam = 300.0_dbl_kind,
00689 alpha = 0.66_dbl_kind ,
00690 m1 = 1.6e-6_dbl_kind ,
00691
00692 m2 = 1.36_dbl_kind
00693
00694
00695 do j = 1, ny_block
00696 do i = 1, nx_block
00697 rside(i,j) = c0
00698 Tbot (i,j) = Tf(i,j)
00699 fbot (i,j) = c0
00700 enddo
00701 enddo
00702
00703
00704
00705
00706
00707 imelt = 0
00708 do j = jlo, jhi
00709 do i = ilo, ihi
00710 if (aice(i,j) > puny .and. frzmlt(i,j) < c0) then
00711 imelt = imelt + 1
00712 indxi(imelt) = i
00713 indxj(imelt) = j
00714 endif
00715 enddo
00716 enddo
00717
00718 allocate(etot (imelt))
00719 allocate(fside(imelt))
00720
00721 do ij = 1, imelt
00722 i = indxi(ij)
00723 j = indxj(ij)
00724
00725 fside(ij) = c0
00726
00727
00728
00729
00730
00731
00732 deltaT = max((sst(i,j)-Tbot(i,j)),c0)
00733
00734
00735 ustar = sqrt (sqrt(strocnxT(i,j)**2+strocnyT(i,j)**2)/rhow)
00736 ustar = max (ustar,ustar_min*ustar_scale)
00737
00738 fbot(i,j) = cpchr * deltaT * ustar
00739 fbot(i,j) = max (fbot(i,j), frzmlt(i,j))
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750 wlat = m1 * deltaT**m2
00751 rside(i,j) = wlat*dt*pi/(alpha*floediam)
00752 rside(i,j) = max(c0,min(rside(i,j),c1))
00753
00754 enddo
00755
00756
00757
00758
00759
00760 do n = 1, ncat
00761
00762 do ij = 1, imelt
00763 etot(ij) = c0
00764 enddo
00765
00766
00767
00768 do k = 1, nslyr
00769
00770
00771
00772 do ij = 1, imelt
00773 i = indxi(ij)
00774 j = indxj(ij)
00775 etot(ij) = etot(ij) + esnon(i,j,slyr1(n)+k-1)
00776 enddo
00777 enddo
00778
00779 do k = 1, nilyr
00780
00781
00782
00783 do ij = 1, imelt
00784 i = indxi(ij)
00785 j = indxj(ij)
00786 etot(ij) = etot(ij) + eicen(i,j,ilyr1(n)+k-1)
00787 enddo
00788 enddo
00789
00790
00791
00792
00793 do ij = 1, imelt
00794 i = indxi(ij)
00795 j = indxj(ij)
00796
00797 fside(ij) = fside(ij) + rside(i,j)*etot(ij)/dt
00798 enddo
00799
00800 enddo
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811 do ij = 1, imelt
00812 i = indxi(ij)
00813 j = indxj(ij)
00814
00815 xtmp = frzmlt(i,j)/(fbot(i,j) + fside(ij) + puny)
00816 xtmp = min(xtmp, c1)
00817 fbot(i,j) = fbot(i,j) * xtmp
00818 rside(i,j) = rside(i,j) * xtmp
00819 enddo
00820
00821 deallocate(etot)
00822 deallocate(fside)
00823
00824 end subroutine frzmlt_bottom_lateral
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844 subroutine init_vertical_profile(nx_block, ny_block, &
00845 my_task, istep1, &
00846 icells, &
00847 indxi, indxj, &
00848 aicen, vicen, &
00849 vsnon, Tsfcn, &
00850 eicen, esnon, &
00851 hin, hilyr, &
00852 hsn, hslyr, &
00853 qin, Tin, &
00854 qsn, Tsn, &
00855 Tsf, einit, &
00856 l_stop, &
00857 istop, jstop)
00858
00859
00860
00861
00862
00863 integer (kind=int_kind), intent(in) ::
00864 nx_block, ny_block,
00865 my_task ,
00866 istep1 ,
00867 icells
00868
00869 integer (kind=int_kind), dimension(nx_block*ny_block),
00870 intent(in) ::
00871 indxi, indxj
00872
00873 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
00874 aicen ,
00875 vicen ,
00876 vsnon ,
00877 Tsfcn
00878
00879 real (kind=dbl_kind), dimension(nx_block,ny_block,nilyr),
00880 intent(in) ::
00881 eicen
00882
00883 real (kind=dbl_kind), dimension(nx_block,ny_block,nslyr),
00884 intent(in) ::
00885 esnon
00886
00887 real (kind=dbl_kind), dimension(icells), intent(out)::
00888 hilyr ,
00889 hslyr ,
00890 Tsf ,
00891 einit
00892
00893 real (kind=dbl_kind), dimension(icells), intent(out)::
00894 hin ,
00895 hsn
00896
00897 real (kind=dbl_kind), dimension (icells,nilyr),
00898 intent(out) ::
00899 qin ,
00900 Tin
00901
00902 real (kind=dbl_kind), dimension (icells,nslyr),
00903 intent(out) ::
00904 qsn ,
00905 Tsn
00906
00907 logical (kind=log_kind), intent(inout) ::
00908 l_stop
00909
00910 integer (kind=int_kind), intent(inout) ::
00911 istop, jstop
00912
00913
00914
00915 real (kind=dbl_kind), parameter ::
00916 Tmin = -100._dbl_kind
00917
00918 integer (kind=int_kind) ::
00919 i, j ,
00920 ij ,
00921 k
00922
00923 real (kind=dbl_kind) ::
00924 rnslyr,
00925 aa1, bb1, cc1,
00926 Tmax
00927
00928 logical (kind=log_kind) ::
00929 tsno_high ,
00930 tice_high ,
00931 tsno_low ,
00932 tice_low
00933
00934
00935
00936
00937
00938 rnslyr = real(nslyr,kind=dbl_kind)
00939
00940 do ij = 1, icells
00941 einit(ij) = c0
00942 enddo
00943
00944 tsno_high = .false.
00945 tice_high = .false.
00946 tsno_low = .false.
00947 tice_low = .false.
00948
00949
00950
00951
00952
00953
00954
00955 do ij = 1, icells
00956 i = indxi(ij)
00957 j = indxj(ij)
00958
00959
00960
00961
00962
00963
00964 Tsf(ij) = Tsfcn(i,j)
00965 hin(ij) = vicen(i,j) / aicen(i,j)
00966 hsn(ij) = vsnon(i,j) / aicen(i,j)
00967 hilyr(ij) = hin(ij) / real(nilyr,kind=dbl_kind)
00968 hslyr(ij) = hsn(ij) / rnslyr
00969
00970 enddo
00971
00972
00973
00974
00975
00976
00977 do k = 1, nslyr
00978
00979
00980
00981 do ij = 1, icells
00982 i = indxi(ij)
00983 j = indxj(ij)
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993 if (hslyr(ij) > hs_min/rnslyr .and. heat_capacity) then
00994
00995 qsn (ij,k) = esnon(i,j,k)*rnslyr/vsnon(i,j)
00996 Tmax = -qsn(ij,k)*puny*rnslyr / &
00997 (rhos*cp_ice*vsnon(i,j))
00998 else
00999 qsn (ij,k) = -rhos * Lfresh
01000 Tmax = puny
01001 endif
01002
01003
01004
01005
01006
01007 Tsn(ij,k) = (Lfresh + qsn(ij,k)/rhos)/cp_ice
01008
01009
01010
01011
01012 if (Tsn(ij,k) > Tmax) then
01013 tsno_high = .true.
01014 elseif (Tsn(ij,k) < Tmin) then
01015 tsno_low = .true.
01016 endif
01017
01018 enddo
01019 enddo
01020
01021
01022
01023
01024
01025 if (tsno_high .and. heat_capacity) then
01026 do k = 1, nslyr
01027 do ij = 1, icells
01028 i = indxi(ij)
01029 j = indxj(ij)
01030
01031 if (hslyr(ij) > hs_min/rnslyr) then
01032 Tmax = -qsn(ij,k)*puny*rnslyr / &
01033 (rhos*cp_ice*vsnon(i,j))
01034 else
01035 Tmax = puny
01036 endif
01037
01038 if (Tsn(ij,k) > Tmax) then
01039 write(nu_diag,*) ' '
01040 write(nu_diag,*) 'Starting thermo, Tsn > Tmax'
01041 write(nu_diag,*) 'Tsn=',Tsn(ij,k)
01042 write(nu_diag,*) 'Tmax=',Tmax
01043 write(nu_diag,*) 'istep1, my_task, i, j:', &
01044 istep1, my_task, i, j
01045 write(nu_diag,*) 'qsn',qsn(ij,k)
01046 l_stop = .true.
01047 istop = i
01048 jstop = j
01049 return
01050 endif
01051
01052 enddo
01053 enddo
01054 endif
01055
01056 if (tsno_low .and. heat_capacity) then
01057 do k = 1, nslyr
01058 do ij = 1, icells
01059 i = indxi(ij)
01060 j = indxj(ij)
01061
01062 if (Tsn(ij,k) < Tmin) then
01063 write(nu_diag,*) ' '
01064 write(nu_diag,*) 'Starting thermo, Tsn < Tmin'
01065 write(nu_diag,*) 'Tsn=', Tsn(ij,k)
01066 write(nu_diag,*) 'Tmin=', Tmin
01067 write(nu_diag,*) 'istep1, my_task, i, j:', &
01068 istep1, my_task, i, j
01069 write(nu_diag,*) 'qsn', qsn(ij,k)
01070 l_stop = .true.
01071 istop = i
01072 jstop = j
01073 return
01074 endif
01075
01076 enddo
01077 enddo
01078 endif
01079
01080 do k = 1, nslyr
01081
01082
01083
01084 do ij = 1, icells
01085
01086 if (Tsn(ij,k) > c0) then
01087 Tsn(ij,k) = c0
01088 qsn(ij,k) = -rhos*Lfresh
01089 endif
01090
01091
01092
01093
01094 einit(ij) = einit(ij) + hslyr(ij)*qsn(ij,k)
01095
01096 enddo
01097 enddo
01098
01099 do k = 1, nilyr
01100
01101
01102
01103 do ij = 1, icells
01104 i = indxi(ij)
01105 j = indxj(ij)
01106
01107
01108
01109
01110
01111
01112 qin(ij,k) = eicen(i,j,k)*real(nilyr,kind=dbl_kind)
01113 /vicen(i,j)
01114
01115
01116
01117
01118
01119 if (l_brine) then
01120 aa1 = cp_ice
01121 bb1 = (cp_ocn-cp_ice)*Tmlt(k) - qin(ij,k)/rhoi - Lfresh
01122 cc1 = Lfresh * Tmlt(k)
01123 Tin(ij,k) = (-bb1 - sqrt(bb1*bb1 - c4*aa1*cc1)) / &
01124 (c2*aa1)
01125 Tmax = Tmlt(k)
01126
01127 else
01128 Tin(ij,k) = (Lfresh + qin(ij,k)/rhoi) / cp_ice
01129 Tmax = -qin(ij,k)*puny/(rhos*cp_ice*vicen(i,j))
01130
01131 endif
01132
01133
01134
01135
01136 if (Tin(ij,k) > Tmax) then
01137 tice_high = .true.
01138 elseif (Tin(ij,k) < Tmin) then
01139 tice_low = .true.
01140 endif
01141
01142 enddo
01143
01144
01145
01146
01147
01148 if (tice_high .and. heat_capacity) then
01149 do ij = 1, icells
01150 i = indxi(ij)
01151 j = indxj(ij)
01152
01153 if (l_brine) then
01154 Tmax = Tmlt(k)
01155 else
01156 Tmax = -qin(ij,k)*puny/(rhos*cp_ice*vicen(i,j))
01157 endif
01158
01159 if (Tin(ij,k) > Tmax) then
01160 write(nu_diag,*) ' '
01161 write(nu_diag,*) 'Starting thermo, T > Tmax, layer', k
01162 write(nu_diag,*) 'Tin=',Tin(ij,k),', Tmax=',Tmax
01163 write(nu_diag,*) 'istep1, my_task, i, j:', &
01164 istep1, my_task, i, j
01165 write(nu_diag,*) 'qin',qin(ij,k)
01166 l_stop = .true.
01167 istop = i
01168 jstop = j
01169 return
01170 endif
01171 enddo
01172 endif
01173
01174 if (tice_low .and. heat_capacity) then
01175 do ij = 1, icells
01176 i = indxi(ij)
01177 j = indxj(ij)
01178
01179 if (Tin(ij,k) < Tmin) then
01180 write(nu_diag,*) ' '
01181 write(nu_diag,*) 'Starting thermo T < Tmin, layer', k
01182 write(nu_diag,*) 'Tin =', Tin(ij,k)
01183 write(nu_diag,*) 'Tmin =', Tmin
01184 write(nu_diag,*) 'istep1, my_task, i, j:', &
01185 istep1, my_task, i, j
01186 l_stop = .true.
01187 istop = i
01188 jstop = j
01189 return
01190 endif
01191 enddo
01192 endif
01193
01194
01195
01196
01197
01198
01199
01200
01201 do ij = 1, icells
01202
01203 if (Tin(ij,k) > c0) then
01204 Tin(ij,k) = c0
01205 qin(ij,k) = -rhoi*Lfresh
01206 endif
01207
01208 einit(ij) = einit(ij) + hilyr(ij)*qin(ij,k)
01209
01210 enddo
01211
01212
01213 enddo
01214
01215 end subroutine init_vertical_profile
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243 subroutine temperature_changes (nx_block, ny_block, &
01244 my_task, istep1, &
01245 dt, icells, &
01246 indxi, indxj, &
01247 rhoa, flw, &
01248 potT, Qa, &
01249 shcoef, lhcoef, &
01250 fswsfc, fswint, &
01251 fswthrun, Sswabs, &
01252 Iswabs, &
01253 hilyr, hslyr, &
01254 qin, Tin, &
01255 qsn, Tsn, &
01256 Tsf, Tbot, &
01257 fsensn, flatn, &
01258 fswabsn, flwoutn, &
01259 fsurfn, &
01260 fcondtopn,fcondbot, &
01261 einit, l_stop, &
01262 istop, jstop)
01263
01264
01265
01266
01267
01268 integer (kind=int_kind), intent(in) ::
01269 nx_block, ny_block,
01270 my_task ,
01271 istep1 ,
01272 icells
01273
01274 real (kind=dbl_kind), intent(in) ::
01275 dt
01276
01277 integer (kind=int_kind), dimension(nx_block*ny_block),
01278 intent(in) ::
01279 indxi, indxj
01280
01281 real (kind=dbl_kind), dimension (nx_block,ny_block),
01282 intent(in) ::
01283 rhoa ,
01284 flw ,
01285 potT ,
01286 Qa ,
01287 shcoef ,
01288 lhcoef ,
01289 Tbot
01290
01291 real (kind=dbl_kind), dimension (nx_block,ny_block),
01292 intent(inout) ::
01293 fswsfc ,
01294 fswint ,
01295 fswthrun
01296
01297 real (kind=dbl_kind), dimension (icells), intent(in) ::
01298 hilyr ,
01299 hslyr ,
01300 einit
01301
01302 real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr),
01303 intent(inout) ::
01304 Sswabs
01305
01306 real (kind=dbl_kind), dimension (nx_block,ny_block,nilyr),
01307 intent(inout) ::
01308 Iswabs
01309
01310 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout)::
01311 fsurfn ,
01312 fcondtopn ,
01313 fsensn ,
01314 flatn ,
01315 fswabsn ,
01316 flwoutn
01317
01318 real (kind=dbl_kind), dimension (icells), intent(out)::
01319 fcondbot
01320
01321 real (kind=dbl_kind), dimension (icells),
01322 intent(inout)::
01323 Tsf
01324
01325 real (kind=dbl_kind), dimension (icells,nilyr),
01326 intent(inout) ::
01327 qin ,
01328 Tin
01329
01330 real (kind=dbl_kind), dimension (icells,nslyr),
01331 intent(inout) ::
01332 qsn ,
01333 Tsn
01334
01335 logical (kind=log_kind), intent(inout) ::
01336 l_stop
01337
01338 integer (kind=int_kind), intent(inout) ::
01339 istop, jstop
01340
01341
01342
01343 integer (kind=int_kind), parameter ::
01344 nitermax = 50,
01345 nmat = nslyr + nilyr + 1
01346
01347 real (kind=dbl_kind), parameter ::
01348 Tsf_errmax = 5.e-4_dbl_kind
01349
01350
01351 integer (kind=int_kind) ::
01352 i, j ,
01353 ij, m ,
01354 k ,
01355 niter
01356
01357 integer (kind=int_kind) ::
01358 isolve
01359
01360 integer (kind=int_kind), dimension (icells) ::
01361 indxii, indxjj
01362
01363 integer (kind=int_kind), dimension (icells) ::
01364 indxij
01365
01366 logical (kind=log_kind), dimension (icells) ::
01367 l_snow ,
01368 l_cold
01369
01370 real (kind=dbl_kind), dimension (:), allocatable ::
01371 Tsf_start ,
01372 dTsf ,
01373 dTi1 ,
01374 dfsurf_dT ,
01375 avg_Tsi ,
01376 enew
01377
01378 real (kind=dbl_kind), dimension (icells) ::
01379 dTsf_prev ,
01380 dTi1_prev ,
01381 dfsens_dT ,
01382 dflat_dT ,
01383 dflwout_dT ,
01384 dt_rhoi_hlyr
01385
01386 real (kind=dbl_kind), dimension (icells,nilyr) ::
01387 Tin_init ,
01388 Tin_start
01389
01390 real (kind=dbl_kind), dimension (icells,nslyr) ::
01391 Tsn_init ,
01392 Tsn_start ,
01393 etas
01394
01395 real (kind=dbl_kind), dimension (:,:), allocatable ::
01396 etai ,
01397 sbdiag ,
01398 diag ,
01399 spdiag ,
01400 rhs ,
01401 Tmat
01402
01403 real (kind=dbl_kind), dimension(icells,nilyr+nslyr+1)::
01404 kh
01405
01406 real (kind=dbl_kind) ::
01407 ci ,
01408 avg_Tsf ,
01409 ferr ,
01410 Iswabs_tmp ,
01411 Sswabs_tmp ,
01412 etai_tmp
01413
01414 logical (kind=log_kind), dimension (icells) ::
01415 converged
01416
01417 logical (kind=log_kind) ::
01418 all_converged
01419
01420
01421
01422
01423
01424 all_converged = .false.
01425
01426 do ij = 1, icells
01427
01428 converged (ij) = .false.
01429 l_snow (ij) = .false.
01430 l_cold (ij) = .true.
01431 fcondbot (ij) = c0
01432 dTsf_prev (ij) = c0
01433 dTi1_prev (ij) = c0
01434 dfsens_dT (ij) = c0
01435 dflat_dT (ij) = c0
01436 dflwout_dT(ij) = c0
01437 dt_rhoi_hlyr(ij) = dt / (rhoi*hilyr(ij))
01438 if (hslyr(ij) > hs_min/real(nslyr,kind=dbl_kind))
01439 l_snow(ij) = .true.
01440 enddo
01441
01442 do k = 1, nslyr
01443 do ij = 1, icells
01444 Tsn_init (ij,k) = Tsn(ij,k)
01445 Tsn_start(ij,k) = Tsn(ij,k)
01446 if (l_snow(ij)) then
01447 etas(ij,k) = dt/(rhos*cp_ice*hslyr(ij))
01448 else
01449 etas(ij,k) = c0
01450 endif
01451 enddo
01452 enddo
01453
01454 do k = 1, nilyr
01455 do ij = 1, icells
01456 Tin_init (ij,k) = Tin(ij,k)
01457 Tin_start(ij,k) = Tin(ij,k)
01458 enddo
01459 enddo
01460
01461
01462
01463
01464
01465
01466
01467
01468 call conductivity (nx_block, ny_block, &
01469 l_snow, icells, &
01470 indxi, indxj, indxij, &
01471 hilyr, hslyr, &
01472 Tin, kh )
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482 do k = 1, nilyr
01483 do ij = 1, icells
01484 i = indxi(ij)
01485 j = indxj(ij)
01486
01487 if (Tin_init(ij,k) <= (Tmlt(k) - p01)) then
01488
01489 if (l_brine) then
01490 ci = cp_ice - Lfresh*Tmlt(k) / &
01491 max( Tin_init(ij,k)**2, Tmlt(k)**2 )
01492 etai_tmp = dt_rhoi_hlyr(ij) / ci
01493 Iswabs_tmp = min(Iswabs(i,j,k), &
01494 (Tmlt(k)-Tin_init(ij,k))/etai_tmp)
01495 else
01496 ci = cp_ice
01497 etai_tmp = dt_rhoi_hlyr(ij) / ci
01498 Iswabs_tmp = min(Iswabs(i,j,k), &
01499 Tin_init(ij,k)/etai_tmp)
01500 endif
01501
01502 fswsfc(i,j) = fswsfc(i,j) + (Iswabs(i,j,k) - Iswabs_tmp)
01503 fswint(i,j) = fswint(i,j) - (Iswabs(i,j,k) - Iswabs_tmp)
01504 Iswabs(i,j,k) = Iswabs_tmp
01505
01506 else
01507
01508 fswsfc(i,j) = fswsfc(i,j) + Iswabs(i,j,k)
01509 fswint(i,j) = fswint(i,j) - Iswabs(i,j,k)
01510 Iswabs(i,j,k) = c0
01511
01512 endif
01513 enddo
01514 enddo
01515
01516 do k = 1, nslyr
01517 do ij = 1, icells
01518 if (l_snow(ij)) then
01519 i = indxi(ij)
01520 j = indxj(ij)
01521
01522 if (Tsn_init(ij,k) <= -p01) then
01523
01524 Sswabs_tmp = min(Sswabs(i,j,k), &
01525 -0.9_dbl_kind*Tsn_init(ij,k)/etas(ij,k))
01526
01527 fswsfc(i,j) = fswsfc(i,j) + (Sswabs(i,j,k) - Sswabs_tmp)
01528 fswint(i,j) = fswint(i,j) - (Sswabs(i,j,k) - Sswabs_tmp)
01529 Sswabs(i,j,k) = Sswabs_tmp
01530
01531 else
01532
01533 fswsfc(i,j) = fswsfc(i,j) + Sswabs(i,j,k)
01534 fswint(i,j) = fswint(i,j) - Sswabs(i,j,k)
01535 Sswabs(i,j,k) = c0
01536
01537 endif
01538 endif
01539 enddo
01540 enddo
01541
01542
01543
01544
01545
01546 do ij = 1, icells
01547 i = indxi(ij)
01548 j = indxj(ij)
01549 fswabsn(i,j) = fswsfc(i,j) + fswint(i,j) + fswthrun(i,j)
01550 enddo
01551
01552
01553
01554
01555
01556
01557 do niter = 1, nitermax
01558
01559
01560
01561
01562
01563 if (all_converged) then
01564 exit
01565 else
01566 isolve = 0
01567 do ij = 1, icells
01568 i = indxi(ij)
01569 j = indxj(ij)
01570 if (.not.converged(ij)) then
01571 isolve = isolve + 1
01572 indxii(isolve) = i
01573 indxjj(isolve) = j
01574 indxij(isolve) = ij
01575 endif
01576 enddo
01577 endif
01578
01579
01580
01581
01582
01583 allocate( sbdiag(isolve,nilyr+nslyr+1))
01584 allocate( diag(isolve,nilyr+nslyr+1))
01585 allocate( spdiag(isolve,nilyr+nslyr+1))
01586 allocate( rhs(isolve,nilyr+nslyr+1))
01587 allocate( Tmat(isolve,nilyr+nslyr+1))
01588 allocate( etai(isolve,nilyr))
01589 allocate(Tsf_start(isolve))
01590 allocate( dTsf(isolve))
01591 allocate(dfsurf_dT(isolve))
01592 allocate( avg_Tsi(isolve))
01593 allocate( enew(isolve))
01594 allocate( dTi1(isolve))
01595
01596 all_converged = .true.
01597
01598 do ij = 1, isolve
01599 m = indxij(ij)
01600 converged(m) = .true.
01601 dfsurf_dT(ij) = c0
01602 avg_Tsi (ij) = c0
01603 enew (ij) = c0
01604 enddo
01605
01606
01607
01608
01609
01610
01611
01612
01613 do k = 1, nilyr
01614 do ij = 1, isolve
01615 m = indxij(ij)
01616 i = indxii(ij)
01617 j = indxjj(ij)
01618
01619 if (l_brine) then
01620 ci = cp_ice - Lfresh*Tmlt(k) / &
01621 (Tin(m,k)*Tin_init(m,k))
01622 else
01623 ci = cp_ice
01624 endif
01625 etai(ij,k) = dt_rhoi_hlyr(m) / ci
01626
01627 enddo
01628 enddo
01629
01630 if (calc_Tsfc) then
01631
01632
01633
01634
01635
01636
01637 call surface_fluxes(nx_block, ny_block, &
01638 isolve, icells, &
01639 indxii, indxjj, indxij, &
01640 Tsf, fswsfc, &
01641 rhoa, flw, &
01642 potT, Qa, &
01643 shcoef, lhcoef, &
01644 flwoutn, fsensn, &
01645 flatn, fsurfn, &
01646 dflwout_dT, dfsens_dT, &
01647 dflat_dT, dfsurf_dT)
01648
01649
01650
01651
01652 do ij = 1, isolve
01653 i = indxii(ij)
01654 j = indxjj(ij)
01655 m = indxij(ij)
01656
01657
01658
01659
01660
01661
01662
01663 if (l_snow(m)) then
01664 fcondtopn(i,j) = kh(m,1) * (Tsf(m) - Tsn(m,1))
01665 else
01666 fcondtopn(i,j) = kh(m,1+nslyr) * (Tsf(m) - Tin(m,1))
01667 endif
01668
01669 if (fsurfn(i,j) < fcondtopn(i,j)) &
01670 Tsf(m) = min (Tsf(m), -puny)
01671
01672
01673
01674
01675
01676 Tsf_start(ij) = Tsf(m)
01677
01678 if (Tsf(m) <= -puny) then
01679 l_cold(m) = .true.
01680 else
01681 l_cold(m) = .false.
01682 endif
01683 enddo
01684
01685
01686
01687
01688
01689 call get_matrix_elements_calc_Tsfc &
01690 (nx_block, ny_block, &
01691 isolve, icells, &
01692 indxii, indxjj, indxij, &
01693 l_snow, l_cold, &
01694 Tsf, Tbot, &
01695 fsurfn, dfsurf_dT, &
01696 Tin_init, Tsn_init, &
01697 kh, Sswabs(:,:,:), &
01698 Iswabs, &
01699 etai, etas, &
01700 sbdiag, diag, &
01701 spdiag, rhs)
01702
01703 else
01704 call get_matrix_elements_know_Tsfc &
01705 (nx_block, ny_block, &
01706 isolve, icells, &
01707 indxii, indxjj, indxij, &
01708 l_snow, Tbot, &
01709 Tin_init, Tsn_init, &
01710 kh, Sswabs(:,:,:), &
01711 Iswabs, &
01712 etai, etas, &
01713 sbdiag, diag, &
01714 spdiag, rhs, &
01715 fcondtopn)
01716 endif
01717
01718
01719
01720
01721
01722 call tridiag_solver (nx_block, ny_block, &
01723 isolve, icells, &
01724 indxii, indxjj, &
01725 nmat, sbdiag, &
01726 diag, spdiag, &
01727 rhs, Tmat)
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755
01756 if (calc_Tsfc) then
01757
01758
01759
01760
01761 do ij = 1, isolve
01762 m = indxij(ij)
01763
01764
01765
01766
01767
01768 if (l_cold(m)) then
01769 if (l_snow(m)) then
01770 Tsf(m) = Tmat(ij,1)
01771 else
01772 Tsf(m) = Tmat(ij,1+nslyr)
01773 endif
01774 else
01775 Tsf(m) = c0
01776 endif
01777
01778
01779
01780
01781
01782
01783
01784
01785 dTsf(ij) = Tsf(m) - Tsf_start(ij)
01786 avg_Tsf = c0
01787
01788
01789
01790
01791
01792
01793
01794 if (Tsf(m) > puny) then
01795 Tsf(m) = c0
01796 dTsf(ij) = -Tsf_start(ij)
01797 if (l_brine) avg_Tsi(ij) = c1
01798 converged(m) = .false.
01799 all_converged = .false.
01800
01801
01802
01803
01804
01805
01806 elseif (niter > 1 &
01807 .and. Tsf_start(ij) <= -puny &
01808 .and. abs(dTsf(ij)) > puny &
01809 .and. abs(dTsf_prev(m)) > puny &
01810 .and. -dTsf(ij)/(dTsf_prev(m)+puny*puny) > p5) then
01811
01812 if (l_brine) then
01813 avg_Tsf = c1
01814 avg_Tsi(ij) = c1
01815 endif
01816 dTsf(ij) = p5 * dTsf(ij)
01817 converged(m) = .false.
01818 all_converged = .false.
01819 endif
01820
01821
01822
01823
01824
01825
01826
01827 Tsf(m) = Tsf(m) &
01828 + avg_Tsf * p5 * (Tsf_start(ij) - Tsf(m))
01829
01830 enddo
01831
01832 endif
01833
01834 do k = 1, nslyr
01835
01836
01837
01838 do ij = 1, isolve
01839 m = indxij(ij)
01840
01841
01842
01843
01844
01845 if (l_snow(m)) then
01846 Tsn(m,k) = Tmat(ij,k+1)
01847 else
01848 Tsn(m,k) = c0
01849 endif
01850 if (l_brine) Tsn(m,k) = min(Tsn(m,k), c0)
01851
01852
01853
01854
01855
01856 Tsn(m,k) = Tsn(m,k) &
01857 + avg_Tsi(ij)*p5*(Tsn_start(m,k)-Tsn(m,k))
01858
01859
01860
01861
01862 qsn(m,k) = -rhos * (Lfresh - cp_ice*Tsn(m,k))
01863 enew(ij) = enew(ij) + hslyr(m) * qsn(m,k)
01864
01865 Tsn_start(m,k) = Tsn(m,k)
01866
01867 enddo
01868 enddo
01869
01870 do k = 1, nilyr
01871
01872
01873
01874 do ij = 1, isolve
01875 m = indxij(ij)
01876
01877
01878
01879
01880 Tin(m,k) = Tmat(ij,k+1+nslyr)
01881 if (l_brine) Tin(m,k) = min(Tin(m,k), Tmlt(k))
01882
01883
01884
01885
01886
01887
01888 if (k==1 .and. .not.calc_Tsfc) then
01889 dTi1(ij) = Tin(m,k) - Tin_start(m,k)
01890
01891 if (niter > 1 &
01892 .and. abs(dTi1(ij)) > puny &
01893 .and. abs(dTi1_prev(m)) > puny &
01894 .and. -dTi1(ij)/(dTi1_prev(m)+puny*puny) > p5) then
01895
01896 if (l_brine) avg_Tsi(ij) = c1
01897 dTi1(ij) = p5 * dTi1(ij)
01898 converged(m) = .false.
01899 all_converged = .false.
01900 endif
01901 dTi1_prev(m) = dTi1(ij)
01902 endif
01903
01904
01905
01906
01907
01908 Tin(m,k) = Tin(m,k) &
01909 + avg_Tsi(ij)*p5*(Tin_start(m,k)-Tin(m,k))
01910
01911
01912
01913
01914 qin(m,k) = -rhoi * (cp_ice*(Tmlt(k)-Tin(m,k)) &
01915 + Lfresh*(c1-Tmlt(k)/Tin(m,k)) &
01916 - cp_ocn*Tmlt(k))
01917 enew(ij) = enew(ij) + hilyr(m) * qin(m,k)
01918
01919 Tin_start(m,k) = Tin(m,k)
01920
01921 enddo
01922 enddo
01923
01924 if (calc_Tsfc) then
01925
01926
01927
01928
01929 do ij = 1, isolve
01930 i = indxii(ij)
01931 j = indxjj(ij)
01932 m = indxij(ij)
01933
01934
01935
01936
01937
01938 if (abs(dTsf(ij)) > Tsf_errmax) then
01939 converged(m) = .false.
01940 all_converged = .false.
01941 endif
01942
01943
01944
01945
01946
01947 fsurfn(i,j) = fsurfn(i,j) + dTsf(ij)*dfsurf_dT(ij)
01948 if (l_snow(m)) then
01949 fcondtopn(i,j) = kh(m,1) * (Tsf(m)-Tsn(m,1))
01950 else
01951 fcondtopn(i,j) = kh(m,1+nslyr) * (Tsf(m)-Tin(m,1))
01952 endif
01953
01954 if (Tsf(m) > -puny .and. fsurfn(i,j) < fcondtopn(i,j)) then
01955 converged(m) = .false.
01956 all_converged = .false.
01957 endif
01958
01959 dTsf_prev(m) = dTsf(ij)
01960
01961 enddo
01962 endif
01963
01964
01965
01966
01967
01968
01969 do ij = 1, isolve
01970 i = indxii(ij)
01971 j = indxjj(ij)
01972 m = indxij(ij)
01973
01974 fcondbot(m) = kh(m,1+nslyr+nilyr) * &
01975 (Tin(m,nilyr) - Tbot(i,j))
01976
01977 ferr = abs( (enew(ij)-einit(m))/dt &
01978 - (fcondtopn(i,j) - fcondbot(m) + fswint(i,j)) )
01979
01980
01981 if (ferr > 0.9_dbl_kind*ferrmax) then
01982 converged(m) = .false.
01983 all_converged = .false.
01984 endif
01985
01986 enddo
01987
01988 deallocate(sbdiag)
01989 deallocate(diag)
01990 deallocate(spdiag)
01991 deallocate(rhs)
01992 deallocate(Tmat)
01993 deallocate(etai)
01994 deallocate(Tsf_start)
01995 deallocate(dTsf)
01996 deallocate(dfsurf_dT)
01997 deallocate(avg_Tsi)
01998 deallocate(enew)
01999 deallocate(dTi1)
02000
02001 enddo
02002
02003 if (.not.all_converged) then
02004
02005 do ij = 1, icells
02006 i = indxi(ij)
02007 j = indxj(ij)
02008
02009
02010
02011
02012 if (.not.converged(ij)) then
02013 write(nu_diag,*) 'Thermo iteration does not converge,', &
02014 'istep1, my_task, i, j:', &
02015 istep1, my_task, i, j
02016 write(nu_diag,*) 'Ice thickness:', hilyr(ij)*nilyr
02017 write(nu_diag,*) 'Snow thickness:', hslyr(ij)*nslyr
02018 write(nu_diag,*) 'dTsf, Tsf_errmax:',dTsf_prev(ij), &
02019 Tsf_errmax
02020 write(nu_diag,*) 'Tsf:', Tsf(ij)
02021 write(nu_diag,*) 'fsurf:', fsurfn(i,j)
02022 write(nu_diag,*) 'fcondtop, fcondbot, fswint', &
02023 fcondtopn(i,j), fcondbot(ij), fswint(i,j)
02024 write(nu_diag,*) 'fswsfc, fswthrun', &
02025 fswsfc(i,j), fswthrun(i,j)
02026 write(nu_diag,*) 'Flux conservation error =', ferr
02027 write(nu_diag,*) 'Internal snow absorption:'
02028 write(nu_diag,*) (Sswabs(i,j,k),k=1,nslyr)
02029 write(nu_diag,*) 'Internal ice absorption:'
02030 write(nu_diag,*) (Iswabs(i,j,k),k=1,nilyr)
02031 write(nu_diag,*) 'Initial snow temperatures:'
02032 write(nu_diag,*) (Tsn_init(ij,k),k=1,nslyr)
02033 write(nu_diag,*) 'Initial ice temperatures:'
02034 write(nu_diag,*) (Tin_init(ij,k),k=1,nilyr)
02035 write(nu_diag,*) 'Final snow temperatures:'
02036 write(nu_diag,*) (Tsn(ij,k),k=1,nslyr)
02037 write(nu_diag,*) 'Final ice temperatures:'
02038 write(nu_diag,*) (Tin(ij,k),k=1,nilyr)
02039 l_stop = .true.
02040 istop = i
02041 jstop = j
02042 return
02043 endif
02044 enddo
02045 endif
02046
02047 if (calc_Tsfc) then
02048
02049
02050
02051 do ij = 1, icells
02052 i = indxi(ij)
02053 j = indxj(ij)
02054
02055
02056 flwoutn(i,j) = flwoutn(i,j) + dTsf_prev(ij) * dflwout_dT(ij)
02057 fsensn(i,j) = fsensn(i,j) + dTsf_prev(ij) * dfsens_dT(ij)
02058 flatn(i,j) = flatn(i,j) + dTsf_prev(ij) * dflat_dT(ij)
02059
02060 enddo
02061 endif
02062
02063 end subroutine temperature_changes
02064
02065
02066
02067
02068
02069
02070
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081
02082
02083
02084 subroutine conductivity (nx_block, ny_block, &
02085 l_snow, icells, &
02086 indxi, indxj, indxij, &
02087 hilyr, hslyr, &
02088 Tin, kh)
02089
02090
02091
02092
02093
02094 integer (kind=int_kind), intent(in) ::
02095 nx_block, ny_block,
02096 icells
02097
02098 logical (kind=log_kind), dimension(icells),
02099 intent(in) ::
02100 l_snow
02101
02102 integer (kind=int_kind), dimension(nx_block*ny_block),
02103 intent(in) ::
02104 indxi, indxj
02105
02106 integer (kind=int_kind), dimension (icells),
02107 intent(in) ::
02108 indxij
02109
02110 real (kind=dbl_kind), dimension (icells), intent(in) ::
02111 hilyr ,
02112 hslyr
02113
02114 real (kind=dbl_kind), dimension (icells,nilyr),
02115 intent(in) ::
02116 Tin
02117
02118 real (kind=dbl_kind), dimension (icells,nilyr+nslyr+1),
02119 intent(out) ::
02120 kh
02121
02122
02123
02124 integer (kind=int_kind) ::
02125 i, j ,
02126 ij, m ,
02127 k
02128
02129 real (kind=dbl_kind), dimension (icells,nilyr) ::
02130 kilyr
02131
02132 real (kind=dbl_kind), dimension (icells,nslyr) ::
02133 kslyr
02134
02135
02136 do k = 1, nslyr
02137 do ij = 1, icells
02138 kslyr(ij,k) = ksno
02139 enddo
02140 enddo
02141
02142
02143 do k = 1, nilyr
02144
02145
02146
02147 do ij = 1, icells
02148 kilyr(ij,k) = kice + betak*salin(k)/min(-puny,Tin(ij,k))
02149 kilyr(ij,k) = max (kilyr(ij,k), kimin)
02150 enddo
02151 enddo
02152
02153
02154 do ij = 1, icells
02155
02156 if (l_snow(ij)) then
02157 kh(ij,1) = c2 * kslyr(ij,1) / hslyr(ij)
02158 kh(ij,1+nslyr) = c2 * kslyr(ij,nslyr) * kilyr(ij,1) / &
02159 ( kslyr(ij,nslyr)*hilyr(ij) + &
02160 kilyr(ij,1 )*hslyr(ij) )
02161 else
02162 kh(ij,1) = c0
02163 kh(ij,1+nslyr) = c2 * kilyr(ij,1) / hilyr(ij)
02164 endif
02165
02166
02167 kh(ij,1+nslyr+nilyr) = c2 * kilyr(ij,nilyr) / hilyr(ij)
02168
02169 enddo
02170
02171
02172
02173 if (nslyr > 1) then
02174 do k = 2, nslyr
02175 do ij = 1, icells
02176 if (l_snow(ij)) then
02177 kh(ij,k) = c2 * kslyr(ij,k-1) * kslyr(ij,k) / &
02178 ((kslyr(ij,k-1) + kslyr(ij,k))*hslyr(ij))
02179 else
02180 kh(ij,k) = c0
02181 endif
02182 enddo
02183 enddo
02184 endif
02185
02186
02187 do k = 2, nilyr
02188 do ij = 1, icells
02189 kh(ij,k+nslyr) = c2 * kilyr(ij,k-1) * kilyr(ij,k) / &
02190 ((kilyr(ij,k-1) + kilyr(ij,k))*hilyr(ij))
02191 enddo
02192 enddo
02193
02194 end subroutine conductivity
02195
02196
02197
02198
02199
02200
02201
02202
02203
02204
02205
02206
02207
02208
02209
02210
02211
02212
02213 subroutine surface_fluxes (nx_block, ny_block, &
02214 isolve, icells, &
02215 indxii, indxjj, indxij, &
02216 Tsf, fswsfc, &
02217 rhoa, flw, &
02218 potT, Qa, &
02219 shcoef, lhcoef, &
02220 flwoutn, fsensn, &
02221 flatn, fsurfn, &
02222 dflwout_dT, dfsens_dT, &
02223 dflat_dT, dfsurf_dT)
02224
02225
02226
02227
02228
02229 integer (kind=int_kind), intent(in) ::
02230 nx_block, ny_block,
02231 isolve ,
02232 icells
02233
02234 integer (kind=int_kind), dimension(icells),
02235 intent(in) ::
02236 indxii, indxjj
02237
02238 integer (kind=int_kind), dimension (icells) ::
02239 indxij
02240
02241 real (kind=dbl_kind), dimension (icells), intent(in) ::
02242 Tsf
02243
02244 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
02245 fswsfc ,
02246 rhoa ,
02247 flw ,
02248 potT ,
02249 Qa ,
02250 shcoef ,
02251 lhcoef
02252
02253 real (kind=dbl_kind), dimension (nx_block,ny_block),
02254 intent(inout) ::
02255 fsensn ,
02256 flatn ,
02257 flwoutn ,
02258 fsurfn
02259
02260 real (kind=dbl_kind), dimension (icells),
02261 intent(inout) ::
02262 dfsens_dT ,
02263 dflat_dT ,
02264 dflwout_dT
02265
02266 real (kind=dbl_kind), dimension (isolve),
02267 intent(inout) ::
02268 dfsurf_dT
02269
02270
02271
02272 integer (kind=int_kind) ::
02273 i, j ,
02274 ij, m ,
02275 k
02276
02277 real (kind=dbl_kind) ::
02278 TsfK ,
02279 Qsfc ,
02280 dQsfcdT ,
02281 qsat ,
02282 flwdabs ,
02283 tmpvar
02284
02285
02286
02287
02288 do ij = 1, isolve
02289 i = indxii(ij)
02290 j = indxjj(ij)
02291 m = indxij(ij)
02292
02293
02294 TsfK = Tsf(m) + Tffresh
02295 tmpvar = c1/TsfK
02296
02297
02298 qsat = qqqice * exp(-TTTice*tmpvar)
02299 Qsfc = qsat / rhoa(i,j)
02300 dQsfcdT = TTTice * tmpvar*tmpvar * Qsfc
02301
02302
02303 flwdabs = emissivity * flw(i,j)
02304 flwoutn(i,j) = -emissivity * stefan_boltzmann * TsfK**4
02305
02306
02307 fsensn(i,j) = shcoef(i,j) * (potT(i,j) - TsfK)
02308 flatn(i,j) = lhcoef(i,j) * (Qa(i,j) - Qsfc)
02309
02310
02311 dflwout_dT(m) = - emissivity*stefan_boltzmann * c4*TsfK**3
02312 dfsens_dT(m) = - shcoef(i,j)
02313 dflat_dT(m) = - lhcoef(i,j) * dQsfcdT
02314
02315 fsurfn(i,j) = fswsfc(i,j) + flwdabs + flwoutn(i,j) &
02316 + fsensn(i,j) + flatn(i,j)
02317 dfsurf_dT(ij) = dflwout_dT(m) &
02318 + dfsens_dT(m) + dflat_dT(m)
02319
02320 enddo
02321
02322 end subroutine surface_fluxes
02323
02324
02325
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338
02339
02340
02341
02342
02343
02344
02345
02346 subroutine get_matrix_elements_calc_Tsfc &
02347 (nx_block, ny_block, &
02348 isolve, icells, &
02349 indxii, indxjj, indxij, &
02350 l_snow, l_cold, &
02351 Tsf, Tbot, &
02352 fsurfn, dfsurf_dT, &
02353 Tin_init, Tsn_init, &
02354 kh, Sswabs, &
02355 Iswabs, &
02356 etai, etas, &
02357 sbdiag, diag, &
02358 spdiag, rhs)
02359
02360
02361
02362
02363
02364 integer (kind=int_kind), intent(in) ::
02365 nx_block, ny_block,
02366 isolve ,
02367 icells
02368
02369 integer (kind=int_kind), dimension(icells),
02370 intent(in) ::
02371 indxii, indxjj
02372
02373 integer (kind=int_kind), dimension (icells),
02374 intent(in) ::
02375 indxij
02376
02377 logical (kind=log_kind), dimension (icells),
02378 intent(in) ::
02379 l_snow ,
02380 l_cold
02381
02382 real (kind=dbl_kind), dimension (icells), intent(in) ::
02383 Tsf
02384
02385 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
02386 fsurfn ,
02387 Tbot
02388
02389 real (kind=dbl_kind), dimension (isolve), intent(in) ::
02390 dfsurf_dT
02391
02392 real (kind=dbl_kind), dimension (isolve,nilyr),
02393 intent(in) ::
02394 etai
02395
02396 real (kind=dbl_kind), dimension (icells,nilyr),
02397 intent(in) ::
02398 Tin_init
02399
02400 real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr),
02401 intent(in) ::
02402 Sswabs
02403
02404 real (kind=dbl_kind), dimension (nx_block,ny_block,nilyr),
02405 intent(in) ::
02406 Iswabs
02407
02408 real (kind=dbl_kind), dimension (icells,nslyr),
02409 intent(in) ::
02410 etas ,
02411 Tsn_init
02412
02413
02414 real (kind=dbl_kind), dimension (icells,nslyr+nilyr+1),
02415 intent(in) ::
02416 kh
02417
02418 real (kind=dbl_kind), dimension (isolve,nslyr+nilyr+1),
02419 intent(inout) ::
02420 sbdiag ,
02421 diag ,
02422 spdiag ,
02423 rhs
02424
02425
02426
02427 integer (kind=int_kind) ::
02428 i, j ,
02429 ij, m ,
02430 k, ks, ki, kr
02431
02432
02433
02434
02435
02436
02437
02438
02439 do k = 1, nslyr+1
02440 do ij = 1, isolve
02441 sbdiag(ij,k) = c0
02442 diag (ij,k) = c1
02443 spdiag(ij,k) = c0
02444 rhs (ij,k) = c0
02445 enddo
02446 enddo
02447
02448
02449
02450
02451
02452
02453
02454
02455
02456
02457
02458 do ij = 1, isolve
02459 i = indxii(ij)
02460 j = indxjj(ij)
02461 m = indxij(ij)
02462
02463
02464
02465
02466 if (l_cold(m)) then
02467 if (l_snow(m)) then
02468 k = 1
02469 else
02470 k = 1 + nslyr
02471 endif
02472 kr = k
02473
02474 sbdiag(ij,kr) = c0
02475 diag (ij,kr) = dfsurf_dT(ij) - kh(m,k)
02476 spdiag(ij,kr) = kh(m,k)
02477 rhs (ij,kr) = dfsurf_dT(ij)*Tsf(m) - fsurfn(i,j)
02478
02479 endif
02480
02481
02482
02483
02484
02485
02486
02487 if (l_snow(m)) then
02488 if (l_cold(m)) then
02489 sbdiag(ij,2) = -etas(m,1) * kh(m,1)
02490 spdiag(ij,2) = -etas(m,1) * kh(m,2)
02491 diag (ij,2) = c1 &
02492 + etas(m,1) * (kh(m,1) + kh(m,2))
02493 rhs (ij,2) = Tsn_init(m,1) &
02494 + etas(m,1) * Sswabs(i,j,1)
02495 else
02496 sbdiag(ij,2) = c0
02497 spdiag(ij,2) = -etas(m,1) * kh(m,2)
02498 diag (ij,2) = c1 &
02499 + etas(m,1) * (kh(m,1) + kh(m,2))
02500 rhs (ij,2) = Tsn_init(m,1) &
02501 + etas(m,1)*kh(m,1)*Tsf(m) &
02502 + etas(m,1) * Sswabs(i,j,1)
02503 endif
02504 endif
02505
02506 enddo
02507
02508
02509
02510
02511
02512 if (nslyr > 1) then
02513
02514 do k = 2, nslyr
02515 kr = k + 1
02516
02517 do ij = 1, isolve
02518 i = indxii(ij)
02519 j = indxjj(ij)
02520 m = indxij(ij)
02521
02522 if (l_snow(m)) then
02523 sbdiag(ij,kr) = -etas(m,k) * kh(m,k)
02524 spdiag(ij,kr) = -etas(m,k) * kh(m,k+1)
02525 diag (ij,kr) = c1 &
02526 + etas(m,k) * (kh(m,k) + kh(m,k+1))
02527 rhs (ij,kr) = Tsn_init(m,k) &
02528 + etas(m,k) * Sswabs(i,j,k)
02529 endif
02530 enddo
02531 enddo
02532
02533 endif
02534
02535
02536 if (nilyr > 1) then
02537
02538
02539
02540
02541
02542 ki = 1
02543 k = ki + nslyr
02544 kr = k + 1
02545
02546 do ij = 1, isolve
02547 i = indxii(ij)
02548 j = indxjj(ij)
02549 m = indxij(ij)
02550
02551 if (l_snow(m) .or. l_cold(m)) then
02552 sbdiag(ij,kr) = -etai(ij,ki) * kh(m,k)
02553 spdiag(ij,kr) = -etai(ij,ki) * kh(m,k+1)
02554 diag (ij,kr) = c1 &
02555 + etai(ij,ki) * (kh(m,k) + kh(m,k+1))
02556 rhs (ij,kr) = Tin_init(m,ki) &
02557 + etai(ij,ki)*Iswabs(i,j,ki)
02558 else
02559 sbdiag(ij,kr) = c0
02560 spdiag(ij,kr) = -etai(ij,ki) * kh(m,k+1)
02561 diag (ij,kr) = c1 &
02562 + etai(ij,ki) * (kh(m,k) + kh(m,k+1))
02563 rhs (ij,kr) = Tin_init(m,ki) &
02564 + etai(ij,ki)*Iswabs(i,j,ki) &
02565 + etai(ij,ki)*kh(m,k)*Tsf(m)
02566 endif
02567
02568 enddo
02569
02570
02571
02572
02573
02574 ki = nilyr
02575 k = ki + nslyr
02576 kr = k + 1
02577
02578 do ij = 1, isolve
02579 i = indxii(ij)
02580 j = indxjj(ij)
02581 m = indxij(ij)
02582 sbdiag(ij,kr) = -etai(ij,ki) * kh(m,k)
02583 spdiag(ij,kr) = c0
02584 diag (ij,kr) = c1 &
02585 + etai(ij,ki) * (kh(m,k) + kh(m,k+1))
02586 rhs (ij,kr) = Tin_init(m,ki) &
02587 + etai(ij,ki)*Iswabs(i,j,ki) &
02588 + etai(ij,ki)*kh(m,k+1)*Tbot(i,j)
02589
02590 enddo
02591
02592 else
02593
02594
02595
02596
02597
02598 ki = 1
02599 k = ki + nslyr
02600 kr = k + 1
02601
02602 do ij = 1, isolve
02603 i = indxii(ij)
02604 j = indxjj(ij)
02605 m = indxij(ij)
02606
02607 if (l_snow(m) .or. l_cold(m)) then
02608 sbdiag(ij,kr) = -etai(ij,ki) * kh(m,k)
02609 spdiag(ij,kr) = c0
02610 diag (ij,kr) = c1 &
02611 + etai(ij,ki) * (kh(m,k) + kh(m,k+1))
02612 rhs (ij,kr) = Tin_init(m,ki) &
02613 + etai(ij,ki) * Iswabs(i,j,ki) &
02614 + etai(ij,ki) * kh(m,k+1)*Tbot(i,j)
02615 else
02616 sbdiag(ij,kr) = c0
02617 spdiag(ij,kr) = c0
02618 diag (ij,kr) = c1 &
02619 + etai(ij,ki) * (kh(m,k) + kh(m,k+1))
02620 rhs (ij,kr) = Tin_init(m,ki) &
02621 + etai(ij,ki) * Iswabs(i,j,ki) &
02622 + etai(ij,ki) * kh(m,k)*Tsf(m) &
02623 + etai(ij,ki) * kh(m,k+1)*Tbot(i,j)
02624 endif
02625 enddo
02626
02627 endif
02628
02629
02630
02631
02632
02633 do ki = 2, nilyr-1
02634
02635 k = ki + nslyr
02636 kr = k + 1
02637 do ij = 1, isolve
02638 i = indxii(ij)
02639 j = indxjj(ij)
02640 m = indxij(ij)
02641
02642 sbdiag(ij,kr) = -etai(ij,ki) * kh(m,k)
02643 spdiag(ij,kr) = -etai(ij,ki) * kh(m,k+1)
02644 diag (ij,kr) = c1 &
02645 + etai(ij,ki) * (kh(m,k) + kh(m,k+1))
02646 rhs (ij,kr) = Tin_init(m,ki) &
02647 + etai(ij,ki)*Iswabs(i,j,ki)
02648
02649 enddo
02650 enddo
02651
02652 end subroutine get_matrix_elements_calc_Tsfc
02653
02654
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664
02665
02666
02667
02668
02669
02670
02671
02672
02673
02674
02675
02676 subroutine get_matrix_elements_know_Tsfc &
02677 (nx_block, ny_block, &
02678 isolve, icells, &
02679 indxii, indxjj, indxij, &
02680 l_snow, Tbot, &
02681 Tin_init, Tsn_init, &
02682 kh, Sswabs, &
02683 Iswabs, &
02684 etai, etas, &
02685 sbdiag, diag, &
02686 spdiag, rhs, &
02687 fcondtopn)
02688
02689
02690
02691
02692
02693 integer (kind=int_kind), intent(in) ::
02694 nx_block, ny_block,
02695 isolve ,
02696 icells
02697
02698 integer (kind=int_kind), dimension(icells),
02699 intent(in) ::
02700 indxii, indxjj
02701
02702 integer (kind=int_kind), dimension (icells),
02703 intent(in) ::
02704 indxij
02705
02706 logical (kind=log_kind), dimension (icells),
02707 intent(in) ::
02708 l_snow
02709
02710 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
02711 Tbot
02712
02713 real (kind=dbl_kind), dimension (isolve,nilyr),
02714 intent(in) ::
02715 etai
02716
02717 real (kind=dbl_kind), dimension (icells,nilyr),
02718 intent(in) ::
02719 Tin_init
02720
02721 real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr),
02722 intent(in) ::
02723 Sswabs
02724
02725 real (kind=dbl_kind), dimension (nx_block,ny_block,nilyr),
02726 intent(in) ::
02727 Iswabs
02728
02729 real (kind=dbl_kind), dimension (icells,nslyr),
02730 intent(in) ::
02731 etas ,
02732 Tsn_init
02733
02734
02735 real (kind=dbl_kind), dimension (icells,nslyr+nilyr+1),
02736 intent(in) ::
02737 kh
02738
02739 real (kind=dbl_kind), dimension (isolve,nslyr+nilyr+1),
02740 intent(inout) ::
02741 sbdiag ,
02742 diag ,
02743 spdiag ,
02744 rhs
02745
02746 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in),
02747 optional ::
02748 fcondtopn
02749
02750
02751
02752 integer (kind=int_kind) ::
02753 i, j ,
02754 ij, m ,
02755 k, ks, ki, kr
02756
02757
02758
02759
02760
02761
02762
02763
02764 do k = 1, nslyr+1
02765 do ij = 1, isolve
02766 sbdiag(ij,k) = c0
02767 diag (ij,k) = c1
02768 spdiag(ij,k) = c0
02769 rhs (ij,k) = c0
02770 enddo
02771 enddo
02772
02773
02774
02775
02776
02777
02778
02779
02780
02781
02782
02783
02784
02785
02786
02787
02788
02789 do ij = 1, isolve
02790 i = indxii(ij)
02791 j = indxjj(ij)
02792 m = indxij(ij)
02793
02794 if (l_snow(m)) then
02795 sbdiag(ij,2) = c0
02796 spdiag(ij,2) = -etas(m,1) * kh(m,2)
02797 diag (ij,2) = c1 &
02798 + etas(m,1) * kh(m,2)
02799 rhs (ij,2) = Tsn_init(m,1) &
02800 + etas(m,1) * Sswabs(i,j,1) &
02801 + etas(m,1) * fcondtopn(i,j)
02802 endif
02803 enddo
02804
02805
02806
02807
02808
02809 if (nslyr > 1) then
02810
02811 do k = 2, nslyr
02812 kr = k + 1
02813
02814 do ij = 1, isolve
02815 i = indxii(ij)
02816 j = indxjj(ij)
02817 m = indxij(ij)
02818
02819 if (l_snow(m)) then
02820 sbdiag(ij,kr) = -etas(m,k) * kh(m,k)
02821 spdiag(ij,kr) = -etas(m,k) * kh(m,k+1)
02822 diag (ij,kr) = c1 &
02823 + etas(m,k) * (kh(m,k) + kh(m,k+1))
02824 rhs (ij,kr) = Tsn_init(m,k) &
02825 + etas(m,k) * Sswabs(i,j,k)
02826 endif
02827 enddo
02828 enddo
02829
02830 endif
02831
02832
02833 if (nilyr > 1) then
02834
02835
02836
02837
02838
02839 ki = 1
02840 k = ki + nslyr
02841 kr = k + 1
02842
02843 do ij = 1, isolve
02844 i = indxii(ij)
02845 j = indxjj(ij)
02846 m = indxij(ij)
02847
02848 if (l_snow(m)) then
02849
02850 sbdiag(ij,kr) = -etai(ij,ki) * kh(m,k)
02851 spdiag(ij,kr) = -etai(ij,ki) * kh(m,k+1)
02852 diag (ij,kr) = c1 &
02853 + etai(ij,ki) * (kh(m,k) + kh(m,k+1))
02854 rhs (ij,kr) = Tin_init(m,ki) &
02855 + etai(ij,ki) * Iswabs(i,j,ki)
02856 else
02857 sbdiag(ij,kr) = c0
02858 spdiag(ij,kr) = -etai(ij,ki) * kh(m,k+1)
02859 diag (ij,kr) = c1 &
02860 + etai(ij,ki) * kh(m,k+1)
02861 rhs (ij,kr) = Tin_init(m,ki) &
02862 + etai(ij,ki) * Iswabs(i,j,ki) &
02863 + etai(ij,ki) * fcondtopn(i,j)
02864 endif
02865 enddo
02866
02867
02868
02869
02870
02871 ki = nilyr
02872 k = ki + nslyr
02873 kr = k + 1
02874
02875 do ij = 1, isolve
02876 i = indxii(ij)
02877 j = indxjj(ij)
02878 m = indxij(ij)
02879 sbdiag(ij,kr) = -etai(ij,ki) * kh(m,k)
02880 spdiag(ij,kr) = c0
02881 diag (ij,kr) = c1 &
02882 + etai(ij,ki) * (kh(m,k) + kh(m,k+1))
02883 rhs (ij,kr) = Tin_init(m,ki) &
02884 + etai(ij,ki)*Iswabs(i,j,ki) &
02885 + etai(ij,ki)*kh(m,k+1)*Tbot(i,j)
02886
02887 enddo
02888
02889 else
02890
02891
02892
02893
02894
02895 ki = 1
02896 k = ki + nslyr
02897 kr = k + 1
02898
02899 do ij = 1, isolve
02900 i = indxii(ij)
02901 j = indxjj(ij)
02902 m = indxij(ij)
02903
02904 if (l_snow(m)) then
02905 sbdiag(ij,kr) = -etai(ij,ki) * kh(m,k)
02906 spdiag(ij,kr) = c0
02907 diag (ij,kr) = c1 &
02908 + etai(ij,ki) * (kh(m,k) + kh(m,k+1))
02909 rhs (ij,kr) = Tin_init(m,ki) &
02910 + etai(ij,ki) * Iswabs(i,j,ki) &
02911 + etai(ij,ki) * kh(m,k+1)*Tbot(i,j)
02912 else
02913 sbdiag(ij,kr) = c0
02914 spdiag(ij,kr) = c0
02915 diag (ij,kr) = c1 &
02916 + etai(ij,ki) * kh(m,k+1)
02917 rhs (ij,kr) = Tin_init(m,ki) &
02918 + etai(ij,ki) * Iswabs(i,j,ki) &
02919 + etai(ij,ki) * fcondtopn(i,j) &
02920 + etai(ij,ki) * kh(m,k+1)*Tbot(i,j)
02921 endif
02922 enddo
02923
02924 endif
02925
02926
02927
02928
02929
02930 do ki = 2, nilyr-1
02931
02932 k = ki + nslyr
02933 kr = k + 1
02934 do ij = 1, isolve
02935 i = indxii(ij)
02936 j = indxjj(ij)
02937 m = indxij(ij)
02938
02939 sbdiag(ij,kr) = -etai(ij,ki) * kh(m,k)
02940 spdiag(ij,kr) = -etai(ij,ki) * kh(m,k+1)
02941 diag (ij,kr) = c1 &
02942 + etai(ij,ki) * (kh(m,k) + kh(m,k+1))
02943 rhs (ij,kr) = Tin_init(m,ki) &
02944 + etai(ij,ki)*Iswabs(i,j,ki)
02945
02946 enddo
02947 enddo
02948
02949 end subroutine get_matrix_elements_know_Tsfc
02950
02951
02952
02953
02954
02955
02956
02957
02958
02959
02960
02961
02962
02963
02964
02965
02966
02967
02968 subroutine tridiag_solver (nx_block, ny_block, &
02969 isolve, icells, &
02970 indxii, indxjj, &
02971 nmat, sbdiag, &
02972 diag, spdiag, &
02973 rhs, xout)
02974
02975
02976
02977
02978
02979 integer (kind=int_kind), intent(in) ::
02980 nx_block, ny_block,
02981 isolve ,
02982 icells
02983
02984 integer (kind=int_kind), dimension(icells),
02985 intent(in) ::
02986 indxii, indxjj
02987
02988 integer (kind=int_kind), intent(in) ::
02989 nmat
02990
02991 real (kind=dbl_kind), dimension (isolve,nmat),
02992 intent(in) ::
02993 sbdiag ,
02994 diag ,
02995 spdiag ,
02996 rhs
02997
02998 real (kind=dbl_kind), dimension (isolve,nmat),
02999 intent(inout) ::
03000 xout
03001
03002
03003
03004 integer (kind=int_kind) ::
03005 i, j ,
03006 ij ,
03007 k
03008
03009 real (kind=dbl_kind), dimension (isolve) ::
03010 wbeta
03011
03012 real (kind=dbl_kind), dimension(isolve,nilyr+nslyr+1)::
03013 wgamma
03014
03015
03016
03017
03018 do ij = 1, isolve
03019 wbeta(ij) = diag(ij,1)
03020 xout(ij,1) = rhs(ij,1) / wbeta(ij)
03021 enddo
03022
03023 do k = 2, nmat
03024
03025
03026
03027 do ij = 1, isolve
03028 wgamma(ij,k) = spdiag(ij,k-1) / wbeta(ij)
03029 wbeta(ij) = diag(ij,k) - sbdiag(ij,k)*wgamma(ij,k)
03030 xout(ij,k) = (rhs(ij,k) - sbdiag(ij,k)*xout(ij,k-1)) &
03031 / wbeta(ij)
03032 enddo
03033 enddo
03034
03035 do k = nmat-1, 1, -1
03036
03037
03038
03039 do ij = 1, isolve
03040 xout(ij,k) = xout(ij,k) - wgamma(ij,k+1)*xout(ij,k+1)
03041 enddo
03042 enddo
03043
03044 end subroutine tridiag_solver
03045
03046
03047
03048
03049
03050
03051
03052
03053
03054
03055
03056
03057
03058
03059
03060
03061
03062
03063
03064
03065
03066
03067 subroutine zerolayer_temperature(nx_block, ny_block, &
03068 my_task, istep1, &
03069 dt, icells, &
03070 indxi, indxj, &
03071 rhoa, flw, &
03072 potT, Qa, &
03073 shcoef, lhcoef, &
03074 fswsfc, fswthrun, &
03075 hilyr, hslyr, &
03076 Tsf, Tbot, &
03077 fsensn, flatn, &
03078 fswabsn, flwoutn, &
03079 fsurfn, &
03080 fcondtopn,fcondbot, &
03081 l_stop, &
03082 istop, jstop)
03083
03084
03085
03086
03087
03088 integer (kind=int_kind), intent(in) ::
03089 nx_block, ny_block,
03090 my_task ,
03091 istep1 ,
03092 icells
03093
03094 real (kind=dbl_kind), intent(in) ::
03095 dt
03096
03097 integer (kind=int_kind), dimension(nx_block*ny_block),
03098 intent(in) ::
03099 indxi, indxj
03100
03101 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
03102 rhoa ,
03103 flw ,
03104 potT ,
03105 Qa ,
03106 shcoef ,
03107 lhcoef ,
03108 Tbot ,
03109 fswsfc ,
03110 fswthrun
03111
03112 real (kind=dbl_kind), dimension (icells), intent(in) ::
03113 hilyr ,
03114 hslyr
03115
03116 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout)::
03117 fsensn ,
03118 fswabsn ,
03119 flatn ,
03120 flwoutn ,
03121 fsurfn ,
03122 fcondtopn
03123
03124 real (kind=dbl_kind), dimension (icells), intent(out)::
03125 fcondbot
03126
03127 real (kind=dbl_kind), dimension (icells),
03128 intent(inout)::
03129 Tsf
03130
03131 logical (kind=log_kind), intent(inout) ::
03132 l_stop
03133
03134 integer (kind=int_kind), intent(inout) ::
03135 istop, jstop
03136
03137
03138
03139 logical (kind=log_kind), parameter ::
03140 l_zerolayerchecks = .true.
03141
03142 integer (kind=int_kind), parameter ::
03143 nitermax = 50
03144
03145 real (kind=dbl_kind), parameter ::
03146 Tsf_errmax = 5.e-4_dbl_kind
03147
03148
03149 integer (kind=int_kind) ::
03150 i, j ,
03151 ij, m ,
03152 niter
03153
03154 integer (kind=int_kind) ::
03155 isolve
03156
03157 integer (kind=int_kind), dimension (icells) ::
03158 indxii, indxjj
03159
03160 integer (kind=int_kind), dimension (icells) ::
03161 indxij
03162
03163 real (kind=dbl_kind), dimension (:), allocatable ::
03164 Tsf_start ,
03165 dTsf ,
03166 dfsurf_dT
03167
03168 real (kind=dbl_kind), dimension (icells) ::
03169 dTsf_prev ,
03170 dfsens_dT ,
03171 dflat_dT ,
03172 dflwout_dT
03173
03174 real (kind=dbl_kind), dimension (:), allocatable ::
03175 heff ,
03176
03177 kh , &
03178 diag , &
03179 rhs
03180
03181 real (kind=dbl_kind) ::
03182 kratio ,
03183 avg_Tsf
03184
03185 logical (kind=log_kind), dimension (icells) ::
03186 converged
03187
03188 logical (kind=log_kind) ::
03189 all_converged
03190
03191
03192
03193
03194
03195 all_converged = .false.
03196
03197 do ij = 1, icells
03198 fcondbot(ij) = c0
03199
03200 converged (ij) = .false.
03201
03202 dTsf_prev (ij) = c0
03203
03204 enddo
03205
03206
03207
03208
03209
03210
03211
03212 do niter = 1, nitermax
03213
03214 if (all_converged) then
03215 exit
03216 else
03217 isolve = 0
03218 do ij = 1, icells
03219 i = indxi(ij)
03220 j = indxj(ij)
03221 if (.not.converged(ij)) then
03222 isolve = isolve + 1
03223 indxii(isolve) = i
03224 indxjj(isolve) = j
03225 indxij(isolve) = ij
03226 endif
03227 enddo
03228 endif
03229
03230 allocate( diag(isolve))
03231 allocate( rhs(isolve))
03232 allocate( kh(isolve))
03233 allocate( heff(isolve))
03234 allocate(Tsf_start(isolve))
03235 allocate( dTsf(isolve))
03236 allocate(dfsurf_dT(isolve))
03237
03238
03239
03240
03241
03242
03243 call surface_fluxes (nx_block, ny_block, &
03244 isolve, icells, &
03245 indxii, indxjj, indxij, &
03246 Tsf, fswsfc, &
03247 rhoa, flw, &
03248 potT, Qa, &
03249 shcoef, lhcoef, &
03250 flwoutn, fsensn, &
03251 flatn, fsurfn, &
03252 dflwout_dT, dfsens_dT, &
03253 dflat_dT, dfsurf_dT)
03254
03255
03256
03257
03258
03259
03260 kratio = kseaice/ksno
03261
03262 do ij = 1, isolve
03263 m = indxij(ij)
03264
03265 heff(ij) = hilyr(m) + kratio * hslyr(m)
03266 kh(ij) = kseaice / heff(ij)
03267 enddo
03268
03269
03270
03271
03272
03273 do ij = 1, isolve
03274 i = indxii(ij)
03275 j = indxjj(ij)
03276 m = indxij(ij)
03277
03278
03279
03280
03281
03282
03283
03284 fcondtopn(i,j) = kh(ij) * (Tsf(m) - Tbot(i,j))
03285
03286 if (fsurfn(i,j) < fcondtopn(i,j)) &
03287 Tsf(m) = min (Tsf(m), -puny)
03288
03289
03290
03291
03292
03293 Tsf_start(ij) = Tsf(m)
03294
03295 enddo
03296
03297
03298
03299
03300
03301 do ij = 1, isolve
03302 i = indxii(ij)
03303 j = indxjj(ij)
03304 m = indxij(ij)
03305
03306 diag(ij) = dfsurf_dT(ij) - kh(ij)
03307 rhs(ij) = dfsurf_dT(ij)*Tsf(m) - fsurfn(i,j) &
03308 - kh(ij)*Tbot(i,j)
03309 Tsf(m) = rhs(ij) / diag(ij)
03310
03311 enddo
03312
03313
03314
03315
03316
03317
03318
03319
03320
03321
03322
03323
03324
03325
03326
03327
03328
03329 all_converged = .true.
03330
03331
03332
03333
03334 do ij = 1, isolve
03335 m = indxij(ij)
03336
03337
03338
03339
03340
03341
03342
03343
03344 converged(m) = .true.
03345 dTsf(ij) = Tsf(m) - Tsf_start(ij)
03346 avg_Tsf = c0
03347
03348
03349
03350
03351
03352
03353 if (Tsf(m) > puny) then
03354 Tsf(m) = c0
03355 dTsf(ij) = -Tsf_start(ij)
03356
03357
03358
03359
03360
03361
03362
03363 elseif (niter > 1 &
03364 .and. Tsf_start(ij) <= -puny &
03365 .and. abs(dTsf(ij)) > puny &
03366 .and. abs(dTsf_prev(m)) > puny &
03367 .and. -dTsf(ij)/(dTsf_prev(m)+puny*puny) > p5) then
03368
03369 avg_Tsf = c1
03370 dTsf(ij) = p5 * dTsf(ij)
03371 converged(m) = .false.
03372 all_converged = .false.
03373 endif
03374
03375
03376
03377
03378
03379 Tsf(m) = Tsf(m) &
03380 + avg_Tsf * p5 * (Tsf_start(ij) - Tsf(m))
03381
03382 enddo
03383
03384
03385
03386
03387 do ij = 1, isolve
03388 i = indxii(ij)
03389 j = indxjj(ij)
03390 m = indxij(ij)
03391
03392
03393
03394
03395
03396 if (abs(dTsf(ij)) > Tsf_errmax) then
03397 converged(m) = .false.
03398 all_converged = .false.
03399 endif
03400
03401
03402
03403
03404
03405 fsurfn(i,j) = fsurfn(i,j) + dTsf(ij)*dfsurf_dT(ij)
03406 fcondtopn(i,j) = kh(ij) * (Tsf(m)-Tbot(i,j))
03407
03408 if (Tsf(m) > -puny .and. fsurfn(i,j) < fcondtopn(i,j)) then
03409 converged(m) = .false.
03410 all_converged = .false.
03411 endif
03412
03413 fcondbot(m) = fcondtopn(i,j)
03414
03415 dTsf_prev(m) = dTsf(ij)
03416
03417 enddo
03418
03419 deallocate(diag)
03420 deallocate(rhs)
03421 deallocate(kh)
03422 deallocate(heff)
03423 deallocate(Tsf_start)
03424 deallocate(dTsf)
03425 deallocate(dfsurf_dT)
03426
03427 enddo
03428
03429 if (.not.all_converged) then
03430
03431 do ij = 1, icells
03432 i = indxi(ij)
03433 j = indxj(ij)
03434
03435
03436
03437
03438 if (.not.converged(ij)) then
03439 write(nu_diag,*) 'Thermo iteration does not converge,', &
03440 'istep1, my_task, i, j:', &
03441 istep1, my_task, i, j
03442 write(nu_diag,*) 'Ice thickness:', hilyr(ij)*nilyr
03443 write(nu_diag,*) 'Snow thickness:', hslyr(ij)*nslyr
03444 write(nu_diag,*) 'dTsf, Tsf_errmax:',dTsf_prev(ij), &
03445 Tsf_errmax
03446 write(nu_diag,*) 'Tsf:', Tsf(ij)
03447 write(nu_diag,*) 'fsurfn:', fsurfn(i,j)
03448 write(nu_diag,*) 'fcondtopn, fcondbot', &
03449 fcondtopn(i,j), fcondbot(ij)
03450 l_stop = .true.
03451 istop = i
03452 jstop = j
03453 return
03454 endif
03455 enddo
03456 endif
03457
03458
03459
03460
03461
03462 if (l_zerolayerchecks) then
03463 do ij = 1, icells
03464 i = indxi(ij)
03465 j = indxj(ij)
03466
03467 if (Tsf(ij) < c0 .and. &
03468 abs(fcondtopn(i,j)-fsurfn(i,j)) > puny) then
03469
03470 write(nu_diag,*) 'fcondtopn does not equal fsurfn,', &
03471 'istep1, my_task, i, j:', &
03472 istep1, my_task, i, j
03473 write(nu_diag,*) 'Tsf=',Tsf(ij)
03474 write(nu_diag,*) 'fcondtopn=',fcondtopn(i,j)
03475 write(nu_diag,*) 'fsurfn=',fsurfn(i,j)
03476 l_stop = .true.
03477 istop = i
03478 jstop = j
03479 return
03480 endif
03481 enddo
03482 endif
03483
03484
03485
03486
03487
03488 do ij = 1, icells
03489 i = indxi(ij)
03490 j = indxj(ij)
03491
03492
03493 flwoutn(i,j) = flwoutn(i,j) + dTsf_prev(ij) * dflwout_dT(ij)
03494 fsensn(i,j) = fsensn(i,j) + dTsf_prev(ij) * dfsens_dT(ij)
03495 flatn(i,j) = flatn(i,j) + dTsf_prev(ij) * dflat_dT(ij)
03496
03497
03498 fswabsn(i,j) = fswsfc(i,j) + fswthrun(i,j)
03499
03500 enddo
03501
03502 end subroutine zerolayer_temperature
03503
03504
03505
03506
03507
03508
03509
03510
03511
03512
03513
03514
03515
03516
03517
03518
03519
03520
03521 subroutine thickness_changes (nx_block, ny_block, &
03522 dt, &
03523 yday, icells, &
03524 indxi, indxj, &
03525 efinal, &
03526 hin, hilyr, &
03527 hsn, hslyr, &
03528 qin, qsn, &
03529 fbot, Tbot, &
03530 flatn, fsurfn, &
03531 fcondtopn, fcondbot, &
03532 fsnow, hsn_new, &
03533 fhocnn, evapn, &
03534 meltt, melts, &
03535 meltb, iage, &
03536 congel, snoice, &
03537 mlt_onset, frz_onset)
03538
03539
03540
03541
03542
03543 integer (kind=int_kind), intent(in) ::
03544 nx_block, ny_block,
03545 icells
03546
03547 integer (kind=int_kind), dimension(nx_block*ny_block),
03548 intent(in) ::
03549 indxi, indxj
03550
03551 real (kind=dbl_kind), intent(in) ::
03552 dt ,
03553 yday
03554
03555 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
03556 fbot ,
03557 Tbot ,
03558 fsnow ,
03559 flatn ,
03560 fsurfn ,
03561 fcondtopn
03562
03563 real (kind=dbl_kind), dimension (icells), intent(in) ::
03564 fcondbot
03565
03566 real (kind=dbl_kind), dimension (icells,nilyr),
03567 intent(inout) ::
03568 qin
03569
03570 real (kind=dbl_kind), dimension (icells,nslyr),
03571 intent(inout) ::
03572 qsn
03573
03574 real (kind=dbl_kind), dimension (icells),
03575 intent(inout) ::
03576 hilyr ,
03577 hslyr
03578
03579 real (kind=dbl_kind), dimension (nx_block,ny_block),
03580 intent(inout) ::
03581 meltt ,
03582 melts ,
03583 meltb ,
03584 congel ,
03585 snoice ,
03586 iage ,
03587 mlt_onset ,
03588 frz_onset
03589
03590 real (kind=dbl_kind), dimension (icells),
03591 intent(inout) ::
03592 hin ,
03593 hsn
03594
03595 real (kind=dbl_kind), dimension (icells), intent(out)::
03596 efinal
03597
03598 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out)::
03599 fhocnn ,
03600 evapn
03601
03602 real (kind=dbl_kind), dimension (icells), intent(out)::
03603 hsn_new
03604
03605
03606
03607 real (kind=dbl_kind), parameter ::
03608 qbotmax = -p5*rhoi*Lfresh
03609
03610 integer (kind=int_kind) ::
03611 i, j ,
03612 ij ,
03613 k
03614
03615 real (kind=dbl_kind), dimension (icells) ::
03616 esub ,
03617 econ ,
03618 etop_mlt ,
03619 ebot_mlt ,
03620 ebot_gro
03621
03622 real (kind=dbl_kind) ::
03623 dhi ,
03624 dhs ,
03625 Ti ,
03626 Ts ,
03627 qbot ,
03628 qsub ,
03629 hqtot ,
03630 wk1 ,
03631 qsnew ,
03632 hstot
03633
03634 real (kind=dbl_kind), dimension (icells,nilyr+1) ::
03635 zi1 ,
03636 zi2
03637
03638 real (kind=dbl_kind), dimension (icells,nslyr+1) ::
03639 zs1 ,
03640 zs2
03641
03642 real (kind=dbl_kind), dimension (icells,nilyr) ::
03643 dzi
03644
03645 real (kind=dbl_kind), dimension (icells,nslyr) ::
03646 dzs
03647
03648
03649
03650
03651
03652 hsn_new (:) = c0
03653
03654 do k = 1, nilyr
03655 do ij = 1, icells
03656 dzi(ij,k) = hilyr(ij)
03657 enddo
03658 enddo
03659
03660 do k = 1, nslyr
03661 do ij = 1, icells
03662 dzs(ij,k) = hslyr(ij)
03663 enddo
03664 enddo
03665
03666
03667
03668
03669
03670
03671
03672 if (.not. l_brine) then
03673
03674 do k = 1, nslyr
03675
03676
03677
03678 do ij = 1, icells
03679
03680 Ts = (Lfresh + qsn(ij,k)/rhos) / cp_ice
03681 if (Ts > c0) then
03682 dhs = cp_ice*Ts*dzs(ij,k) / Lfresh
03683 dzs(ij,k) = dzs(ij,k) - dhs
03684 qsn(ij,k) = -rhos*Lfresh
03685 endif
03686 enddo
03687 enddo
03688
03689 do k = 1, nilyr
03690
03691
03692
03693 do ij = 1, icells
03694
03695 Ti = (Lfresh + qin(ij,k)/rhoi) / cp_ice
03696 if (Ti > c0) then
03697 dhi = cp_ice*Ti*dzi(ij,k) / Lfresh
03698 dzi(ij,k) = dzi(ij,k) - dhi
03699 qin(ij,k) = -rhoi*Lfresh
03700 endif
03701 enddo
03702 enddo
03703
03704 endif
03705
03706
03707
03708
03709
03710
03711 do ij = 1, icells
03712 i = indxi(ij)
03713 j = indxj(ij)
03714
03715 wk1 = -flatn(i,j) * dt
03716 esub(ij) = max(wk1, c0)
03717 econ(ij) = min(wk1, c0)
03718
03719 wk1 = (fsurfn(i,j) - fcondtopn(i,j)) * dt
03720 etop_mlt(ij) = max(wk1, c0)
03721
03722 wk1 = (fcondbot(ij) - fbot(i,j)) * dt
03723 ebot_mlt(ij) = max(wk1, c0)
03724 ebot_gro(ij) = min(wk1, c0)
03725
03726
03727
03728
03729
03730
03731 evapn (i,j) = c0
03732
03733 if (hsn(ij) > puny) then
03734 dhs = econ(ij) / (qsn(ij,1) - rhos*Lvap)
03735 dzs(ij,1) = dzs(ij,1) + dhs
03736 evapn(i,j) = evapn(i,j) + dhs*rhos
03737 else
03738 dhi = econ(ij) / (qin(ij,1) - rhoi*Lvap)
03739 dzi(ij,1) = dzi(ij,1) + dhi
03740 evapn(i,j) = evapn(i,j) + dhi*rhoi
03741 endif
03742
03743
03744
03745
03746
03747
03748 if (heat_capacity) then
03749 qbot = -rhoi * (cp_ice * (Tmlt(nilyr+1)-Tbot(i,j)) &
03750 + Lfresh * (c1-Tmlt(nilyr+1)/Tbot(i,j)) &
03751 - cp_ocn * Tmlt(nilyr+1))
03752 qbot = min (qbot, qbotmax)
03753 else
03754 qbot = -rhoi * Lfresh
03755 endif
03756
03757 dhi = ebot_gro(ij) / qbot
03758
03759 hqtot = dzi(ij,nilyr)*qin(ij,nilyr) + dhi*qbot
03760 dzi(ij,nilyr) = dzi(ij,nilyr) + dhi
03761
03762 if (dzi(ij,nilyr) > puny) &
03763 qin(ij,nilyr) = hqtot / dzi(ij,nilyr)
03764
03765
03766
03767
03768
03769
03770 congel(i,j) = congel(i,j) + dhi
03771 if (dhi > puny .and. frz_onset(i,j) < puny) &
03772 frz_onset(i,j) = yday
03773
03774 enddo
03775
03776 do k = 1, nslyr
03777
03778
03779
03780 do ij = 1, icells
03781 i = indxi(ij)
03782 j = indxj(ij)
03783
03784
03785
03786
03787
03788 qsub = qsn(ij,k) - rhos*Lvap
03789 dhs = max (-dzs(ij,k), esub(ij)/qsub)
03790 dzs(ij,k) = dzs(ij,k) + dhs
03791 esub(ij) = esub(ij) - dhs*qsub
03792 esub(ij) = max(esub(ij), c0)
03793 evapn(i,j) = evapn(i,j) + dhs*rhos
03794
03795
03796
03797
03798
03799 dhs = max(-dzs(ij,k), etop_mlt(ij)/qsn(ij,k))
03800 dzs(ij,k) = dzs(ij,k) + dhs
03801 etop_mlt(ij) = etop_mlt(ij) - dhs*qsn(ij,k)
03802 etop_mlt(ij) = max(etop_mlt(ij), c0)
03803
03804
03805 if (dhs < -puny .and. mlt_onset(i,j) < puny) &
03806 mlt_onset(i,j) = yday
03807 melts(i,j) = melts(i,j) - dhs
03808
03809 enddo
03810 enddo
03811
03812 do k = 1, nilyr
03813
03814
03815
03816 do ij = 1, icells
03817 i = indxi(ij)
03818 j = indxj(ij)
03819
03820
03821
03822
03823
03824 qsub = qin(ij,k) - rhoi*Lvap
03825 dhi = max (-dzi(ij,k), esub(ij)/qsub)
03826 dzi(ij,k) = dzi(ij,k) + dhi
03827 esub(ij) = esub(ij) - dhi*qsub
03828 esub(ij) = max(esub(ij), c0)
03829 evapn(i,j) = evapn(i,j) + dhi*rhoi
03830
03831
03832
03833
03834
03835 dhi = max(-dzi(ij,k), etop_mlt(ij)/qin(ij,k))
03836 dzi(ij,k) = dzi(ij,k) + dhi
03837 etop_mlt(ij) = etop_mlt(ij) - dhi*qin(ij,k)
03838 etop_mlt(ij) = max(etop_mlt(ij), c0)
03839
03840
03841 if (dhi < -puny .and. mlt_onset(i,j) < puny) &
03842 mlt_onset(i,j) = yday
03843 meltt(i,j) = meltt(i,j) - dhi
03844
03845 enddo
03846 enddo
03847
03848 do k = nilyr, 1, -1
03849
03850
03851
03852 do ij = 1, icells
03853 i = indxi(ij)
03854 j = indxj(ij)
03855
03856
03857
03858
03859
03860 dhi = max(-dzi(ij,k), ebot_mlt(ij)/qin(ij,k))
03861 dzi(ij,k) = dzi(ij,k) + dhi
03862 ebot_mlt(ij) = ebot_mlt(ij) - dhi*qin(ij,k)
03863 ebot_mlt(ij) = max(ebot_mlt(ij), c0)
03864
03865
03866 meltb(i,j) = meltb(i,j) - dhi
03867
03868 enddo
03869 enddo
03870
03871 do k = nslyr, 1, -1
03872
03873
03874
03875 do ij = 1, icells
03876
03877
03878
03879
03880
03881 dhs = max(-dzs(ij,k), ebot_mlt(ij)/qsn(ij,k))
03882 dzs(ij,k) = dzs(ij,k) + dhs
03883 ebot_mlt(ij) = ebot_mlt(ij) - dhs*qsn(ij,k)
03884 ebot_mlt(ij) = max(ebot_mlt(ij), c0)
03885
03886 enddo
03887 enddo
03888
03889
03890
03891
03892
03893
03894 do ij = 1, icells
03895 i = indxi(ij)
03896 j = indxj(ij)
03897 fhocnn(i,j) = fbot(i,j) &
03898 + (esub(ij) + etop_mlt(ij) + ebot_mlt(ij))/dt
03899 enddo
03900
03901
03902
03903
03904
03905
03906
03907
03908 do ij = 1, icells
03909 i = indxi(ij)
03910 j = indxj(ij)
03911
03912
03913
03914
03915
03916
03917
03918 if (fsnow(i,j) > c0) then
03919
03920 hsn_new(ij) = fsnow(i,j)/rhos * dt
03921 qsnew = -rhos*Lfresh
03922 hstot = dzs(ij,1) + hsn_new(ij)
03923
03924 if (hstot > c0) then
03925 qsn(ij,1) = (dzs(ij,1) * qsn(ij,1) &
03926 + hsn_new(ij) * qsnew) / hstot
03927
03928 qsn(ij,1) = min(qsn(ij,1), -rhos*Lfresh)
03929
03930 dzs(ij,1) = hstot
03931 endif
03932 endif
03933
03934
03935
03936
03937
03938 hin(ij) = c0
03939 hsn(ij) = c0
03940 enddo
03941
03942 do k = 1, nilyr
03943
03944
03945
03946 do ij = 1, icells
03947 hin(ij) = hin(ij) + dzi(ij,k)
03948 enddo
03949 enddo
03950
03951 do k = 1, nslyr
03952
03953
03954
03955 do ij = 1, icells
03956 hsn(ij) = hsn(ij) + dzs(ij,k)
03957 enddo
03958 enddo
03959
03960
03961
03962
03963
03964 call freeboard (nx_block, ny_block, &
03965 icells, &
03966 indxi, indxj, &
03967 dt, &
03968 snoice, &
03969 iage, &
03970 hin, hsn, &
03971 qin, qsn, &
03972 dzi, dzs)
03973
03974
03975
03976
03977
03978
03979
03980
03981
03982
03983 do ij = 1, icells
03984
03985 if (hin(ij) > c0) then
03986 hilyr(ij) = hin(ij) / real(nilyr,kind=dbl_kind)
03987 else
03988 hin(ij) = c0
03989 hilyr(ij) = c0
03990 endif
03991 if (hsn(ij) > c0) then
03992 hslyr(ij) = hsn(ij) / real(nslyr,kind=dbl_kind)
03993 else
03994 hsn(ij) = c0
03995 hslyr(ij) = c0
03996 endif
03997
03998
03999
04000
04001
04002
04003 zi1(ij,1) = c0
04004 zi1(ij,1+nilyr) = hin(ij)
04005
04006 zi2(ij,1) = c0
04007 zi2(ij,1+nilyr) = hin(ij)
04008
04009 enddo
04010
04011 if (heat_capacity) then
04012
04013 do k = 1, nilyr-1
04014
04015
04016
04017 do ij = 1, icells
04018 zi1(ij,k+1) = zi1(ij,k) + dzi(ij,k)
04019 zi2(ij,k+1) = zi2(ij,k) + hilyr(ij)
04020 end do
04021 enddo
04022
04023
04024
04025
04026
04027 call adjust_enthalpy (nx_block, ny_block, &
04028 nilyr, icells, &
04029 indxi, indxj, &
04030 zi1, zi2, &
04031 hilyr, hin, &
04032 qin)
04033
04034 else
04035
04036 do ij = 1, icells
04037 qin(ij,1) = -rhoi * Lfresh
04038 qsn(ij,1) = -rhos * Lfresh
04039 end do
04040
04041 endif
04042
04043 if (nslyr > 1) then
04044
04045
04046
04047
04048
04049
04050 do ij = 1, icells
04051 zs1(ij,1) = c0
04052 zs1(ij,1+nslyr) = hsn(ij)
04053
04054 zs2(ij,1) = c0
04055 zs2(ij,1+nslyr) = hsn(ij)
04056 enddo
04057
04058 do k = 1, nslyr-1
04059
04060
04061
04062 do ij = 1, icells
04063 zs1(ij,k+1) = zs1(ij,k) + dzs(ij,k)
04064 zs2(ij,k+1) = zs2(ij,k) + hslyr(ij)
04065 end do
04066 enddo
04067
04068
04069
04070
04071
04072 call adjust_enthalpy (nx_block, ny_block, &
04073 nslyr, icells, &
04074 indxi, indxj, &
04075 zs1, zs2, &
04076 hslyr, hsn, &
04077 qsn)
04078
04079 endif
04080
04081
04082
04083
04084
04085
04086 do ij = 1, icells
04087 i = indxi(ij)
04088 j = indxj(ij)
04089 efinal(ij) = -evapn(i,j)*Lvap
04090 evapn(i,j) = evapn(i,j)/dt
04091 enddo
04092
04093 do k = 1, nslyr
04094
04095
04096
04097 do ij = 1, icells
04098 efinal(ij) = efinal(ij) + hslyr(ij)*qsn(ij,k)
04099 enddo
04100 enddo
04101
04102 do k = 1, nilyr
04103
04104
04105
04106 do ij = 1, icells
04107 efinal(ij) = efinal(ij) + hilyr(ij)*qin(ij,k)
04108 enddo
04109 enddo
04110
04111 end subroutine thickness_changes
04112
04113
04114
04115
04116
04117
04118
04119
04120
04121
04122
04123
04124
04125
04126
04127
04128
04129
04130
04131 subroutine freeboard (nx_block, ny_block, &
04132 icells, &
04133 indxi, indxj, &
04134 dt, &
04135 snoice, &
04136 iage, &
04137 hin, hsn, &
04138 qin, qsn, &
04139 dzi, dzs)
04140
04141
04142
04143
04144
04145 integer (kind=int_kind), intent(in) ::
04146 nx_block, ny_block,
04147 icells
04148
04149 integer (kind=int_kind), dimension(nx_block*ny_block),
04150 intent(in) ::
04151 indxi, indxj
04152
04153 real (kind=dbl_kind), intent(in) ::
04154 dt
04155
04156 real (kind=dbl_kind), dimension (nx_block,ny_block),
04157 intent(inout) ::
04158 snoice ,
04159 iage
04160
04161 real (kind=dbl_kind), dimension (icells),
04162 intent(inout) ::
04163 hin ,
04164 hsn
04165
04166 real (kind=dbl_kind), dimension (icells,nilyr),
04167 intent(inout) ::
04168 qin
04169
04170 real (kind=dbl_kind), dimension (icells,nilyr),
04171 intent(inout) ::
04172 dzi
04173
04174 real (kind=dbl_kind), dimension (icells,nslyr),
04175 intent(in) ::
04176 qsn
04177
04178 real (kind=dbl_kind), dimension (icells,nslyr),
04179 intent(inout) ::
04180 dzs
04181
04182
04183
04184 integer (kind=int_kind) ::
04185 i, j ,
04186 ij ,
04187 k
04188
04189 real (kind=dbl_kind), dimension (icells) ::
04190 dhin ,
04191 dhsn ,
04192 hqs
04193
04194 real (kind=dbl_kind) ::
04195 wk1 ,
04196 dhs
04197
04198
04199
04200
04201
04202 do ij = 1, icells
04203
04204 dhin(ij) = c0
04205 dhsn(ij) = c0
04206 hqs (ij) = c0
04207
04208 wk1 = hsn(ij) - hin(ij)*(rhow-rhoi)/rhos
04209
04210 if (wk1 > puny .and. hsn(ij) > puny) then
04211 dhsn(ij) = min(wk1*rhoi/rhow, hsn(ij))
04212 dhin(ij) = dhsn(ij) * rhos/rhoi
04213 endif
04214 enddo
04215
04216
04217
04218
04219
04220
04221 do k = nslyr, 1, -1
04222
04223
04224
04225 do ij = 1, icells
04226 if (dhin(ij) > puny) then
04227 dhs = min(dhsn(ij), dzs(ij,k))
04228 hsn(ij) = hsn(ij) - dhs
04229 dzs(ij,k) = dzs(ij,k) - dhs
04230 dhsn(ij) = dhsn(ij) - dhs
04231 dhsn(ij) = max(dhsn(ij),c0)
04232 hqs(ij) = hqs(ij) + dhs * qsn(ij,k)
04233 endif
04234 enddo
04235 enddo
04236
04237
04238
04239
04240
04241
04242
04243
04244 do ij = 1, icells
04245 i = indxi(ij)
04246 j = indxj(ij)
04247
04248 if (dhin(ij) > puny) then
04249
04250
04251
04252
04253 wk1 = dzi(ij,1) + dhin(ij)
04254 hin(ij) = hin(ij) + dhin(ij)
04255 qin(ij,1) = (dzi(ij,1)*qin(ij,1) + hqs(ij)) / wk1
04256 dzi(ij,1) = wk1
04257
04258
04259 snoice(i,j) = snoice(i,j) + dhin(ij)
04260 endif
04261
04262 enddo
04263
04264 end subroutine freeboard
04265
04266
04267
04268
04269
04270
04271
04272
04273
04274
04275
04276
04277
04278
04279
04280
04281
04282
04283 subroutine adjust_enthalpy (nx_block, ny_block, &
04284 nlyr, icells, &
04285 indxi, indxj, &
04286 z1, z2, &
04287 hlyr, hn, &
04288 qn)
04289
04290
04291
04292
04293
04294 integer (kind=int_kind), intent(in) ::
04295 nx_block, ny_block,
04296 nlyr ,
04297 icells
04298
04299 integer (kind=int_kind), dimension (nx_block*ny_block),
04300 intent(in) ::
04301 indxi, indxj
04302
04303 real (kind=dbl_kind), dimension (icells,nlyr+1),
04304 intent(in) ::
04305 z1 ,
04306 z2
04307
04308 real (kind=dbl_kind), dimension (icells), intent(in) ::
04309 hlyr
04310
04311 real (kind=dbl_kind), dimension (icells), intent(in) ::
04312 hn
04313
04314 real (kind=dbl_kind), dimension (icells,nlyr),
04315 intent(inout) ::
04316 qn
04317
04318
04319
04320 integer (kind=int_kind) ::
04321 i, j ,
04322 ij ,
04323 k, k1, k2
04324
04325 real (kind=dbl_kind) ::
04326 hovlp
04327
04328 real (kind=dbl_kind), dimension (icells) ::
04329 rhlyr
04330
04331 real (kind=dbl_kind), dimension (icells,nlyr) ::
04332 hq
04333
04334
04335
04336
04337
04338 do ij = 1, icells
04339 rhlyr(ij) = c0
04340 if (hn(ij) > puny) rhlyr(ij) = c1 / hlyr(ij)
04341 enddo
04342
04343
04344
04345
04346
04347 do k2 = 1, nlyr
04348
04349 do ij = 1, icells
04350 hq(ij,k2) = c0
04351 enddo
04352
04353 do k1 = 1, nlyr
04354
04355
04356
04357 do ij = 1, icells
04358 hovlp = min (z1(ij,k1+1), z2(ij,k2+1)) &
04359 - max (z1(ij,k1), z2(ij,k2))
04360 hovlp = max (hovlp, c0)
04361
04362 hq(ij,k2) = hq(ij,k2) + hovlp*qn(ij,k1)
04363 enddo
04364 enddo
04365 enddo
04366
04367
04368
04369
04370
04371 do k = 1, nlyr
04372 do ij = 1, icells
04373 i = indxi(ij)
04374 j = indxj(ij)
04375 qn(ij,k) = hq(ij,k) * rhlyr(ij)
04376 enddo
04377 enddo
04378
04379 end subroutine adjust_enthalpy
04380
04381
04382
04383
04384
04385
04386
04387
04388
04389
04390
04391
04392
04393
04394
04395
04396
04397
04398 subroutine conservation_check_vthermo(nx_block, ny_block, &
04399 my_task, istep1, &
04400 dt, icells, &
04401 indxi, indxj, &
04402 fsurfn, flatn, &
04403 fhocnn, fswint, &
04404 fsnow, &
04405 einit, efinal, &
04406 l_stop, &
04407 istop, jstop)
04408
04409
04410
04411
04412
04413 integer (kind=int_kind), intent(in) ::
04414 nx_block, ny_block,
04415 my_task ,
04416 istep1 ,
04417 icells
04418
04419 integer (kind=int_kind), dimension(nx_block*ny_block),
04420 intent(in) ::
04421 indxi, indxj
04422
04423 real (kind=dbl_kind), intent(in) ::
04424 dt
04425
04426 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
04427 fsurfn ,
04428 flatn ,
04429 fhocnn ,
04430 fswint ,
04431 fsnow
04432
04433
04434 real (kind=dbl_kind), dimension (icells), intent(in) ::
04435 einit ,
04436 efinal
04437
04438 logical (kind=log_kind), intent(inout) ::
04439 l_stop
04440
04441 integer (kind=int_kind), intent(inout) ::
04442 istop, jstop
04443
04444
04445
04446 integer (kind=int_kind) ::
04447 i, j ,
04448 ij
04449
04450 real (kind=dbl_kind) ::
04451 einp ,
04452 ferr
04453
04454
04455
04456
04457
04458
04459
04460 do ij = 1, icells
04461 i = indxi(ij)
04462 j = indxj(ij)
04463
04464 einp = (fsurfn(i,j) - flatn(i,j) + fswint(i,j) - fhocnn(i,j) &
04465 - fsnow(i,j)*Lfresh) * dt
04466 ferr = abs(efinal(ij)-einit(ij)-einp) / dt
04467 if (ferr > ferrmax) then
04468 l_stop = .true.
04469 istop = i
04470 jstop = j
04471
04472
04473
04474
04475
04476
04477
04478
04479
04480
04481
04482
04483 write(nu_diag,*) 'Thermo energy conservation error'
04484 write(nu_diag,*) 'istep1, my_task, i, j:', &
04485 istep1, my_task, i, j
04486 write(nu_diag,*) 'Flux error (W/m^2) =', ferr
04487 write(nu_diag,*) 'Energy error (J) =', ferr*dt
04488 write(nu_diag,*) 'Initial energy =', einit(ij)
04489 write(nu_diag,*) 'Final energy =', efinal(ij)
04490 write(nu_diag,*) 'efinal - einit =', &
04491 efinal(ij)-einit(ij)
04492 write(nu_diag,*) 'Input energy =', einp
04493 return
04494 endif
04495 enddo
04496
04497 end subroutine conservation_check_vthermo
04498
04499
04500
04501
04502
04503
04504
04505
04506
04507
04508
04509
04510
04511
04512
04513
04514
04515
04516
04517
04518 subroutine update_state_vthermo (nx_block, ny_block, &
04519 icells, &
04520 indxi, indxj, &
04521 Tf, Tsf, &
04522 hin, hsn, &
04523 qin, qsn, &
04524 aicen, vicen, &
04525 vsnon, Tsfcn, &
04526 eicen, esnon)
04527
04528
04529
04530
04531
04532 integer (kind=int_kind), intent(in) ::
04533 nx_block, ny_block,
04534 icells
04535
04536 integer (kind=int_kind), dimension(nx_block*ny_block),
04537 intent(in) ::
04538 indxi, indxj
04539
04540 real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) ::
04541 Tf
04542
04543 real (kind=dbl_kind), dimension(icells), intent(in) ::
04544 Tsf
04545
04546 real (kind=dbl_kind), dimension(icells), intent(in) ::
04547 hin ,
04548 hsn
04549
04550 real (kind=dbl_kind), dimension (icells,nilyr),
04551 intent(in) ::
04552 qin
04553
04554 real (kind=dbl_kind), dimension (icells,nslyr),
04555 intent(in) ::
04556 qsn
04557
04558 real (kind=dbl_kind), dimension (nx_block,ny_block),
04559 intent(inout) ::
04560 aicen ,
04561 vicen ,
04562 vsnon ,
04563 Tsfcn
04564
04565 real (kind=dbl_kind), dimension (nx_block,ny_block,nilyr),
04566 intent(inout) ::
04567 eicen
04568
04569 real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr),
04570 intent(inout) ::
04571 esnon
04572
04573
04574
04575 integer (kind=int_kind) ::
04576 i, j ,
04577 ij ,
04578 k
04579
04580
04581
04582
04583 do ij = 1, icells
04584 i = indxi(ij)
04585 j = indxj(ij)
04586
04587 if (hin(ij) > c0) then
04588
04589 vicen(i,j) = aicen(i,j) * hin(ij)
04590 vsnon(i,j) = aicen(i,j) * hsn(ij)
04591 Tsfcn(i,j) = Tsf(ij)
04592 else
04593 aicen(i,j) = c0
04594 vicen(i,j) = c0
04595 vsnon(i,j) = c0
04596 Tsfcn(i,j) = Tf(i,j)
04597 endif
04598
04599 enddo
04600
04601 do k = 1, nilyr
04602 do ij = 1, icells
04603 i = indxi(ij)
04604 j = indxj(ij)
04605
04606 if (hin(ij) > c0) then
04607 eicen(i,j,k) = qin(ij,k) * vicen(i,j) &
04608 /real(nilyr,kind=dbl_kind)
04609 else
04610 eicen(i,j,k) = c0
04611 endif
04612
04613 enddo
04614 enddo
04615
04616 do k = 1, nslyr
04617 do ij = 1, icells
04618 i = indxi(ij)
04619 j = indxj(ij)
04620
04621 if (hin(ij) > c0) then
04622 esnon(i,j,k) = qsn(ij,k) * vsnon(i,j) &
04623 /real(nslyr,kind=dbl_kind)
04624 else
04625 esnon(i,j,k) = c0
04626 endif
04627
04628 enddo
04629 enddo
04630
04631 end subroutine update_state_vthermo
04632
04633
04634
04635
04636
04637 end module ice_therm_vertical
04638
04639