- Timestamp:
- 2017-11-24T17:56:51+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
r8637 r8813 15 15 !! 'key_lim3' ESIM sea-ice model 16 16 !!---------------------------------------------------------------------- 17 !! ice_thd_zdf : select the appropriate routine for vertical heat diffusion calculation 18 !! ice_thd_zdf_BL99 : 19 !! ice_thd_enmelt : 20 !! ice_thd_zdf_init : 21 !!---------------------------------------------------------------------- 17 22 USE dom_oce ! ocean space and time domain 18 23 USE phycst ! physical constants (ocean directory) … … 34 39 LOGICAL :: ln_cndi_U64 ! thermal conductivity: Untersteiner (1964) 35 40 LOGICAL :: ln_cndi_P07 ! thermal conductivity: Pringle et al (2007) 36 REAL(wp) :: rn_cnd_s ! thermal conductivity of the snow [W/m/K]37 41 REAL(wp) :: rn_kappa_i ! coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 38 LOGICAL :: ln_dqns_i ! change non-solar surface flux with changing surface temperature (T) or not (F) 39 42 REAL(wp), PUBLIC :: rn_cnd_s ! thermal conductivity of the snow [W/m/K] 43 44 INTEGER :: nice_zdf ! Choice of the type of vertical heat diffusion formulation 45 ! ! associated indices: 46 INTEGER, PARAMETER :: np_BL99 = 1 ! Bitz and Lipscomb (1999) 47 48 INTEGER, PARAMETER :: np_zdf_jules_OFF = 0 ! compute all temperatures from qsr and qns 49 INTEGER, PARAMETER :: np_zdf_jules_SND = 1 ! compute conductive heat flux and surface temperature from qsr and qns 50 INTEGER, PARAMETER :: np_zdf_jules_RCV = 2 ! compute snow and inner ice temperatures from qcnd 51 40 52 !!---------------------------------------------------------------------- 41 53 !! NEMO/ICE 4.0 , NEMO Consortium (2017) … … 48 60 !!------------------------------------------------------------------- 49 61 !! *** ROUTINE ice_thd_zdf *** 50 !! ** Purpose : 51 !! This routine determines the time evolution of snow and sea-ice 52 !! temperature profiles. 53 !! ** Method : 54 !! This is done by solving the heat equation diffusion with 55 !! a Neumann boundary condition at the surface and a Dirichlet one 56 !! at the bottom. Solar radiation is partially absorbed into the ice. 57 !! The specific heat and thermal conductivities depend on ice salinity 58 !! and temperature to take into account brine pocket melting. The 59 !! numerical 60 !! scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid 61 !! in the ice and snow system. 62 !! 63 !! ** Purpose : select the appropriate routine for the computation 64 !! of vertical diffusion 65 !!------------------------------------------------------------------- 66 67 SELECT CASE ( nice_zdf ) ! Choose the vertical heat diffusion solver 68 ! 69 ! !------------- 70 CASE( np_BL99 ) ! BL99 solver 71 ! !------------- 72 SELECT CASE( nice_jules ) 73 CASE( np_jules_OFF ) ; CALL ice_thd_zdf_BL99 ( np_zdf_jules_OFF ) ! No Jules coupler 74 CASE( np_jules_EMULE ) ; CALL ice_thd_zdf_BL99 ( np_zdf_jules_SND ) ! Jules coupler is emulated (send/recieve) 75 CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) 76 CASE( np_jules_ACTIVE ) ; CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) ! Jules coupler is active (receive only) 77 END SELECT 78 END SELECT 79 80 END SUBROUTINE ice_thd_zdf 81 82 83 SUBROUTINE ice_thd_zdf_BL99( k_jules ) 84 !!------------------------------------------------------------------- 85 !! *** ROUTINE ice_thd_zdf_BL99 *** 86 !! 87 !! ** Purpose : computes the time evolution of snow and sea-ice temperature 88 !! profiles, using the original Bitz and Lipscomb (1999) algorithm 89 !! 90 !! ** Method : solves the heat equation diffusion with a Neumann boundary 91 !! condition at the surface and a Dirichlet one at the bottom. 92 !! Solar radiation is partially absorbed into the ice. 93 !! The specific heat and thermal conductivities depend on ice 94 !! salinity and temperature to take into account brine pocket 95 !! melting. The numerical scheme is an iterative Crank-Nicolson 96 !! on a non-uniform multilayer grid in the ice and snow system. 62 97 !! 63 98 !! The successive steps of this routine are … … 83 118 !! total ice/snow thickness : h_i_1d, h_s_1d 84 119 !!------------------------------------------------------------------- 120 INTEGER, INTENT(in) :: k_jules ! Jules coupling (0=OFF, 1=RECEIVE, 2=SEND) 121 ! 85 122 INTEGER :: ji, jk ! spatial loop index 86 123 INTEGER :: jm ! current reference number of equation … … 88 125 INTEGER :: iconv ! number of iterations in iterative procedure 89 126 INTEGER :: iconv_max = 50 ! max number of iterations in iterative procedure 90 127 ! 91 128 INTEGER, DIMENSION(jpij) :: jm_min ! reference number of top equation 92 129 INTEGER, DIMENSION(jpij) :: jm_max ! reference number of bottom equation 93 130 ! 94 131 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system 95 132 REAL(wp) :: zg1 = 2._wp ! … … 101 138 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature 102 139 REAL(wp) :: ztmelt_i ! ice melting temperature 103 REAL(wp) :: z1_hsu104 140 REAL(wp) :: zdti_max ! current maximal error on temperature 105 141 REAL(wp) :: zcpi ! Ice specific heat 106 142 REAL(wp) :: zhfx_err, zdq ! diag errors on heat 107 143 REAL(wp) :: zfac ! dummy factor 108 144 ! 109 145 REAL(wp), DIMENSION(jpij) :: isnow ! switch for presence (1) or absence (0) of snow 110 146 REAL(wp), DIMENSION(jpij) :: ztsub ! surface temperature at previous iteration 111 147 REAL(wp), DIMENSION(jpij) :: zh_i, z1_h_i ! ice layer thickness 112 148 REAL(wp), DIMENSION(jpij) :: zh_s, z1_h_s ! snow layer thickness 113 REAL(wp), DIMENSION(jpij) :: zfsw ! solar radiation absorbed at the surface114 149 REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface 115 150 REAL(wp), DIMENSION(jpij) :: zfnet ! surface flux function 116 151 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function 117 REAL(wp), DIMENSION(jpij) :: zftrice ! solar radiation transmitted through the ice118 152 ! 153 REAL(wp), DIMENSION(jpij ) :: ztsuold ! Old surface temperature in the ice 119 154 REAL(wp), DIMENSION(jpij,nlay_i) :: ztiold ! Old temperature in the ice 120 155 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsold ! Old temperature in the snow … … 136 171 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat 137 172 REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 138 173 ! 174 REAL(wp) :: zfr1, zfr2, zfrqsr_tr_i ! dummy factor 175 ! 139 176 ! Mono-category 140 177 REAL(wp) :: zepsilon ! determines thres. above which computation of G(h) is done … … 162 199 ELSEWHERE ; z1_h_i(1:npti) = 0._wp 163 200 END WHERE 164 201 ! 165 202 WHERE( zh_s(1:npti) >= epsi10 ) ; z1_h_s(1:npti) = 1._wp / zh_s(1:npti) 166 203 ELSEWHERE ; z1_h_s(1:npti) = 0._wp 167 204 END WHERE 168 205 ! 169 ! temperatures 170 ztsub (1:npti) = t_su_1d(1:npti) ! temperature at the previous iteration 206 ! Store initial temperatures and non solar heat fluxes 207 IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN ! OFF or SND mode 208 ! 209 ztsub (1:npti) = t_su_1d(1:npti) ! surface temperature at iteration n-1 210 ztsuold(1:npti) = t_su_1d(1:npti) ! surface temperature initial value 211 zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti) ! derivative of incoming nonsolar flux 212 zqns_ice_b (1:npti) = qns_ice_1d(1:npti) ! store previous qns_ice_1d value 213 ! 214 t_su_1d(1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err ) ! required to leave the choice between melting or not 215 ! 216 ENDIF 217 ! 171 218 ztsold (1:npti,:) = t_s_1d(1:npti,:) ! Old snow temperature 172 219 ztiold (1:npti,:) = t_i_1d(1:npti,:) ! Old ice temperature 173 t_su_1d(1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err ) ! necessary 174 ! 220 175 221 !------------- 176 222 ! 2) Radiation 177 223 !------------- 178 z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0179 DO ji = 1, npti180 ! --- Computation of i0 --- !181 ! i0 describes the fraction of solar radiation which does not contribute182 ! to the surface energy budget but rather penetrates inside the ice.183 ! We assume that no radiation is transmitted through the snow184 ! If there is no no snow185 ! zfsw = (1-i0).qsr_ice is absorbed at the surface186 ! zftrice = io.qsr_ice is below the surface187 ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice188 ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover189 zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * z1_hsu ) )190 i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zfac * fr2_i0_1d(ji) )191 192 ! --- Solar radiation absorbed / transmitted at the surface --- !193 ! Derivative of the non solar flux194 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface195 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer196 zdqns_ice_b(ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux197 zqns_ice_b (ji) = qns_ice_1d(ji) ! store previous qns_ice_1d value198 END DO199 200 224 ! --- Transmission/absorption of solar radiation in the ice --- ! 201 zradtr_s(1:npti,0) = zftrice(1:npti) 225 ! zfr1 = ( 0.18 * ( 1.0 - 0.81 ) + 0.35 * 0.81 ) ! standard value 226 ! zfr2 = ( 0.82 * ( 1.0 - 0.81 ) + 0.65 * 0.81 ) ! zfr2 such that zfr1 + zfr2 to equal 1 227 ! 228 ! DO ji = 1, npti 229 ! 230 ! zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * 10._wp ) ) 231 ! 232 ! zfrqsr_tr_i = zfr1 + zfac * zfr2 ! below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness 233 ! IF ( h_s_1d(ji) >= 0.0_wp ) zfrqsr_tr_i = 0._wp ! snow fully opaque 234 ! 235 ! qsr_ice_tr_1d(ji) = zfrqsr_tr_i * qsr_ice_1d(ji) ! transmitted solar radiation 236 ! 237 ! zfsw(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) 238 ! zftrice(ji) = qsr_ice_tr_1d(ji) 239 ! i0(ji) = zfrqsr_tr_i 240 ! 241 ! END DO 242 243 zradtr_s(1:npti,0) = qsr_ice_tr_1d(1:npti) 202 244 DO jk = 1, nlay_s 203 245 DO ji = 1, npti … … 208 250 END DO 209 251 END DO 210 211 zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + zftrice(1:npti) * ( 1._wp - isnow(1:npti) )252 ! 253 zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qsr_ice_tr_1d(1:npti) * ( 1._wp - isnow(1:npti) ) 212 254 DO jk = 1, nlay_i 213 255 DO ji = 1, npti … … 218 260 END DO 219 261 END DO 220 262 ! 221 263 ftr_ice_1d(1:npti) = zradtr_i(1:npti,nlay_i) ! record radiation transmitted below the ice 222 264 ! … … 273 315 ! 274 316 SELECT CASE ( nn_monocat ) 275 317 ! 276 318 CASE ( 1 , 3 ) 277 319 ! 278 320 zepsilon = 0.1_wp 279 321 DO ji = 1, npti … … 284 326 ENDIF 285 327 END DO 286 328 ! 287 329 END SELECT 288 330 ! … … 330 372 END DO 331 373 END DO 332 ! 333 !---------------------------- 334 ! 6) surface flux computation 335 !---------------------------- 336 IF ( ln_dqns_i ) THEN 337 DO ji = 1, npti 338 ! update of the non solar flux according to the update in T_su 374 375 ! 376 !----------------------------------------! 377 ! ! 378 ! JULES IF (OFF or SND MODE) ! 379 ! ! 380 !----------------------------------------! 381 ! 382 383 IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN ! OFF or SND mode 384 385 ! 386 ! In OFF mode the original BL99 temperature computation is used 387 ! (with qsr_ice, qns_ice and dqns_ice as inputs) 388 ! 389 ! In SND mode, the computation is required to compute the conduction fluxes 390 ! 391 392 !---------------------------- 393 ! 6) surface flux computation 394 !---------------------------- 395 396 DO ji = 1, npti 397 ! update of the non solar flux according to the update in T_su 339 398 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 340 399 END DO 341 ENDIF 342 343 DO ji = 1, npti 344 zfnet(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH) 345 END DO 346 ! 347 !---------------------------- 348 ! 7) tridiagonal system terms 349 !---------------------------- 350 !!layer denotes the number of the layer in the snow or in the ice 351 !!jm denotes the reference number of the equation in the tridiagonal 352 !!system, terms of tridiagonal system are indexed as following : 353 !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 354 355 !!ice interior terms (top equation has the same form as the others) 356 ztrid (1:npti,:,:) = 0._wp 357 zindterm(1:npti,:) = 0._wp 358 zindtbis(1:npti,:) = 0._wp 359 zdiagbis(1:npti,:) = 0._wp 360 361 DO jm = nlay_s + 2, nlay_s + nlay_i 400 401 DO ji = 1, npti 402 zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net solar - transmitted solar + non solar 403 END DO 404 ! 405 !---------------------------- 406 ! 7) tridiagonal system terms 407 !---------------------------- 408 !!layer denotes the number of the layer in the snow or in the ice 409 !!jm denotes the reference number of the equation in the tridiagonal 410 !!system, terms of tridiagonal system are indexed as following : 411 !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 412 413 !!ice interior terms (top equation has the same form as the others) 414 ztrid (1:npti,:,:) = 0._wp 415 zindterm(1:npti,:) = 0._wp 416 zindtbis(1:npti,:) = 0._wp 417 zdiagbis(1:npti,:) = 0._wp 418 419 DO jm = nlay_s + 2, nlay_s + nlay_i 362 420 DO ji = 1, npti 363 421 jk = jm - nlay_s - 1 … … 367 425 zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 368 426 END DO 369 ENDDO370 371 jm = nlay_s + nlay_i + 1372 DO ji = 1, npti427 END DO 428 429 jm = nlay_s + nlay_i + 1 430 DO ji = 1, npti 373 431 !!ice bottom term 374 432 ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) … … 377 435 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 378 436 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 379 ENDDO380 381 382 DO ji = 1, npti437 END DO 438 439 440 DO ji = 1, npti 383 441 ! !---------------------! 384 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells !442 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! 385 443 ! !---------------------! 386 444 ! snow interior terms (bottom equation has the same form as the others) … … 428 486 & ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 429 487 ENDIF 430 !!---------------------!488 ! !---------------------! 431 489 ELSE ! cells without snow ! 432 490 ! !---------------------! … … 453 511 ztrid(ji,jm_min(ji),1) = 0.0 454 512 ztrid(ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 455 ztrid(ji,jm_min(ji),3) = zkappa_i(ji,0) * 2.0513 ztrid(ji,jm_min(ji),3) = zkappa_i(ji,0) * 2.0 456 514 ztrid(ji,jm_min(ji)+1,1) = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 457 515 ztrid(ji,jm_min(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) … … 488 546 zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 489 547 ! 490 END DO491 !492 !------------------------------493 ! 8) tridiagonal system solving494 !------------------------------495 ! Solve the tridiagonal system with Gauss elimination method.496 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984497 jm_maxt = 0498 jm_mint = nlay_i+5499 DO ji = 1, npti500 jm_mint = MIN(jm_min(ji),jm_mint)501 jm_maxt = MAX(jm_max(ji),jm_maxt)502 END DO503 504 DO jk = jm_mint+1, jm_maxt548 END DO 549 ! 550 !------------------------------ 551 ! 8) tridiagonal system solving 552 !------------------------------ 553 ! Solve the tridiagonal system with Gauss elimination method. 554 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 555 jm_maxt = 0 556 jm_mint = nlay_i+5 557 DO ji = 1, npti 558 jm_mint = MIN(jm_min(ji),jm_mint) 559 jm_maxt = MAX(jm_max(ji),jm_maxt) 560 END DO 561 562 DO jk = jm_mint+1, jm_maxt 505 563 DO ji = 1, npti 506 564 jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) … … 508 566 zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 509 567 END DO 510 END DO511 512 DO ji = 1, npti568 END DO 569 570 DO ji = 1, npti 513 571 ! ice temperatures 514 572 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 515 END DO516 517 DO jm = nlay_i + nlay_s, nlay_s + 2, -1573 END DO 574 575 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 518 576 DO ji = 1, npti 519 577 jk = jm - nlay_s - 1 520 578 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 521 579 END DO 522 END DO523 524 DO ji = 1, npti580 END DO 581 582 DO ji = 1, npti 525 583 ! snow temperatures 526 584 IF( h_s_1d(ji) > 0._wp ) THEN 527 585 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) & 528 586 & / zdiagbis(ji,nlay_s+1) 529 587 ENDIF 530 588 ! surface temperature … … 534 592 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 535 593 ENDIF 536 END DO537 !538 !--------------------------------------------------------------539 ! 9) Has the scheme converged ?, end of the iterative procedure540 !--------------------------------------------------------------541 ! check that nowhere it has started to melt542 ! zdti_max is a measure of error, it has to be under zdti_bnd543 zdti_max = 0._wp544 DO ji = 1, npti594 END DO 595 ! 596 !-------------------------------------------------------------- 597 ! 9) Has the scheme converged ?, end of the iterative procedure 598 !-------------------------------------------------------------- 599 ! check that nowhere it has started to melt 600 ! zdti_max is a measure of error, it has to be under zdti_bnd 601 zdti_max = 0._wp 602 DO ji = 1, npti 545 603 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 546 604 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 547 END DO548 549 DO jk = 1, nlay_s605 END DO 606 607 DO jk = 1, nlay_s 550 608 DO ji = 1, npti 551 609 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 552 610 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 553 611 END DO 554 END DO555 556 DO jk = 1, nlay_i612 END DO 613 614 DO jk = 1, nlay_i 557 615 DO ji = 1, npti 558 616 ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0 … … 560 618 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 561 619 END DO 562 END DO 563 564 ! Compute spatial maximum over all errors 565 ! note that this could be optimized substantially by iterating only the non-converging points 566 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 567 620 END DO 621 622 ! Compute spatial maximum over all errors 623 ! note that this could be optimized substantially by iterating only the non-converging points 624 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 625 ! 626 !----------------------------------------! 627 ! ! 628 ! JULES IF (RCV MODE) ! 629 ! ! 630 !----------------------------------------! 631 ! 632 633 ELSE IF ( k_jules == np_zdf_jules_RCV ) THEN ! RCV mode 634 635 ! 636 ! In RCV mode, we use a modified BL99 solver 637 ! with conduction flux (qcn_ice) as forcing term 638 ! 639 !---------------------------- 640 ! 7) tridiagonal system terms 641 !---------------------------- 642 !!layer denotes the number of the layer in the snow or in the ice 643 !!jm denotes the reference number of the equation in the tridiagonal 644 !!system, terms of tridiagonal system are indexed as following : 645 !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 646 647 !!ice interior terms (top equation has the same form as the others) 648 ztrid (1:npti,:,:) = 0._wp 649 zindterm(1:npti,:) = 0._wp 650 zindtbis(1:npti,:) = 0._wp 651 zdiagbis(1:npti,:) = 0._wp 652 653 DO jm = nlay_s + 2, nlay_s + nlay_i 654 DO ji = 1, npti 655 jk = jm - nlay_s - 1 656 ztrid(ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 657 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 658 ztrid(ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk) 659 zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 660 END DO 661 ENDDO 662 663 jm = nlay_s + nlay_i + 1 664 DO ji = 1, npti 665 !!ice bottom term 666 ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) 667 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 668 ztrid(ji,jm,3) = 0.0 669 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 670 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 671 ENDDO 672 673 674 DO ji = 1, npti 675 ! !---------------------! 676 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! 677 ! !---------------------! 678 ! snow interior terms (bottom equation has the same form as the others) 679 DO jm = 3, nlay_s + 1 680 jk = jm - 1 681 ztrid(ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 682 ztrid(ji,jm,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 683 ztrid(ji,jm,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 684 zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 685 END DO 686 687 ! case of only one layer in the ice (ice equation is altered) 688 IF ( nlay_i == 1 ) THEN 689 ztrid(ji,nlay_s+2,3) = 0.0 690 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 691 ENDIF 692 693 jm_min(ji) = 2 694 jm_max(ji) = nlay_i + nlay_s + 1 695 696 ! first layer of snow equation 697 ztrid(ji,2,1) = 0.0 698 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1) 699 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 700 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * & 701 & ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 702 703 ! !---------------------! 704 ELSE ! cells without snow ! 705 ! !---------------------! 706 707 jm_min(ji) = nlay_s + 2 708 jm_max(ji) = nlay_i + nlay_s + 1 709 710 ! first layer of ice equation 711 ztrid(ji,jm_min(ji),1) = 0.0 712 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 713 ztrid(ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 714 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * & 715 & ( zradab_i(ji,1) + qcn_ice_1d(ji) ) 716 717 ! case of only one layer in the ice (surface & ice equations are altered) 718 IF ( nlay_i == 1 ) THEN 719 ztrid(ji,jm_min(ji),1) = 0.0 720 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 721 ztrid(ji,jm_min(ji),3) = 0.0 722 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) & 723 & + qcn_ice_1d(ji) ) 724 725 ENDIF 726 727 ENDIF 728 ! 729 zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) 730 zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 731 ! 732 END DO 733 ! 734 !------------------------------ 735 ! 8) tridiagonal system solving 736 !------------------------------ 737 ! Solve the tridiagonal system with Gauss elimination method. 738 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 739 jm_maxt = 0 740 jm_mint = nlay_i+5 741 DO ji = 1, npti 742 jm_mint = MIN(jm_min(ji),jm_mint) 743 jm_maxt = MAX(jm_max(ji),jm_maxt) 744 END DO 745 746 DO jk = jm_mint+1, jm_maxt 747 DO ji = 1, npti 748 jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 749 zdiagbis(ji,jm) = ztrid(ji,jm,2) - ztrid(ji,jm,1) * ztrid(ji,jm-1,3) / zdiagbis(ji,jm-1) 750 zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 751 END DO 752 END DO 753 754 DO ji = 1, npti 755 ! ice temperatures 756 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 757 END DO 758 759 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 760 DO ji = 1, npti 761 jk = jm - nlay_s - 1 762 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 763 END DO 764 END DO 765 766 DO ji = 1, npti 767 ! snow temperatures 768 IF( h_s_1d(ji) > 0._wp ) THEN 769 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) & 770 & / zdiagbis(ji,nlay_s+1) 771 ENDIF 772 END DO 773 ! 774 !-------------------------------------------------------------- 775 ! 9) Has the scheme converged ?, end of the iterative procedure 776 !-------------------------------------------------------------- 777 ! check that nowhere it has started to melt 778 ! zdti_max is a measure of error, it has to be under zdti_bnd 779 zdti_max = 0._wp 780 781 DO jk = 1, nlay_s 782 DO ji = 1, npti 783 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 784 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 785 END DO 786 END DO 787 788 DO jk = 1, nlay_i 789 DO ji = 1, npti 790 ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0 791 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) 792 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 793 END DO 794 END DO 795 796 ! Compute spatial maximum over all errors 797 ! note that this could be optimized substantially by iterating only the non-converging points 798 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 799 800 ENDIF ! k_jules 801 568 802 END DO ! End of the do while iterative procedure 569 803 570 804 IF( ln_icectl .AND. lwp ) THEN 571 805 WRITE(numout,*) ' zdti_max : ', zdti_max 572 806 WRITE(numout,*) ' iconv : ', iconv 573 807 ENDIF 808 574 809 ! 575 810 !----------------------------- 576 811 ! 10) Fluxes at the interfaces 577 812 !----------------------------- 813 ! 814 ! --- update conduction fluxes 815 ! 578 816 DO ji = 1, npti 579 817 ! ! surface ice conduction flux 580 818 fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & 581 819 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji)) 582 820 ! ! bottom ice conduction flux 583 821 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 584 822 END DO 585 586 ! --- computes sea ice energy of melting compulsory for icethd_dh --- ! 587 CALL ice_thd_enmelt 588 589 ! --- diagnose the change in non-solar flux due to surface temperature change --- ! 590 IF ( ln_dqns_i ) THEN 823 824 ! 825 ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- ! 826 ! 827 IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN ! OFF or SND mode 591 828 DO ji = 1, npti 592 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 593 END DO 594 END IF 595 596 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 597 ! hfx_dif = Heat flux used to warm/cool ice in W.m-2 598 ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 599 DO ji = 1, npti 600 zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & 601 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 602 603 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 604 zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 605 ELSE ! case T_su = 0degC 606 zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 607 ENDIF 608 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 609 610 ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation) 611 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 612 613 END DO 614 615 ! --- Diagnostics SIMIP --- ! 616 DO ji = 1, npti 617 !--- Conduction fluxes (positive downwards) 618 diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 619 diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji) * a_i_1d(ji) / at_i_1d(ji) 620 621 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 622 zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 623 IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 624 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + & 625 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 626 ELSE 627 t_si_1d(ji) = t_su_1d(ji) 628 ENDIF 629 END DO 630 ! 631 END SUBROUTINE ice_thd_zdf 829 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 830 END DO 831 ELSE ! RCV mode 832 DO ji = 1, npti 833 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji) 834 END DO 835 ENDIF 836 837 ! 838 ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- ! 839 ! 840 841 IF ( ( k_jules == np_zdf_jules_OFF ) .OR. ( k_jules == np_zdf_jules_RCV ) ) THEN ! OFF 842 843 CALL ice_thd_enmelt 844 845 ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 846 DO ji = 1, npti 847 zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & 848 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 849 850 IF ( ( k_jules == np_zdf_jules_OFF ) ) THEN 851 852 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 853 zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 854 ELSE ! case T_su = 0degC 855 zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 856 ENDIF 857 858 ELSE ! RCV CASE 859 860 zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 861 862 ENDIF 863 ! 864 ! total heat sink to be sent to the ocean 865 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 866 ! 867 ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2 868 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 869 ! 870 END DO 871 ! 872 ! --- SIMIP diagnostics 873 ! 874 DO ji = 1, npti 875 !--- Conduction fluxes (positive downwards) 876 diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 877 diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji) * a_i_1d(ji) / at_i_1d(ji) 878 879 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 880 zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 881 IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 882 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + & 883 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 884 ELSE 885 t_si_1d(ji) = t_su_1d(ji) 886 ENDIF 887 END DO 888 ! 889 ENDIF 890 ! 891 !--------------------------------------------------------------------------------------- 892 ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes 893 !--------------------------------------------------------------------------------------- 894 ! 895 IF ( k_jules == np_zdf_jules_SND ) THEN ! --- Jules coupling in "SND" mode 896 ! 897 ! Restore temperatures to their initial values 898 t_s_1d(1:npti,:) = ztsold (1:npti,:) 899 t_i_1d(1:npti,:) = ztiold (1:npti,:) 900 qcn_ice_1d(1:npti) = fc_su(1:npti) 901 ! 902 ENDIF 903 ! 904 END SUBROUTINE ice_thd_zdf_BL99 632 905 633 906 … … 675 948 !! ** input : Namelist namthd_zdf 676 949 !!------------------------------------------------------------------- 677 INTEGER :: ios ! Local integer output status for namelist read950 INTEGER :: ios, ioptio ! Local integer 678 951 !! 679 NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i , ln_dqns_i952 NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i 680 953 !!------------------------------------------------------------------- 681 954 ! … … 694 967 WRITE(numout,*) '~~~~~~~~~~~~~~~~' 695 968 WRITE(numout,*) ' Namelist namthd_zdf:' 696 WRITE(numout,*) ' Diffusion follows a Bitz and Lipscomb (1999)ln_zdf_BL99 = ', ln_zdf_BL99969 WRITE(numout,*) ' Bitz and Lipscomb (1999) formulation ln_zdf_BL99 = ', ln_zdf_BL99 697 970 WRITE(numout,*) ' thermal conductivity in the ice (Untersteiner 1964) ln_cndi_U64 = ', ln_cndi_U64 698 971 WRITE(numout,*) ' thermal conductivity in the ice (Pringle et al 2007) ln_cndi_P07 = ', ln_cndi_P07 699 972 WRITE(numout,*) ' thermal conductivity in the snow rn_cnd_s = ', rn_cnd_s 700 973 WRITE(numout,*) ' extinction radiation parameter in sea ice rn_kappa_i = ', rn_kappa_i 701 WRITE(numout,*) ' change the surface non-solar flux with Tsu or not ln_dqns_i = ', ln_dqns_i702 974 ENDIF 975 703 976 ! 704 977 IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN 705 978 CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' ) 706 979 ENDIF 980 981 ! !== set the choice of ice vertical thermodynamic formulation ==! 982 ioptio = 0 983 ! !--- BL99 thermo dynamics (linear liquidus + constant thermal properties) 984 IF( ln_zdf_BL99 ) THEN ; ioptio = ioptio + 1 ; nice_zdf = np_BL99 ; ENDIF 985 IF( ioptio /= 1 ) CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' ) 707 986 ! 708 987 END SUBROUTINE ice_thd_zdf_init
Note: See TracChangeset
for help on using the changeset viewer.
