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
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 module ice_shortwave
00051
00052
00053
00054 use ice_kinds_mod
00055 use ice_domain_size
00056 use ice_constants
00057 use ice_blocks
00058
00059
00060
00061 implicit none
00062 save
00063
00064 character (len=char_len) ::
00065 shortwave,
00066 albedo_type
00067
00068
00069
00070 real (kind=dbl_kind) ::
00071 albicev ,
00072 albicei ,
00073 albsnowv,
00074 albsnowi
00075
00076
00077 real (kind=dbl_kind),
00078 dimension (nx_block,ny_block,ncat,max_blocks) ::
00079 alvdrn ,
00080 alidrn ,
00081 alvdfn ,
00082 alidfn
00083
00084
00085 real (kind=dbl_kind),
00086 dimension (nx_block,ny_block,ncat,max_blocks) ::
00087 albicen ,
00088 albsnon ,
00089 albpndn
00090
00091
00092 real (kind=dbl_kind),
00093 dimension (nx_block,ny_block,ntilyr,max_blocks) ::
00094 Iswabsn
00095
00096 real (kind=dbl_kind),
00097 dimension (nx_block,ny_block,ntslyr,max_blocks) ::
00098 Sswabsn
00099
00100 real (kind=dbl_kind), dimension (nx_block,ny_block,ncat,max_blocks) ::
00101 fswsfcn ,
00102 fswthrun ,
00103 fswintn
00104
00105 real (kind=dbl_kind) ::
00106 rnilyr ,
00107 rnslyr
00108
00109
00110 real (kind=dbl_kind) ::
00111 R_ice ,
00112 R_pnd ,
00113 R_snw
00114
00115 real (kind=dbl_kind) ::
00116 dT_mlt_in ,
00117 rsnw_melt_in
00118
00119
00120 real (kind=dbl_kind) ::
00121 exp_min
00122
00123
00124
00125 contains
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140 subroutine init_shortwave
00141
00142
00143
00144 use ice_domain_size
00145 use ice_blocks
00146 use ice_calendar
00147 use ice_domain
00148 use ice_flux
00149 use ice_grid
00150 use ice_itd
00151 use ice_orbital
00152 use ice_state
00153
00154
00155
00156
00157
00158
00159
00160 integer (kind=int_kind) ::
00161 icells
00162
00163 integer (kind=int_kind), dimension(nx_block*ny_block) ::
00164 indxi, indxj
00165
00166 integer (kind=int_kind) ::
00167 i, j, ij ,
00168 iblk ,
00169 ilo,ihi,jlo,jhi,
00170 n ,
00171 il1, il2 ,
00172 sl1, sl2
00173
00174 real (kind=dbl_kind) :: cszn
00175
00176 type (block) ::
00177 this_block
00178
00179 do iblk=1,nblocks
00180 do j = 1, ny_block
00181 do i = 1, nx_block
00182 alvdf(i,j,iblk) = c0
00183 alidf(i,j,iblk) = c0
00184 alvdr(i,j,iblk) = c0
00185 alidr(i,j,iblk) = c0
00186 alvdr_gbm(i,j,iblk) = c0
00187 alidr_gbm(i,j,iblk) = c0
00188 alvdf_gbm(i,j,iblk) = c0
00189 alidf_gbm(i,j,iblk) = c0
00190 enddo
00191 enddo
00192 enddo
00193
00194 if (trim(shortwave) == 'dEdd') then
00195
00196 call init_orbit
00197 call init_dEdd
00198
00199 else
00200
00201 coszen(:,:,:) = p5
00202
00203 do iblk=1,nblocks
00204 this_block = get_block(blocks_ice(iblk),iblk)
00205 ilo = this_block%ilo
00206 ihi = this_block%ihi
00207 jlo = this_block%jlo
00208 jhi = this_block%jhi
00209
00210 do n = 1, ncat
00211
00212 icells = 0
00213 do j = jlo, jhi
00214 do i = ilo, ihi
00215 if (aicen(i,j,n,iblk) > puny) then
00216 icells = icells + 1
00217 indxi(icells) = i
00218 indxj(icells) = j
00219 endif
00220 enddo
00221 enddo
00222
00223 il1 = ilyr1(n)
00224 il2 = ilyrn(n)
00225 sl1 = slyr1(n)
00226 sl2 = slyrn(n)
00227
00228 Sswabsn(:,:,sl1:sl2,iblk) = c0
00229
00230 call shortwave_ccsm3(nx_block, ny_block, &
00231 icells, &
00232 indxi, indxj, &
00233 aicen(:,:,n,iblk), vicen(:,:,n,iblk), &
00234 vsnon(:,:,n,iblk), &
00235 trcrn(:,:,nt_Tsfc,n,iblk), &
00236 swvdr(:,:, iblk), swvdf(:,:, iblk), &
00237 swidr(:,:, iblk), swidf(:,:, iblk), &
00238 alvdrn(:,:,n,iblk),alidrn(:,:,n,iblk), &
00239 alvdfn(:,:,n,iblk),alidfn(:,:,n,iblk), &
00240 fswsfcn(:,:,n,iblk),fswintn(:,:,n,iblk),&
00241 fswthrun(:,:,n,iblk), &
00242 Iswabsn(:,:,il1:il2,iblk), &
00243 albicen(:,:,n,iblk),albsnon(:,:,n,iblk))
00244
00245
00246 enddo
00247 enddo
00248
00249 endif
00250
00251
00252
00253
00254
00255 do iblk=1,nblocks
00256 this_block = get_block(blocks_ice(iblk),iblk)
00257 ilo = this_block%ilo
00258 ihi = this_block%ihi
00259 jlo = this_block%jlo
00260 jhi = this_block%jhi
00261
00262 do n = 1, ncat
00263
00264 icells = 0
00265 do j = jlo, jhi
00266 do i = ilo, ihi
00267 if (aicen(i,j,n,iblk) > puny) then
00268 icells = icells + 1
00269 indxi(icells) = i
00270 indxj(icells) = j
00271 endif
00272 enddo
00273 enddo
00274
00275 do ij = 1, icells
00276 i = indxi(ij)
00277 j = indxj(ij)
00278
00279 alvdf(i,j,iblk) = alvdf(i,j,iblk) &
00280 + alvdfn(i,j,n,iblk)*aicen(i,j,n,iblk)
00281 alidf(i,j,iblk) = alidf(i,j,iblk) &
00282 + alidfn(i,j,n,iblk)*aicen(i,j,n,iblk)
00283 alvdr(i,j,iblk) = alvdr(i,j,iblk) &
00284 + alvdrn(i,j,n,iblk)*aicen(i,j,n,iblk)
00285 alidr(i,j,iblk) = alidr(i,j,iblk) &
00286 + alidrn(i,j,n,iblk)*aicen(i,j,n,iblk)
00287
00288 if (coszen(i,j,iblk) > puny) then
00289 albice(i,j,iblk) = albice(i,j,iblk) &
00290 + albicen(i,j,n,iblk)*aicen(i,j,n,iblk)
00291 albsno(i,j,iblk) = albsno(i,j,iblk) &
00292 + albsnon(i,j,n,iblk)*aicen(i,j,n,iblk)
00293 albpnd(i,j,iblk) = albpnd(i,j,iblk) &
00294 + albpndn(i,j,n,iblk)*aicen(i,j,n,iblk)
00295 endif
00296 enddo
00297
00298 enddo
00299
00300
00301
00302
00303
00304 do j = 1, ny_block
00305 do i = 1, nx_block
00306 alvdf_gbm (i,j,iblk) = alvdf (i,j,iblk)
00307 alidf_gbm (i,j,iblk) = alidf (i,j,iblk)
00308 alvdr_gbm (i,j,iblk) = alvdr (i,j,iblk)
00309 alidr_gbm (i,j,iblk) = alidr (i,j,iblk)
00310
00311
00312 cszn = c0
00313 if (coszen(i,j,iblk) > puny) cszn = c1
00314 do n = 1, nstreams
00315 albcnt(i,j,iblk,n) = albcnt(i,j,iblk,n) + cszn
00316 enddo
00317 enddo
00318 enddo
00319
00320 enddo
00321
00322 end subroutine init_shortwave
00323
00324
00325
00326
00327
00328
00329
00330
00331 subroutine shortwave_ccsm3 (nx_block, ny_block, &
00332 icells, &
00333 indxi, indxj, &
00334 aicen, vicen, &
00335 vsnon, Tsfcn, &
00336 swvdr, swvdf, &
00337 swidr, swidf, &
00338 alvdrn, alidrn, &
00339 alvdfn, alidfn, &
00340 fswsfc, fswint, &
00341 fswthru, Iswabs, &
00342 albin, albsn)
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356 integer (kind=int_kind), intent(in) ::
00357 nx_block, ny_block,
00358 icells
00359
00360 integer (kind=int_kind), dimension (nx_block*ny_block),
00361 intent(in) ::
00362 indxi ,
00363 indxj
00364
00365 real (kind=dbl_kind), dimension (nx_block,ny_block),
00366 intent(in) ::
00367 aicen ,
00368 vicen ,
00369 vsnon ,
00370 Tsfcn ,
00371 swvdr ,
00372 swvdf ,
00373 swidr ,
00374 swidf
00375
00376 real (kind=dbl_kind), dimension (nx_block,ny_block),
00377 intent(out) ::
00378 alvdrn ,
00379 alidrn ,
00380 alvdfn ,
00381 alidfn ,
00382 fswsfc ,
00383 fswint ,
00384 fswthru ,
00385 albin ,
00386 albsn
00387
00388 real (kind=dbl_kind), dimension (nx_block,ny_block,nilyr),
00389 intent(out) ::
00390 Iswabs
00391
00392
00393
00394
00395
00396 real (kind=dbl_kind), dimension (nx_block,ny_block)::
00397 alvdrni,
00398 alidrni,
00399 alvdfni,
00400 alidfni,
00401 alvdrns,
00402 alidrns,
00403 alvdfns,
00404 alidfns
00405
00406
00407
00408
00409
00410 if (trim(albedo_type) == 'constant') then
00411 call constant_albedos (nx_block, ny_block, &
00412 icells, &
00413 indxi, indxj, &
00414 aicen, &
00415 vsnon, Tsfcn, &
00416 alvdrni, alidrni, &
00417 alvdfni, alidfni, &
00418 alvdrns, alidrns, &
00419 alvdfns, alidfns, &
00420 alvdrn, alidrn, &
00421 alvdfn, alidfn, &
00422 albin, albsn)
00423 else
00424 call compute_albedos (nx_block, ny_block, &
00425 icells, &
00426 indxi, indxj, &
00427 aicen, vicen, &
00428 vsnon, Tsfcn, &
00429 alvdrni, alidrni, &
00430 alvdfni, alidfni, &
00431 alvdrns, alidrns, &
00432 alvdfns, alidfns, &
00433 alvdrn, alidrn, &
00434 alvdfn, alidfn, &
00435 albin, albsn)
00436 endif
00437
00438
00439
00440
00441
00442 call absorbed_solar (nx_block, ny_block, &
00443 icells, &
00444 indxi, indxj, &
00445 aicen, &
00446 vicen, vsnon, &
00447 swvdr, swvdf, &
00448 swidr, swidf, &
00449 alvdrni, alvdfni, &
00450 alidrni, alidfni, &
00451 alvdrns, alvdfns, &
00452 alidrns, alidfns, &
00453 fswsfc, fswint, &
00454 fswthru, Iswabs)
00455
00456 end subroutine shortwave_ccsm3
00457
00458
00459
00460
00461
00462
00463
00464
00465 subroutine compute_albedos (nx_block, ny_block, &
00466 icells, &
00467 indxi, indxj, &
00468 aicen, vicen, &
00469 vsnon, Tsfcn, &
00470 alvdrni, alidrni, &
00471 alvdfni, alidfni, &
00472 alvdrns, alidrns, &
00473 alvdfns, alidfns, &
00474 alvdrn, alidrn, &
00475 alvdfn, alidfn, &
00476 albin, albsn)
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490 integer (kind=int_kind), intent(in) ::
00491 nx_block, ny_block,
00492 icells
00493
00494 integer (kind=int_kind), dimension (nx_block*ny_block),
00495 intent(in) ::
00496 indxi ,
00497 indxj
00498
00499 real (kind=dbl_kind), dimension (nx_block,ny_block),
00500 intent(in) ::
00501 aicen ,
00502 vicen ,
00503 vsnon ,
00504 Tsfcn
00505
00506 real (kind=dbl_kind), dimension (nx_block,ny_block),
00507 intent(out) ::
00508 alvdrni ,
00509 alidrni ,
00510 alvdfni ,
00511 alidfni ,
00512 alvdrns ,
00513 alidrns ,
00514 alvdfns ,
00515 alidfns ,
00516 alvdrn ,
00517 alidrn ,
00518 alvdfn ,
00519 alidfn ,
00520 albin ,
00521 albsn
00522
00523
00524
00525 real (kind=dbl_kind), parameter ::
00526 ahmax = p5 ,
00527
00528 dT_mlt = c1 , &
00529
00530 dalb_mlt = -0.075_dbl_kind, &
00531
00532 dalb_mltv = -p1 , &
00533
00534 dalb_mlti = -p15
00535
00536
00537 integer (kind=int_kind) ::
00538 i, j, n
00539
00540 real (kind=dbl_kind) ::
00541 hi ,
00542 hs ,
00543 albo,
00544 fh ,
00545 fT ,
00546 dTs ,
00547 fhtan,
00548 asnow
00549
00550 integer (kind=int_kind) ::
00551 ij
00552
00553 fhtan = atan(ahmax*c4)
00554
00555 do j = 1, ny_block
00556 do i = 1, nx_block
00557 alvdrni(i,j) = albocn
00558 alidrni(i,j) = albocn
00559 alvdfni(i,j) = albocn
00560 alidfni(i,j) = albocn
00561
00562 alvdrns(i,j) = albocn
00563 alidrns(i,j) = albocn
00564 alvdfns(i,j) = albocn
00565 alidfns(i,j) = albocn
00566
00567 alvdrn(i,j) = albocn
00568 alidrn(i,j) = albocn
00569 alvdfn(i,j) = albocn
00570 alidfn(i,j) = albocn
00571
00572 albin(i,j) = c0
00573 albsn(i,j) = c0
00574 enddo
00575 enddo
00576
00577
00578
00579
00580
00581
00582
00583
00584 do ij = 1, icells
00585 i = indxi(ij)
00586 j = indxj(ij)
00587
00588 hi = vicen(i,j) / aicen(i,j)
00589 hs = vsnon(i,j) / aicen(i,j)
00590
00591
00592 fh = min(atan(hi*c4)/fhtan,c1)
00593 albo = albocn*(c1-fh)
00594 alvdfni(i,j) = albicev*fh + albo
00595 alidfni(i,j) = albicei*fh + albo
00596
00597
00598 dTs = Timelt - Tsfcn(i,j)
00599 fT = min(dTs/dT_mlt-c1,c0)
00600 alvdfni(i,j) = alvdfni(i,j) - dalb_mlt*fT
00601 alidfni(i,j) = alidfni(i,j) - dalb_mlt*fT
00602
00603
00604 alvdfni(i,j) = max (alvdfni(i,j), albocn)
00605 alidfni(i,j) = max (alidfni(i,j), albocn)
00606
00607 if (hs > puny) then
00608
00609 alvdfns(i,j) = albsnowv
00610 alidfns(i,j) = albsnowi
00611
00612
00613 alvdfns(i,j) = alvdfns(i,j) - dalb_mltv*fT
00614 alidfns(i,j) = alidfns(i,j) - dalb_mlti*fT
00615
00616 endif
00617
00618
00619 alvdrni(i,j) = alvdfni(i,j)
00620 alidrni(i,j) = alidfni(i,j)
00621 alvdrns(i,j) = alvdfns(i,j)
00622 alidrns(i,j) = alidfns(i,j)
00623
00624
00625
00626 if (hs > puny) then
00627 asnow = hs / (hs + snowpatch)
00628 else
00629 asnow = c0
00630 endif
00631
00632
00633 alvdfn(i,j) = alvdfni(i,j)*(c1-asnow) + &
00634 alvdfns(i,j)*asnow
00635 alidfn(i,j) = alidfni(i,j)*(c1-asnow) + &
00636 alidfns(i,j)*asnow
00637 alvdrn(i,j) = alvdrni(i,j)*(c1-asnow) + &
00638 alvdrns(i,j)*asnow
00639 alidrn(i,j) = alidrni(i,j)*(c1-asnow) + &
00640 alidrns(i,j)*asnow
00641
00642
00643 albin(i,j) = awtvdr*alvdrni(i,j) + awtidr*alidrni(i,j) &
00644 + awtvdf*alvdfni(i,j) + awtidf*alidfni(i,j)
00645 albsn(i,j) = awtvdr*alvdrns(i,j) + awtidr*alidrns(i,j) &
00646 + awtvdf*alvdfns(i,j) + awtidf*alidfns(i,j)
00647
00648 enddo
00649
00650 end subroutine compute_albedos
00651
00652
00653
00654
00655
00656
00657
00658
00659 subroutine constant_albedos (nx_block, ny_block, &
00660 icells, &
00661 indxi, indxj, &
00662 aicen, &
00663 vsnon, Tsfcn, &
00664 alvdrni, alidrni, &
00665 alvdfni, alidfni, &
00666 alvdrns, alidrns, &
00667 alvdfns, alidfns, &
00668 alvdrn, alidrn, &
00669 alvdfn, alidfn, &
00670 albin, albsn)
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684 integer (kind=int_kind), intent(in) ::
00685 nx_block, ny_block,
00686 icells
00687
00688 integer (kind=int_kind), dimension (nx_block*ny_block),
00689 intent(in) ::
00690 indxi ,
00691 indxj
00692
00693 real (kind=dbl_kind), dimension (nx_block,ny_block),
00694 intent(in) ::
00695 aicen ,
00696 vsnon ,
00697 Tsfcn
00698
00699 real (kind=dbl_kind), dimension (nx_block,ny_block),
00700 intent(out) ::
00701 alvdrni ,
00702 alidrni ,
00703 alvdfni ,
00704 alidfni ,
00705 alvdrns ,
00706 alidrns ,
00707 alvdfns ,
00708 alidfns ,
00709 alvdrn ,
00710 alidrn ,
00711 alvdfn ,
00712 alidfn ,
00713 albin ,
00714 albsn
00715
00716
00717
00718 real (kind=dbl_kind), parameter ::
00719 warmice = 0.68_dbl_kind,
00720 coldice = 0.70_dbl_kind,
00721 warmsnow = 0.77_dbl_kind,
00722 coldsnow = 0.81_dbl_kind
00723
00724 integer (kind=int_kind) ::
00725 i, j, n
00726
00727 real (kind=dbl_kind) ::
00728 hs
00729
00730 integer (kind=int_kind) ::
00731 ij
00732
00733 do j = 1, ny_block
00734 do i = 1, nx_block
00735 alvdrn(i,j) = albocn
00736 alidrn(i,j) = albocn
00737 alvdfn(i,j) = albocn
00738 alidfn(i,j) = albocn
00739
00740 albin(i,j) = c0
00741 albsn(i,j) = c0
00742 enddo
00743 enddo
00744
00745
00746
00747
00748
00749
00750
00751
00752 do ij = 1, icells
00753 i = indxi(ij)
00754 j = indxj(ij)
00755
00756 hs = vsnon(i,j) / aicen(i,j)
00757
00758 if (hs > puny) then
00759
00760 if (Tsfcn(i,j) >= -c2*puny) then
00761 alvdfn(i,j) = warmsnow
00762 alidfn(i,j) = warmsnow
00763 else
00764 alvdfn(i,j) = coldsnow
00765 alidfn(i,j) = coldsnow
00766 endif
00767 else
00768
00769 if (Tsfcn(i,j) >= -c2*puny) then
00770 alvdfn(i,j) = warmice
00771 alidfn(i,j) = warmice
00772 else
00773 alvdfn(i,j) = coldice
00774 alidfn(i,j) = coldice
00775 endif
00776 endif
00777
00778
00779 alvdrn (i,j) = alvdfn(i,j)
00780 alidrn (i,j) = alidfn(i,j)
00781
00782 alvdrni(i,j) = alvdrn(i,j)
00783 alidrni(i,j) = alidrn(i,j)
00784 alvdrns(i,j) = alvdrn(i,j)
00785 alidrns(i,j) = alidrn(i,j)
00786 alvdfni(i,j) = alvdfn(i,j)
00787 alidfni(i,j) = alidfn(i,j)
00788 alvdfns(i,j) = alvdfn(i,j)
00789 alidfns(i,j) = alidfn(i,j)
00790
00791
00792 albin(i,j) = awtvdr*alvdrni(i,j) + awtidr*alidrni(i,j) &
00793 + awtvdf*alvdfni(i,j) + awtidf*alidfni(i,j)
00794 albsn(i,j) = awtvdr*alvdrns(i,j) + awtidr*alidrns(i,j) &
00795 + awtvdf*alvdfns(i,j) + awtidf*alidfns(i,j)
00796
00797 enddo
00798
00799 end subroutine constant_albedos
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817 subroutine absorbed_solar (nx_block, ny_block, &
00818 icells, &
00819 indxi, indxj, &
00820 aicen, &
00821 vicen, vsnon, &
00822 swvdr, swvdf, &
00823 swidr, swidf, &
00824 alvdrni, alvdfni, &
00825 alidrni, alidfni, &
00826 alvdrns, alvdfns, &
00827 alidrns, alidfns, &
00828 fswsfc, fswint, &
00829 fswthru, Iswabs)
00830
00831
00832
00833 use ice_therm_vertical, only: heat_capacity
00834
00835
00836
00837 integer (kind=int_kind), intent(in) ::
00838 nx_block, ny_block,
00839 icells
00840
00841 integer (kind=int_kind), dimension(nx_block*ny_block),
00842 intent(in) ::
00843 indxi, indxj
00844
00845 real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) ::
00846 aicen ,
00847 vicen ,
00848 vsnon ,
00849 swvdr ,
00850 swvdf ,
00851 swidr ,
00852 swidf ,
00853 alvdrni ,
00854 alidrni ,
00855 alvdfni ,
00856 alidfni ,
00857 alvdrns ,
00858 alidrns ,
00859 alvdfns ,
00860 alidfns
00861
00862 real (kind=dbl_kind), dimension (nx_block,ny_block),
00863 intent(out)::
00864 fswsfc ,
00865 fswint ,
00866 fswthru
00867
00868 real (kind=dbl_kind), dimension (nx_block,ny_block,nilyr),
00869 intent(out) ::
00870 Iswabs
00871
00872
00873
00874 real (kind=dbl_kind), parameter ::
00875 i0vis = 0.70_dbl_kind
00876
00877 integer (kind=int_kind) ::
00878 i, j ,
00879 ij ,
00880 k
00881
00882 real (kind=dbl_kind), dimension (nx_block,ny_block) ::
00883 fswpen ,
00884 trantop ,
00885 tranbot
00886
00887 real (kind=dbl_kind) ::
00888 swabs ,
00889 swabsv ,
00890 swabsi ,
00891 fswpenvdr ,
00892 fswpenvdf ,
00893 hi ,
00894 hs ,
00895 hilyr ,
00896 asnow
00897
00898
00899
00900
00901
00902 rnilyr = real(nilyr,kind=dbl_kind)
00903
00904 do j = 1, ny_block
00905 do i = 1, nx_block
00906 fswsfc (i,j) = c0
00907 fswint (i,j) = c0
00908 fswthru(i,j) = c0
00909 fswpen (i,j) = c0
00910 trantop(i,j) = c0
00911 tranbot(i,j) = c0
00912 enddo
00913 enddo
00914 Iswabs (:,:,:) = c0
00915
00916 do ij = 1, icells
00917 i = indxi(ij)
00918 j = indxj(ij)
00919
00920 hs = vsnon(i,j) / aicen(i,j)
00921
00922
00923
00924
00925 if (hs > puny) then
00926 asnow = hs / (hs + snowpatch)
00927 else
00928 asnow = c0
00929 endif
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939 swabsv = swvdr(i,j) * ( (c1-alvdrni(i,j))*(c1-asnow) &
00940 + (c1-alvdrns(i,j))*asnow ) &
00941 + swvdf(i,j) * ( (c1-alvdfni(i,j))*(c1-asnow) &
00942 + (c1-alvdfns(i,j))*asnow )
00943
00944 swabsi = swidr(i,j) * ( (c1-alidrni(i,j))*(c1-asnow) &
00945 + (c1-alidrns(i,j))*asnow ) &
00946 + swidf(i,j) * ( (c1-alidfni(i,j))*(c1-asnow) &
00947 + (c1-alidfns(i,j))*asnow )
00948
00949 swabs = swabsv + swabsi
00950
00951 fswpenvdr = swvdr(i,j) * (c1-alvdrni(i,j)) * (c1-asnow) * i0vis
00952 fswpenvdf = swvdf(i,j) * (c1-alvdfni(i,j)) * (c1-asnow) * i0vis
00953
00954
00955
00956
00957
00958 fswpen(i,j) = fswpenvdr + fswpenvdf
00959
00960 fswsfc(i,j) = swabs - fswpen(i,j)
00961
00962 trantop(i,j) = c1
00963
00964 enddo
00965
00966
00967
00968
00969
00970 do k = 1, nilyr
00971
00972
00973
00974 do ij = 1, icells
00975 i = indxi(ij)
00976 j = indxj(ij)
00977
00978 hi = vicen(i,j) / aicen(i,j)
00979 hilyr = hi / rnilyr
00980
00981 tranbot(i,j) = exp (-kappav * hilyr * real(k,kind=dbl_kind))
00982 Iswabs(i,j,k) = fswpen(i,j) * (trantop(i,j)-tranbot(i,j))
00983
00984
00985 trantop(i,j) = tranbot(i,j)
00986
00987 enddo
00988 enddo
00989
00990
00991
00992
00993 do ij = 1, icells
00994 i = indxi(ij)
00995 j = indxj(ij)
00996
00997
00998 fswthru(i,j) = fswpen(i,j) * tranbot(i,j)
00999
01000
01001 fswint(i,j) = fswpen(i,j) - fswthru(i,j)
01002
01003 enddo
01004
01005
01006
01007
01008
01009
01010 if (.not. heat_capacity) then
01011
01012
01013
01014
01015 do ij = 1, icells
01016 i = indxi(ij)
01017 j = indxj(ij)
01018
01019
01020 fswsfc(i,j) = fswsfc(i,j) + fswint(i,j)
01021
01022
01023 fswint(i,j) = c0
01024 Iswabs(i,j,1) = c0
01025
01026 enddo
01027
01028 endif
01029
01030 end subroutine absorbed_solar
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042 subroutine init_dEdd
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055 use ice_domain_size
01056 use ice_blocks
01057 use ice_calendar
01058 use ice_domain
01059 use ice_flux
01060 use ice_grid
01061 use ice_itd
01062 use ice_meltpond
01063 use ice_orbital
01064 use ice_state
01065
01066
01067
01068
01069
01070
01071
01072 integer (kind=int_kind) ::
01073 icells
01074
01075 integer (kind=int_kind), dimension(nx_block*ny_block) ::
01076 indxi, indxj
01077
01078
01079
01080 real (kind=dbl_kind), dimension (nx_block,ny_block) ::
01081 fsn
01082 real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr) ::
01083 rhosnwn ,
01084 rsnwn
01085
01086
01087 real (kind=dbl_kind), dimension (nx_block,ny_block) ::
01088 fpn ,
01089 hpn
01090
01091 integer (kind=int_kind) ::
01092 i, j, ij ,
01093 iblk ,
01094 ilo,ihi,jlo,jhi,
01095 n ,
01096 il1, il2 ,
01097 sl1, sl2
01098
01099 type (block) ::
01100 this_block
01101
01102 exp_min = exp(-c10)
01103
01104 do iblk=1,nblocks
01105 this_block = get_block(blocks_ice(iblk),iblk)
01106 ilo = this_block%ilo
01107 ihi = this_block%ihi
01108 jlo = this_block%jlo
01109 jhi = this_block%jhi
01110
01111
01112 icells = 0
01113 do j = 1, ny_block
01114 do i = 1, nx_block
01115 if (tmask(i,j,iblk)) then
01116 icells = icells + 1
01117 indxi(icells) = i
01118 indxj(icells) = j
01119 endif
01120 enddo
01121 enddo
01122
01123 call compute_coszen (nx_block, ny_block, &
01124 icells, &
01125 indxi, indxj, &
01126 tlat (:,:,iblk), tlon(:,:,iblk), &
01127 coszen(:,:,iblk), dt)
01128
01129 do n = 1, ncat
01130
01131 icells = 0
01132 do j = jlo, jhi
01133 do i = ilo, ihi
01134 if (aicen(i,j,n,iblk) > puny) then
01135 icells = icells + 1
01136 indxi(icells) = i
01137 indxj(icells) = j
01138 endif
01139 enddo
01140 enddo
01141
01142 il1 = ilyr1(n)
01143 il2 = ilyrn(n)
01144 sl1 = slyr1(n)
01145 sl2 = slyrn(n)
01146
01147
01148
01149
01150
01151 call shortwave_dEdd_set_snow(nx_block, ny_block, &
01152 icells, &
01153 indxi, indxj, &
01154 aicen(:,:,n,iblk), vsnon(:,:,n,iblk), &
01155 trcrn(:,:,1,n,iblk), fsn, &
01156 rhosnwn, rsnwn)
01157
01158
01159 if (.not. tr_pond) then
01160
01161
01162 call shortwave_dEdd_set_pond(nx_block, ny_block, &
01163 icells, &
01164 indxi, indxj, &
01165 aicen(:,:,n,iblk), trcrn(:,:,1,n,iblk), &
01166 fsn, fpn, &
01167 hpn)
01168
01169 else
01170
01171 fpn(:,:) = apondn(:,:,n,iblk)
01172 hpn(:,:) = hpondn(:,:,n,iblk)
01173
01174 endif
01175
01176 call shortwave_dEdd(nx_block, ny_block, &
01177 icells, &
01178 indxi, indxj, &
01179 coszen(:,:, iblk), &
01180 aicen(:,:,n,iblk), vicen(:,:,n,iblk), &
01181 vsnon(:,:,n,iblk), fsn, &
01182 rhosnwn, rsnwn, &
01183 fpn, hpn, &
01184 trcrn(:,:,:,n,iblk),tarea(:,:,iblk), &
01185 swvdr(:,:, iblk), swvdf(:,:, iblk), &
01186 swidr(:,:, iblk), swidf(:,:, iblk), &
01187 alvdrn(:,:,n,iblk),alvdfn(:,:,n,iblk), &
01188 alidrn(:,:,n,iblk),alidfn(:,:,n,iblk), &
01189 fswsfcn(:,:,n,iblk),fswintn(:,:,n,iblk),&
01190 fswthrun(:,:,n,iblk), &
01191 Sswabsn(:,:,sl1:sl2,iblk), &
01192 Iswabsn(:,:,il1:il2,iblk), &
01193 albicen(:,:,n,iblk),albsnon(:,:,n,iblk),&
01194 albpndn(:,:,n,iblk))
01195
01196 enddo
01197 enddo
01198
01199 end subroutine init_dEdd
01200
01201
01202
01203
01204
01205
01206
01207
01208 subroutine shortwave_dEdd (nx_block, ny_block, &
01209 icells, indxi, &
01210 indxj, coszen, &
01211 aice, vice, &
01212 vsno, fs, &
01213 rhosnw, rsnw, &
01214 fp, hp, &
01215 trcr, tarea, &
01216 swvdr, swvdf, &
01217 swidr, swidf, &
01218 alvdr, alvdf, &
01219 alidr, alidf, &
01220 fswsfc, fswint, &
01221 fswthru, Sswabs, &
01222 Iswabs, albice, &
01223 albsno, albpnd)
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258 use ice_calendar
01259 use ice_state, only: nt_aero, tr_aero
01260
01261 use ice_diagnostics
01262
01263
01264
01265 integer (kind=int_kind),
01266 intent(in) ::
01267 nx_block, ny_block,
01268 icells
01269
01270 integer (kind=int_kind), dimension (nx_block*ny_block),
01271 intent(in) ::
01272 indxi ,
01273 indxj
01274
01275 real (kind=dbl_kind), dimension (nx_block,ny_block),
01276 intent(in) ::
01277 coszen ,
01278 aice ,
01279 vice ,
01280 vsno ,
01281 fs
01282
01283 real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr),
01284 intent(in) ::
01285 rhosnw ,
01286 rsnw
01287
01288 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr),
01289 intent(in) ::
01290 trcr
01291
01292 real (kind=dbl_kind), dimension (nx_block,ny_block),
01293 intent(in) ::
01294 tarea
01295
01296 real (kind=dbl_kind), dimension (nx_block,ny_block),
01297 intent(in) ::
01298 fp ,
01299 hp ,
01300 swvdr ,
01301 swvdf ,
01302 swidr ,
01303 swidf
01304
01305 real (kind=dbl_kind), dimension (nx_block,ny_block),
01306 intent(out) ::
01307 alvdr ,
01308 alvdf ,
01309 alidr ,
01310 alidf ,
01311 fswsfc ,
01312 fswint ,
01313 fswthru
01314
01315 real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr),
01316 intent(out) ::
01317 Sswabs
01318
01319 real (kind=dbl_kind), dimension (nx_block,ny_block,nilyr),
01320 intent(out) ::
01321 Iswabs
01322
01323 real (kind=dbl_kind), dimension (nx_block,ny_block),
01324 intent(out) ::
01325 albice ,
01326 albsno ,
01327 albpnd
01328
01329
01330
01331
01332
01333 real (kind=dbl_kind),dimension (nx_block,ny_block) ::
01334 fnidr
01335
01336 real (kind=dbl_kind), dimension(nx_block,ny_block) ::
01337 hs ,
01338 hi ,
01339 fi
01340
01341 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr) ::
01342 aero_mp
01343
01344 integer (kind=int_kind), dimension(nx_block,ny_block) ::
01345 srftyp
01346
01347 integer (kind=int_kind) ::
01348 i ,
01349 j ,
01350 ij ,
01351 k ,
01352 na ,
01353 icells_DE
01354
01355 integer (kind=int_kind), dimension (nx_block*ny_block) ::
01356 indxi_DE ,
01357 indxj_DE
01358
01359 real (kind=dbl_kind) ::
01360 hpmin ,
01361 hsmax ,
01362 hs_ssl ,
01363 frcadj
01364 data hpmin / .005_dbl_kind /
01365 data hs_ssl / .040_dbl_kind /
01366
01367
01368 integer (kind=int_kind) ::
01369 n
01370
01371 logical (kind=log_kind) ::
01372 dbug
01373
01374 real (kind=dbl_kind) ::
01375 swdn ,
01376 swab ,
01377 swalb
01378
01379
01380 real (kind=dbl_kind), dimension (nx_block,ny_block) ::
01381 avdrl ,
01382 avdfl ,
01383 aidrl ,
01384 aidfl
01385
01386
01387
01388 do j = 1, ny_block
01389 do i = 1, nx_block
01390
01391 hs(i,j) = c0
01392 hi(i,j) = c0
01393 fi(i,j) = c0
01394 srftyp(i,j) = 0
01395 alvdr(i,j) = c0
01396 alvdf(i,j) = c0
01397 alidr(i,j) = c0
01398 alidf(i,j) = c0
01399 avdrl(i,j) = c0
01400 avdfl(i,j) = c0
01401 aidrl(i,j) = c0
01402 aidfl(i,j) = c0
01403 fswsfc(i,j) = c0
01404 fswint(i,j) = c0
01405 fswthru(i,j) = c0
01406
01407 fnidr(i,j) = c0
01408 if( swidr(i,j) + swidf(i,j) > puny ) then
01409 fnidr(i,j) = swidr(i,j)/(swidr(i,j)+swidf(i,j))
01410 endif
01411 albice(i,j) = c0
01412 albsno(i,j) = c0
01413 albpnd(i,j) = c0
01414 enddo
01415 enddo
01416 Sswabs(:,:,:) = c0
01417 Iswabs(:,:,:) = c0
01418
01419
01420
01421 aero_mp(:,:,:) = c0
01422 if( tr_aero ) then
01423
01424
01425
01426
01427
01428 do na=1,4*n_aero,4
01429 do ij = 1, icells
01430 i = indxi(ij)
01431 j = indxj(ij)
01432
01433 if (aice(i,j) > puny .and. coszen(i,j) > puny) then
01434 aero_mp(i,j,na ) = trcr(i,j,nt_aero-1+na )*vsno(i,j)
01435 aero_mp(i,j,na+1) = trcr(i,j,nt_aero-1+na+1)*vsno(i,j)
01436 aero_mp(i,j,na+2) = trcr(i,j,nt_aero-1+na+2)*vice(i,j)
01437 aero_mp(i,j,na+3) = trcr(i,j,nt_aero-1+na+3)*vice(i,j)
01438 endif
01439 enddo
01440 enddo
01441 endif
01442
01443
01444
01445
01446
01447
01448
01449
01450 icells_DE = 0
01451 do ij = 1, icells
01452 i = indxi(ij)
01453 j = indxj(ij)
01454
01455 if (aice(i,j) > puny .and. coszen(i,j) > puny) then
01456
01457 hi(i,j) = vice(i,j) / aice(i,j)
01458 fi(i,j) = c1 - fs(i,j) - fp(i,j)
01459
01460 if(fi(i,j) > c0) then
01461 icells_DE = icells_DE + 1
01462 indxi_DE(icells_DE) = i
01463 indxj_DE(icells_DE) = j
01464
01465 srftyp(i,j) = 0
01466 endif
01467 endif
01468 enddo
01469
01470
01471 call compute_dEdd &
01472 (nx_block,ny_block, &
01473 icells_DE, indxi_DE, indxj_DE, fnidr, coszen, &
01474 swvdr, swvdf, swidr, swidf, srftyp, &
01475 hs, rhosnw, rsnw, hi, hp, &
01476 fi, aero_mp, avdrl, avdfl, &
01477 aidrl, aidfl, &
01478 fswsfc, fswint, &
01479 fswthru, Sswabs(:,:,:), &
01480 Iswabs)
01481
01482
01483
01484
01485 do ij = 1, icells_DE
01486 i = indxi_DE(ij)
01487 j = indxj_DE(ij)
01488 alvdr(i,j) = alvdr(i,j) + avdrl(i,j) *fi(i,j)
01489 alvdf(i,j) = alvdf(i,j) + avdfl(i,j) *fi(i,j)
01490 alidr(i,j) = alidr(i,j) + aidrl(i,j) *fi(i,j)
01491 alidf(i,j) = alidf(i,j) + aidfl(i,j) *fi(i,j)
01492
01493 albice(i,j) = albice(i,j) &
01494 + awtvdr*avdrl(i,j) + awtidr*aidrl(i,j) &
01495 + awtvdf*avdfl(i,j) + awtidf*aidfl(i,j)
01496 enddo
01497
01498
01499
01500
01501
01502 icells_DE = 0
01503 do ij = 1, icells
01504 i = indxi(ij)
01505 j = indxj(ij)
01506
01507 if (aice(i,j) > puny .and. coszen(i,j) > puny) then
01508
01509 hs(i,j) = vsno(i,j) / aice(i,j)
01510
01511 if(fs(i,j) > c0) then
01512 icells_DE = icells_DE + 1
01513 indxi_DE(icells_DE) = i
01514 indxj_DE(icells_DE) = j
01515
01516 srftyp(i,j) = 1
01517 endif
01518 endif
01519 enddo
01520
01521
01522 call compute_dEdd &
01523 (nx_block,ny_block, &
01524 icells_DE, indxi_DE, indxj_DE, fnidr, coszen, &
01525 swvdr, swvdf, swidr, swidf, srftyp, &
01526 hs, rhosnw, rsnw, hi, hp, &
01527 fs, aero_mp, avdrl, avdfl, &
01528 aidrl, aidfl, &
01529 fswsfc, fswint, &
01530 fswthru, Sswabs(:,:,:), &
01531 Iswabs)
01532
01533
01534
01535
01536 do ij = 1, icells_DE
01537 i = indxi_DE(ij)
01538 j = indxj_DE(ij)
01539 alvdr(i,j) = alvdr(i,j) + avdrl(i,j) *fs(i,j)
01540 alvdf(i,j) = alvdf(i,j) + avdfl(i,j) *fs(i,j)
01541 alidr(i,j) = alidr(i,j) + aidrl(i,j) *fs(i,j)
01542 alidf(i,j) = alidf(i,j) + aidfl(i,j) *fs(i,j)
01543
01544 albsno(i,j) = albsno(i,j) &
01545 + awtvdr*avdrl(i,j) + awtidr*aidrl(i,j) &
01546 + awtvdf*avdfl(i,j) + awtidf*aidfl(i,j)
01547 enddo
01548
01549
01550
01551
01552
01553 icells_DE = 0
01554 do ij = 1, icells
01555 i = indxi(ij)
01556 j = indxj(ij)
01557 hi(i,j) = c0
01558
01559 if (aice(i,j) > puny .and. coszen(i,j) > puny) then
01560 hi(i,j) = vice(i,j) / aice(i,j)
01561
01562 if( fp(i,j) > puny .and. hp(i,j) > hpmin ) then
01563 icells_DE = icells_DE + 1
01564 indxi_DE(icells_DE) = i
01565 indxj_DE(icells_DE) = j
01566
01567 srftyp(i,j) = 2
01568 endif
01569 endif
01570 enddo
01571
01572
01573 call compute_dEdd &
01574 (nx_block,ny_block, &
01575 icells_DE, indxi_DE, indxj_DE, fnidr, coszen, &
01576 swvdr, swvdf, swidr, swidf, srftyp, &
01577 hs, rhosnw, rsnw, hi, hp, &
01578 fp, aero_mp, avdrl, avdfl, &
01579 aidrl, aidfl, &
01580 fswsfc, fswint, &
01581 fswthru, Sswabs(:,:,:), &
01582 Iswabs)
01583
01584
01585
01586
01587 do ij = 1, icells_DE
01588 i = indxi_DE(ij)
01589 j = indxj_DE(ij)
01590 alvdr(i,j) = alvdr(i,j) + avdrl(i,j) *fp(i,j)
01591 alvdf(i,j) = alvdf(i,j) + avdfl(i,j) *fp(i,j)
01592 alidr(i,j) = alidr(i,j) + aidrl(i,j) *fp(i,j)
01593 alidf(i,j) = alidf(i,j) + aidfl(i,j) *fp(i,j)
01594
01595 albpnd(i,j) = albpnd(i,j) &
01596 + awtvdr*avdrl(i,j) + awtidr*aidrl(i,j) &
01597 + awtvdf*avdfl(i,j) + awtidf*aidfl(i,j)
01598 enddo
01599
01600 dbug = .false.
01601 if (dbug .and. print_points) then
01602 do n = 1, npnt
01603 if (my_task == pmloc(n)) then
01604 i = piloc(n)
01605 j = pjloc(n)
01606 if( coszen(i,j) > .01_dbl_kind ) then
01607 write(nu_diag,*) ' my_task = ',my_task &
01608 ,' printing point = ',n &
01609 ,' i and j = ',i,j
01610 write(nu_diag,*) ' coszen = ', &
01611 coszen(i,j)
01612 write(nu_diag,*) ' swvdr swvdf = ', &
01613 swvdr(i,j),swvdf(i,j)
01614 write(nu_diag,*) ' swidr swidf = ', &
01615 swidr(i,j),swidf(i,j)
01616 write(nu_diag,*) ' aice = ', &
01617 aice(i,j)
01618 write(nu_diag,*) ' hs = ', &
01619 hs(i,j)
01620 write(nu_diag,*) ' hp = ', &
01621 hp(i,j)
01622 write(nu_diag,*) ' fs = ', &
01623 fs(i,j)
01624 write(nu_diag,*) ' fi = ', &
01625 fi(i,j)
01626 write(nu_diag,*) ' fp = ', &
01627 fp(i,j)
01628 write(nu_diag,*) ' hi = ', &
01629 hi(i,j)
01630 write(nu_diag,*) ' srftyp = ', &
01631 srftyp(i,j)
01632 write(nu_diag,*) ' alvdr alvdf = ', &
01633 alvdr(i,j),alvdf(i,j)
01634 write(nu_diag,*) ' alidr alidf = ', &
01635 alidr(i,j),alidf(i,j)
01636 write(nu_diag,*) ' fswsfc fswint fswthru = ', &
01637 fswsfc(i,j),fswint(i,j),fswthru(i,j)
01638 swdn = swvdr(i,j)+swvdf(i,j)+swidr(i,j)+swidf(i,j)
01639 swab = fswsfc(i,j)+fswint(i,j)+fswthru(i,j)
01640 swalb = (1.-swab/(swdn+.0001))
01641 write(nu_diag,*) ' swdn swab swalb = ',swdn,swab,swalb
01642 do k = 1, nslyr
01643 write(nu_diag,*) ' snow layer k = ', k, &
01644 ' rhosnw = ', &
01645 rhosnw(i,j,k), &
01646 ' rsnw = ', &
01647 rsnw(i,j,k)
01648 enddo
01649 do k = 1, nslyr
01650 write(nu_diag,*) ' snow layer k = ', k, &
01651 ' Sswabs(k) = ', Sswabs(i,j,k)
01652 enddo
01653 do k = 1, nilyr
01654 write(nu_diag,*) ' sea ice layer k = ', k, &
01655 ' Iswabs(k) = ', Iswabs(i,j,k)
01656 enddo
01657 endif
01658 endif
01659 enddo
01660 endif
01661
01662 end subroutine shortwave_dEdd
01663
01664
01665
01666
01667
01668
01669
01670
01671 subroutine compute_dEdd &
01672 (nx_block,ny_block, &
01673 icells_DE, indxi_DE, indxj_DE, fnidr, coszen, &
01674 swvdr, swvdf, swidr, swidf, srftyp, &
01675 hs, rhosnw, rsnw, hi, hp, &
01676 fi, aero_mp, alvdr, alvdf, &
01677 alidr, alidf, &
01678 fswsfc, fswint, &
01679 fswthru, Sswabs, &
01680 Iswabs)
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695 use ice_therm_vertical, only: heat_capacity
01696
01697
01698
01699 integer (kind=int_kind),
01700 intent(in) ::
01701 nx_block, ny_block,
01702 icells_DE
01703
01704 integer (kind=int_kind), dimension(nx_block*ny_block),
01705 intent(in) ::
01706 indxi_DE,
01707 indxj_DE
01708
01709 real (kind=dbl_kind), dimension (nx_block,ny_block),
01710 intent(in) ::
01711 fnidr ,
01712 coszen ,
01713 swvdr ,
01714 swvdf ,
01715 swidr ,
01716 swidf
01717
01718 integer (kind=int_kind), dimension(nx_block,ny_block),
01719 intent(in) ::
01720 srftyp
01721
01722 real (kind=dbl_kind), dimension(nx_block,ny_block),
01723 intent(in) ::
01724 hs
01725
01726 real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr),
01727 intent(in) ::
01728 rhosnw ,
01729 rsnw
01730
01731 real (kind=dbl_kind), dimension(nx_block,ny_block),
01732 intent(in) ::
01733 hi ,
01734 hp ,
01735 fi
01736
01737 real (kind=dbl_kind), dimension (nx_block,ny_block,ntrcr),
01738 intent(in) ::
01739 aero_mp
01740
01741 real (kind=dbl_kind), dimension (nx_block,ny_block),
01742 intent(inout) ::
01743 alvdr ,
01744 alvdf ,
01745 alidr ,
01746 alidf ,
01747 fswsfc ,
01748 fswint ,
01749 fswthru
01750
01751 real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr),
01752 intent(inout) ::
01753 Sswabs
01754
01755 real (kind=dbl_kind), dimension (nx_block,ny_block,nilyr),
01756 intent(inout) ::
01757 Iswabs
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864 integer (kind=int_kind) ::
01865 i ,
01866 j ,
01867 k ,
01868 ij ,
01869 ns ,
01870 nr ,
01871 ksa ,
01872 ki ,
01873 km ,
01874 kp ,
01875 ksrf ,
01876 ksnow ,
01877 kii
01878
01879 integer (kind=int_kind), parameter ::
01880 klev = nslyr + nilyr + 1 ,
01881 klevp = klev + 1
01882
01883
01884 integer (kind=int_kind), parameter ::
01885 nspint = 3 ,
01886 nmbrad = 32
01887
01888 real (kind=dbl_kind), dimension(icells_DE) ::
01889 avdr ,
01890 avdf ,
01891 aidr ,
01892 aidf
01893
01894 real (kind=dbl_kind), dimension(icells_DE) ::
01895 fsfc ,
01896 fint ,
01897 fthru
01898
01899 real (kind=dbl_kind), dimension(icells_DE,nslyr) ::
01900 Sabs
01901
01902 real (kind=dbl_kind), dimension(icells_DE,nilyr) ::
01903 Iabs
01904
01905 real (kind=dbl_kind), dimension (icells_DE,nspint) ::
01906 wghtns
01907
01908 real (kind=dbl_kind), parameter ::
01909 cp67 = 0.67_dbl_kind ,
01910 cp33 = 0.33_dbl_kind ,
01911 cp78 = 0.78_dbl_kind ,
01912 cp22 = 0.22_dbl_kind ,
01913 cp01 = 0.01_dbl_kind
01914
01915 real (kind=dbl_kind), dimension (0:klev,icells_DE) ::
01916 tau ,
01917 w0 ,
01918 g
01919
01920
01921
01922 real (kind=dbl_kind), dimension (0:klevp,icells_DE) ::
01923 trndir ,
01924 trntdr ,
01925 trndif ,
01926 rupdir ,
01927 rupdif ,
01928 rdndif
01929
01930 real (kind=dbl_kind) ::
01931 refk
01932
01933 real (kind=dbl_kind), dimension (0:klevp,icells_DE) ::
01934 fdirup ,
01935 fdirdn ,
01936 fdifup ,
01937 fdifdn
01938
01939
01940 real (kind=dbl_kind), dimension (nspint) ::
01941 Qs ,
01942 ks ,
01943 ws ,
01944 gs
01945
01946 real (kind=dbl_kind), dimension (nmbrad) ::
01947 rsnw_tab
01948
01949 real (kind=dbl_kind), dimension (nspint,nmbrad) ::
01950 Qs_tab ,
01951 ws_tab ,
01952 gs_tab
01953 real (kind=dbl_kind) ::
01954 delr ,
01955 rhoi ,
01956 fr ,
01957 fr_max ,
01958 fr_min
01959
01960
01961
01962 real (kind=dbl_kind), dimension (nspint) ::
01963 ki_ssl_mn ,
01964 wi_ssl_mn ,
01965 gi_ssl_mn ,
01966 ki_dl_mn ,
01967 wi_dl_mn ,
01968 gi_dl_mn ,
01969 ki_int_mn ,
01970 wi_int_mn ,
01971 gi_int_mn ,
01972 ki_p_ssl_mn ,
01973 wi_p_ssl_mn ,
01974 gi_p_ssl_mn ,
01975 ki_p_int_mn ,
01976 wi_p_int_mn ,
01977 gi_p_int_mn
01978
01979
01980
01981 real (kind=dbl_kind), dimension (nspint) ::
01982 ki_ssl ,
01983 wi_ssl ,
01984 gi_ssl ,
01985 ki_dl ,
01986 wi_dl ,
01987 gi_dl ,
01988 ki_int ,
01989 wi_int ,
01990 gi_int ,
01991 ki_p_ssl ,
01992 wi_p_ssl ,
01993 gi_p_ssl ,
01994 ki_p_int ,
01995 wi_p_int ,
01996 gi_p_int
01997
01998 real (kind=dbl_kind) ::
01999 hi_ssl ,
02000 hs_ssl ,
02001 dz ,
02002 dz_ssl ,
02003 fs
02004
02005
02006
02007
02008 real (kind=dbl_kind) ::
02009 kalg ,
02010 sig ,
02011 kabs ,
02012 sigp
02013
02014
02015 real (kind=dbl_kind), dimension (nspint) ::
02016 kw ,
02017 ww ,
02018 gw
02019 real (kind=dbl_kind), dimension (icells_DE) ::
02020 albodr ,
02021 albodf
02022
02023
02024 real (kind=dbl_kind) ::
02025 fp_ice ,
02026 fm_ice ,
02027 fp_pnd ,
02028 fm_pnd
02029
02030
02031 real (kind=dbl_kind) ::
02032 hpmin ,
02033 hp0 ,
02034 sig_i ,
02035 sig_p ,
02036 kext
02037
02038
02039
02040
02041
02042
02043
02044
02045 integer (kind=int_kind), parameter ::
02046 nmbaer = 6
02047
02048 integer (kind=int_kind) ::
02049 nmbaer_actual,
02050 na
02051
02052 real (kind=dbl_kind) ::
02053 kaer_tab(nspint,nmbaer),
02054 waer_tab(nspint,nmbaer),
02055 gaer_tab(nspint,nmbaer),
02056 taer ,
02057 waer ,
02058 gaer
02059
02060
02061 data rsnw_tab/ &
02062 5._dbl_kind, 7._dbl_kind, 10._dbl_kind, 15._dbl_kind, &
02063 20._dbl_kind, 30._dbl_kind, 40._dbl_kind, 50._dbl_kind, &
02064 65._dbl_kind, 80._dbl_kind, 100._dbl_kind, 120._dbl_kind, &
02065 140._dbl_kind, 170._dbl_kind, 200._dbl_kind, 240._dbl_kind, &
02066 290._dbl_kind, 350._dbl_kind, 420._dbl_kind, 500._dbl_kind, &
02067 570._dbl_kind, 660._dbl_kind, 760._dbl_kind, 870._dbl_kind, &
02068 1000._dbl_kind, 1100._dbl_kind, 1250._dbl_kind, 1400._dbl_kind, &
02069 1600._dbl_kind, 1800._dbl_kind, 2000._dbl_kind, 2500._dbl_kind/
02070
02071
02072 data Qs_tab/ &
02073 2.131798_dbl_kind, 2.187756_dbl_kind, 2.267358_dbl_kind, &
02074 2.104499_dbl_kind, 2.148345_dbl_kind, 2.236078_dbl_kind, &
02075 2.081580_dbl_kind, 2.116885_dbl_kind, 2.175067_dbl_kind, &
02076 2.062595_dbl_kind, 2.088937_dbl_kind, 2.130242_dbl_kind, &
02077 2.051403_dbl_kind, 2.072422_dbl_kind, 2.106610_dbl_kind, &
02078 2.039223_dbl_kind, 2.055389_dbl_kind, 2.080586_dbl_kind, &
02079 2.032383_dbl_kind, 2.045751_dbl_kind, 2.066394_dbl_kind, &
02080 2.027920_dbl_kind, 2.039388_dbl_kind, 2.057224_dbl_kind, &
02081 2.023444_dbl_kind, 2.033137_dbl_kind, 2.048055_dbl_kind, &
02082 2.020412_dbl_kind, 2.028840_dbl_kind, 2.041874_dbl_kind, &
02083 2.017608_dbl_kind, 2.024863_dbl_kind, 2.036046_dbl_kind, &
02084 2.015592_dbl_kind, 2.022021_dbl_kind, 2.031954_dbl_kind, &
02085 2.014083_dbl_kind, 2.019887_dbl_kind, 2.028853_dbl_kind, &
02086 2.012368_dbl_kind, 2.017471_dbl_kind, 2.025353_dbl_kind, &
02087 2.011092_dbl_kind, 2.015675_dbl_kind, 2.022759_dbl_kind, &
02088 2.009837_dbl_kind, 2.013897_dbl_kind, 2.020168_dbl_kind, &
02089 2.008668_dbl_kind, 2.012252_dbl_kind, 2.017781_dbl_kind, &
02090 2.007627_dbl_kind, 2.010813_dbl_kind, 2.015678_dbl_kind, &
02091 2.006764_dbl_kind, 2.009577_dbl_kind, 2.013880_dbl_kind, &
02092 2.006037_dbl_kind, 2.008520_dbl_kind, 2.012382_dbl_kind, &
02093 2.005528_dbl_kind, 2.007807_dbl_kind, 2.011307_dbl_kind, &
02094 2.005025_dbl_kind, 2.007079_dbl_kind, 2.010280_dbl_kind, &
02095 2.004562_dbl_kind, 2.006440_dbl_kind, 2.009333_dbl_kind, &
02096 2.004155_dbl_kind, 2.005898_dbl_kind, 2.008523_dbl_kind, &
02097 2.003794_dbl_kind, 2.005379_dbl_kind, 2.007795_dbl_kind, &
02098 2.003555_dbl_kind, 2.005041_dbl_kind, 2.007329_dbl_kind, &
02099 2.003264_dbl_kind, 2.004624_dbl_kind, 2.006729_dbl_kind, &
02100 2.003037_dbl_kind, 2.004291_dbl_kind, 2.006230_dbl_kind, &
02101 2.002776_dbl_kind, 2.003929_dbl_kind, 2.005700_dbl_kind, &
02102 2.002590_dbl_kind, 2.003627_dbl_kind, 2.005276_dbl_kind, &
02103 2.002395_dbl_kind, 2.003391_dbl_kind, 2.004904_dbl_kind, &
02104 2.002071_dbl_kind, 2.002922_dbl_kind, 2.004241_dbl_kind/
02105
02106
02107 data ws_tab/ &
02108 0.9999994_dbl_kind, 0.9999673_dbl_kind, 0.9954589_dbl_kind, &
02109 0.9999992_dbl_kind, 0.9999547_dbl_kind, 0.9938576_dbl_kind, &
02110 0.9999990_dbl_kind, 0.9999382_dbl_kind, 0.9917989_dbl_kind, &
02111 0.9999985_dbl_kind, 0.9999123_dbl_kind, 0.9889724_dbl_kind, &
02112 0.9999979_dbl_kind, 0.9998844_dbl_kind, 0.9866190_dbl_kind, &
02113 0.9999970_dbl_kind, 0.9998317_dbl_kind, 0.9823021_dbl_kind, &
02114 0.9999960_dbl_kind, 0.9997800_dbl_kind, 0.9785269_dbl_kind, &
02115 0.9999951_dbl_kind, 0.9997288_dbl_kind, 0.9751601_dbl_kind, &
02116 0.9999936_dbl_kind, 0.9996531_dbl_kind, 0.9706974_dbl_kind, &
02117 0.9999922_dbl_kind, 0.9995783_dbl_kind, 0.9667577_dbl_kind, &
02118 0.9999903_dbl_kind, 0.9994798_dbl_kind, 0.9621007_dbl_kind, &
02119 0.9999885_dbl_kind, 0.9993825_dbl_kind, 0.9579541_dbl_kind, &
02120 0.9999866_dbl_kind, 0.9992862_dbl_kind, 0.9541924_dbl_kind, &
02121 0.9999838_dbl_kind, 0.9991434_dbl_kind, 0.9490959_dbl_kind, &
02122 0.9999810_dbl_kind, 0.9990025_dbl_kind, 0.9444940_dbl_kind, &
02123 0.9999772_dbl_kind, 0.9988171_dbl_kind, 0.9389141_dbl_kind, &
02124 0.9999726_dbl_kind, 0.9985890_dbl_kind, 0.9325819_dbl_kind, &
02125 0.9999670_dbl_kind, 0.9983199_dbl_kind, 0.9256405_dbl_kind, &
02126 0.9999605_dbl_kind, 0.9980117_dbl_kind, 0.9181533_dbl_kind, &
02127 0.9999530_dbl_kind, 0.9976663_dbl_kind, 0.9101540_dbl_kind, &
02128 0.9999465_dbl_kind, 0.9973693_dbl_kind, 0.9035031_dbl_kind, &
02129 0.9999382_dbl_kind, 0.9969939_dbl_kind, 0.8953134_dbl_kind, &
02130 0.9999289_dbl_kind, 0.9965848_dbl_kind, 0.8865789_dbl_kind, &
02131 0.9999188_dbl_kind, 0.9961434_dbl_kind, 0.8773350_dbl_kind, &
02132 0.9999068_dbl_kind, 0.9956323_dbl_kind, 0.8668233_dbl_kind, &
02133 0.9998975_dbl_kind, 0.9952464_dbl_kind, 0.8589990_dbl_kind, &
02134 0.9998837_dbl_kind, 0.9946782_dbl_kind, 0.8476493_dbl_kind, &
02135 0.9998699_dbl_kind, 0.9941218_dbl_kind, 0.8367318_dbl_kind, &
02136 0.9998515_dbl_kind, 0.9933966_dbl_kind, 0.8227881_dbl_kind, &
02137 0.9998332_dbl_kind, 0.9926888_dbl_kind, 0.8095131_dbl_kind, &
02138 0.9998148_dbl_kind, 0.9919968_dbl_kind, 0.7968620_dbl_kind, &
02139 0.9997691_dbl_kind, 0.9903277_dbl_kind, 0.7677887_dbl_kind/
02140
02141
02142 data gs_tab / &
02143 0.859913_dbl_kind, 0.848003_dbl_kind, 0.824415_dbl_kind, &
02144 0.867130_dbl_kind, 0.858150_dbl_kind, 0.848445_dbl_kind, &
02145 0.873381_dbl_kind, 0.867221_dbl_kind, 0.861714_dbl_kind, &
02146 0.878368_dbl_kind, 0.874879_dbl_kind, 0.874036_dbl_kind, &
02147 0.881462_dbl_kind, 0.879661_dbl_kind, 0.881299_dbl_kind, &
02148 0.884361_dbl_kind, 0.883903_dbl_kind, 0.890184_dbl_kind, &
02149 0.885937_dbl_kind, 0.886256_dbl_kind, 0.895393_dbl_kind, &
02150 0.886931_dbl_kind, 0.887769_dbl_kind, 0.899072_dbl_kind, &
02151 0.887894_dbl_kind, 0.889255_dbl_kind, 0.903285_dbl_kind, &
02152 0.888515_dbl_kind, 0.890236_dbl_kind, 0.906588_dbl_kind, &
02153 0.889073_dbl_kind, 0.891127_dbl_kind, 0.910152_dbl_kind, &
02154 0.889452_dbl_kind, 0.891750_dbl_kind, 0.913100_dbl_kind, &
02155 0.889730_dbl_kind, 0.892213_dbl_kind, 0.915621_dbl_kind, &
02156 0.890026_dbl_kind, 0.892723_dbl_kind, 0.918831_dbl_kind, &
02157 0.890238_dbl_kind, 0.893099_dbl_kind, 0.921540_dbl_kind, &
02158 0.890441_dbl_kind, 0.893474_dbl_kind, 0.924581_dbl_kind, &
02159 0.890618_dbl_kind, 0.893816_dbl_kind, 0.927701_dbl_kind, &
02160 0.890762_dbl_kind, 0.894123_dbl_kind, 0.930737_dbl_kind, &
02161 0.890881_dbl_kind, 0.894397_dbl_kind, 0.933568_dbl_kind, &
02162 0.890975_dbl_kind, 0.894645_dbl_kind, 0.936148_dbl_kind, &
02163 0.891035_dbl_kind, 0.894822_dbl_kind, 0.937989_dbl_kind, &
02164 0.891097_dbl_kind, 0.895020_dbl_kind, 0.939949_dbl_kind, &
02165 0.891147_dbl_kind, 0.895212_dbl_kind, 0.941727_dbl_kind, &
02166 0.891189_dbl_kind, 0.895399_dbl_kind, 0.943339_dbl_kind, &
02167 0.891225_dbl_kind, 0.895601_dbl_kind, 0.944915_dbl_kind, &
02168 0.891248_dbl_kind, 0.895745_dbl_kind, 0.945950_dbl_kind, &
02169 0.891277_dbl_kind, 0.895951_dbl_kind, 0.947288_dbl_kind, &
02170 0.891299_dbl_kind, 0.896142_dbl_kind, 0.948438_dbl_kind, &
02171 0.891323_dbl_kind, 0.896388_dbl_kind, 0.949762_dbl_kind, &
02172 0.891340_dbl_kind, 0.896623_dbl_kind, 0.950916_dbl_kind, &
02173 0.891356_dbl_kind, 0.896851_dbl_kind, 0.951945_dbl_kind, &
02174 0.891386_dbl_kind, 0.897399_dbl_kind, 0.954156_dbl_kind/
02175
02176
02177 data ki_ssl_mn / 1000.1_dbl_kind, 1003.7_dbl_kind, 7042._dbl_kind/
02178 data wi_ssl_mn / .9999_dbl_kind, .9963_dbl_kind, .9088_dbl_kind/
02179 data gi_ssl_mn / .94_dbl_kind, .94_dbl_kind, .94_dbl_kind/
02180
02181
02182 data ki_dl_mn / 100.2_dbl_kind, 107.7_dbl_kind, 1309._dbl_kind /
02183 data wi_dl_mn / .9980_dbl_kind, .9287_dbl_kind, .0305_dbl_kind /
02184 data gi_dl_mn / .94_dbl_kind, .94_dbl_kind, .94_dbl_kind /
02185
02186
02187 data ki_int_mn / 20.2_dbl_kind, 27.7_dbl_kind, 1445._dbl_kind /
02188 data wi_int_mn / .9901_dbl_kind, .7223_dbl_kind, .0277_dbl_kind /
02189 data gi_int_mn / .94_dbl_kind, .94_dbl_kind, .94_dbl_kind /
02190
02191
02192 data ki_p_ssl_mn / 70.2_dbl_kind, 77.7_dbl_kind, 1309._dbl_kind/
02193 data wi_p_ssl_mn / .9972_dbl_kind, .9009_dbl_kind, .0305_dbl_kind/
02194 data gi_p_ssl_mn / .94_dbl_kind, .94_dbl_kind, .94_dbl_kind /
02195
02196
02197 data ki_p_int_mn / 20.2_dbl_kind, 27.7_dbl_kind, 1445._dbl_kind/
02198 data wi_p_int_mn / .9901_dbl_kind, .7223_dbl_kind, .0277_dbl_kind/
02199 data gi_p_int_mn / .94_dbl_kind, .94_dbl_kind, .94_dbl_kind /
02200
02201
02202 data kw / 0.20_dbl_kind, 12.0_dbl_kind, 729._dbl_kind /
02203 data ww / 0.00_dbl_kind, 0.00_dbl_kind, 0.00_dbl_kind /
02204 data gw / 0.00_dbl_kind, 0.00_dbl_kind, 0.00_dbl_kind /
02205
02206
02207 data hs_ssl / 0.040_dbl_kind /
02208 data rhoi /917.0_dbl_kind /
02209 data fr_max / 1.00_dbl_kind /
02210 data fr_min / 0.80_dbl_kind /
02211
02212
02213 data hi_ssl / 0.050_dbl_kind /
02214
02215
02216 data kalg / 0.00_dbl_kind /
02217
02218
02219 data fp_ice / 0.15_dbl_kind /
02220 data fm_ice / 0.15_dbl_kind /
02221 data fp_pnd / 2.00_dbl_kind /
02222 data fm_pnd / 0.50_dbl_kind /
02223
02224
02225 data hpmin / .005_dbl_kind /
02226 data hp0 / .200_dbl_kind /
02227
02228
02229
02230
02231 data kaer_tab/ &
02232 11580.61872, 5535.41835, 2793.79690, &
02233 25798.96479, 11536.03871, 4688.24207, &
02234 196.49772, 204.14078, 214.42287, &
02235 2665.85867, 2256.71027, 820.36024, &
02236 840.78295, 1028.24656, 1163.03298, &
02237 387.51211, 414.68808, 450.29814/
02238 data waer_tab/ &
02239 0.29003, 0.17349, 0.06613, &
02240 0.51731, 0.41609, 0.21324, &
02241 0.84467, 0.94216, 0.95666, &
02242 0.97764, 0.99402, 0.98552, &
02243 0.94146, 0.98527, 0.99093, &
02244 0.90034, 0.96543, 0.97678/
02245 data gaer_tab/ &
02246 0.35445, 0.19838, 0.08857, &
02247 0.52581, 0.32384, 0.14970, &
02248 0.83162, 0.78306, 0.74375, &
02249 0.68861, 0.70836, 0.54171, &
02250 0.70239, 0.66115, 0.71983, &
02251 0.78734, 0.73580, 0.64411/
02252
02253
02254
02255
02256
02257
02258
02259
02260
02261
02262
02263
02264
02265
02266
02267
02268
02269
02270
02271
02272
02273
02274
02275
02276
02277 rnilyr = real(nilyr,kind=dbl_kind)
02278 rnslyr = real(nslyr,kind=dbl_kind)
02279
02280
02281 do ij = 1, icells_DE
02282 avdr(ij) = c0
02283 avdf(ij) = c0
02284 aidr(ij) = c0
02285 aidf(ij) = c0
02286 fsfc(ij) = c0
02287 fint(ij) = c0
02288 fthru(ij) = c0
02289 enddo
02290 Sabs(:,:) = c0
02291 Iabs(:,:) = c0
02292
02293
02294
02295
02296
02297 do ij = 1, icells_DE
02298 i = indxi_DE(ij)
02299 j = indxj_DE(ij)
02300 wghtns(ij,1) = c1
02301 wghtns(ij,2) = cp67 + (cp78-cp67)*(c1-fnidr(i,j))
02302 wghtns(ij,3) = cp33 + (cp22-cp33)*(c1-fnidr(i,j))
02303 enddo
02304
02305
02306
02307
02308
02309
02310 if( R_ice >= c0 ) then
02311 do ns = 1, nspint
02312 sigp = ki_ssl_mn(ns)*wi_ssl_mn(ns)*(c1+fp_ice*R_ice)
02313 ki_ssl(ns) = sigp+ki_ssl_mn(ns)*(c1-wi_ssl_mn(ns))
02314 wi_ssl(ns) = sigp/ki_ssl(ns)
02315 gi_ssl(ns) = gi_ssl_mn(ns)
02316
02317 sigp = ki_dl_mn(ns)*wi_dl_mn(ns)*(c1+fp_ice*R_ice)
02318 ki_dl(ns) = sigp+ki_dl_mn(ns)*(c1-wi_dl_mn(ns))
02319 wi_dl(ns) = sigp/ki_dl(ns)
02320 gi_dl(ns) = gi_dl_mn(ns)
02321
02322 sigp = ki_int_mn(ns)*wi_int_mn(ns)*(c1+fp_ice*R_ice)
02323 ki_int(ns) = sigp+ki_int_mn(ns)*(c1-wi_int_mn(ns))
02324 wi_int(ns) = sigp/ki_int(ns)
02325 gi_int(ns) = gi_int_mn(ns)
02326 enddo
02327 else
02328 do ns = 1, nspint
02329 sigp = ki_ssl_mn(ns)*wi_ssl_mn(ns)*(c1+fm_ice*R_ice)
02330 sigp = max(sigp, c0)
02331 ki_ssl(ns) = sigp+ki_ssl_mn(ns)*(c1-wi_ssl_mn(ns))
02332 wi_ssl(ns) = sigp/ki_ssl(ns)
02333 gi_ssl(ns) = gi_ssl_mn(ns)
02334
02335 sigp = ki_dl_mn(ns)*wi_dl_mn(ns)*(c1+fm_ice*R_ice)
02336 sigp = max(sigp, c0)
02337 ki_dl(ns) = sigp+ki_dl_mn(ns)*(c1-wi_dl_mn(ns))
02338 wi_dl(ns) = sigp/ki_dl(ns)
02339 gi_dl(ns) = gi_dl_mn(ns)
02340
02341 sigp = ki_int_mn(ns)*wi_int_mn(ns)*(c1+fm_ice*R_ice)
02342 sigp = max(sigp, c0)
02343 ki_int(ns) = sigp+ki_int_mn(ns)*(c1-wi_int_mn(ns))
02344 wi_int(ns) = sigp/ki_int(ns)
02345 gi_int(ns) = gi_int_mn(ns)
02346 enddo
02347 endif
02348
02349
02350 if( R_pnd >= c0 ) then
02351 do ns = 1, nspint
02352 sigp = ki_p_ssl_mn(ns)*wi_p_ssl_mn(ns)*(c1+fp_pnd*R_pnd)
02353 ki_p_ssl(ns) = sigp+ki_p_ssl_mn(ns)*(c1-wi_p_ssl_mn(ns))
02354 wi_p_ssl(ns) = sigp/ki_p_ssl(ns)
02355 gi_p_ssl(ns) = gi_p_ssl_mn(ns)
02356
02357 sigp = ki_p_int_mn(ns)*wi_p_int_mn(ns)*(c1+fp_pnd*R_pnd)
02358 ki_p_int(ns) = sigp+ki_p_int_mn(ns)*(c1-wi_p_int_mn(ns))
02359 wi_p_int(ns) = sigp/ki_p_int(ns)
02360 gi_p_int(ns) = gi_p_int_mn(ns)
02361 enddo
02362 else
02363 do ns = 1, nspint
02364 sigp = ki_p_ssl_mn(ns)*wi_p_ssl_mn(ns)*(c1+fm_pnd*R_pnd)
02365 sigp = max(sigp, c0)
02366 ki_p_ssl(ns) = sigp+ki_p_ssl_mn(ns)*(c1-wi_p_ssl_mn(ns))
02367 wi_p_ssl(ns) = sigp/ki_p_ssl(ns)
02368 gi_p_ssl(ns) = gi_p_ssl_mn(ns)
02369
02370 sigp = ki_p_int_mn(ns)*wi_p_int_mn(ns)*(c1+fm_pnd*R_pnd)
02371 sigp = max(sigp, c0)
02372 ki_p_int(ns) = sigp+ki_p_int_mn(ns)*(c1-wi_p_int_mn(ns))
02373 wi_p_int(ns) = sigp/ki_p_int(ns)
02374 gi_p_int(ns) = gi_p_int_mn(ns)
02375 enddo
02376 endif
02377
02378
02379
02380
02381 do ns = 1, nspint
02382
02383
02384 do ij = 1, icells_DE
02385 i = indxi_DE(ij)
02386 j = indxj_DE(ij)
02387
02388 if( srftyp(i,j) == 0 ) then
02389 do k=0,nslyr
02390 tau(k,ij) = c0
02391 w0(k,ij) = c0
02392 g(k,ij) = c0
02393 enddo
02394
02395 else if( srftyp(i,j) == 1 ) then
02396 dz_ssl = hs_ssl
02397 dz = hs(i,j)/rnslyr
02398
02399 dz_ssl = min(dz_ssl, dz/c2)
02400
02401
02402
02403
02404 fr = fr_max*fnidr(i,j) + fr_min*(c1-fnidr(i,j))
02405
02406
02407 ksnow = 1
02408 do k=0,nslyr
02409
02410 if( k > 1 ) ksnow = k
02411
02412 if( fr*rsnw(i,j,ksnow) < rsnw_tab(1) ) then
02413 Qs(ns) = Qs_tab(ns,1)
02414 ws(ns) = ws_tab(ns,1)
02415 gs(ns) = gs_tab(ns,1)
02416 else if( fr*rsnw(i,j,ksnow) >= rsnw_tab(nmbrad) ) then
02417 Qs(ns) = Qs_tab(ns,nmbrad)
02418 ws(ns) = ws_tab(ns,nmbrad)
02419 gs(ns) = gs_tab(ns,nmbrad)
02420 else
02421
02422 do nr=2,nmbrad
02423 if( rsnw_tab(nr-1) <= fr*rsnw(i,j,ksnow) .and. &
02424 fr*rsnw(i,j,ksnow) < rsnw_tab(nr)) then
02425 delr = (fr*rsnw(i,j,ksnow) - rsnw_tab(nr-1)) / &
02426 (rsnw_tab(nr) - rsnw_tab(nr-1))
02427 Qs(ns) = Qs_tab(ns,nr-1)*(c1-delr) + &
02428 Qs_tab(ns,nr)*delr
02429 ws(ns) = ws_tab(ns,nr-1)*(c1-delr) + &
02430 ws_tab(ns,nr)*delr
02431 gs(ns) = gs_tab(ns,nr-1)*(c1-delr) + &
02432 gs_tab(ns,nr)*delr
02433 endif
02434 enddo
02435 endif
02436 ks(ns) = Qs(ns)*((rhosnw(i,j,ksnow)/rhoi)*3._dbl_kind / &
02437 (4._dbl_kind*fr*rsnw(i,j,ksnow)*1.0e-6_dbl_kind))
02438 if( k == 0 ) then
02439 tau(k,ij) = ks(ns)*dz_ssl
02440 else if( k == 1 ) then
02441 tau(k,ij) = ks(ns)*(dz-dz_ssl)
02442 else
02443 tau(k,ij) = ks(ns)*dz
02444 endif
02445 w0(k,ij) = ws(ns)
02446 g(k,ij) = gs(ns)
02447
02448 nmbaer_actual = min(n_aero,nmbaer)
02449 if( k == 0 ) then
02450 taer = c0
02451 waer = c0
02452 gaer = c0
02453 do na=1,4*nmbaer_actual,4
02454 taer = taer + &
02455 aero_mp(i,j,na)*kaer_tab(ns,(1+(na-1)/4))
02456 waer = waer + &
02457 aero_mp(i,j,na)*kaer_tab(ns,(1+(na-1)/4))* &
02458 waer_tab(ns,(1+(na-1)/4))
02459 gaer = gaer + &
02460 aero_mp(i,j,na)*kaer_tab(ns,(1+(na-1)/4))* &
02461 waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
02462 enddo
02463 gaer = gaer/(waer+puny)
02464 waer = waer/(taer+puny)
02465 else if ( k > 0 ) then
02466 taer = c0
02467 waer = c0
02468 gaer = c0
02469 do na=1,4*nmbaer_actual,4
02470 taer = taer + &
02471 (aero_mp(i,j,na+1)/rnslyr)*kaer_tab(ns,(1+(na-1)/4))
02472 waer = waer + &
02473 (aero_mp(i,j,na+1)/rnslyr)*kaer_tab(ns,(1+(na-1)/4))* &
02474 waer_tab(ns,(1+(na-1)/4))
02475 gaer = gaer + &
02476 (aero_mp(i,j,na+1)/rnslyr)*kaer_tab(ns,(1+(na-1)/4))* &
02477 waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
02478 enddo
02479 gaer = gaer/(waer+puny)
02480 waer = waer/(taer+puny)
02481 endif
02482 g(k,ij) = (g(k,ij)*w0(k,ij)*tau(k,ij) + gaer*waer*taer) / &
02483 (w0(k,ij)*tau(k,ij) + waer*taer)
02484 w0(k,ij) = (w0(k,ij)*tau(k,ij) + waer*taer) / &
02485 (tau(k,ij) + taer)
02486 tau(k,ij) = tau(k,ij) + taer
02487 enddo
02488
02489 else
02490
02491 dz = hp(i,j)/(rnslyr+c1)
02492 do k=0,nslyr
02493 tau(k,ij) = kw(ns)*dz
02494 w0(k,ij) = ww(ns)
02495 g(k,ij) = gw(ns)
02496
02497 enddo
02498 endif
02499 enddo
02500
02501
02502 kii = nslyr + 1
02503 do ij = 1, icells_DE
02504 i = indxi_DE(ij)
02505 j = indxj_DE(ij)
02506 dz_ssl = hi_ssl
02507 dz = hi(i,j)/rnilyr
02508
02509
02510 if( hi(i,j) < 1.5_dbl_kind ) dz_ssl = hi(i,j)/30._dbl_kind
02511
02512 dz_ssl = min(dz_ssl, dz/c2)
02513
02514 if( srftyp(i,j) <= 1 ) then
02515
02516 k = kii
02517 tau(k,ij) = ki_ssl(ns)*dz_ssl
02518 w0(k,ij) = wi_ssl(ns)
02519 g(k,ij) = gi_ssl(ns)
02520
02521 k = kii + 1
02522
02523 fs = rnilyr/c4
02524 tau(k,ij) = ki_dl(ns)*(dz-dz_ssl)*fs
02525 w0(k,ij) = wi_dl(ns)
02526 g(k,ij) = gi_dl(ns)
02527
02528 if (kii+2 <= klev-1) then
02529 do k = kii+2, klev-1
02530 tau(k,ij) = ki_int(ns)*dz
02531 w0(k,ij) = wi_int(ns)
02532 g(k,ij) = gi_int(ns)
02533 enddo
02534 endif
02535
02536 k = klev
02537
02538 kabs = ki_int(ns)*(c1-wi_int(ns))
02539 if( ns == 1 ) then
02540
02541
02542 kabs = kabs + kalg*(0.50_dbl_kind/dz)
02543 endif
02544 sig = ki_int(ns)*wi_int(ns)
02545 tau(k,ij) = (kabs+sig)*dz
02546 w0(k,ij) = (sig/(sig+kabs))
02547 g(k,ij) = gi_int(ns)
02548
02549 nmbaer_actual = min(n_aero,nmbaer)
02550 do k = kii, klev
02551 if( k == kii ) then
02552 taer = c0
02553 waer = c0
02554 gaer = c0
02555 do na=1,4*nmbaer_actual,4
02556 taer = taer + &
02557 aero_mp(i,j,na+2)*kaer_tab(ns,(1+(na-1)/4))
02558 waer = waer + &
02559 aero_mp(i,j,na+2)*kaer_tab(ns,(1+(na-1)/4))* &
02560 waer_tab(ns,(1+(na-1)/4))
02561 gaer = gaer + &
02562 aero_mp(i,j,na+2)*kaer_tab(ns,(1+(na-1)/4))* &
02563 waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
02564 enddo
02565 gaer = gaer/(waer+puny)
02566 waer = waer/(taer+puny)
02567 else if ( k > kii ) then
02568 taer = c0
02569 waer = c0
02570 gaer = c0
02571 do na=1,4*nmbaer_actual,4
02572 taer = taer + &
02573 (aero_mp(i,j,na+3)/rnilyr)*kaer_tab(ns,(1+(na-1)/4))
02574 waer = waer + &
02575 (aero_mp(i,j,na+3)/rnilyr)*kaer_tab(ns,(1+(na-1)/4))* &
02576 waer_tab(ns,(1+(na-1)/4))
02577 gaer = gaer + &
02578 (aero_mp(i,j,na+3)/rnilyr)*kaer_tab(ns,(1+(na-1)/4))* &
02579 waer_tab(ns,(1+(na-1)/4))*gaer_tab(ns,(1+(na-1)/4))
02580 enddo
02581 gaer = gaer/(waer+puny)
02582 waer = waer/(taer+puny)
02583 endif
02584 g(k,ij) = (g(k,ij)*w0(k,ij)*tau(k,ij) + gaer*waer*taer) / &
02585 (w0(k,ij)*tau(k,ij) + waer*taer)
02586 w0(k,ij) = (w0(k,ij)*tau(k,ij) + waer*taer) / &
02587 (tau(k,ij) + taer)
02588 tau(k,ij) = tau(k,ij) + taer
02589 enddo
02590
02591 else
02592 k = kii
02593 tau(k,ij) = ki_p_ssl(ns)*dz_ssl
02594 w0(k,ij) = wi_p_ssl(ns)
02595 g(k,ij) = gi_p_ssl(ns)
02596 k = kii + 1
02597 tau(k,ij) = ki_p_int(ns)*(dz-dz_ssl)
02598 w0(k,ij) = wi_p_int(ns)
02599 g(k,ij) = gi_p_int(ns)
02600 if (kii+2 <= klev) then
02601 do k = kii+2, klev
02602 tau(k,ij) = ki_p_int(ns)*dz
02603 w0(k,ij) = wi_p_int(ns)
02604 g(k,ij) = gi_p_int(ns)
02605 enddo
02606 endif
02607
02608 if( hpmin <= hp(i,j) .and. hp(i,j) <= hp0 ) then
02609 k = kii
02610 sig_i = ki_ssl(ns)*wi_ssl(ns)
02611 sig_p = ki_p_ssl(ns)*wi_p_ssl(ns)
02612 sig = sig_i + (sig_p-sig_i)*(hp(i,j)/hp0)
02613 kext = sig + ki_p_ssl(ns)*(c1-wi_p_ssl(ns))
02614 tau(k,ij) = kext*dz_ssl
02615 w0(k,ij) = sig/kext
02616 g(k,ij) = gi_p_int(ns)
02617 k = kii + 1
02618
02619 fs = rnilyr/c4
02620 sig_i = ki_dl(ns)*wi_dl(ns)*fs
02621 sig_p = ki_p_int(ns)*wi_p_int(ns)
02622 sig = sig_i + (sig_p-sig_i)*(hp(i,j)/hp0)
02623 kext = sig + ki_p_int(ns)*(c1-wi_p_int(ns))
02624 tau(k,ij) = kext*(dz-dz_ssl)
02625 w0(k,ij) = sig/kext
02626 g(k,ij) = gi_p_int(ns)
02627 if (kii+2 <= klev) then
02628 do k = kii+2, klev
02629 sig_i = ki_int(ns)*wi_int(ns)
02630 sig_p = ki_p_int(ns)*wi_p_int(ns)
02631 sig = sig_i + (sig_p-sig_i)*(hp(i,j)/hp0)
02632 kext = sig + ki_p_int(ns)*(c1-wi_p_int(ns))
02633 tau(k,ij) = kext*dz
02634 w0(k,ij) = sig/kext
02635 g(k,ij) = gi_p_int(ns)
02636 enddo
02637 endif
02638 endif
02639 endif
02640 enddo
02641
02642
02643 if(ns == 1) then
02644 do ij = 1, icells_DE
02645 i = indxi_DE(ij)
02646 j = indxj_DE(ij)
02647 albodr(ij) = cp01
02648 albodf(ij) = cp01
02649 enddo
02650 else
02651 do ij = 1, icells_DE
02652 i = indxi_DE(ij)
02653 j = indxj_DE(ij)
02654 albodr(ij) = c0
02655 albodf(ij) = c0
02656 enddo
02657 endif
02658
02659
02660
02661
02662
02663
02664
02665
02666
02667 call solution_dEdd &
02668 (nx_block, ny_block, &
02669 icells_DE, indxi_DE, indxj_DE, coszen, srftyp, &
02670 tau, w0, g, albodr, albodf, &
02671 trndir, trntdr, trndif, rupdir, rupdif, &
02672 rdndif)
02673
02674
02675
02676
02677
02678
02679
02680
02681
02682
02683
02684
02685 do ij = 1, icells_DE
02686 i = indxi_DE(ij)
02687 j = indxj_DE(ij)
02688 do k=0,klevp
02689
02690 refk = c1/(c1 - rdndif(k,ij)*rupdif(k,ij))
02691
02692
02693 fdirup(k,ij) = (trndir(k,ij)*rupdir(k,ij) + &
02694 (trntdr(k,ij)-trndir(k,ij)) &
02695 *rupdif(k,ij))*refk
02696
02697
02698 fdirdn(k,ij) = trndir(k,ij) + (trntdr(k,ij) &
02699 - trndir(k,ij) + trndir(k,ij) &
02700 *rupdir(k,ij)*rdndif(k,ij))*refk
02701
02702 fdifup(k,ij) = trndif(k,ij)*rupdif(k,ij)*refk
02703
02704 fdifdn(k,ij) = trndif(k,ij)*refk
02705 enddo
02706 enddo
02707
02708
02709
02710
02711 if( ns == 1) then
02712
02713 do ij = 1, icells_DE
02714 i = indxi_DE(ij)
02715 j = indxj_DE(ij)
02716 avdr(ij) = rupdir(0,ij)
02717 avdf(ij) = rupdif(0,ij)
02718
02719
02720 if( srftyp(i,j) == 1 ) then
02721
02722 ksrf = 1
02723 else
02724
02725 ksrf = nslyr + 2
02726 endif
02727
02728 fsfc(ij) = fsfc(ij) + &
02729 ((fdirdn(0,ij)-fdirup(0,ij))*swvdr(i,j) + &
02730 (fdifdn(0,ij)-fdifup(0,ij))*swvdf(i,j)) - &
02731 ((fdirdn(ksrf,ij)-fdirup(ksrf,ij))*swvdr(i,j) + &
02732 (fdifdn(ksrf,ij)-fdifup(ksrf,ij))*swvdf(i,j))
02733
02734 fint(ij) = fint(ij) + &
02735 ((fdirdn(ksrf,ij)-fdirup(ksrf,ij))*swvdr(i,j) + &
02736 (fdifdn(ksrf,ij)-fdifup(ksrf,ij))*swvdf(i,j)) - &
02737 ((fdirdn(klevp,ij)-fdirup(klevp,ij))*swvdr(i,j) + &
02738 (fdifdn(klevp,ij)-fdifup(klevp,ij))*swvdf(i,j))
02739
02740 fthru(ij) = fthru(ij) + &
02741 (fdirdn(klevp,ij)-fdirup(klevp,ij))*swvdr(i,j) + &
02742 (fdifdn(klevp,ij)-fdifup(klevp,ij))*swvdf(i,j)
02743
02744
02745 if( srftyp(i,j) == 1 ) then
02746 ksa = 0
02747 do k=1,nslyr
02748
02749
02750 km = k
02751 kp = km + 1
02752 ksa = ksa + 1
02753 Sabs(ij,ksa) = Sabs(ij,ksa) + &
02754 ((fdirdn(km,ij)-fdirup(km,ij))*swvdr(i,j) + &
02755 (fdifdn(km,ij)-fdifup(km,ij))*swvdf(i,j)) - &
02756 ((fdirdn(kp,ij)-fdirup(kp,ij))*swvdr(i,j) + &
02757 (fdifdn(kp,ij)-fdifup(kp,ij))*swvdf(i,j))
02758 enddo
02759 endif
02760
02761
02762 ki = 0
02763 do k=nslyr+2,nslyr+1+nilyr
02764
02765 km = k
02766 kp = km + 1
02767
02768 if( srftyp(i,j) == 1 ) then
02769
02770 if( k == nslyr+2 ) then
02771 km = k - 1
02772 kp = km + 2
02773 endif
02774 endif
02775 ki = ki + 1
02776 Iabs(ij,ki) = Iabs(ij,ki) + &
02777 ((fdirdn(km,ij)-fdirup(km,ij))*swvdr(i,j) + &
02778 (fdifdn(km,ij)-fdifup(km,ij))*swvdf(i,j)) - &
02779 ((fdirdn(kp,ij)-fdirup(kp,ij))*swvdr(i,j) + &
02780 (fdifdn(kp,ij)-fdifup(kp,ij))*swvdf(i,j))
02781 enddo
02782 enddo
02783
02784 else
02785
02786 do ij = 1, icells_DE
02787 i = indxi_DE(ij)
02788 j = indxj_DE(ij)
02789
02790
02791
02792
02793
02794
02795
02796
02797 aidr(ij) = aidr(ij) + rupdir(0,ij)*wghtns(ij,ns)
02798 aidf(ij) = aidf(ij) + rupdif(0,ij)*wghtns(ij,ns)
02799
02800
02801 if( srftyp(i,j) == 1 ) then
02802
02803 ksrf = 1
02804 else
02805
02806 ksrf = nslyr + 2
02807 endif
02808
02809 fsfc(ij) = fsfc(ij) + &
02810 ( ((fdirdn(0,ij)-fdirup(0,ij))*swidr(i,j) + &
02811 (fdifdn(0,ij)-fdifup(0,ij))*swidf(i,j)) - &
02812 ((fdirdn(ksrf,ij)-fdirup(ksrf,ij))*swidr(i,j) + &
02813 (fdifdn(ksrf,ij)-fdifup(ksrf,ij))*swidf(i,j)) ) &
02814 *wghtns(ij,ns)
02815
02816 fint(ij) = fint(ij) + &
02817 ( ((fdirdn(ksrf,ij)-fdirup(ksrf,ij))*swidr(i,j) + &
02818 (fdifdn(ksrf,ij)-fdifup(ksrf,ij))*swidf(i,j)) - &
02819 ((fdirdn(klevp,ij)-fdirup(klevp,ij))*swidr(i,j) + &
02820 (fdifdn(klevp,ij)-fdifup(klevp,ij))*swidf(i,j)) ) &
02821 *wghtns(ij,ns)
02822
02823 fthru(ij) = fthru(ij) + &
02824 ((fdirdn(klevp,ij)-fdirup(klevp,ij))*swidr(i,j) + &
02825 (fdifdn(klevp,ij)-fdifup(klevp,ij))*swidf(i,j)) &
02826 *wghtns(ij,ns)
02827
02828
02829 if( srftyp(i,j) == 1 ) then
02830 ksa = 0
02831 do k=1,nslyr
02832
02833
02834 km = k
02835 kp = km + 1
02836 ksa = ksa + 1
02837 Sabs(ij,ksa) = Sabs(ij,ksa) + &
02838 ( ((fdirdn(km,ij)-fdirup(km,ij))*swidr(i,j) + &
02839 (fdifdn(km,ij)-fdifup(km,ij))*swidf(i,j)) - &
02840 ((fdirdn(kp,ij)-fdirup(kp,ij))*swidr(i,j) + &
02841 (fdifdn(kp,ij)-fdifup(kp,ij))*swidf(i,j)) ) &
02842 *wghtns(ij,ns)
02843 enddo
02844 endif
02845
02846
02847 ki = 0
02848 do k=nslyr+2,nslyr+1+nilyr
02849
02850 km = k
02851 kp = km + 1
02852
02853 if( srftyp(i,j) == 1 ) then
02854
02855 if( k == nslyr+2 ) then
02856 km = k - 1
02857 kp = km + 2
02858 endif
02859 endif
02860 ki = ki + 1
02861 Iabs(ij,ki) = Iabs(ij,ki) + &
02862 ( ((fdirdn(km,ij)-fdirup(km,ij))*swidr(i,j) + &
02863 (fdifdn(km,ij)-fdifup(km,ij))*swidf(i,j)) - &
02864 ((fdirdn(kp,ij)-fdirup(kp,ij))*swidr(i,j) + &
02865 (fdifdn(kp,ij)-fdifup(kp,ij))*swidf(i,j)) ) &
02866 *wghtns(ij,ns)
02867 enddo
02868 enddo
02869
02870 endif
02871
02872 enddo
02873
02874
02875
02876
02877
02878 do ij = 1, icells_DE
02879 i = indxi_DE(ij)
02880 j = indxj_DE(ij)
02881 alvdr(i,j) = avdr(ij)
02882 alvdf(i,j) = avdf(ij)
02883 alidr(i,j) = aidr(ij)
02884 alidf(i,j) = aidf(ij)
02885 fswsfc(i,j) = fswsfc(i,j) + fsfc(ij) *fi(i,j)
02886 fswint(i,j) = fswint(i,j) + fint(ij) *fi(i,j)
02887 fswthru(i,j) = fswthru(i,j) + fthru(ij)*fi(i,j)
02888 enddo
02889
02890 do k = 1, nslyr
02891
02892
02893
02894 do ij = 1, icells_DE
02895 i = indxi_DE(ij)
02896 j = indxj_DE(ij)
02897 Sswabs(i,j,k) = Sswabs(i,j,k) + Sabs(ij,k)*fi(i,j)
02898 enddo
02899 enddo
02900
02901 do k = 1, nilyr
02902
02903
02904
02905 do ij = 1, icells_DE
02906 i = indxi_DE(ij)
02907 j = indxj_DE(ij)
02908 Iswabs(i,j,k) = Iswabs(i,j,k) + Iabs(ij,k)*fi(i,j)
02909 enddo
02910 enddo
02911
02912
02913
02914
02915
02916
02917
02918 if (.not. heat_capacity) then
02919
02920
02921
02922
02923 do ij = 1, icells_DE
02924 i = indxi_DE(ij)
02925 j = indxj_DE(ij)
02926
02927
02928 fswsfc(i,j) = fswsfc(i,j) + Iswabs(i,j,1) + Sswabs(i,j,1)
02929
02930
02931 fswint(i,j) = c0
02932 Iswabs(i,j,1) = c0
02933 Sswabs(i,j,1) = c0
02934
02935 enddo
02936
02937 endif
02938
02939 end subroutine compute_dEdd
02940
02941
02942
02943
02944
02945
02946
02947
02948 subroutine solution_dEdd &
02949 (nx_block, ny_block, &
02950 icells_DE, indxi_DE, indxj_DE, coszen, srftyp, &
02951 tau, w0, g, albodr, albodf, &
02952 trndir, trntdr, trndif, rupdir, rupdif, &
02953 rdndif)
02954
02955
02956
02957
02958
02959
02960
02961
02962
02963
02964
02965
02966
02967
02968
02969 integer (kind=int_kind),
02970 intent(in) ::
02971 nx_block, ny_block,
02972 icells_DE
02973
02974 integer (kind=int_kind), dimension(nx_block*ny_block),
02975 intent(in) ::
02976 indxi_DE,
02977 indxj_DE
02978
02979 real (kind=dbl_kind), dimension (nx_block,ny_block),
02980 intent(in) ::
02981 coszen
02982
02983 integer (kind=int_kind), dimension(nx_block,ny_block),
02984 intent(in) ::
02985 srftyp
02986
02987 integer (kind=int_kind), parameter ::
02988 klev = nslyr + nilyr + 1 ,
02989 klevp = klev + 1
02990
02991 real (kind=dbl_kind), dimension(0:klev,icells_DE),
02992 intent(in) ::
02993 tau ,
02994 w0 ,
02995 g
02996
02997 real (kind=dbl_kind), dimension(icells_DE),
02998 intent(in) ::
02999 albodr ,
03000 albodf
03001
03002
03003
03004 real (kind=dbl_kind), dimension (0:klevp,icells_DE),
03005 intent(out) ::
03006 trndir ,
03007 trntdr ,
03008 trndif ,
03009 rupdir ,
03010 rupdif ,
03011 rdndif
03012
03013
03014
03015
03016
03017
03018
03019
03020
03021
03022
03023
03024
03025
03026
03027
03028
03029
03030
03031
03032
03033
03034
03035
03036
03037
03038
03039
03040
03041
03042
03043
03044
03045
03046
03047
03048
03049
03050
03051
03052
03053
03054
03055
03056
03057
03058
03059
03060
03061
03062
03063
03064
03065
03066
03067 integer (kind=int_kind) ::
03068 kfrsnl
03069
03070
03071
03072
03073
03074
03075 real (kind=dbl_kind), dimension (0:klev,icells_DE) ::
03076 rdir ,
03077 rdif_a ,
03078 rdif_b ,
03079 tdir ,
03080 tdif_a ,
03081 tdif_b ,
03082 trnlay
03083
03084 integer (kind=int_kind) ::
03085 i ,
03086 j ,
03087 ij ,
03088 k
03089
03090 real (kind=dbl_kind), parameter ::
03091 trmin = 0.001_dbl_kind
03092
03093
03094
03095 real (kind=dbl_kind) ::
03096 tautot ,
03097 wtot ,
03098 gtot ,
03099 ftot ,
03100 ts ,
03101 ws ,
03102 gs ,
03103 rintfc ,
03104 refkp1 ,
03105 refkm1 ,
03106 tdrrdir ,
03107 tdndif
03108
03109
03110 real (kind=dbl_kind) ::
03111 R1 ,
03112 R2 ,
03113 T1 ,
03114 T2 ,
03115 Rf_dir_a ,
03116 Tf_dir_a ,
03117 Rf_dif_a ,
03118 Rf_dif_b ,
03119 Tf_dif_a ,
03120 Tf_dif_b
03121
03122
03123
03124 real (kind=dbl_kind), parameter ::
03125 refindx = 1.310_dbl_kind ,
03126 cp063 = 0.063_dbl_kind ,
03127 cp455 = 0.455_dbl_kind
03128
03129 real (kind=dbl_kind) ::
03130 mu0 ,
03131 mu0n
03132
03133 real (kind=dbl_kind) ::
03134 alpha ,
03135 gamma ,
03136 el ,
03137 taus ,
03138 omgs ,
03139 asys ,
03140 u ,
03141 n ,
03142 lm ,
03143 mu ,
03144 ne
03145
03146 real (kind=dbl_kind) ::
03147 w ,
03148 uu ,
03149 gg ,
03150 e ,
03151 f ,
03152 t ,
03153 et
03154
03155 real (kind=dbl_kind) ::
03156 alp ,
03157 gam ,
03158 ue ,
03159 arg ,
03160 extins ,
03161 amg ,
03162 apg
03163
03164 integer (kind=int_kind), parameter ::
03165 ngmax = 8
03166
03167 real (kind=dbl_kind), dimension (ngmax) ::
03168 gauspt ,
03169 gauswt
03170
03171 data gauspt/ &
03172 .9894009_dbl_kind, .9445750_dbl_kind, &
03173 .8656312_dbl_kind, .7554044_dbl_kind, &
03174 .6178762_dbl_kind, .4580168_dbl_kind, &
03175 .2816036_dbl_kind, .0950125_dbl_kind /
03176 data gauswt/ &
03177 .0271525_dbl_kind, .0622535_dbl_kind, &
03178 .0951585_dbl_kind, .1246290_dbl_kind, &
03179 .1495960_dbl_kind, .1691565_dbl_kind, &
03180 .1826034_dbl_kind, .1894506_dbl_kind /
03181
03182 integer (kind=int_kind) ::
03183 ng
03184
03185 real (kind=dbl_kind) ::
03186 gwt ,
03187 swt ,
03188 trn ,
03189 rdr ,
03190 tdr ,
03191 smr ,
03192 smt
03193
03194
03195 alpha(w,uu,gg,e) = p75*w*uu*((c1 + gg*(c1-w))/(c1 - e*e*uu*uu))
03196 gamma(w,uu,gg,e) = p5*w*((c1 + c3*gg*(c1-w)*uu*uu) &
03197 / (c1-e*e*uu*uu))
03198 n(uu,et) = ((uu+c1)*(uu+c1)/et ) - ((uu-c1)*(uu-c1)*et)
03199 u(w,gg,e) = c1p5*(c1 - w*gg)/e
03200 el(w,gg) = sqrt(c3*(c1-w)*(c1 - w*gg))
03201 taus(w,f,t) = (c1 - w*f)*t
03202 omgs(w,f) = (c1 - f)*w/(c1 - w*f)
03203 asys(gg,f) = (gg - f)/(c1 - f)
03204
03205
03206
03207
03208 do ij = 1, icells_DE
03209 do k = 0, klevp
03210 trndir(k,ij) = c0
03211 trntdr(k,ij) = c0
03212 trndif(k,ij) = c0
03213 rupdir(k,ij) = c0
03214 rupdif(k,ij) = c0
03215 rdndif(k,ij) = c0
03216 enddo
03217
03218
03219 do k = 0, klev
03220 rdir(k,ij) = c0
03221 rdif_a(k,ij) = c0
03222 rdif_b(k,ij) = c0
03223 tdir(k,ij) = c0
03224 tdif_a(k,ij) = c0
03225 tdif_b(k,ij) = c0
03226 trnlay(k,ij) = c0
03227 enddo
03228
03229
03230 trndir(0,ij) = c1
03231 trntdr(0,ij) = c1
03232 trndif(0,ij) = c1
03233 rdndif(0,ij) = c0
03234 enddo
03235
03236
03237
03238
03239
03240 do ij = 1, icells_DE
03241 i = indxi_DE(ij)
03242 j = indxj_DE(ij)
03243
03244
03245 do k=0,klev
03246
03247
03248
03249
03250 if ( k > 0 ) then
03251
03252
03253
03254
03255
03256
03257
03258
03259
03260
03261
03262
03263 trndir(k,ij) = trndir(k-1,ij)*trnlay(k-1,ij)
03264 refkm1 = c1/(c1 - rdndif(k-1,ij)*rdif_a(k-1,ij))
03265 tdrrdir = trndir(k-1,ij)*rdir(k-1,ij)
03266 tdndif = trntdr(k-1,ij) - trndir(k-1,ij)
03267 trntdr(k,ij) = trndir(k-1,ij)*tdir(k-1,ij) + &
03268 (tdndif + tdrrdir*rdndif(k-1,ij))*refkm1*tdif_a(k-1,ij)
03269 rdndif(k,ij) = rdif_b(k-1,ij) + &
03270 (tdif_b(k-1,ij)*rdndif(k-1,ij)*refkm1*tdif_a(k-1,ij))
03271 trndif(k,ij) = trndif(k-1,ij)*refkm1*tdif_a(k-1,ij)
03272 endif
03273
03274
03275
03276 if (trntdr(k,ij) > trmin ) then
03277
03278
03279
03280 tautot = tau(k,ij)
03281 wtot = w0(k,ij)
03282 gtot = g(k,ij)
03283 ftot = gtot*gtot
03284
03285 ts = taus(wtot,ftot,tautot)
03286 ws = omgs(wtot,ftot)
03287 gs = asys(gtot,ftot)
03288 lm = el(ws,gs)
03289 ue = u(ws,gs,lm)
03290
03291
03292 if( srftyp(i,j) < 2 ) then
03293
03294
03295
03296
03297 kfrsnl = nslyr + 2
03298 else
03299
03300 kfrsnl = 0
03301 endif
03302
03303
03304
03305
03306
03307 mu0 = max(coszen(i,j),p01)
03308
03309
03310
03311
03312
03313
03314 mu0n = sqrt(c1-((c1-mu0*mu0)/(refindx*refindx)))
03315
03316
03317
03318
03319 if( srftyp(i,j) < 2 .and. k < kfrsnl ) mu0n = mu0
03320
03321 extins = max(exp_min, exp(-lm*ts))
03322 ne = n(ue,extins)
03323
03324
03325
03326 rdif_a(k,ij) = (ue+c1)*(ue-c1)*(c1/extins - extins)/ne
03327 tdif_a(k,ij) = c4*ue/ne
03328
03329
03330 trnlay(k,ij) = max(exp_min, exp(-ts/mu0n))
03331 alp = alpha(ws,mu0n,gs,lm)
03332 gam = gamma(ws,mu0n,gs,lm)
03333 apg = alp + gam
03334 amg = alp - gam
03335 rdir(k,ij) = amg*(tdif_a(k,ij)*trnlay(k,ij) - c1) + &
03336 apg*rdif_a(k,ij)
03337 tdir(k,ij) = apg*tdif_a(k,ij) + &
03338 (amg*rdif_a(k,ij) - (apg-c1))*trnlay(k,ij)
03339
03340
03341
03342
03343
03344 swt = c0
03345 smr = c0
03346 smt = c0
03347 do ng=1,ngmax
03348 mu = gauspt(ng)
03349 gwt = gauswt(ng)
03350 swt = swt + mu*gwt
03351 trn = max(exp_min, exp(-ts/mu))
03352 alp = alpha(ws,mu,gs,lm)
03353 gam = gamma(ws,mu,gs,lm)
03354 apg = alp + gam
03355 amg = alp - gam
03356 rdr = amg*(tdif_a(k,ij)*trn-c1) + &
03357 apg*rdif_a(k,ij)
03358 tdr = apg*tdif_a(k,ij) + &
03359 (amg*rdif_a(k,ij)-(apg-c1))*trn
03360 smr = smr + mu*rdr*gwt
03361 smt = smt + mu*tdr*gwt
03362 enddo
03363 rdif_a(k,ij) = smr/swt
03364 tdif_a(k,ij) = smt/swt
03365
03366
03367 rdif_b(k,ij) = rdif_a(k,ij)
03368 tdif_b(k,ij) = tdif_a(k,ij)
03369
03370
03371
03372
03373
03374 if( k == kfrsnl ) then
03375
03376
03377
03378 R1 = (mu0 - refindx*mu0n) / &
03379 (mu0 + refindx*mu0n)
03380 R2 = (refindx*mu0 - mu0n) / &
03381 (refindx*mu0 + mu0n)
03382 T1 = c2*mu0 / &
03383 (mu0 + refindx*mu0n)
03384 T2 = c2*mu0 / &
03385 (refindx*mu0 + mu0n)
03386
03387
03388 Rf_dir_a = p5 * (R1*R1 + R2*R2)
03389 Tf_dir_a = p5 * (T1*T1 + T2*T2)*refindx*mu0n/mu0
03390
03391
03392
03393
03394
03395
03396
03397
03398 Rf_dif_a = cp063
03399 Tf_dif_a = c1 - Rf_dif_a
03400
03401 Rf_dif_b = cp455
03402 Tf_dif_b = c1 - Rf_dif_b
03403
03404
03405
03406
03407
03408 rintfc = c1 / (c1-Rf_dif_b*rdif_a(kfrsnl,ij))
03409 tdir(kfrsnl,ij) = Tf_dir_a*tdir(kfrsnl,ij) + &
03410 Tf_dir_a*rdir(kfrsnl,ij) * &
03411 Rf_dif_b*rintfc*tdif_a(kfrsnl,ij)
03412 rdir(kfrsnl,ij) = Rf_dir_a + &
03413 Tf_dir_a*rdir(kfrsnl,ij) * &
03414 rintfc*Tf_dif_b
03415 rdif_a(kfrsnl,ij) = Rf_dif_a + &
03416 Tf_dif_a*rdif_a(kfrsnl,ij) * &
03417 rintfc*Tf_dif_b
03418 rdif_b(kfrsnl,ij) = rdif_b(kfrsnl,ij) + &
03419 tdif_b(kfrsnl,ij)*Rf_dif_b * &
03420 rintfc*tdif_a(kfrsnl,ij)
03421 tdif_a(kfrsnl,ij) = Tf_dif_a*rintfc*tdif_a(kfrsnl,ij)
03422 tdif_b(kfrsnl,ij) = tdif_b(kfrsnl,ij)*rintfc*Tf_dif_b
03423
03424
03425 trnlay(kfrsnl,ij) = Tf_dir_a*trnlay(kfrsnl,ij)
03426 endif
03427
03428 endif
03429
03430 enddo
03431
03432
03433
03434
03435
03436
03437
03438
03439
03440
03441
03442
03443
03444
03445
03446 k = klevp
03447 trndir(k,ij) = trndir(k-1,ij)*trnlay(k-1,ij)
03448 refkm1 = c1/(c1 - rdndif(k-1,ij)*rdif_a(k-1,ij))
03449 tdrrdir = trndir(k-1,ij)*rdir(k-1,ij)
03450 tdndif = trntdr(k-1,ij) - trndir(k-1,ij)
03451 trntdr(k,ij) = trndir(k-1,ij)*tdir(k-1,ij) + &
03452 (tdndif + tdrrdir*rdndif(k-1,ij))*refkm1*tdif_a(k-1,ij)
03453 rdndif(k,ij) = rdif_b(k-1,ij) + &
03454 (tdif_b(k-1,ij)*rdndif(k-1,ij)*refkm1*tdif_a(k-1,ij))
03455 trndif(k,ij) = trndif(k-1,ij)*refkm1*tdif_a(k-1,ij)
03456
03457
03458
03459
03460
03461
03462
03463
03464
03465
03466
03467
03468
03469 rupdir(klevp,ij) = albodr(ij)
03470 rupdif(klevp,ij) = albodf(ij)
03471 enddo
03472 do k=klev,0,-1
03473 do ij = 1, icells_DE
03474 i = indxi_DE(ij)
03475 j = indxj_DE(ij)
03476
03477 refkp1 = c1/( c1 - rdif_b(k,ij)*rupdif(k+1,ij))
03478
03479
03480
03481 rupdir(k,ij) = rdir(k,ij) + &
03482 ( trnlay(k,ij)*rupdir(k+1,ij) + &
03483 (tdir(k,ij)-trnlay(k,ij))*rupdif(k+1,ij) ) * &
03484 refkp1*tdif_b(k,ij)
03485
03486
03487 rupdif(k,ij) = rdif_a(k,ij) + &
03488 tdif_a(k,ij)*rupdif(k+1,ij)* &
03489 refkp1*tdif_b(k,ij)
03490 enddo
03491 enddo
03492
03493 end subroutine solution_dEdd
03494
03495
03496
03497
03498
03499
03500
03501
03502
03503 subroutine shortwave_dEdd_set_snow(nx_block, ny_block, &
03504 icells, &
03505 indxi, indxj, &
03506 aice, vsno, &
03507 Tsfc, fs, &
03508 rhosnw, rsnw)
03509
03510
03511
03512
03513
03514
03515
03516
03517
03518
03519
03520
03521
03522 integer (kind=int_kind),
03523 intent(in) ::
03524 nx_block, ny_block,
03525 icells
03526
03527 integer (kind=int_kind), dimension (nx_block*ny_block),
03528 intent(in) ::
03529 indxi ,
03530 indxj
03531
03532 real (kind=dbl_kind), dimension (nx_block,ny_block),
03533 intent(in) ::
03534 aice ,
03535 vsno ,
03536 Tsfc
03537
03538 real (kind=dbl_kind), dimension (nx_block,ny_block),
03539 intent(out) ::
03540 fs
03541
03542 real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr),
03543 intent(out) ::
03544 rhosnw ,
03545 rsnw
03546
03547
03548
03549
03550
03551
03552 integer (kind=int_kind) ::
03553 i ,
03554 j ,
03555 ij ,
03556 ks
03557
03558 real (kind=dbl_kind) ::
03559 hs ,
03560 fT ,
03561 dTs ,
03562 rsnw_nm
03563
03564 real (kind=dbl_kind), parameter ::
03565
03566
03567
03568
03569 rsnw_fresh = 100._dbl_kind, &
03570 rsnw_nonmelt = 500._dbl_kind, &
03571 rsnw_sig = 250._dbl_kind
03572
03573 real (kind=dbl_kind) ::
03574 dT_mlt ,
03575
03576 rsnw_melt
03577
03578
03579
03580 dT_mlt = dT_mlt_in
03581 rsnw_melt = rsnw_melt_in
03582
03583 fs(:,:) = c0
03584 do ks = 1, nslyr
03585 do j = 1, ny_block
03586 do i = 1, nx_block
03587 rhosnw(i,j,ks) = c0
03588 rsnw(i,j,ks) = c0
03589 enddo
03590 enddo
03591 enddo
03592
03593
03594
03595
03596 do ij = 1, icells
03597 i = indxi(ij)
03598 j = indxj(ij)
03599
03600 if( aice(i,j) > puny ) then
03601 hs = vsno(i,j) / aice(i,j)
03602 if( hs < hsmin ) then
03603 fs(i,j) = c0
03604 else if( hs <= hs0 ) then
03605 fs(i,j) = hs/hs0
03606 else
03607 fs(i,j) = c1
03608 endif
03609
03610 dTs = Timelt - Tsfc(i,j)
03611 fT = -min(dTs/dT_mlt-c1,c0)
03612
03613
03614
03615 rsnw_nm = rsnw_nonmelt - R_snw*rsnw_sig
03616 rsnw_nm = max(rsnw_nm, rsnw_fresh)
03617 rsnw_nm = min(rsnw_nm, rsnw_melt)
03618 do ks = 1, nslyr
03619
03620 rhosnw(i,j,ks) = rhos
03621
03622 rsnw(i,j,ks) = rsnw_nm + (rsnw_melt-rsnw_nm)*fT
03623 rsnw(i,j,ks) = max(rsnw(i,j,ks), rsnw_fresh)
03624 rsnw(i,j,ks) = min(rsnw(i,j,ks), rsnw_melt)
03625 enddo
03626 endif
03627 enddo
03628
03629 end subroutine shortwave_dEdd_set_snow
03630
03631
03632
03633
03634
03635
03636
03637
03638
03639 subroutine shortwave_dEdd_set_pond(nx_block, ny_block, &
03640 icells, &
03641 indxi, indxj, &
03642 aice, Tsfc, &
03643 fs, fp, &
03644 hp)
03645
03646
03647
03648
03649
03650
03651
03652
03653
03654
03655
03656
03657
03658 integer (kind=int_kind),
03659 intent(in) ::
03660 nx_block, ny_block,
03661 icells
03662
03663 integer (kind=int_kind), dimension (nx_block*ny_block),
03664 intent(in) ::
03665 indxi ,
03666 indxj
03667
03668 real (kind=dbl_kind), dimension (nx_block,ny_block),
03669 intent(in) ::
03670 aice ,
03671 Tsfc ,
03672 fs
03673
03674 real (kind=dbl_kind), dimension (nx_block,ny_block),
03675 intent(out) ::
03676 fp ,
03677 hp
03678
03679
03680
03681
03682
03683
03684
03685 integer (kind=int_kind) ::
03686 i ,
03687 j ,
03688 ij
03689
03690 real (kind=dbl_kind) ::
03691 fT ,
03692 dTs
03693
03694 real (kind=dbl_kind), parameter ::
03695 dT_mlt = c1
03696
03697
03698
03699 do j = 1, ny_block
03700 do i = 1, nx_block
03701 fp(i,j) = c0
03702 hp(i,j) = c0
03703 enddo
03704 enddo
03705
03706
03707
03708
03709
03710 do ij = 1, icells
03711 i = indxi(ij)
03712 j = indxj(ij)
03713 if (aice(i,j) > puny) then
03714
03715 dTs = Timelt - Tsfc(i,j)
03716 fT = -min(dTs/dT_mlt-c1,c0)
03717
03718 fp(i,j) = 0.3_dbl_kind*fT*(c1-fs(i,j))
03719 hp(i,j) = 0.3_dbl_kind*fT*(c1-fs(i,j))
03720 endif
03721 enddo
03722
03723 end subroutine shortwave_dEdd_set_pond
03724
03725
03726
03727
03728 end module ice_shortwave
03729
03730