Changeset 8752 for branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
- Timestamp:
- 2017-11-20T13:54:32+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
r8738 r8752 34 34 LOGICAL :: ln_cndi_U64 ! thermal conductivity: Untersteiner (1964) 35 35 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 36 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 37 REAL(wp), PUBLIC :: rn_cnd_s ! thermal conductivity of the snow [W/m/K] 38 39 INTEGER :: nice_zdf ! Choice of the type of vertical heat diffusion formulation 40 ! ! associated indices: 41 INTEGER, PARAMETER :: np_BL99 = 1 ! Bitz and Lipscomb (1999) 42 43 INTEGER , PARAMETER :: np_zdf_jules_OFF = 0 ! compute all temperatures from qsr and qns 44 INTEGER , PARAMETER :: np_zdf_jules_SND = 1 ! compute conductive heat flux and surface temperature from qsr and qns 45 INTEGER , PARAMETER :: np_zdf_jules_RCV = 2 ! compute snow and inner ice temperatures from qcnd 46 40 47 !!---------------------------------------------------------------------- 41 48 !! NEMO/ICE 4.0 , NEMO Consortium (2017) … … 45 52 CONTAINS 46 53 47 SUBROUTINE ice_thd_zdf 54 SUBROUTINE ice_thd_zdf 55 56 !! 48 57 !!------------------------------------------------------------------- 49 58 !! *** ROUTINE ice_thd_zdf *** 50 59 !! ** Purpose : 60 !! This chooses between the appropriate routine for the 61 !! computation of vertical diffusion 62 !! 63 !!------------------------------------------------------------------- 64 !! 65 66 SELECT CASE ( nice_zdf ) ! Choose the vertical heat diffusion solver 67 68 !------------- 69 CASE( np_BL99 ) ! BL99 solver 70 !------------- 71 72 IF ( nice_jules == np_jules_OFF ) THEN ! No Jules coupler 73 74 CALL ice_thd_zdf_BL99 ( np_zdf_jules_OFF ) 75 76 ENDIF 77 78 IF ( nice_jules == np_jules_EMULE ) THEN ! Jules coupler is emulated 79 80 CALL ice_thd_zdf_BL99 ( np_zdf_jules_SND ) 81 CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) 82 83 ENDIF 84 85 IF ( nice_jules == np_jules_ACTIVE ) THEN ! Jules coupler is emulated 86 87 CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) 88 89 ENDIF 90 91 END SELECT 92 93 END SUBROUTINE ice_thd_zdf 94 95 96 97 SUBROUTINE ice_thd_zdf_BL99(k_jules) 98 !!------------------------------------------------------------------- 99 !! *** ROUTINE ice_thd_zdf_BL99 *** 100 !! ** Purpose : 51 101 !! This routine determines the time evolution of snow and sea-ice 52 !! temperature profiles. 102 !! temperature profiles, using the original Bitz and Lipscomb (1999) algorithm 103 !! 53 104 !! ** Method : 54 105 !! This is done by solving the heat equation diffusion with … … 83 134 !! total ice/snow thickness : h_i_1d, h_s_1d 84 135 !!------------------------------------------------------------------- 136 INTEGER, INTENT(in) :: k_jules ! Jules coupling (0=OFF, 1=RECEIVE, 2=SEND) 137 85 138 INTEGER :: ji, jk ! spatial loop index 86 139 INTEGER :: jm ! current reference number of equation … … 101 154 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature 102 155 REAL(wp) :: ztmelt_i ! ice melting temperature 103 REAL(wp) :: z1_hsu104 156 REAL(wp) :: zdti_max ! current maximal error on temperature 105 157 REAL(wp) :: zcpi ! Ice specific heat … … 111 163 REAL(wp), DIMENSION(jpij) :: zh_i, z1_h_i ! ice layer thickness 112 164 REAL(wp), DIMENSION(jpij) :: zh_s, z1_h_s ! snow layer thickness 113 REAL(wp), DIMENSION(jpij) :: zfsw ! solar radiation absorbed at the surface114 165 REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface 115 166 REAL(wp), DIMENSION(jpij) :: zfnet ! surface flux function 116 167 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 168 169 REAL(wp), DIMENSION(jpij ) :: ztsuold ! Old surface temperature in the ice 119 170 REAL(wp), DIMENSION(jpij,nlay_i) :: ztiold ! Old temperature in the ice 120 171 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsold ! Old temperature in the snow … … 136 187 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat 137 188 REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 189 190 REAL(wp) :: zfr1, zfr2, zfrqsr_tr_i ! dummy factor 138 191 139 192 ! Mono-category … … 166 219 ELSEWHERE ; z1_h_s(1:npti) = 0._wp 167 220 END WHERE 168 ! 169 ! temperatures 170 ztsub (1:npti) = t_su_1d(1:npti) ! temperature at the previous iteration 221 222 ! Store initial temperatures and non solar heat fluxes 223 IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN ! OFF or SND mode 224 225 ztsub (1:npti) = t_su_1d(1:npti) ! surface temperature at iteration n-1 226 ztsuold(1:npti) = t_su_1d(1:npti) ! surface temperature initial value 227 zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti) ! derivative of incoming nonsolar flux 228 zqns_ice_b (1:npti) = qns_ice_1d(1:npti) ! store previous qns_ice_1d value 229 230 t_su_1d(1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err ) ! required to leave the choice between melting or not 231 232 ENDIF 233 171 234 ztsold (1:npti,:) = t_s_1d(1:npti,:) ! Old snow temperature 172 235 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 ! 236 175 237 !------------- 176 238 ! 2) Radiation 177 239 !------------- 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 240 ! --- Transmission/absorption of solar radiation in the ice --- ! 201 zradtr_s(1:npti,0) = zftrice(1:npti) 241 ! zfr1 = ( 0.18 * ( 1.0 - 0.81 ) + 0.35 * 0.81 ) ! standard value 242 ! zfr2 = ( 0.82 * ( 1.0 - 0.81 ) + 0.65 * 0.81 ) ! zfr2 such that zfr1 + zfr2 to equal 1 243 244 ! DO ji = 1, npti 245 246 ! zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * 10._wp ) ) 247 248 ! zfrqsr_tr_i = zfr1 + zfac * zfr2 ! below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness 249 ! IF ( h_s_1d(ji) >= 0.0_wp ) zfrqsr_tr_i = 0._wp ! snow fully opaque 250 251 ! qsr_ice_tr_1d(ji) = zfrqsr_tr_i * qsr_ice_1d(ji) ! transmitted solar radiation 252 253 ! zfsw(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) 254 ! zftrice(ji) = qsr_ice_tr_1d(ji) 255 ! i0(ji) = zfrqsr_tr_i 256 257 ! END DO 258 259 zradtr_s(1:npti,0) = qsr_ice_tr_1d(1:npti) 202 260 DO jk = 1, nlay_s 203 261 DO ji = 1, npti … … 209 267 END DO 210 268 211 zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + zftrice(1:npti) * ( 1._wp - isnow(1:npti) )269 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 270 DO jk = 1, nlay_i 213 271 DO ji = 1, npti … … 330 388 END DO 331 389 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 390 391 ! 392 !----------------------------------------! 393 ! ! 394 ! JULES IF (OFF or SND MODE) ! 395 ! ! 396 !----------------------------------------! 397 ! 398 399 IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN ! OFF or SND mode 400 401 ! 402 ! In OFF mode the original BL99 temperature computation is used 403 ! (with qsr_ice, qns_ice and dqns_ice as inputs) 404 ! 405 ! In SND mode, the computation is required to compute the conduction fluxes 406 ! 407 408 !---------------------------- 409 ! 6) surface flux computation 410 !---------------------------- 411 412 DO ji = 1, npti 413 ! update of the non solar flux according to the update in T_su 339 414 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 340 415 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 416 417 DO ji = 1, npti 418 zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net solar - transmitted solar + non solar 419 END DO 420 ! 421 !---------------------------- 422 ! 7) tridiagonal system terms 423 !---------------------------- 424 !!layer denotes the number of the layer in the snow or in the ice 425 !!jm denotes the reference number of the equation in the tridiagonal 426 !!system, terms of tridiagonal system are indexed as following : 427 !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 428 429 !!ice interior terms (top equation has the same form as the others) 430 ztrid (1:npti,:,:) = 0._wp 431 zindterm(1:npti,:) = 0._wp 432 zindtbis(1:npti,:) = 0._wp 433 zdiagbis(1:npti,:) = 0._wp 434 435 DO jm = nlay_s + 2, nlay_s + nlay_i 362 436 DO ji = 1, npti 363 437 jk = jm - nlay_s - 1 … … 367 441 zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 368 442 END DO 369 ENDDO370 371 jm = nlay_s + nlay_i + 1372 DO ji = 1, npti443 ENDDO 444 445 jm = nlay_s + nlay_i + 1 446 DO ji = 1, npti 373 447 !!ice bottom term 374 448 ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) … … 377 451 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 378 452 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 379 ENDDO380 381 382 DO ji = 1, npti453 ENDDO 454 455 456 DO ji = 1, npti 383 457 ! !---------------------! 384 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells !458 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! 385 459 ! !---------------------! 386 460 ! snow interior terms (bottom equation has the same form as the others) 387 461 DO jm = 3, nlay_s + 1 388 389 390 391 392 462 jk = jm - 1 463 ztrid(ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 464 ztrid(ji,jm,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 465 ztrid(ji,jm,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 466 zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 393 467 END DO 394 468 395 469 ! case of only one layer in the ice (ice equation is altered) 396 470 IF ( nlay_i == 1 ) THEN 397 398 471 ztrid(ji,nlay_s+2,3) = 0.0 472 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 399 473 ENDIF 400 474 401 475 IF ( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting 402 476 403 404 405 406 407 408 409 410 411 412 413 414 415 416 477 jm_min(ji) = 1 478 jm_max(ji) = nlay_i + nlay_s + 1 479 480 ! surface equation 481 ztrid(ji,1,1) = 0.0 482 ztrid(ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 483 ztrid(ji,1,3) = zg1s * zkappa_s(ji,0) 484 zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 485 486 ! first layer of snow equation 487 ztrid(ji,2,1) = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1) 488 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 489 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 490 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 417 491 418 492 ELSE !-- case 2 : surface is melting 419 420 421 422 423 424 425 426 427 428 493 ! 494 jm_min(ji) = 2 495 jm_max(ji) = nlay_i + nlay_s + 1 496 497 ! first layer of snow equation 498 ztrid(ji,2,1) = 0.0 499 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 500 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 501 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * & 502 & ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 429 503 ENDIF 430 504 ! !---------------------! … … 433 507 ! 434 508 IF ( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 509 ! 510 jm_min(ji) = nlay_s + 1 511 jm_max(ji) = nlay_i + nlay_s + 1 512 513 ! surface equation 514 ztrid(ji,jm_min(ji),1) = 0.0 515 ztrid(ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1 516 ztrid(ji,jm_min(ji),3) = zkappa_i(ji,0)*zg1 517 zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji) 518 519 ! first layer of ice equation 520 ztrid(ji,jm_min(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 521 ztrid(ji,jm_min(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 522 ztrid(ji,jm_min(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1) 523 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 524 525 ! case of only one layer in the ice (surface & ice equations are altered) 526 IF ( nlay_i == 1 ) THEN 527 ztrid(ji,jm_min(ji),1) = 0.0 528 ztrid(ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 529 ztrid(ji,jm_min(ji),3) = zkappa_i(ji,0) * 2.0 530 ztrid(ji,jm_min(ji)+1,1) = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 531 ztrid(ji,jm_min(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 532 ztrid(ji,jm_min(ji)+1,3) = 0.0 533 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * & 534 & ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 535 ENDIF 462 536 463 537 ELSE !-- case 2 : surface is melting 464 538 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 539 jm_min(ji) = nlay_s + 2 540 jm_max(ji) = nlay_i + nlay_s + 1 541 542 ! first layer of ice equation 543 ztrid(ji,jm_min(ji),1) = 0.0 544 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 545 ztrid(ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 546 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * & 547 & ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 548 549 ! case of only one layer in the ice (surface & ice equations are altered) 550 IF ( nlay_i == 1 ) THEN 551 ztrid(ji,jm_min(ji),1) = 0.0 552 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 553 ztrid(ji,jm_min(ji),3) = 0.0 554 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) & 555 & + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 556 ENDIF 483 557 484 558 ENDIF … … 488 562 zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 489 563 ! 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, npti564 END DO 565 ! 566 !------------------------------ 567 ! 8) tridiagonal system solving 568 !------------------------------ 569 ! Solve the tridiagonal system with Gauss elimination method. 570 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 571 jm_maxt = 0 572 jm_mint = nlay_i+5 573 DO ji = 1, npti 500 574 jm_mint = MIN(jm_min(ji),jm_mint) 501 575 jm_maxt = MAX(jm_max(ji),jm_maxt) 502 END DO503 504 DO jk = jm_mint+1, jm_maxt576 END DO 577 578 DO jk = jm_mint+1, jm_maxt 505 579 DO ji = 1, npti 506 580 jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) … … 508 582 zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 509 583 END DO 510 END DO511 512 DO ji = 1, npti584 END DO 585 586 DO ji = 1, npti 513 587 ! ice temperatures 514 588 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, -1589 END DO 590 591 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 518 592 DO ji = 1, npti 519 593 jk = jm - nlay_s - 1 520 594 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 521 595 END DO 522 END DO523 524 DO ji = 1, npti596 END DO 597 598 DO ji = 1, npti 525 599 ! snow temperatures 526 600 IF( h_s_1d(ji) > 0._wp ) THEN 527 601 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) & 528 602 & / zdiagbis(ji,nlay_s+1) 529 603 ENDIF 530 604 ! surface temperature … … 532 606 IF( t_su_1d(ji) < rt0 ) THEN 533 607 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 534 535 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, npti608 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 609 ENDIF 610 END DO 611 ! 612 !-------------------------------------------------------------- 613 ! 9) Has the scheme converged ?, end of the iterative procedure 614 !-------------------------------------------------------------- 615 ! check that nowhere it has started to melt 616 ! zdti_max is a measure of error, it has to be under zdti_bnd 617 zdti_max = 0._wp 618 DO ji = 1, npti 545 619 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 546 620 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 547 END DO548 549 DO jk = 1, nlay_s621 END DO 622 623 DO jk = 1, nlay_s 550 624 DO ji = 1, npti 551 625 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 552 626 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 553 627 END DO 554 END DO555 556 DO jk = 1, nlay_i628 END DO 629 630 DO jk = 1, nlay_i 557 631 DO ji = 1, npti 558 632 ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0 … … 560 634 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 561 635 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 636 END DO 637 638 ! Compute spatial maximum over all errors 639 ! note that this could be optimized substantially by iterating only the non-converging points 640 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 641 ! 642 !----------------------------------------! 643 ! ! 644 ! JULES IF (RCV MODE) ! 645 ! ! 646 !----------------------------------------! 647 ! 648 649 ELSE IF ( k_jules == np_zdf_jules_RCV ) THEN ! RCV mode 650 651 ! 652 ! In RCV mode, we use a modified BL99 solver 653 ! with conduction flux (qcn_ice) as forcing term 654 ! 655 !---------------------------- 656 ! 7) tridiagonal system terms 657 !---------------------------- 658 !!layer denotes the number of the layer in the snow or in the ice 659 !!jm denotes the reference number of the equation in the tridiagonal 660 !!system, terms of tridiagonal system are indexed as following : 661 !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 662 663 !!ice interior terms (top equation has the same form as the others) 664 ztrid (1:npti,:,:) = 0._wp 665 zindterm(1:npti,:) = 0._wp 666 zindtbis(1:npti,:) = 0._wp 667 zdiagbis(1:npti,:) = 0._wp 668 669 DO jm = nlay_s + 2, nlay_s + nlay_i 670 DO ji = 1, npti 671 jk = jm - nlay_s - 1 672 ztrid(ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 673 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 674 ztrid(ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk) 675 zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 676 END DO 677 ENDDO 678 679 jm = nlay_s + nlay_i + 1 680 DO ji = 1, npti 681 !!ice bottom term 682 ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) 683 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 684 ztrid(ji,jm,3) = 0.0 685 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 686 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 687 ENDDO 688 689 690 DO ji = 1, npti 691 ! !---------------------! 692 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! 693 ! !---------------------! 694 ! snow interior terms (bottom equation has the same form as the others) 695 DO jm = 3, nlay_s + 1 696 jk = jm - 1 697 ztrid(ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 698 ztrid(ji,jm,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 699 ztrid(ji,jm,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 700 zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 701 END DO 702 703 ! case of only one layer in the ice (ice equation is altered) 704 IF ( nlay_i == 1 ) THEN 705 ztrid(ji,nlay_s+2,3) = 0.0 706 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 707 ENDIF 708 709 jm_min(ji) = 2 710 jm_max(ji) = nlay_i + nlay_s + 1 711 712 ! first layer of snow equation 713 ztrid(ji,2,1) = 0.0 714 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1) 715 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 716 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * & 717 & ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 718 719 ! !---------------------! 720 ELSE ! cells without snow ! 721 ! !---------------------! 722 723 jm_min(ji) = nlay_s + 2 724 jm_max(ji) = nlay_i + nlay_s + 1 725 726 ! first layer of ice equation 727 ztrid(ji,jm_min(ji),1) = 0.0 728 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 729 ztrid(ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 730 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * & 731 & ( zradab_i(ji,1) + qcn_ice_1d(ji) ) 732 733 ! case of only one layer in the ice (surface & ice equations are altered) 734 IF ( nlay_i == 1 ) THEN 735 ztrid(ji,jm_min(ji),1) = 0.0 736 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 737 ztrid(ji,jm_min(ji),3) = 0.0 738 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) & 739 & + qcn_ice_1d(ji) ) 740 741 ENDIF 742 743 ENDIF 744 ! 745 zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) 746 zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 747 ! 748 END DO 749 ! 750 !------------------------------ 751 ! 8) tridiagonal system solving 752 !------------------------------ 753 ! Solve the tridiagonal system with Gauss elimination method. 754 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 755 jm_maxt = 0 756 jm_mint = nlay_i+5 757 DO ji = 1, npti 758 jm_mint = MIN(jm_min(ji),jm_mint) 759 jm_maxt = MAX(jm_max(ji),jm_maxt) 760 END DO 761 762 DO jk = jm_mint+1, jm_maxt 763 DO ji = 1, npti 764 jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 765 zdiagbis(ji,jm) = ztrid(ji,jm,2) - ztrid(ji,jm,1) * ztrid(ji,jm-1,3) / zdiagbis(ji,jm-1) 766 zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 767 END DO 768 END DO 769 770 DO ji = 1, npti 771 ! ice temperatures 772 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 773 END DO 774 775 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 776 DO ji = 1, npti 777 jk = jm - nlay_s - 1 778 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 779 END DO 780 END DO 781 782 DO ji = 1, npti 783 ! snow temperatures 784 IF( h_s_1d(ji) > 0._wp ) THEN 785 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) & 786 & / zdiagbis(ji,nlay_s+1) 787 ENDIF 788 END DO 789 ! 790 !-------------------------------------------------------------- 791 ! 9) Has the scheme converged ?, end of the iterative procedure 792 !-------------------------------------------------------------- 793 ! check that nowhere it has started to melt 794 ! zdti_max is a measure of error, it has to be under zdti_bnd 795 zdti_max = 0._wp 796 797 DO jk = 1, nlay_s 798 DO ji = 1, npti 799 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 800 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 801 END DO 802 END DO 803 804 DO jk = 1, nlay_i 805 DO ji = 1, npti 806 ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0 807 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) 808 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 809 END DO 810 END DO 811 812 ! Compute spatial maximum over all errors 813 ! note that this could be optimized substantially by iterating only the non-converging points 814 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 815 816 ENDIF ! k_jules 817 568 818 END DO ! End of the do while iterative procedure 569 819 570 820 IF( ln_icectl .AND. lwp ) THEN 571 821 WRITE(numout,*) ' zdti_max : ', zdti_max 572 822 WRITE(numout,*) ' iconv : ', iconv 573 823 ENDIF 824 574 825 ! 575 826 !----------------------------- 576 827 ! 10) Fluxes at the interfaces 577 828 !----------------------------- 829 ! 830 ! --- update conduction fluxes 831 ! 578 832 DO ji = 1, npti 579 833 ! ! surface ice conduction flux 580 834 fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & 581 835 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji)) 582 836 ! ! bottom ice conduction flux 583 837 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 584 838 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 839 840 ! 841 ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- ! 842 ! 843 DO ji = 1, npti 844 IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND ) THEN 845 ! OFF or SND mode 846 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 847 ELSE ! RCV mode 848 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji) 849 ENDIF 850 END DO 851 852 ! 853 ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- ! 854 ! 855 856 IF ( ( k_jules == np_zdf_jules_OFF ) .OR. ( k_jules == np_zdf_jules_RCV ) ) THEN ! OFF 857 858 CALL ice_thd_enmelt 859 860 ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 591 861 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) 862 zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & 863 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 864 865 IF ( ( k_jules == np_zdf_jules_OFF ) ) THEN 866 867 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 868 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) 869 ELSE ! case T_su = 0degC 870 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) 871 ENDIF 872 873 ELSE ! RCV CASE 874 875 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) 876 877 ENDIF 878 879 ! total heat sink to be sent to the ocean 880 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 881 882 ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2 883 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 884 885 END DO 886 887 ! 888 ! --- SIMIP diagnostics 889 ! 890 891 DO ji = 1, npti 892 !--- Conduction fluxes (positive downwards) 893 diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 894 diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji) * a_i_1d(ji) / at_i_1d(ji) 895 896 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 897 zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 898 IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 899 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + & 900 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 901 ELSE 902 t_si_1d(ji) = t_su_1d(ji) 903 ENDIF 593 904 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 905 906 ENDIF 907 ! 908 !--------------------------------------------------------------------------------------- 909 ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes 910 !--------------------------------------------------------------------------------------- 911 ! 912 IF ( k_jules == np_zdf_jules_SND ) THEN ! --- Jules coupling in "SND" mode 913 914 ! Restore temperatures to their initial values 915 t_s_1d(1:npti,:) = ztsold (1:npti,:) 916 t_i_1d(1:npti,:) = ztiold (1:npti,:) 917 qcn_ice_1d(1:npti) = fc_su(1:npti) 918 919 ENDIF 920 921 END SUBROUTINE ice_thd_zdf_BL99 922 632 923 633 924 … … 663 954 664 955 956 665 957 SUBROUTINE ice_thd_zdf_init 666 958 !!----------------------------------------------------------------------- … … 675 967 !! ** input : Namelist namthd_zdf 676 968 !!------------------------------------------------------------------- 677 INTEGER :: ios ! Local integer output status for namelist read969 INTEGER :: ios, ioptio ! Local integer output status for namelist read 678 970 !! 679 NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i , ln_dqns_i971 NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i 680 972 !!------------------------------------------------------------------- 681 973 ! … … 694 986 WRITE(numout,*) '~~~~~~~~~~~~~~~~' 695 987 WRITE(numout,*) ' Namelist namthd_zdf:' 696 WRITE(numout,*) ' Diffusion follows a Bitz and Lipscomb (1999)ln_zdf_BL99 = ', ln_zdf_BL99988 WRITE(numout,*) ' Bitz and Lipscomb (1999) formulation ln_zdf_BL99 = ', ln_zdf_BL99 697 989 WRITE(numout,*) ' thermal conductivity in the ice (Untersteiner 1964) ln_cndi_U64 = ', ln_cndi_U64 698 990 WRITE(numout,*) ' thermal conductivity in the ice (Pringle et al 2007) ln_cndi_P07 = ', ln_cndi_P07 699 991 WRITE(numout,*) ' thermal conductivity in the snow rn_cnd_s = ', rn_cnd_s 700 992 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 993 ENDIF 994 703 995 ! 704 996 IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN 705 997 CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' ) 706 998 ENDIF 999 1000 ! !== set the choice of ice vertical thermodynamic formulation ==! 1001 ioptio = 0 1002 ! !--- BL99 thermo dynamics (linear liquidus + constant thermal properties) 1003 IF( ln_zdf_BL99 ) THEN ; ioptio = ioptio + 1 ; nice_zdf = np_BL99 ; ENDIF 1004 IF( ioptio /= 1 ) CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' ) 707 1005 ! 708 1006 END SUBROUTINE ice_thd_zdf_init … … 711 1009 !!---------------------------------------------------------------------- 712 1010 !! Default option Dummy Module No ESIM sea-ice model 713 !!--------------------------------------------------------------------- -1011 !!--------------------------------------------------------------------- 714 1012 #endif 715 1013
Note: See TracChangeset
for help on using the changeset viewer.