Changeset 14789 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/ICE/icethd_zdf_bl99.F90
- Timestamp:
- 2021-05-05T13:18:04+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icethd_zdf_bl99.F90 (modified) (36 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev _r12970_AGRIF_CMEMSext/AGRIF5 ^/vendors/AGRIF/dev@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 ^/vendors/PPR@HEAD ext/PPR 8 9 9 10 # SETTE 10 ^/utils/CI/sette@1 3559sette11 ^/utils/CI/sette@14244 sette
-
- Property svn:externals
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/ICE/icethd_zdf_bl99.F90
r13472 r14789 2 2 !!====================================================================== 3 3 !! *** MODULE icethd_zdf_BL99 *** 4 !! sea-ice: vertical heat diffusion in sea ice (computation of temperatures) 4 !! sea-ice: vertical heat diffusion in sea ice (computation of temperatures) 5 5 !!====================================================================== 6 6 !! History : ! 2003-02 (M. Vancoppenolle) original 1D code … … 15 15 !!---------------------------------------------------------------------- 16 16 USE dom_oce ! ocean space and time domain 17 USE phycst ! physical constants (ocean directory) 17 USE phycst ! physical constants (ocean directory) 18 18 USE ice ! sea-ice: variables 19 19 USE ice1D ! sea-ice: thermodynamics variables … … 44 44 !! 45 45 !! ** Method : solves the heat equation diffusion with a Neumann boundary 46 !! condition at the surface and a Dirichlet one at the bottom. 46 !! condition at the surface and a Dirichlet one at the bottom. 47 47 !! Solar radiation is partially absorbed into the ice. 48 !! The specific heat and thermal conductivities depend on ice 49 !! salinity and temperature to take into account brine pocket 48 !! The specific heat and thermal conductivities depend on ice 49 !! salinity and temperature to take into account brine pocket 50 50 !! melting. The numerical scheme is an iterative Crank-Nicolson 51 51 !! on a non-uniform multilayer grid in the ice and snow system. … … 91 91 REAL(wp) :: zbeta = 0.117_wp ! for thermal conductivity (could be 0.13) 92 92 REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity 93 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C 94 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature 95 REAL(wp) :: zhs_ssl = 0.03_wp ! surface scattering layer in the snow 93 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C 94 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature 95 REAL(wp) :: zhs_ssl = 0.03_wp ! surface scattering layer in the snow 96 96 REAL(wp) :: zhi_ssl = 0.10_wp ! surface scattering layer in the ice 97 97 REAL(wp) :: zh_min = 1.e-3_wp ! minimum ice/snow thickness for conduction 98 98 REAL(wp) :: ztmelts ! ice melting temperature 99 REAL(wp) :: zdti_max ! current maximal error on temperature 99 REAL(wp) :: zdti_max ! current maximal error on temperature 100 100 REAL(wp) :: zcpi ! Ice specific heat 101 101 REAL(wp) :: zhfx_err, zdq ! diag errors on heat … … 109 109 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function 110 110 ! 111 REAL(wp), DIMENSION(jpij ) 112 REAL(wp), DIMENSION(jpij,nlay_i) 113 REAL(wp), DIMENSION(jpij,nlay_s) 114 REAL(wp), DIMENSION(jpij,nlay_i) 115 REAL(wp), DIMENSION(jpij,nlay_s) 116 REAL(wp), DIMENSION(jpij,0:nlay_i) 117 REAL(wp), DIMENSION(jpij,0:nlay_i) 118 REAL(wp), DIMENSION(jpij,0:nlay_i) 119 REAL(wp), DIMENSION(jpij,0:nlay_i) 120 REAL(wp), DIMENSION(jpij,0:nlay_i) 121 REAL(wp), DIMENSION(jpij,0:nlay_i) 122 REAL(wp), DIMENSION(jpij,0:nlay_s) 123 REAL(wp), DIMENSION(jpij,0:nlay_s) 124 REAL(wp), DIMENSION(jpij,0:nlay_s) 125 REAL(wp), DIMENSION(jpij,0:nlay_s) 126 REAL(wp), DIMENSION(jpij) 127 REAL(wp), DIMENSION(jpij ,nlay_i+3) :: zindterm ! 'Ind'ependent term128 REAL(wp), DIMENSION(jpij ,nlay_i+3) :: zindtbis ! Temporary 'ind'ependent term129 REAL(wp), DIMENSION(jpij ,nlay_i+3) :: zdiagbis ! Temporary 'dia'gonal term130 REAL(wp), DIMENSION(jpij ,nlay_i+3,3) :: ztrid ! Tridiagonal system terms131 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat132 REAL(wp), DIMENSION(jpij ) :: zghe ! G(he), th. conduct enhancement factor, mono-cat133 REAL(wp), DIMENSION(jpij ) :: za_s_fra ! ice fraction covered by snow134 REAL(wp), DIMENSION(jpij ) :: isnow ! snow presence (1) or not (0)135 REAL(wp), DIMENSION(jpij ) :: isnow_comb ! snow presence for met-office111 REAL(wp), DIMENSION(jpij ) :: ztsuold ! Old surface temperature in the ice 112 REAL(wp), DIMENSION(jpij,nlay_i) :: ztiold ! Old temperature in the ice 113 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsold ! Old temperature in the snow 114 REAL(wp), DIMENSION(jpij,nlay_i) :: ztib ! Temporary temperature in the ice to check the convergence 115 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow to check the convergence 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 118 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice 119 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice 120 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zkappa_i ! Kappa factor in the ice 121 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zeta_i ! Eta factor in the ice 122 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow 123 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow 124 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zkappa_s ! Kappa factor in the snow 125 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeta_s ! Eta factor in the snow 126 REAL(wp), DIMENSION(jpij) :: zkappa_comb ! Combined snow and ice surface conductivity 127 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat 128 REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 129 REAL(wp), DIMENSION(jpij) :: za_s_fra ! ice fraction covered by snow 130 REAL(wp), DIMENSION(jpij) :: isnow ! snow presence (1) or not (0) 131 REAL(wp), DIMENSION(jpij) :: isnow_comb ! snow presence for met-office 132 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zindterm ! 'Ind'ependent term 133 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zindtbis ! Temporary 'ind'ependent term 134 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zdiagbis ! Temporary 'dia'gonal term 135 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1,3) :: ztrid ! Tridiagonal system terms 136 136 ! 137 137 ! Mono-category … … 139 139 REAL(wp) :: zhe ! dummy factor 140 140 REAL(wp) :: zcnd_i ! mean sea ice thermal conductivity 141 !!------------------------------------------------------------------ 141 !!------------------------------------------------------------------ 142 142 143 143 ! --- diag error on heat diffusion - PART 1 --- ! 144 144 DO ji = 1, npti 145 145 zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & 146 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 146 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 147 147 END DO 148 148 149 149 ! calculate ice fraction covered by snow for radiation 150 150 CALL ice_var_snwfra( h_s_1d(1:npti), za_s_fra(1:npti) ) 151 151 152 152 !------------------ 153 153 ! 1) Initialization … … 155 155 ! 156 156 ! extinction radiation in the snow 157 IF ( nn_qtrice == 0 ) THEN ! constant 157 IF ( nn_qtrice == 0 ) THEN ! constant 158 158 zraext_s(1:npti) = rn_kappa_s 159 159 ELSEIF( nn_qtrice == 1 ) THEN ! depends on melting/freezing conditions … … 166 166 DO ji = 1, npti 167 167 ! ice thickness 168 IF( h_i_1d(ji) > 0._wp ) THEN 168 IF( h_i_1d(ji) > 0._wp ) THEN 169 169 zh_i (ji) = MAX( zh_min , h_i_1d(ji) ) * r1_nlay_i ! set a minimum thickness for conduction 170 170 z1_h_i(ji) = 1._wp / zh_i(ji) ! it must be very small … … 198 198 ztsuold (1:npti) = t_su_1d(1:npti) ! surface temperature initial value 199 199 t_su_1d (1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err ) ! required to leave the choice between melting or not 200 zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti) ! derivative of incoming nonsolar flux 200 zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti) ! derivative of incoming nonsolar flux 201 201 zqns_ice_b (1:npti) = qns_ice_1d(1:npti) ! store previous qns_ice_1d value 202 202 ! … … 221 221 ! 222 222 zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * za_s_fra(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - za_s_fra(1:npti) ) 223 DO jk = 1, nlay_i 223 DO jk = 1, nlay_i 224 224 DO ji = 1, npti 225 225 ! ! radiation transmitted below the layer-th ice layer … … 227 227 & * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zh_min ) ) & 228 228 & + ( 1._wp - za_s_fra(ji) ) * qtr_ice_top_1d(ji) & ! part snow free 229 & * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) ) 229 & * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) ) 230 230 ! ! radiation absorbed by the layer-th ice layer 231 231 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) … … 288 288 DO ji = 1, npti 289 289 IF ( .NOT. l_T_converged(ji) ) & 290 ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) ) 290 ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) ) 291 291 END DO 292 292 ! … … 401 401 zdiagbis(1:npti,:) = 0._wp 402 402 403 DO jm = nlay_s + 2, nlay_s + nlay_i 403 DO jm = nlay_s + 2, nlay_s + nlay_i 404 404 DO ji = 1, npti 405 405 jk = jm - nlay_s - 1 … … 414 414 DO ji = 1, npti 415 415 ! ice bottom term 416 ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1) 416 ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1) 417 417 ztrid (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 ) 418 418 ztrid (ji,jm,3) = 0._wp 419 419 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 420 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 420 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 421 421 END DO 422 422 … … 433 433 zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 434 434 END DO 435 435 436 436 ! case of only one layer in the ice (ice equation is altered) 437 437 IF( nlay_i == 1 ) THEN 438 438 ztrid (ji,nlay_s+2,3) = 0._wp 439 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 440 ENDIF 441 439 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 440 ENDIF 441 442 442 IF( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting 443 443 444 444 jm_min(ji) = 1 445 445 jm_max(ji) = nlay_i + nlay_s + 1 446 446 447 447 ! surface equation 448 448 ztrid (ji,1,1) = 0._wp … … 450 450 ztrid (ji,1,3) = zg1s * zkappa_s(ji,0) 451 451 zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 452 452 453 453 ! first layer of snow equation 454 454 ztrid (ji,2,1) = - zeta_s(ji,1) * zkappa_s(ji,0) * zg1s … … 456 456 ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1) 457 457 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 458 458 459 459 ELSE !-- case 2 : surface is melting 460 460 ! 461 461 jm_min(ji) = 2 462 462 jm_max(ji) = nlay_i + nlay_s + 1 463 463 464 464 ! first layer of snow equation 465 465 ztrid (ji,2,1) = 0._wp 466 466 ztrid (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 467 ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1) 468 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 467 ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1) 468 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 469 469 ENDIF 470 470 ! !---------------------! … … 476 476 jm_min(ji) = nlay_s + 1 477 477 jm_max(ji) = nlay_i + nlay_s + 1 478 479 ! surface equation 478 479 ! surface equation 480 480 ztrid (ji,jm_min(ji),1) = 0._wp 481 ztrid (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * zg1 481 ztrid (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * zg1 482 482 ztrid (ji,jm_min(ji),3) = zkappa_i(ji,0) * zg1 483 483 zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 484 484 485 485 ! first layer of ice equation 486 486 ztrid (ji,jm_min(ji)+1,1) = - zeta_i(ji,1) * zkappa_i(ji,0) * zg1 487 487 ztrid (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 488 ztrid (ji,jm_min(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1) 489 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 490 488 ztrid (ji,jm_min(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1) 489 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 490 491 491 ! case of only one layer in the ice (surface & ice equations are altered) 492 492 IF( nlay_i == 1 ) THEN … … 499 499 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji)) 500 500 ENDIF 501 501 502 502 ELSE !-- case 2 : surface is melting 503 503 504 504 jm_min(ji) = nlay_s + 2 505 505 jm_max(ji) = nlay_i + nlay_s + 1 506 506 507 507 ! first layer of ice equation 508 508 ztrid (ji,jm_min(ji),1) = 0._wp 509 ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 509 ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 510 510 ztrid (ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 511 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji)) 512 511 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji)) 512 513 513 ! case of only one layer in the ice (surface & ice equations are altered) 514 514 IF( nlay_i == 1 ) THEN … … 519 519 & + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2._wp 520 520 ENDIF 521 521 522 522 ENDIF 523 523 ENDIF … … 533 533 ! Solve the tridiagonal system with Gauss elimination method. 534 534 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 535 jm_maxt = 0 536 jm_mint = nlay_i+5 537 DO ji = 1, npti 538 jm_mint = MIN(jm_min(ji),jm_mint) 539 jm_maxt = MAX(jm_max(ji),jm_maxt) 540 END DO 541 542 DO jk = jm_mint+1, jm_maxt 543 DO ji = 1, npti 544 jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 535 !!$ jm_maxt = 0 536 !!$ jm_mint = nlay_i+5 537 !!$ DO ji = 1, npti 538 !!$ jm_mint = MIN(jm_min(ji),jm_mint) 539 !!$ jm_maxt = MAX(jm_max(ji),jm_maxt) 540 !!$ END DO 541 !!$ !!clem SNWLAY => check why LIM1D does not get this loop. Is nlay_i+5 correct? 542 !!$ 543 !!$ DO jk = jm_mint+1, jm_maxt 544 !!$ DO ji = 1, npti 545 !!$ jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 546 !!$ zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 547 !!$ zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1) 548 !!$ END DO 549 !!$ END DO 550 ! clem: maybe one should find a way to reverse this loop for mpi performance 551 DO ji = 1, npti 552 jm_mint = jm_min(ji) 553 jm_maxt = jm_max(ji) 554 DO jm = jm_mint+1, jm_maxt 545 555 zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 546 556 zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1) … … 564 574 END DO 565 575 576 ! snow temperatures 566 577 DO ji = 1, npti 567 578 ! Variables used after iterations 568 579 ! Value must be frozen after convergence for MPP independance reason 569 IF ( .NOT. l_T_converged(ji) ) THEN 570 ! snow temperatures 571 IF( h_s_1d(ji) > 0._wp ) THEN 572 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) 573 ENDIF 574 ! surface temperature 580 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 581 & 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) 582 END DO 583 !!clem SNWLAY 584 DO jm = nlay_s, 2, -1 585 DO ji = 1, npti 586 jk = jm - 1 587 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 588 & t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm) 589 END DO 590 END DO 591 592 ! surface temperature 593 DO ji = 1, npti 594 IF( .NOT. l_T_converged(ji) ) THEN 575 595 ztsub(ji) = t_su_1d(ji) 576 596 IF( t_su_1d(ji) < rt0 ) THEN 577 t_su_1d(ji) = ( 578 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) *t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))597 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 598 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 579 599 ENDIF 580 600 ENDIF 581 601 END DO 582 !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1)583 602 ! 584 603 !-------------------------------------------------------------- … … 609 628 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 610 629 END DO 611 630 612 631 ! convergence test 613 632 IF( ln_zdf_chkcvg ) THEN … … 646 665 zdiagbis(1:npti,:) = 0._wp 647 666 648 DO jm = nlay_s + 2, nlay_s + nlay_i 667 DO jm = nlay_s + 2, nlay_s + nlay_i 649 668 DO ji = 1, npti 650 669 jk = jm - nlay_s - 1 … … 659 678 DO ji = 1, npti 660 679 ! ice bottom term 661 ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1) 680 ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1) 662 681 ztrid (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 ) 663 682 ztrid (ji,jm,3) = 0._wp 664 683 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 665 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 684 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 666 685 ENDDO 667 686 … … 678 697 zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 679 698 END DO 680 699 681 700 ! case of only one layer in the ice (ice equation is altered) 682 701 IF ( nlay_i == 1 ) THEN 683 702 ztrid (ji,nlay_s+2,3) = 0._wp 684 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 685 ENDIF 686 703 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 704 ENDIF 705 687 706 jm_min(ji) = 2 688 707 jm_max(ji) = nlay_i + nlay_s + 1 689 708 690 709 ! first layer of snow equation 691 710 ztrid (ji,2,1) = 0._wp 692 711 ztrid (ji,2,2) = 1._wp + zeta_s(ji,1) * zkappa_s(ji,1) 693 ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1) 694 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 695 712 ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1) 713 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 714 696 715 ! !---------------------! 697 716 ELSE ! cells without snow ! … … 699 718 jm_min(ji) = nlay_s + 2 700 719 jm_max(ji) = nlay_i + nlay_s + 1 701 720 702 721 ! first layer of ice equation 703 722 ztrid (ji,jm_min(ji),1) = 0._wp 704 ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1) 723 ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1) 705 724 ztrid (ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 706 725 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) ) 707 726 708 727 ! case of only one layer in the ice (surface & ice equations are altered) 709 728 IF( nlay_i == 1 ) THEN … … 714 733 & ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) + qcn_ice_1d(ji) ) 715 734 ENDIF 716 735 717 736 ENDIF 718 737 ! … … 727 746 ! Solve the tridiagonal system with Gauss elimination method. 728 747 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 729 jm_maxt = 0 730 jm_mint = nlay_i+5 731 DO ji = 1, npti 732 jm_mint = MIN(jm_min(ji),jm_mint) 733 jm_maxt = MAX(jm_max(ji),jm_maxt) 734 END DO 735 736 DO jk = jm_mint+1, jm_maxt 737 DO ji = 1, npti 738 jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 748 !!$ jm_maxt = 0 749 !!$ jm_mint = nlay_i+5 750 !!$ DO ji = 1, npti 751 !!$ jm_mint = MIN(jm_min(ji),jm_mint) 752 !!$ jm_maxt = MAX(jm_max(ji),jm_maxt) 753 !!$ END DO 754 !!$ 755 !!$ DO jk = jm_mint+1, jm_maxt 756 !!$ DO ji = 1, npti 757 !!$ jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 758 !!$ zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 759 !!$ zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 760 !!$ END DO 761 !!$ END DO 762 ! clem: maybe one should find a way to reverse this loop for mpi performance 763 DO ji = 1, npti 764 jm_mint = jm_min(ji) 765 jm_maxt = jm_max(ji) 766 DO jm = jm_mint+1, jm_maxt 739 767 zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 740 zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1)/ zdiagbis(ji,jm-1)768 zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1) 741 769 END DO 742 770 END DO 743 771 744 772 ! ice temperatures 745 773 DO ji = 1, npti … … 758 786 END DO 759 787 END DO 760 761 ! snow temperatures 762 DO ji = 1, npti 763 ! Variable used after iterations788 789 ! snow temperatures 790 DO ji = 1, npti 791 ! Variables used after iterations 764 792 ! Value must be frozen after convergence for MPP independance reason 765 IF ( .NOT. l_T_converged(ji) ) THEN 766 IF( h_s_1d(ji) > 0._wp ) THEN 767 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) 768 ENDIF 769 ENDIF 770 END DO 771 !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1) 793 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 794 & 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) 795 END DO 796 !!clem SNWLAY 797 DO jm = nlay_s, 2, -1 798 DO ji = 1, npti 799 jk = jm - 1 800 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 801 & t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm) 802 END DO 803 END DO 772 804 ! 773 805 !-------------------------------------------------------------- … … 791 823 792 824 DO jk = 1, nlay_i 793 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 825 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 794 826 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 795 827 zdti_max = MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) … … 853 885 ! 854 886 DO ji = 1, npti 855 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 887 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 856 888 END DO 857 889 ! … … 861 893 ! 862 894 IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_ON ) THEN 863 864 CALL ice_var_enthalpy 865 895 896 CALL ice_var_enthalpy 897 866 898 ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 867 899 DO ji = 1, npti 868 900 zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & 869 901 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 870 902 871 903 IF( k_cnd == np_cnd_OFF ) THEN 872 904 873 905 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 874 906 zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji) & … … 878 910 & + zdq * r1_Dt_ice ) * a_i_1d(ji) 879 911 ENDIF 880 912 881 913 ELSEIF( k_cnd == np_cnd_ON ) THEN 882 914 883 915 zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji) & 884 916 & + zdq * r1_Dt_ice ) * a_i_1d(ji) 885 917 886 918 ENDIF 887 919 ! … … 889 921 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 890 922 ! 891 ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2 923 ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2 892 924 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_Dt_ice * a_i_1d(ji) 893 925 ! … … 920 952 ! --- SIMIP diagnostics 921 953 ! 922 DO ji = 1, npti 954 DO ji = 1, npti 923 955 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 924 956 IF( h_s_1d(ji) >= zhs_ssl ) THEN 925 t_si_1d(ji) = ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji, 1) &926 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) &957 t_si_1d(ji) = ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,nlay_s) & 958 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) & 927 959 & / ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i & 928 960 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s )
Note: See TracChangeset
for help on using the changeset viewer.
