Changeset 14037 for NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icethd_zdf_bl99.F90
- Timestamp:
- 2020-12-03T12:20:38+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icethd_zdf_bl99.F90 (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette @13292sette10 ^/utils/CI/sette_wave@13990 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icethd_zdf_bl99.F90
r12489 r14037 85 85 86 86 LOGICAL, DIMENSION(jpij) :: l_T_converged ! true when T converges (per grid point) 87 !87 ! 88 88 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system 89 89 REAL(wp) :: zg1 = 2._wp ! 90 90 REAL(wp) :: zgamma = 18009._wp ! for specific heat 91 91 REAL(wp) :: zbeta = 0.117_wp ! for thermal conductivity (could be 0.13) 92 REAL(wp) :: zraext_s = 10._wp ! extinction coefficient of radiation in the snow93 92 REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity 94 93 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C 95 94 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature 96 REAL(wp) :: zhs_min = 0.01_wp ! minimum snow thickness for conductivity calculation 95 REAL(wp) :: zhs_ssl = 0.03_wp ! surface scattering layer in the snow 96 REAL(wp) :: zhi_ssl = 0.10_wp ! surface scattering layer in the ice 97 REAL(wp) :: zh_min = 1.e-3_wp ! minimum ice/snow thickness for conduction 97 98 REAL(wp) :: ztmelts ! ice melting temperature 98 99 REAL(wp) :: zdti_max ! current maximal error on temperature 99 100 REAL(wp) :: zcpi ! Ice specific heat 100 101 REAL(wp) :: zhfx_err, zdq ! diag errors on heat 101 REAL(wp) :: zfac ! dummy factor 102 ! 103 REAL(wp), DIMENSION(jpij) :: isnow ! switch for presence (1) or absence (0) of snow 102 ! 103 REAL(wp), DIMENSION(jpij) :: zraext_s ! extinction coefficient of radiation in the snow 104 104 REAL(wp), DIMENSION(jpij) :: ztsub ! surface temperature at previous iteration 105 105 REAL(wp), DIMENSION(jpij) :: zh_i, z1_h_i ! ice layer thickness … … 109 109 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function 110 110 ! 111 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,nlay_i+3) :: zindterm ! 'Ind'ependent term 127 REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindtbis ! Temporary 'ind'ependent term 128 REAL(wp), DIMENSION(jpij,nlay_i+3) :: zdiagbis ! Temporary 'dia'gonal term 129 REAL(wp), DIMENSION(jpij,nlay_i+3,3) :: ztrid ! Tridiagonal system terms 130 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat 131 REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 111 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 132 136 ! 133 137 ! Mono-category … … 143 147 END DO 144 148 149 ! calculate ice fraction covered by snow for radiation 150 CALL ice_var_snwfra( h_s_1d(1:npti), za_s_fra(1:npti) ) 151 145 152 !------------------ 146 153 ! 1) Initialization 147 154 !------------------ 155 ! 156 ! extinction radiation in the snow 157 IF ( nn_qtrice == 0 ) THEN ! constant 158 zraext_s(1:npti) = rn_kappa_s 159 ELSEIF( nn_qtrice == 1 ) THEN ! depends on melting/freezing conditions 160 WHERE( t_su_1d(1:npti) < rt0 ) ; zraext_s(1:npti) = rn_kappa_sdry ! no surface melting 161 ELSEWHERE ; zraext_s(1:npti) = rn_kappa_smlt ! surface melting 162 END WHERE 163 ENDIF 164 ! 165 ! thicknesses 148 166 DO ji = 1, npti 149 isnow(ji) = 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) ) ! is there snow or not 150 ! layer thickness 151 zh_i(ji) = h_i_1d(ji) * r1_nlay_i 152 zh_s(ji) = h_s_1d(ji) * r1_nlay_s 167 ! ice thickness 168 IF( h_i_1d(ji) > 0._wp ) THEN 169 zh_i (ji) = MAX( zh_min , h_i_1d(ji) ) * r1_nlay_i ! set a minimum thickness for conduction 170 z1_h_i(ji) = 1._wp / zh_i(ji) ! it must be very small 171 ELSE 172 zh_i (ji) = 0._wp 173 z1_h_i(ji) = 0._wp 174 ENDIF 175 ! snow thickness 176 IF( h_s_1d(ji) > 0._wp ) THEN 177 zh_s (ji) = MAX( zh_min , h_s_1d(ji) ) * r1_nlay_s ! set a minimum thickness for conduction 178 z1_h_s(ji) = 1._wp / zh_s(ji) ! it must be very small 179 isnow (ji) = 1._wp 180 ELSE 181 zh_s (ji) = 0._wp 182 z1_h_s(ji) = 0._wp 183 isnow (ji) = 0._wp 184 ENDIF 185 ! for Met-Office 186 IF( h_s_1d(ji) < zh_min ) THEN 187 isnow_comb(ji) = h_s_1d(ji) / zh_min 188 ELSE 189 isnow_comb(ji) = 1._wp 190 ENDIF 153 191 END DO 154 ! 155 WHERE( zh_i(1:npti) >= epsi10 ) ; z1_h_i(1:npti) = 1._wp / zh_i(1:npti) 156 ELSEWHERE ; z1_h_i(1:npti) = 0._wp 157 END WHERE 158 ! 159 WHERE( zh_s(1:npti) > 0._wp ) zh_s(1:npti) = MAX( zhs_min * r1_nlay_s, zh_s(1:npti) ) 160 ! 161 WHERE( zh_s(1:npti) > 0._wp ) ; z1_h_s(1:npti) = 1._wp / zh_s(1:npti) 162 ELSEWHERE ; z1_h_s(1:npti) = 0._wp 163 END WHERE 192 ! clem: we should apply correction on snow thickness to take into account snow fraction 193 ! it must be a distribution, so it is a bit complicated 164 194 ! 165 195 ! Store initial temperatures and non solar heat fluxes 166 196 IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN 167 !168 197 ztsub (1:npti) = t_su_1d(1:npti) ! surface temperature at iteration n-1 169 198 ztsuold (1:npti) = t_su_1d(1:npti) ! surface temperature initial value … … 185 214 DO ji = 1, npti 186 215 ! ! radiation transmitted below the layer-th snow layer 187 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * h_s_1d(ji) * r1_nlay_s * REAL(jk) )216 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s(ji) * MAX( 0._wp, zh_s(ji) * REAL(jk) - zhs_ssl ) ) 188 217 ! ! radiation absorbed by the layer-th snow layer 189 218 zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) … … 191 220 END DO 192 221 ! 193 zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - isnow(1:npti) )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) ) 194 223 DO jk = 1, nlay_i 195 224 DO ji = 1, npti 196 225 ! ! radiation transmitted below the layer-th ice layer 197 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) 226 zradtr_i(ji,jk) = za_s_fra(ji) * zradtr_s(ji,nlay_s) & ! part covered by snow 227 & * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zh_min ) ) & 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 ) ) 198 230 ! ! radiation absorbed by the layer-th ice layer 199 231 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) … … 203 235 qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti,nlay_i) ! record radiation transmitted below the ice 204 236 ! 205 iconv 237 iconv = 0 ! number of iterations 206 238 ! 207 239 l_T_converged(:) = .FALSE. … … 230 262 DO ji = 1, npti 231 263 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 )264 & MIN( -epsi10, 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 ) 233 265 END DO 234 266 END DO … … 238 270 DO ji = 1, npti 239 271 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 )272 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 241 273 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 )274 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 243 275 END DO 244 276 DO jk = 1, nlay_i-1 245 277 DO ji = 1, npti 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 )278 ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 279 & MIN( -epsi10, 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 ) & 280 & - 0.011_wp * ( 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 ) 249 281 END DO 250 282 END DO … … 290 322 END DO 291 323 DO ji = 1, npti ! Snow-ice interface 292 IF ( .NOT. l_T_converged(ji) ) THEN 293 zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 294 IF( zfac > epsi10 ) THEN 295 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 296 ELSE 297 zkappa_s(ji,nlay_s) = 0._wp 298 ENDIF 299 ENDIF 324 IF ( .NOT. l_T_converged(ji) ) & 325 zkappa_s(ji,nlay_s) = isnow(ji) * zghe(ji) * rn_cnd_s * ztcond_i(ji,0) & 326 & / ( 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) ) 300 327 END DO 301 328 … … 310 337 END DO 311 338 DO ji = 1, npti ! Snow-ice interface 312 IF ( .NOT. l_T_converged(ji) ) & 313 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 339 IF ( .NOT. l_T_converged(ji) ) THEN 340 ! Calculate combined surface snow and ice conductivity to pass through the coupler (met-office) 341 zkappa_comb(ji) = isnow_comb(ji) * zkappa_s(ji,0) + ( 1._wp - isnow_comb(ji) ) * zkappa_i(ji,0) 342 ! If there is snow then use the same snow-ice interface conductivity for the top layer of ice 343 IF( h_s_1d(ji) > 0._wp ) zkappa_i(ji,0) = zkappa_s(ji,nlay_s) 344 ENDIF 314 345 END DO 315 346 ! … … 320 351 DO ji = 1, npti 321 352 zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 322 zeta_i(ji,jk) = rDt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi )353 zeta_i(ji,jk) = rDt_ice * r1_rhoi * z1_h_i(ji) / zcpi 323 354 END DO 324 355 END DO … … 502 533 ! Solve the tridiagonal system with Gauss elimination method. 503 534 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 504 jm_maxt = 0 505 jm_mint = nlay_i+5 506 DO ji = 1, npti 507 jm_mint = MIN(jm_min(ji),jm_mint) 508 jm_maxt = MAX(jm_max(ji),jm_maxt) 509 END DO 510 511 DO jk = jm_mint+1, jm_maxt 512 DO ji = 1, npti 513 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 514 555 zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 515 556 zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1) … … 533 574 END DO 534 575 576 ! snow temperatures 535 577 DO ji = 1, npti 536 578 ! Variables used after iterations 537 579 ! Value must be frozen after convergence for MPP independance reason 538 IF ( .NOT. l_T_converged(ji) ) THEN 539 ! snow temperatures 540 IF( h_s_1d(ji) > 0._wp ) THEN 541 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) 542 ENDIF 543 ! 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 544 595 ztsub(ji) = t_su_1d(ji) 545 596 IF( t_su_1d(ji) < rt0 ) THEN … … 561 612 562 613 IF ( .NOT. l_T_converged(ji) ) THEN 614 563 615 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 564 616 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 565 617 566 t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 567 zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 618 IF( h_s_1d(ji) > 0._wp ) THEN 619 DO jk = 1, nlay_s 620 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 621 zdti_max = MAX ( zdti_max , ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 622 END DO 623 ENDIF 568 624 569 625 DO jk = 1, nlay_i … … 572 628 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 573 629 END DO 574 575 IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE. 630 631 ! convergence test 632 IF( ln_zdf_chkcvg ) THEN 633 tice_cvgerr_1d(ji) = zdti_max 634 tice_cvgstp_1d(ji) = REAL(iconv) 635 ENDIF 636 637 IF( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE. 576 638 577 639 ENDIF … … 684 746 ! Solve the tridiagonal system with Gauss elimination method. 685 747 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 686 jm_maxt = 0 687 jm_mint = nlay_i+5 688 DO ji = 1, npti 689 jm_mint = MIN(jm_min(ji),jm_mint) 690 jm_maxt = MAX(jm_max(ji),jm_maxt) 691 END DO 692 693 DO jk = jm_mint+1, jm_maxt 694 DO ji = 1, npti 695 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 696 767 zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 697 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) 698 769 END DO 699 770 END DO 700 771 701 772 ! ice temperatures 702 773 DO ji = 1, npti … … 718 789 ! snow temperatures 719 790 DO ji = 1, npti 720 ! Variable used after iterations791 ! Variables used after iterations 721 792 ! Value must be frozen after convergence for MPP independance reason 722 IF ( .NOT. l_T_converged(ji) ) THEN 723 IF( h_s_1d(ji) > 0._wp ) THEN 724 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) 725 ENDIF 726 ENDIF 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 727 803 END DO 728 804 ! … … 738 814 739 815 IF ( .NOT. l_T_converged(ji) ) THEN 740 ! t_s 741 t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 742 zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 743 ! t_i 816 817 IF( h_s_1d(ji) > 0._wp ) THEN 818 DO jk = 1, nlay_s 819 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 820 zdti_max = MAX ( zdti_max , ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 821 END DO 822 ENDIF 823 744 824 DO jk = 1, nlay_i 745 825 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 … … 748 828 END DO 749 829 750 IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE. 830 ! convergence test 831 IF( ln_zdf_chkcvg ) THEN 832 tice_cvgerr_1d(ji) = zdti_max 833 tice_cvgstp_1d(ji) = REAL(iconv) 834 ENDIF 835 836 IF( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE. 751 837 752 838 ENDIF … … 755 841 756 842 ENDIF ! k_cnd 757 843 758 844 END DO ! End of the do while iterative procedure 759 760 IF( ln_icectl .AND. lwp ) THEN761 WRITE(numout,*) ' zdti_max : ', zdti_max762 WRITE(numout,*) ' iconv : ', iconv763 ENDIF764 765 845 ! 766 846 !----------------------------- … … 771 851 ! bottom ice conduction flux 772 852 DO ji = 1, npti 773 qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1 853 qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1 * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 774 854 END DO 775 855 ! surface ice conduction flux … … 777 857 ! 778 858 DO ji = 1, npti 779 qcn_ice_top_1d(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )&780 & 859 qcn_ice_top_1d(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) & 860 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * ( t_i_1d(ji,1) - t_su_1d(ji) ) 781 861 END DO 782 862 ! … … 792 872 ! 793 873 DO ji = 1, npti 794 t_su_1d(ji) = ( qcn_ice_top_1d(ji) & ! calculate surface temperature 795 & + isnow(ji) * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) & 796 & + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * t_i_1d(ji,1) & 797 & ) / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 ) 874 t_su_1d(ji) = ( qcn_ice_top_1d(ji) + isnow(ji) * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) + & 875 & ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * t_i_1d(ji,1) ) & 876 & / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 ) 798 877 t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp ) ! cap t_su 799 878 END DO … … 853 932 !-------------------------------------------------------------------- 854 933 ! effective conductivity and 1st layer temperature (needed by Met Office) 934 ! this is a conductivity at mid-layer, hence the factor 2 855 935 DO ji = 1, npti 856 IF( h_s_1d(ji) > 0.1_wp ) THEN 857 cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0) 936 IF( h_i_1d(ji) >= zhi_ssl ) THEN 937 cnd_ice_1d(ji) = 2._wp * zkappa_comb(ji) 938 !!cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0) 858 939 ELSE 859 IF( h_i_1d(ji) > 0.1_wp ) THEN 860 cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0) 861 ELSE 862 cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp 863 ENDIF 940 cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / zhi_ssl ! cnd_ice is capped by: cond_i/zhi_ssl 864 941 ENDIF 865 942 t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) … … 877 954 DO ji = 1, npti 878 955 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 879 zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 880 IF( h_s_1d(ji) >= zhs_min ) THEN 881 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + & 882 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac ) 956 IF( h_s_1d(ji) >= zhs_ssl ) THEN 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) ) & 959 & / ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i & 960 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s ) 883 961 ELSE 884 962 t_si_1d(ji) = t_su_1d(ji)
Note: See TracChangeset
for help on using the changeset viewer.
