Changeset 10425 for NEMO/trunk/src/ICE/icethd_zdf_bl99.F90
- Timestamp:
- 2018-12-19T22:54:16+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icethd_zdf_bl99.F90
r10069 r10425 83 83 INTEGER, DIMENSION(jpij) :: jm_min ! reference number of top equation 84 84 INTEGER, DIMENSION(jpij) :: jm_max ! reference number of bottom equation 85 ! 85 86 LOGICAL, DIMENSION(jpij) :: l_T_converged ! true when T converges (per grid point) 87 ! 86 88 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system 87 89 REAL(wp) :: zg1 = 2._wp ! … … 113 115 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow to check the convergence 114 116 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity 117 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i_cp ! copy 115 118 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice 116 119 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice … … 201 204 ! 202 205 iconv = 0 ! number of iterations 203 zdti_max = 1000._wp ! maximal value of error on all points204 !206 ! 207 l_T_converged(:) = .FALSE. 205 208 ! !============================! 206 DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max ) ! Iterative procedure begins ! 209 ! Convergence calculated until all sub-domain grid points have converged 210 ! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation) 211 ! but values are not taken into account (results independant of MPI partitioning) 212 ! 213 DO WHILE ( ( .NOT. ALL (l_T_converged(1:npti)) ) .AND. iconv < iconv_max ) ! Iterative procedure begins ! 207 214 ! !============================! 208 215 iconv = iconv + 1 … … 217 224 ! 218 225 DO ji = 1, npti 219 ztcond_i (ji,0) = rcnd_i + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )220 ztcond_i (ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )226 ztcond_i_cp(ji,0) = rcnd_i + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 227 ztcond_i_cp(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) 221 228 END DO 222 229 DO jk = 1, nlay_i-1 223 230 DO ji = 1, npti 224 ztcond_i (ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &225 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )231 ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 232 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 226 233 END DO 227 234 END DO … … 230 237 ! 231 238 DO ji = 1, npti 232 ztcond_i (ji,0) = rcnd_i + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) &233 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 )234 ztcond_i (ji,nlay_i) = rcnd_i + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) &235 & - 0.011_wp * ( t_bo_1d(ji) - rt0 )239 ztcond_i_cp(ji,0) = rcnd_i + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 240 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 241 ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 242 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 236 243 END DO 237 244 DO jk = 1, nlay_i-1 238 245 DO ji = 1, npti 239 ztcond_i (ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &240 & MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) &241 & - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )246 ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 247 & MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) & 248 & - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 242 249 END DO 243 250 END DO 244 251 ! 245 252 ENDIF 246 ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) ) 253 254 ! Variable used after iterations 255 ! Value must be frozen after convergence for MPP independance reason 256 DO ji = 1, npti 257 IF ( .NOT. l_T_converged(ji) ) & 258 ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) ) 259 END DO 247 260 ! 248 261 !--- G(he) : enhancement of thermal conductivity in mono-category case … … 270 283 !----------------- 271 284 !--- Snow 285 ! Variable used after iterations 286 ! Value must be frozen after convergence for MPP independance reason 272 287 DO jk = 0, nlay_s-1 273 288 DO ji = 1, npti 274 zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 289 IF ( .NOT. l_T_converged(ji) ) & 290 zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 275 291 END DO 276 292 END DO 277 293 DO ji = 1, npti ! Snow-ice interface 278 zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 279 IF( zfac > epsi10 ) THEN 280 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 281 ELSE 282 zkappa_s(ji,nlay_s) = 0._wp 294 IF ( .NOT. l_T_converged(ji) ) THEN 295 zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 296 IF( zfac > epsi10 ) THEN 297 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 298 ELSE 299 zkappa_s(ji,nlay_s) = 0._wp 300 ENDIF 283 301 ENDIF 284 302 END DO 285 303 286 304 !--- Ice 305 ! Variable used after iterations 306 ! Value must be frozen after convergence for MPP independance reason 287 307 DO jk = 0, nlay_i 288 308 DO ji = 1, npti 289 zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 309 IF ( .NOT. l_T_converged(ji) ) & 310 zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 290 311 END DO 291 312 END DO 292 313 DO ji = 1, npti ! Snow-ice interface 293 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 314 IF ( .NOT. l_T_converged(ji) ) & 315 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 294 316 END DO 295 317 ! … … 326 348 ! update of the non solar flux according to the update in T_su 327 349 DO ji = 1, npti 328 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 350 ! Variable used after iterations 351 ! Value must be frozen after convergence for MPP independance reason 352 IF ( .NOT. l_T_converged(ji) ) & 353 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 329 354 END DO 330 355 … … 496 521 ! ice temperatures 497 522 DO ji = 1, npti 498 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 523 ! Variable used after iterations 524 ! Value must be frozen after convergence for MPP independance reason 525 IF ( .NOT. l_T_converged(ji) ) & 526 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 499 527 END DO 500 528 … … 502 530 DO ji = 1, npti 503 531 jk = jm - nlay_s - 1 504 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 532 IF ( .NOT. l_T_converged(ji) ) & 533 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 505 534 END DO 506 535 END DO 507 536 508 537 DO ji = 1, npti 509 ! snow temperatures 510 IF( h_s_1d(ji) > 0._wp ) THEN 511 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 512 ENDIF 513 ! surface temperature 514 ztsub(ji) = t_su_1d(ji) 515 IF( t_su_1d(ji) < rt0 ) THEN 516 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 517 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 538 ! Variables used after iterations 539 ! Value must be frozen after convergence for MPP independance reason 540 IF ( .NOT. l_T_converged(ji) ) THEN 541 ! snow temperatures 542 IF( h_s_1d(ji) > 0._wp ) THEN 543 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 544 ENDIF 545 ! surface temperature 546 ztsub(ji) = t_su_1d(ji) 547 IF( t_su_1d(ji) < rt0 ) THEN 548 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 549 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 550 ENDIF 518 551 ENDIF 519 552 END DO … … 524 557 ! check that nowhere it has started to melt 525 558 ! zdti_max is a measure of error, it has to be under zdti_bnd 526 zdti_max = 0._wp 527 DO ji = 1, npti 528 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 529 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 530 END DO 531 532 DO jk = 1, nlay_s 533 DO ji = 1, npti 534 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 535 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 536 END DO 537 END DO 538 539 DO jk = 1, nlay_i 540 DO ji = 1, npti 541 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 542 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 543 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 544 END DO 545 END DO 546 547 ! Compute spatial maximum over all errors 548 ! note that this could be optimized substantially by iterating only the non-converging points 549 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 550 ! 559 560 DO ji = 1, npti 561 562 zdti_max = 0._wp 563 564 IF ( .NOT. l_T_converged(ji) ) THEN 565 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 566 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 567 568 t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 569 zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 570 571 DO jk = 1, nlay_i 572 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 573 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 574 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 575 END DO 576 577 IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE. 578 579 ENDIF 580 581 END DO 582 551 583 !----------------------------------------! 552 584 ! ! … … 670 702 671 703 ! ice temperatures 672 DO ji = 1, npti 673 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 704 DO ji = 1, npti 705 ! Variable used after iterations 706 ! Value must be frozen after convergence for MPP independance reason 707 IF ( .NOT. l_T_converged(ji) ) & 708 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 674 709 END DO 675 710 676 711 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 677 712 DO ji = 1, npti 678 jk = jm - nlay_s - 1 679 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 713 IF ( .NOT. l_T_converged(ji) ) THEN 714 jk = jm - nlay_s - 1 715 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 716 ENDIF 680 717 END DO 681 718 END DO … … 683 720 ! snow temperatures 684 721 DO ji = 1, npti 685 IF( h_s_1d(ji) > 0._wp ) THEN 686 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 722 ! Variable used after iterations 723 ! Value must be frozen after convergence for MPP independance reason 724 IF ( .NOT. l_T_converged(ji) ) THEN 725 IF( h_s_1d(ji) > 0._wp ) THEN 726 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 727 ENDIF 687 728 ENDIF 688 729 END DO … … 693 734 ! check that nowhere it has started to melt 694 735 ! zdti_max is a measure of error, it has to be under zdti_bnd 695 zdti_max = 0._wp 696 697 DO jk = 1, nlay_s 698 DO ji = 1, npti 699 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 700 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 701 END DO 702 END DO 703 704 DO jk = 1, nlay_i 705 DO ji = 1, npti 706 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 707 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 708 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 709 END DO 710 END DO 711 712 ! Compute spatial maximum over all errors 713 ! note that this could be optimized substantially by iterating only the non-converging points 714 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 736 737 DO ji = 1, npti 738 739 zdti_max = 0._wp 740 741 IF ( .NOT. l_T_converged(ji) ) THEN 742 ! t_s 743 t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 744 zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 745 ! t_i 746 DO jk = 0, nlay_i 747 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 748 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 749 zdti_max = MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 750 END DO 751 752 IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE. 753 754 ENDIF 755 756 END DO 715 757 716 758 ENDIF ! k_jules
Note: See TracChangeset
for help on using the changeset viewer.