- Timestamp:
- 2017-07-15T17:27:14+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r8325 r8342 37 37 CONTAINS 38 38 39 SUBROUTINE lim_thd_dif ( kideb , kiut )39 SUBROUTINE lim_thd_dif 40 40 !!------------------------------------------------------------------ 41 41 !! *** ROUTINE lim_thd_dif *** … … 67 67 !! of temperature 68 68 !! 69 !! ** Arguments :70 !! kideb , kiut : Starting and ending points on which the71 !! the computation is applied72 69 !! 73 70 !! ** Inputs / Ouputs : (global commons) … … 89 86 !! (04-2007) Energy conservation tested by M. Vancoppenolle 90 87 !!------------------------------------------------------------------ 91 INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied92 93 88 !! * Local variables 94 89 INTEGER :: ji ! spatial loop index … … 180 175 ! --- diag error on heat diffusion - PART 1 --- ! 181 176 zdq(:) = 0._wp ; zq_ini(:) = 0._wp 182 DO ji = kideb, kiut177 DO ji = 1, nidx 183 178 zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & 184 179 & SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) … … 188 183 ! 1) Initialization ! 189 184 !------------------------------------------------------------------------------! 190 DO ji = kideb , kiut185 DO ji = 1 , nidx 191 186 isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) ) ! is there snow or not 192 187 ! layer thickness … … 203 198 204 199 DO jk = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 205 DO ji = kideb , kiut200 DO ji = 1 , nidx 206 201 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s 207 202 END DO … … 209 204 210 205 DO jk = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 211 DO ji = kideb , kiut206 DO ji = 1 , nidx 212 207 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i 213 208 END DO … … 230 225 ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 231 226 zhsu = 0.1_wp ! threshold for the computation of i0 232 DO ji = kideb , kiut227 DO ji = 1 , nidx 233 228 ! switches 234 229 isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) … … 243 238 ! Derivative of the non solar flux 244 239 !------------------------------------------------------- 245 DO ji = kideb , kiut240 DO ji = 1 , nidx 246 241 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface 247 242 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer … … 254 249 !--------------------------------------------------------- 255 250 256 DO ji = kideb, kiut! snow initialization251 DO ji = 1, nidx ! snow initialization 257 252 zradtr_s(ji,0) = zftrice(ji) ! radiation penetrating through snow 258 253 END DO 259 254 260 255 DO jk = 1, nlay_s ! Radiation through snow 261 DO ji = kideb, kiut256 DO ji = 1, nidx 262 257 ! ! radiation transmitted below the layer-th snow layer 263 258 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) ) … … 267 262 END DO 268 263 269 DO ji = kideb, kiut! ice initialization264 DO ji = 1, nidx ! ice initialization 270 265 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) ) 271 266 END DO 272 267 273 268 DO jk = 1, nlay_i ! Radiation through ice 274 DO ji = kideb, kiut269 DO ji = 1, nidx 275 270 ! ! radiation transmitted below the layer-th ice layer 276 271 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) ) … … 280 275 END DO 281 276 282 DO ji = kideb, kiut! Radiation transmitted below the ice277 DO ji = 1, nidx ! Radiation transmitted below the ice 283 278 ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 284 279 END DO … … 288 283 !------------------------------------------------------------------------------| 289 284 ! 290 DO ji = kideb, kiut! Old surface temperature285 DO ji = 1, nidx ! Old surface temperature 291 286 ztsub (ji) = t_su_1d(ji) ! temperature at the beg of iter pr. 292 287 ztsubit(ji) = t_su_1d(ji) ! temperature at the previous iter … … 296 291 297 292 DO jk = 1, nlay_s ! Old snow temperature 298 DO ji = kideb , kiut293 DO ji = 1 , nidx 299 294 ztsb(ji,jk) = t_s_1d(ji,jk) 300 295 END DO … … 302 297 303 298 DO jk = 1, nlay_i ! Old ice temperature 304 DO ji = kideb , kiut299 DO ji = 1 , nidx 305 300 ztib(ji,jk) = t_i_1d(ji,jk) 306 301 END DO … … 319 314 ! 320 315 IF( nn_ice_thcon == 0 ) THEN ! Untersteiner (1964) formula 321 DO ji = kideb , kiut316 DO ji = 1 , nidx 322 317 ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 323 318 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 324 319 END DO 325 320 DO jk = 1, nlay_i-1 326 DO ji = kideb , kiut321 DO ji = 1 , nidx 327 322 ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / & 328 323 MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) … … 333 328 334 329 IF( nn_ice_thcon == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 335 DO ji = kideb , kiut330 DO ji = 1 , nidx 336 331 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 337 332 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) … … 339 334 END DO 340 335 DO jk = 1, nlay_i-1 341 DO ji = kideb , kiut336 DO ji = 1 , nidx 342 337 ztcond_i(ji,jk) = rcdic + & 343 338 & 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) & … … 347 342 END DO 348 343 END DO 349 DO ji = kideb , kiut344 DO ji = 1 , nidx 350 345 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 351 346 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) … … 372 367 zh_thres = EXP( 1._wp ) * zepsilon * 0.5_wp 373 368 374 DO ji = kideb, kiut369 DO ji = 1, nidx 375 370 376 371 ! Mean sea ice thermal conductivity … … 400 395 ! 401 396 !--- Snow 402 DO ji = kideb, kiut397 DO ji = 1, nidx 403 398 zfac = 1. / MAX( epsi10 , zh_s(ji) ) 404 399 zkappa_s(ji,0) = zghe(ji) * rn_cdsn * zfac … … 407 402 408 403 DO jk = 1, nlay_s-1 409 DO ji = kideb , kiut404 DO ji = 1 , nidx 410 405 zkappa_s(ji,jk) = zghe(ji) * 2.0 * rn_cdsn / MAX( epsi10, 2.0 * zh_s(ji) ) 411 406 END DO … … 414 409 !--- Ice 415 410 DO jk = 1, nlay_i-1 416 DO ji = kideb , kiut411 DO ji = 1 , nidx 417 412 zkappa_i(ji,jk) = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0 * zh_i(ji) ) 418 413 END DO … … 420 415 421 416 !--- Snow-ice interface 422 DO ji = kideb , kiut417 DO ji = 1 , nidx 423 418 zfac = 1./ MAX( epsi10 , zh_i(ji) ) 424 419 zkappa_i(ji,0) = zghe(ji) * ztcond_i(ji,0) * zfac … … 435 430 ! 436 431 DO jk = 1, nlay_i 437 DO ji = kideb , kiut432 DO ji = 1 , nidx 438 433 ztitemp(ji,jk) = t_i_1d(ji,jk) 439 434 zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 ) … … 443 438 444 439 DO jk = 1, nlay_s 445 DO ji = kideb , kiut440 DO ji = 1 , nidx 446 441 ztstemp(ji,jk) = t_s_1d(ji,jk) 447 442 zeta_s(ji,jk) = rdt_ice / MAX( rhosn * cpic * zh_s(ji), epsi10 ) … … 455 450 ! 456 451 IF ( ln_dqnsice ) THEN 457 DO ji = kideb , kiut452 DO ji = 1 , nidx 458 453 ! update of the non solar flux according to the update in T_su 459 454 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) ) … … 462 457 463 458 ! Update incoming flux 464 DO ji = kideb , kiut459 DO ji = 1 , nidx 465 460 ! update incoming flux 466 461 zf(ji) = zfsw(ji) & ! net absorbed solar radiation … … 481 476 482 477 DO numeq=1,nlay_i+3 483 DO ji = kideb , kiut478 DO ji = 1 , nidx 484 479 ztrid(ji,numeq,1) = 0. 485 480 ztrid(ji,numeq,2) = 0. … … 492 487 493 488 DO numeq = nlay_s + 2, nlay_s + nlay_i 494 DO ji = kideb , kiut489 DO ji = 1 , nidx 495 490 jk = numeq - nlay_s - 1 496 491 ztrid(ji,numeq,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) … … 502 497 503 498 numeq = nlay_s + nlay_i + 1 504 DO ji = kideb , kiut499 DO ji = 1 , nidx 505 500 !!ice bottom term 506 501 ztrid(ji,numeq,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) … … 512 507 513 508 514 DO ji = kideb , kiut509 DO ji = 1 , nidx 515 510 IF ( ht_s_1d(ji) > 0.0 ) THEN 516 511 ! … … 659 654 minnumeqmin = nlay_i+5 660 655 661 DO ji = kideb , kiut656 DO ji = 1 , nidx 662 657 zindtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji)) 663 658 zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2) … … 667 662 668 663 DO jk = minnumeqmin+1, maxnumeqmax 669 DO ji = kideb , kiut664 DO ji = 1 , nidx 670 665 numeq = min(max(numeqmin(ji)+1,jk),numeqmax(ji)) 671 666 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3) / zdiagbis(ji,numeq-1) … … 674 669 END DO 675 670 676 DO ji = kideb , kiut671 DO ji = 1 , nidx 677 672 ! ice temperatures 678 673 t_i_1d(ji,nlay_i) = zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji)) … … 680 675 681 676 DO numeq = nlay_i + nlay_s, nlay_s + 2, -1 682 DO ji = kideb , kiut677 DO ji = 1 , nidx 683 678 jk = numeq - nlay_s - 1 684 679 t_i_1d(ji,jk) = ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq) … … 686 681 END DO 687 682 688 DO ji = kideb , kiut683 DO ji = 1 , nidx 689 684 ! snow temperatures 690 685 IF (ht_s_1d(ji) > 0._wp) & … … 706 701 ! check that nowhere it has started to melt 707 702 ! zdti(ji) is a measure of error, it has to be under zdti_bnd 708 DO ji = kideb , kiut703 DO ji = 1 , nidx 709 704 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , 190._wp ) 710 705 zdti (ji) = ABS( t_su_1d(ji) - ztsubit(ji) ) … … 712 707 713 708 DO jk = 1, nlay_s 714 DO ji = kideb , kiut709 DO ji = 1 , nidx 715 710 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), 190._wp ) 716 711 zdti (ji) = MAX( zdti(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) ) … … 719 714 720 715 DO jk = 1, nlay_i 721 DO ji = kideb , kiut716 DO ji = 1 , nidx 722 717 ztmelt_i = -tmut * s_i_1d(ji,jk) + rt0 723 718 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp ) … … 729 724 ! note that this could be optimized substantially by iterating only the non-converging points 730 725 zdti_max = 0._wp 731 DO ji = kideb, kiut726 DO ji = 1, nidx 732 727 zdti_max = MAX( zdti_max, zdti(ji) ) 733 728 END DO … … 738 733 ! MV SIMIP 2016 739 734 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 740 DO ji = kideb, kiut735 DO ji = 1, nidx 741 736 zfac = 1. / MAX( epsi10 , rn_cdsn * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) ) 742 737 t_si_1d(ji) = ( rn_cdsn * zh_i(ji) * t_s_1d(ji,1) + & … … 755 750 ! 12) Fluxes at the interfaces ! 756 751 !-------------------------------------------------------------------------! 757 DO ji = kideb, kiut752 DO ji = 1, nidx 758 753 ! ! surface ice conduction flux 759 754 isnow(ji) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) … … 771 766 772 767 ! --- computes sea ice energy of melting compulsory for limthd_dh --- ! 773 CALL lim_thd_enmelt ( kideb, kiut )768 CALL lim_thd_enmelt 774 769 775 770 ! --- diagnose the change in non-solar flux due to surface temperature change --- ! 776 771 IF ( ln_dqnsice ) THEN 777 DO ji = kideb, kiut772 DO ji = 1, nidx 778 773 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 779 774 END DO … … 781 776 782 777 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 783 DO ji = kideb, kiut778 DO ji = 1, nidx 784 779 zdq(ji) = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & 785 780 & SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) … … 798 793 ! Heat flux used to warm/cool ice in W.m-2 799 794 !----------------------------------------- 800 DO ji = kideb, kiut795 DO ji = 1, nidx 801 796 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 802 797 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & … … 821 816 END SUBROUTINE lim_thd_dif 822 817 823 SUBROUTINE lim_thd_enmelt ( kideb, kiut )818 SUBROUTINE lim_thd_enmelt 824 819 !!----------------------------------------------------------------------- 825 820 !! *** ROUTINE lim_thd_enmelt *** … … 829 824 !! ** Method : Formula (Bitz and Lipscomb, 1999) 830 825 !!------------------------------------------------------------------- 831 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop832 !833 826 INTEGER :: ji, jk ! dummy loop indices 834 827 REAL(wp) :: ztmelts ! local scalar … … 836 829 ! 837 830 DO jk = 1, nlay_i ! Sea ice energy of melting 838 DO ji = kideb, kiut831 DO ji = 1, nidx 839 832 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 840 833 t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point … … 846 839 END DO 847 840 DO jk = 1, nlay_s ! Snow energy of melting 848 DO ji = kideb, kiut841 DO ji = 1, nidx 849 842 e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) 850 843 END DO
Note: See TracChangeset
for help on using the changeset viewer.
