- Timestamp:
- 2018-10-15T12:22:07+02:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r10190_ice_thd_zdf/src/ICE/icethd_zdf_bl99.F90
r10069 r10192 108 108 ! 109 109 REAL(wp), DIMENSION(jpij ) :: ztsuold ! Old surface temperature in the ice 110 REAL(wp), DIMENSION( jpij,nlay_i) :: ztiold ! Old temperature in the ice111 REAL(wp), DIMENSION( jpij,nlay_s) :: ztsold ! Old temperature in the snow112 REAL(wp), DIMENSION( jpij,nlay_i) :: ztib ! Temporary temperature in the ice to check the convergence113 REAL(wp), DIMENSION( jpij,nlay_s) :: ztsb ! Temporary temperature in the snow to check the convergence114 REAL(wp), DIMENSION( jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity115 REAL(wp), DIMENSION(jpij, 0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice116 REAL(wp), DIMENSION( jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice117 REAL(wp), DIMENSION( jpij,0:nlay_i) :: zkappa_i ! Kappa factor in the ice118 REAL(wp), DIMENSION( jpij,0:nlay_i) :: zeta_i ! Eta factor in the ice110 REAL(wp), DIMENSION(nlay_i, jpij) :: ztiold ! Old temperature in the ice 111 REAL(wp), DIMENSION(nlay_s, jpij) :: ztsold ! Old temperature in the snow 112 REAL(wp), DIMENSION(nlay_i, jpij) :: ztib ! Temporary temperature in the ice to check the convergence 113 REAL(wp), DIMENSION(nlay_s, jpij) :: ztsb ! Temporary temperature in the snow to check the convergence 114 REAL(wp), DIMENSION(0:nlay_i, jpij) :: ztcond_i ! Ice thermal conductivity 115 REAL(wp), DIMENSION(jpij, 0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice 116 REAL(wp), DIMENSION(0:nlay_i, jpij) :: zradab_i ! Radiation absorbed in the ice 117 REAL(wp), DIMENSION(0:nlay_i, jpij) :: zkappa_i ! Kappa factor in the ice 118 REAL(wp), DIMENSION(0:nlay_i, jpij) :: zeta_i ! Eta factor in the ice 119 119 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow 120 REAL(wp), DIMENSION( jpij,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow121 REAL(wp), DIMENSION( jpij,0:nlay_s) :: zkappa_s ! Kappa factor in the snow122 REAL(wp), DIMENSION( jpij,0:nlay_s) :: zeta_s ! Eta factor in the snow123 REAL(wp), DIMENSION( jpij,nlay_i+3) :: zindterm ! 'Ind'ependent term124 REAL(wp), DIMENSION( jpij,nlay_i+3) :: zindtbis ! Temporary 'ind'ependent term125 REAL(wp), DIMENSION( jpij,nlay_i+3) :: zdiagbis ! Temporary 'dia'gonal term126 REAL(wp), DIMENSION( jpij,nlay_i+3,3):: ztrid ! Tridiagonal system terms120 REAL(wp), DIMENSION(0:nlay_s, jpij) :: zradab_s ! Radiation absorbed in the snow 121 REAL(wp), DIMENSION(0:nlay_s, jpij) :: zkappa_s ! Kappa factor in the snow 122 REAL(wp), DIMENSION(0:nlay_s, jpij) :: zeta_s ! Eta factor in the snow 123 REAL(wp), DIMENSION(nlay_i+3, jpij) :: zindterm ! 'Ind'ependent term 124 REAL(wp), DIMENSION(nlay_i+3, jpij) :: zindtbis ! Temporary 'ind'ependent term 125 REAL(wp), DIMENSION(nlay_i+3, jpij) :: zdiagbis ! Temporary 'dia'gonal term 126 REAL(wp), DIMENSION(3, nlay_i+3, jpij):: ztrid ! Tridiagonal system terms 127 127 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat 128 128 REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 129 REAL(wp), DIMENSION(nlay_i, jpij) :: tt_i_1d, tsz_i_1d 130 REAL(wp), DIMENSION(nlay_s, jpij) :: tt_s_1d 129 131 ! 130 132 ! Mono-category … … 132 134 REAL(wp) :: zhe ! dummy factor 133 135 REAL(wp) :: zcnd_i ! mean sea ice thermal conductivity 136 INTEGER :: itot 137 LOGICAL, DIMENSION(jpij) :: liter 134 138 !!------------------------------------------------------------------ 135 139 itot = npti 136 140 ! --- diag error on heat diffusion - PART 1 --- ! 137 141 DO ji = 1, npti … … 171 175 ENDIF 172 176 ! 173 ztsold (1:npti,:) = t_s_1d(1:npti,:) ! Old snow temperature 174 ztiold (1:npti,:) = t_i_1d(1:npti,:) ! Old ice temperature 177 DO jk = 1, nlay_i 178 DO ji = 1, npti 179 tt_i_1d(jk, ji) = t_i_1d(ji, jk) 180 tsz_i_1d(jk, ji) = sz_i_1d(ji, jk) 181 END DO 182 END DO 183 184 DO jk = 1, nlay_s 185 DO ji = 1, npti 186 tt_s_1d(jk, ji) = t_s_1d(ji, jk) 187 END DO 188 END DO 189 190 191 ztsold (:, :) = tt_s_1d(:,:) ! Old snow temperature 192 ztiold (:, :) = tt_i_1d(:,:) ! Old ice temperature 175 193 176 194 !------------- … … 184 202 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * h_s_1d(ji) * r1_nlay_s * REAL(jk) ) 185 203 ! ! radiation absorbed by the layer-th snow layer 186 zradab_s(j i,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)204 zradab_s(jk,ji) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 187 205 END DO 188 206 END DO 189 207 ! 190 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) )208 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) ) 191 209 DO jk = 1, nlay_i 192 210 DO ji = 1, npti 193 211 ! ! radiation transmitted below the layer-th ice layer 194 zradtr_i(ji, jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )212 zradtr_i(ji, jk) = zradtr_i(ji, 0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) 195 213 ! ! radiation absorbed by the layer-th ice layer 196 zradab_i(j i,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)214 zradab_i(jk,ji) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 197 215 END DO 198 216 END DO 199 217 ! 200 qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti, nlay_i) ! record radiation transmitted below the ice218 qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti, nlay_i) ! record radiation transmitted below the ice 201 219 ! 202 220 iconv = 0 ! number of iterations 203 221 zdti_max = 1000._wp ! maximal value of error on all points 204 222 ! 205 ! !============================! 206 DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max ) ! Iterative procedure begins ! 207 ! !============================! 223 liter(:) = .TRUE. !============================! 224 ! DO iconv = 1, iconv_max ! Iterative procedure begins ! 225 DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max) 226 zdti_max = 0._wp 227 ! IF(itot<1) EXIT 208 228 iconv = iconv + 1 209 ! 210 ztib(1:npti,:) = t_i_1d(1:npti,:) 211 ztsb(1:npti,:) = t_s_1d(1:npti,:) 212 ! 213 !-------------------------------- 214 ! 3) Sea ice thermal conductivity 215 !-------------------------------- 216 IF( ln_cndi_U64 ) THEN !-- Untersteiner (1964) formula: k = k0 + beta.S/T 217 ! 218 DO ji = 1, npti 219 ztcond_i(ji,0) = rcnd_i + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 220 ztcond_i(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) 221 END DO 222 DO jk = 1, nlay_i-1 223 DO ji = 1, npti 224 ztcond_i(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 225 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 229 DO ji = 1, npti 230 IF(liter(ji)) THEN 231 ! zdti_max = 0._wp 232 ! 233 ztib(:, ji) = tt_i_1d(:, ji) 234 ztsb(:, ji) = tt_s_1d(:, ji) 235 ! 236 !-------------------------------- 237 ! 3) Sea ice thermal conductivity 238 !-------------------------------- 239 IF( ln_cndi_U64 ) THEN !-- Untersteiner (1964) formula: k = k0 + beta.S/T 240 ! 241 ztcond_i(0, ji) = rcnd_i + zbeta * tsz_i_1d(1, ji) / MIN( -epsi10, tt_i_1d(1, ji) - rt0 ) 242 ztcond_i(nlay_i, ji) = rcnd_i + zbeta * tsz_i_1d(nlay_i, ji) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) 243 DO jk = 1, nlay_i-1 244 ztcond_i(jk, ji) = rcnd_i + zbeta * 0.5_wp * ( tsz_i_1d(jk, ji) + tsz_i_1d(jk+1, ji) ) / & 245 & MIN( -epsi10, 0.5_wp * (tt_i_1d(jk, ji) + tt_i_1d(jk+1, ji)) - rt0 ) 246 END DO 247 ! 248 ELSEIF( ln_cndi_P07 ) THEN !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T 249 ! 250 ztcond_i(0, ji) = rcnd_i + 0.09_wp * tsz_i_1d(1, ji) / MIN( -epsi10, tt_i_1d(1, ji) - rt0 ) & 251 & - 0.011_wp * ( tt_i_1d(1, ji) - rt0 ) 252 ztcond_i(nlay_i, ji) = rcnd_i + 0.09_wp * tsz_i_1d(nlay_i, ji) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 253 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 254 DO jk = 1, nlay_i-1 255 ztcond_i(jk,ji) = rcnd_i + 0.09_wp * 0.5_wp * ( tsz_i_1d(jk,ji) + tsz_i_1d(jk+1,ji) ) / & 256 & MIN( -epsi10, 0.5_wp * ( tt_i_1d (jk,ji) + tt_i_1d (jk+1,ji) ) - rt0 ) & 257 & - 0.011_wp * ( 0.5_wp * ( tt_i_1d (jk,ji) + tt_i_1d (jk+1,ji) ) - rt0 ) 226 258 END DO 227 END DO 228 ! 229 ELSEIF( ln_cndi_P07 ) THEN !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T 230 ! 231 DO ji = 1, npti 232 ztcond_i(ji,0) = rcnd_i + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 233 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 234 ztcond_i(ji,nlay_i) = rcnd_i + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 235 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 236 END DO 237 DO jk = 1, nlay_i-1 238 DO ji = 1, npti 239 ztcond_i(ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 240 & MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) & 241 & - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 242 END DO 243 END DO 244 ! 245 ENDIF 246 ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) ) 247 ! 248 !--- G(he) : enhancement of thermal conductivity in mono-category case 249 ! Computation of effective thermal conductivity G(h) 250 ! Used in mono-category case only to simulate an ITD implicitly 251 ! Fichefet and Morales Maqueda, JGR 1997 252 zghe(1:npti) = 1._wp 253 ! 254 SELECT CASE ( nn_virtual_itd ) 255 ! 256 CASE ( 1 , 2 ) 259 ! 260 ENDIF 261 ztcond_i(:, ji) = MAX( zkimin, ztcond_i(:, ji) ) 262 ! 263 !--- G(he) : enhancement of thermal conductivity in mono-category case 264 ! Computation of effective thermal conductivity G(h) 265 ! Used in mono-category case only to simulate an ITD implicitly 266 ! Fichefet and Morales Maqueda, JGR 1997 267 zghe(ji) = 1._wp 268 ! 269 SELECT CASE ( nn_virtual_itd ) 270 ! 271 CASE ( 1 , 2 ) 257 272 ! 258 273 zepsilon = 0.1_wp 259 DO ji = 1, npti 260 zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp ) ! Mean sea ice thermal conductivity 274 zcnd_i = SUM( ztcond_i(:,ji) ) / REAL( nlay_i+1, wp ) ! Mean sea ice thermal conductivity 261 275 zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i ) ! Effective thickness he (zhe) 262 276 IF( zhe >= zepsilon * 0.5_wp * EXP(1._wp) ) & 263 277 & zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) ) ! G(he) 278 ! 279 END SELECT 280 ! 281 !----------------- 282 ! 4) kappa factors 283 !----------------- 284 !--- Snow 285 DO jk = 0, nlay_s-1 286 zkappa_s(jk,ji) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 264 287 END DO 265 ! 266 END SELECT 267 ! 268 !----------------- 269 ! 4) kappa factors 270 !----------------- 271 !--- Snow 272 DO jk = 0, nlay_s-1 273 DO ji = 1, npti 274 zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 288 ! Snow-ice interface 289 zfac = 0.5_wp * ( ztcond_i(0, ji) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 290 IF( zfac > epsi10 ) THEN 291 zkappa_s(nlay_s, ji) = zghe(ji) * rn_cnd_s * ztcond_i(0, ji) / zfac 292 ELSE 293 zkappa_s(nlay_s, ji) = 0._wp 294 ENDIF 295 296 !--- Ice 297 DO jk = 0, nlay_i 298 zkappa_i(jk,ji) = zghe(ji) * ztcond_i(jk,ji) * z1_h_i(ji) 275 299 END DO 276 END DO 277 DO ji = 1, npti ! Snow-ice interface 278 zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 279 IF( zfac > epsi10 ) THEN 280 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 281 ELSE 282 zkappa_s(ji,nlay_s) = 0._wp 283 ENDIF 284 END DO 285 286 !--- Ice 287 DO jk = 0, nlay_i 288 DO ji = 1, npti 289 zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 300 zkappa_i(0, ji) = zkappa_s(nlay_s, ji) * isnow(ji) + zkappa_i(0, ji) * ( 1._wp - isnow(ji) ) 301 ! 302 !-------------------------------------- 303 ! 5) Sea ice specific heat, eta factors 304 !-------------------------------------- 305 DO jk = 1, nlay_i 306 zcpi = rcpi + zgamma * tsz_i_1d(jk,ji) / MAX( ( tt_i_1d(jk,ji) - rt0 ) * ( ztiold(jk,ji) - rt0 ), epsi10 ) 307 zeta_i(jk,ji) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi ) 290 308 END DO 291 END DO 292 DO ji = 1, npti ! Snow-ice interface 293 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 294 END DO 295 ! 296 !-------------------------------------- 297 ! 5) Sea ice specific heat, eta factors 298 !-------------------------------------- 299 DO jk = 1, nlay_i 300 DO ji = 1, npti 301 zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 302 zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi ) 309 310 DO jk = 1, nlay_s 311 zeta_s(jk,ji) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji) 303 312 END DO 304 END DO305 306 DO jk = 1, nlay_s307 DO ji = 1, npti308 zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji)309 END DO310 END DO311 313 ! 312 314 !----------------------------------------! … … 325 327 !---------------------------- 326 328 ! update of the non solar flux according to the update in T_su 327 DO ji = 1, npti328 329 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 329 END DO 330 331 DO ji = 1, npti 330 332 331 zfnet(ji) = qsr_ice_1d(ji) - qtr_ice_top_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar 333 END DO334 332 ! 335 333 !---------------------------- … … 342 340 343 341 ! ice interior terms (top equation has the same form as the others) 344 ztrid ( 1:npti,:,:) = 0._wp345 zindterm( 1:npti,:) = 0._wp346 zindtbis( 1:npti,:) = 0._wp347 zdiagbis( 1:npti,:) = 0._wp342 ztrid (:, :, ji) = 0._wp 343 zindterm(:, ji) = 0._wp 344 zindtbis(:, ji) = 0._wp 345 zdiagbis(:, ji) = 0._wp 348 346 349 347 DO jm = nlay_s + 2, nlay_s + nlay_i 350 DO ji = 1, npti351 348 jk = jm - nlay_s - 1 352 ztrid (ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 353 ztrid (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 354 ztrid (ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk) 355 zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 356 END DO 349 ztrid (1, jm, ji) = - zeta_i(jk,ji) * zkappa_i(jk-1, ji) 350 ztrid (2, jm, ji) = 1._wp + zeta_i(jk,ji) * ( zkappa_i(jk-1, ji) + zkappa_i(jk, ji) ) 351 ztrid (3, jm, ji) = - zeta_i(jk,ji) * zkappa_i(jk, ji) 352 zindterm(jm,ji) = ztiold(jk,ji) + zeta_i(jk,ji) * zradab_i(jk, ji) 357 353 END DO 358 354 359 355 jm = nlay_s + nlay_i + 1 360 DO ji = 1, npti361 356 ! ice bottom term 362 ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1) 363 ztrid (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 ) 364 ztrid (ji,jm,3) = 0._wp 365 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 366 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 367 END DO 368 369 DO ji = 1, npti 357 ztrid (1, jm, ji) = - zeta_i(nlay_i, ji) * zkappa_i(nlay_i-1, ji) 358 ztrid (2, jm, ji) = 1._wp + zeta_i(nlay_i, ji) * ( zkappa_i(nlay_i-1, ji) + zkappa_i(nlay_i, ji) * zg1 ) 359 ztrid (3, jm, ji) = 0._wp 360 zindterm(jm,ji) = ztiold(nlay_i, ji) + zeta_i(nlay_i, ji) * & 361 & ( zradab_i(nlay_i, ji) + zkappa_i(nlay_i, ji) * zg1 * t_bo_1d(ji) ) 362 370 363 ! !---------------------! 371 364 IF( h_s_1d(ji) > 0._wp ) THEN ! snow-covered cells ! … … 374 367 DO jm = 3, nlay_s + 1 375 368 jk = jm - 1 376 ztrid ( ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)377 ztrid ( ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )378 ztrid ( ji,jm,3) = - zeta_s(ji,jk) * zkappa_s(ji,jk)379 zindterm(j i,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)369 ztrid (1, jm, ji) = - zeta_s(jk,ji) * zkappa_s(jk-1,ji) 370 ztrid (2, jm, ji) = 1._wp + zeta_s(jk,ji) * ( zkappa_s(jk-1,ji) + zkappa_s(jk,ji) ) 371 ztrid (3, jm, ji) = - zeta_s(jk,ji) * zkappa_s(jk,ji) 372 zindterm(jm,ji) = ztsold(jk,ji) + zeta_s(jk,ji) * zradab_s(jk,ji) 380 373 END DO 381 374 382 375 ! case of only one layer in the ice (ice equation is altered) 383 376 IF( nlay_i == 1 ) THEN 384 ztrid ( ji,nlay_s+2,3) = 0._wp385 zindterm( ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji)377 ztrid (3, nlay_s+2, ji) = 0._wp 378 zindterm(nlay_s+2, ji) = zindterm(nlay_s+2, ji) + zeta_i(1, ji) * zkappa_i(1, ji) * t_bo_1d(ji) 386 379 ENDIF 387 380 … … 392 385 393 386 ! surface equation 394 ztrid ( ji,1,1) = 0._wp395 ztrid ( ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)396 ztrid ( ji,1,3) = zg1s * zkappa_s(ji,0)397 zindterm( ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)387 ztrid (1, 1, ji) = 0._wp 388 ztrid (2, 1, ji) = zdqns_ice_b(ji) - zg1s * zkappa_s(0, ji) 389 ztrid (3, 1, ji) = zg1s * zkappa_s(0, ji) 390 zindterm(1, ji) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 398 391 399 392 ! first layer of snow equation 400 ztrid ( ji,2,1) = - zeta_s(ji,1) * zkappa_s(ji,0) * zg1s401 ztrid ( ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )402 ztrid ( ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1)403 zindterm( ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)393 ztrid (1, 2, ji) = - zeta_s(1, ji) * zkappa_s(0, ji) * zg1s 394 ztrid (2, 2, ji) = 1._wp + zeta_s(1, ji) * ( zkappa_s(1, ji) + zkappa_s(0, ji) * zg1s ) 395 ztrid (3, 2, ji) = - zeta_s(1, ji) * zkappa_s(1, ji) 396 zindterm(2, ji) = ztsold(1, ji) + zeta_s(1, ji) * zradab_s(1, ji) 404 397 405 398 ELSE !-- case 2 : surface is melting … … 409 402 410 403 ! first layer of snow equation 411 ztrid ( ji,2,1) = 0._wp412 ztrid ( ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )413 ztrid ( ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1)414 zindterm( ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )404 ztrid (1, 2, ji) = 0._wp 405 ztrid (2, 2, ji) = 1._wp + zeta_s(1, ji) * ( zkappa_s(1, ji) + zkappa_s(0, ji) * zg1s ) 406 ztrid (3, 2, ji) = - zeta_s(1, ji) * zkappa_s(1, ji) 407 zindterm(2, ji) = ztsold(1, ji) + zeta_s(1, ji) * ( zradab_s(1, ji) + zkappa_s(0, ji) * zg1s * t_su_1d(ji) ) 415 408 ENDIF 416 409 ! !---------------------! … … 424 417 425 418 ! surface equation 426 ztrid ( ji,jm_min(ji),1) = 0._wp427 ztrid ( ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * zg1428 ztrid ( ji,jm_min(ji),3) = zkappa_i(ji,0) * zg1429 zindterm(j i,jm_min(ji)) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)419 ztrid (1, jm_min(ji), ji) = 0._wp 420 ztrid (2, jm_min(ji), ji) = zdqns_ice_b(ji) - zkappa_i(0, ji) * zg1 421 ztrid (3, jm_min(ji), ji) = zkappa_i(0, ji) * zg1 422 zindterm(jm_min(ji), ji) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 430 423 431 424 ! first layer of ice equation 432 ztrid ( ji,jm_min(ji)+1,1) = - zeta_i(ji,1) * zkappa_i(ji,0) * zg1433 ztrid ( ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )434 ztrid ( ji,jm_min(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1)435 zindterm(j i,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)425 ztrid (1, jm_min(ji)+1, ji) = - zeta_i(1, ji) * zkappa_i(0, ji) * zg1 426 ztrid (2, jm_min(ji)+1, ji) = 1._wp + zeta_i(1, ji) * ( zkappa_i(1, ji) + zkappa_i(0, ji) * zg1 ) 427 ztrid (3, jm_min(ji)+1, ji) = - zeta_i(1, ji) * zkappa_i(1, ji) 428 zindterm(jm_min(ji)+1, ji) = ztiold(1, ji) + zeta_i(1, ji) * zradab_i(1, ji) 436 429 437 430 ! case of only one layer in the ice (surface & ice equations are altered) 438 431 IF( nlay_i == 1 ) THEN 439 ztrid ( ji,jm_min(ji),1) = 0._wp440 ztrid ( ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2._wp441 ztrid ( ji,jm_min(ji),3) = zkappa_i(ji,0) * 2._wp442 ztrid ( ji,jm_min(ji)+1,1) = - zeta_i(ji,1) * zkappa_i(ji,0) * 2._wp443 ztrid ( ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )444 ztrid ( ji,jm_min(ji)+1,3) = 0._wp445 zindterm(j i,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji))432 ztrid (1, jm_min(ji), ji) = 0._wp 433 ztrid (2, jm_min(ji), ji) = zdqns_ice_b(ji) - zkappa_i(0, ji) * 2._wp 434 ztrid (3, jm_min(ji), ji) = zkappa_i(0, ji) * 2._wp 435 ztrid (1, jm_min(ji)+1, ji) = - zeta_i(1, ji) * zkappa_i(0, ji) * 2._wp 436 ztrid (2, jm_min(ji)+1, ji) = 1._wp + zeta_i(1, ji) * ( zkappa_i(0, ji) * 2._wp + zkappa_i(1, ji) ) 437 ztrid (3, jm_min(ji)+1, ji) = 0._wp 438 zindterm(jm_min(ji)+1, ji) = ztiold(1, ji) + zeta_i(1, ji) * (zradab_i(1, ji) + zkappa_i(1, ji) * t_bo_1d(ji)) 446 439 ENDIF 447 440 … … 452 445 453 446 ! first layer of ice equation 454 ztrid ( ji,jm_min(ji),1) = 0._wp455 ztrid ( ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )456 ztrid ( ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)457 zindterm(j i,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji))447 ztrid (1, jm_min(ji), ji) = 0._wp 448 ztrid (2, jm_min(ji), ji) = 1._wp + zeta_i(1, ji) * ( zkappa_i(1, ji) + zkappa_i(0, ji) * zg1 ) 449 ztrid (3, jm_min(ji), ji) = - zeta_i(1, ji) * zkappa_i(1, ji) 450 zindterm(jm_min(ji), ji) = ztiold(1, ji) + zeta_i(1, ji) * (zradab_i(1, ji) + zkappa_i(0, ji) * zg1 * t_su_1d(ji)) 458 451 459 452 ! case of only one layer in the ice (surface & ice equations are altered) 460 453 IF( nlay_i == 1 ) THEN 461 ztrid ( ji,jm_min(ji),1) = 0._wp462 ztrid ( ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )463 ztrid ( ji,jm_min(ji),3) = 0._wp464 zindterm(j i,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &465 & + t_su_1d(ji) * zeta_i( ji,1) * zkappa_i(ji,0) * 2._wp454 ztrid (1, jm_min(ji), ji) = 0._wp 455 ztrid (2, jm_min(ji), ji) = 1._wp + zeta_i(1, ji) * ( zkappa_i(0, ji) * 2._wp + zkappa_i(1, ji) ) 456 ztrid (3, jm_min(ji), ji) = 0._wp 457 zindterm(jm_min(ji), ji) = ztiold(1, ji) + zeta_i(1, ji) * ( zradab_i(1, ji) + zkappa_i(1, ji) * t_bo_1d(ji) ) & 458 & + t_su_1d(ji) * zeta_i(1, ji) * zkappa_i(0, ji) * 2._wp 466 459 ENDIF 467 460 … … 469 462 ENDIF 470 463 ! 471 zindtbis(j i,jm_min(ji)) = zindterm(ji,jm_min(ji))472 zdiagbis(j i,jm_min(ji)) = ztrid (ji,jm_min(ji),2)464 zindtbis(jm_min(ji), ji) = zindterm(jm_min(ji), ji) 465 zdiagbis(jm_min(ji), ji) = ztrid (2, jm_min(ji), ji) 473 466 ! 474 END DO475 467 ! 476 468 !------------------------------ … … 479 471 ! Solve the tridiagonal system with Gauss elimination method. 480 472 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 481 jm_maxt = 0 482 jm_mint = nlay_i+5 483 DO ji = 1, npti 484 jm_mint = MIN(jm_min(ji),jm_mint) 485 jm_maxt = MAX(jm_max(ji),jm_maxt) 486 END DO 487 488 DO jk = jm_mint+1, jm_maxt 489 DO ji = 1, npti 490 jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 491 zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 492 zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1) 473 474 475 476 477 478 DO jm = jm_min(ji)+1, jm_max(ji) 479 zdiagbis(jm,ji) = ztrid (2, jm, ji) - ztrid(1, jm, ji) * ztrid (3, jm-1, ji) / zdiagbis(jm-1,ji) 480 zindtbis(jm,ji) = zindterm(jm,ji ) - ztrid(1, jm, ji) * zindtbis(jm-1,ji ) / zdiagbis(jm-1,ji) 493 481 END DO 494 END DO495 482 496 483 ! ice temperatures 497 DO ji = 1, npti 498 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 499 END DO 500 501 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 502 DO ji = 1, npti 484 tt_i_1d(nlay_i, ji) = zindtbis(jm_max(ji), ji) / zdiagbis(jm_max(ji), ji) 485 486 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 503 487 jk = jm - nlay_s - 1 504 t _i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)488 tt_i_1d(jk, ji) = ( zindtbis(jm,ji) - ztrid(3, jm,ji) * tt_i_1d(jk+1, ji) ) / zdiagbis(jm,ji) 505 489 END DO 506 END DO 507 508 DO ji = 1, npti 490 509 491 ! snow temperatures 510 492 IF( h_s_1d(ji) > 0._wp ) THEN 511 t _s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)493 tt_s_1d(nlay_s, ji) = ( zindtbis(nlay_s+1, ji) - ztrid(3, nlay_s+1, ji) * tt_i_1d(1, ji) ) / zdiagbis(nlay_s+1, ji) 512 494 ENDIF 513 495 ! surface temperature 514 496 ztsub(ji) = t_su_1d(ji) 515 497 IF( t_su_1d(ji) < rt0 ) THEN 516 t_su_1d(ji) = ( zindtbis(j i,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * &517 & ( isnow(ji) * t _s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))498 t_su_1d(ji) = ( zindtbis(jm_min(ji), ji) - ztrid(3,jm_min(ji),ji) * & 499 & ( isnow(ji) * tt_s_1d(1, ji) + ( 1._wp - isnow(ji) ) * tt_i_1d(1, ji) ) ) / zdiagbis(jm_min(ji), ji) 518 500 ENDIF 519 END DO520 501 ! 521 502 !-------------------------------------------------------------- … … 524 505 ! check that nowhere it has started to melt 525 506 ! zdti_max is a measure of error, it has to be under zdti_bnd 526 zdti_max = 0._wp 527 DO ji = 1, npti 507 528 508 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 529 509 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 530 END DO 531 532 DO jk = 1, nlay_s 533 DO ji = 1, npti 534 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 535 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 510 511 DO jk = 1, nlay_s 512 tt_s_1d(jk,ji) = MAX( MIN( tt_s_1d(jk,ji), rt0 ), rt0 - 100._wp ) 513 zdti_max = MAX( zdti_max, ABS( tt_s_1d(jk,ji) - ztsb(jk,ji) ) ) 536 514 END DO 537 END DO 538 539 DO jk = 1, nlay_i 540 DO ji = 1, npti 541 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 542 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 543 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 515 516 DO jk = 1, nlay_i 517 ztmelts = -rTmlt * tsz_i_1d(jk, ji) + rt0 518 tt_i_1d(jk, ji) = MAX( MIN( tt_i_1d(jk, ji), ztmelts ), rt0 - 100._wp ) 519 zdti_max = MAX( zdti_max, ABS( tt_i_1d(jk, ji) - ztib(jk,ji) ) ) 544 520 END DO 545 END DO546 547 ! Compute spatial maximum over all errors548 ! note that this could be optimized substantially by iterating only the non-converging points549 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )550 521 ! 551 522 !----------------------------------------! … … 568 539 569 540 ! ice interior terms (top equation has the same form as the others) 570 ztrid (1:npti,:,:) = 0._wp 571 zindterm(1:npti,:) = 0._wp 572 zindtbis(1:npti,:) = 0._wp 573 zdiagbis(1:npti,:) = 0._wp 574 575 DO jm = nlay_s + 2, nlay_s + nlay_i 576 DO ji = 1, npti 541 ztrid (:, :, ji) = 0._wp 542 zindterm(:, ji) = 0._wp 543 zindtbis(:, ji) = 0._wp 544 zdiagbis(:, ji) = 0._wp 545 546 DO jm = nlay_s + 2, nlay_s + nlay_i 577 547 jk = jm - nlay_s - 1 578 ztrid (ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 579 ztrid (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 580 ztrid (ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk) 581 zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 582 END DO 583 ENDDO 584 585 jm = nlay_s + nlay_i + 1 586 DO ji = 1, npti 548 ztrid (1, jm, ji) = - zeta_i(jk,ji) * zkappa_i(jk-1, ji) 549 ztrid (2, jm, ji) = 1._wp + zeta_i(jk,ji) * ( zkappa_i(jk-1, ji) + zkappa_i(jk,ji) ) 550 ztrid (3, jm, ji) = - zeta_i(jk,ji) * zkappa_i(jk,ji) 551 zindterm(jm,ji) = ztiold(jk,ji) + zeta_i(jk,ji) * zradab_i(jk,ji) 552 ENDDO 553 554 jm = nlay_s + nlay_i + 1 587 555 ! ice bottom term 588 ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1) 589 ztrid (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 ) 590 ztrid (ji,jm,3) = 0._wp 591 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 592 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 593 ENDDO 594 595 DO ji = 1, npti 556 ztrid (1, jm, ji) = - zeta_i(nlay_i, ji) * zkappa_i(nlay_i-1, ji) 557 ztrid (2, jm, ji) = 1._wp + zeta_i(nlay_i, ji) * ( zkappa_i(nlay_i-1, ji) + zkappa_i(nlay_i, ji) * zg1 ) 558 ztrid (3, jm, ji) = 0._wp 559 zindterm(jm,ji) = ztiold(nlay_i, ji) + zeta_i(nlay_i, ji) * & 560 & ( zradab_i(nlay_i, ji) + zkappa_i(nlay_i, ji) * zg1 * t_bo_1d(ji) ) 561 596 562 ! !---------------------! 597 563 IF( h_s_1d(ji) > 0._wp ) THEN ! snow-covered cells ! … … 600 566 DO jm = 3, nlay_s + 1 601 567 jk = jm - 1 602 ztrid ( ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)603 ztrid ( ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )604 ztrid ( ji,jm,3) = - zeta_s(ji,jk) * zkappa_s(ji,jk)605 zindterm(j i,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)568 ztrid (1, jm, ji) = - zeta_s(jk,ji) * zkappa_s(jk-1,ji) 569 ztrid (2, jm, ji) = 1._wp + zeta_s(jk,ji) * ( zkappa_s(jk-1,ji) + zkappa_s(jk,ji) ) 570 ztrid (3, jm, ji) = - zeta_s(jk,ji) * zkappa_s(jk,ji) 571 zindterm(jm,ji) = ztsold(jk,ji) + zeta_s(jk,ji) * zradab_s(jk,ji) 606 572 END DO 607 573 608 574 ! case of only one layer in the ice (ice equation is altered) 609 575 IF ( nlay_i == 1 ) THEN 610 ztrid ( ji,nlay_s+2,3) = 0._wp611 zindterm( ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji)576 ztrid (3, nlay_s+2, ji) = 0._wp 577 zindterm(nlay_s+2, ji) = zindterm(nlay_s+2, ji) + zeta_i(1, ji) * zkappa_i(1, ji) * t_bo_1d(ji) 612 578 ENDIF 613 579 … … 616 582 617 583 ! first layer of snow equation 618 ztrid ( ji,2,1) = 0._wp619 ztrid ( ji,2,2) = 1._wp + zeta_s(ji,1) * zkappa_s(ji,1)620 ztrid ( ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1)621 zindterm( ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) )584 ztrid (1, 2, ji) = 0._wp 585 ztrid (2, 2, ji) = 1._wp + zeta_s(1, ji) * zkappa_s(1, ji) 586 ztrid (3, 2, ji) = - zeta_s(1, ji) * zkappa_s(1, ji) 587 zindterm(2, ji) = ztsold(1, ji) + zeta_s(1, ji) * ( zradab_s(1, ji) + qcn_ice_1d(ji) ) 622 588 623 589 ! !---------------------! … … 628 594 629 595 ! first layer of ice equation 630 ztrid ( ji,jm_min(ji),1) = 0._wp631 ztrid ( ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)632 ztrid ( ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)633 zindterm(j i,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) )596 ztrid (1, jm_min(ji), ji) = 0._wp 597 ztrid (2, jm_min(ji), ji) = 1._wp + zeta_i(1, ji) * zkappa_i(1, ji) 598 ztrid (3, jm_min(ji), ji) = - zeta_i(1, ji) * zkappa_i(1, ji) 599 zindterm(jm_min(ji), ji) = ztiold(1, ji) + zeta_i(1, ji) * ( zradab_i(1, ji) + qcn_ice_1d(ji) ) 634 600 635 601 ! case of only one layer in the ice (surface & ice equations are altered) 636 602 IF( nlay_i == 1 ) THEN 637 ztrid ( ji,jm_min(ji),1) = 0._wp638 ztrid ( ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)639 ztrid ( ji,jm_min(ji),3) = 0._wp640 zindterm(j i,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * &641 & ( zradab_i( ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) + qcn_ice_1d(ji) )603 ztrid (1, jm_min(ji), ji) = 0._wp 604 ztrid (2, jm_min(ji), ji) = 1._wp + zeta_i(1, ji) * zkappa_i(1, ji) 605 ztrid (3, jm_min(ji), ji) = 0._wp 606 zindterm(jm_min(ji), ji) = ztiold(1, ji) + zeta_i(1, ji) * & 607 & ( zradab_i(1, ji) + zkappa_i(1, ji) * t_bo_1d(ji) + qcn_ice_1d(ji) ) 642 608 ENDIF 643 609 644 610 ENDIF 645 611 ! 646 zindtbis(j i,jm_min(ji)) = zindterm(ji,jm_min(ji))647 zdiagbis(j i,jm_min(ji)) = ztrid (ji,jm_min(ji),2)612 zindtbis(jm_min(ji), ji) = zindterm(jm_min(ji), ji) 613 zdiagbis(jm_min(ji), ji) = ztrid (2, jm_min(ji), ji) 648 614 ! 649 END DO650 615 ! 651 616 !------------------------------ … … 654 619 ! Solve the tridiagonal system with Gauss elimination method. 655 620 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 656 jm_maxt = 0 657 jm_mint = nlay_i+5 658 DO ji = 1, npti 659 jm_mint = MIN(jm_min(ji),jm_mint) 660 jm_maxt = MAX(jm_max(ji),jm_maxt) 661 END DO 662 663 DO jk = jm_mint+1, jm_maxt 664 DO ji = 1, npti 665 jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 666 zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 667 zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 668 END DO 621 622 623 624 625 626 DO jm = jm_min(ji)+1, jm_max(ji) 627 zdiagbis(jm,ji) = ztrid (2, jm, ji) - ztrid(1, jm, ji) * ztrid (3, jm-1, ji) / zdiagbis(jm-1,ji) 628 zindtbis(jm,ji) = zindterm(jm,ji) - ztrid(1, jm,ji) * zindtbis(jm-1,ji) / zdiagbis(jm-1,ji) 669 629 END DO 670 630 671 631 ! ice temperatures 672 DO ji = 1, npti 673 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 674 END DO 675 676 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 677 DO ji = 1, npti 632 tt_i_1d(nlay_i, ji) = zindtbis(jm_max(ji), ji) / zdiagbis(jm_max(ji), ji) 633 634 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 678 635 jk = jm - nlay_s - 1 679 t _i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)636 tt_i_1d(jk, ji) = ( zindtbis(jm,ji) - ztrid(3, jm, ji) * tt_i_1d(jk+1, ji) ) / zdiagbis(jm,ji) 680 637 END DO 681 END DO682 638 683 639 ! snow temperatures 684 DO ji = 1, npti685 640 IF( h_s_1d(ji) > 0._wp ) THEN 686 t _s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)641 tt_s_1d(nlay_s, ji) = ( zindtbis(nlay_s+1, ji) - ztrid(3, nlay_s+1, ji) * tt_i_1d(1, ji) ) / zdiagbis(nlay_s+1, ji) 687 642 ENDIF 688 END DO689 643 ! 690 644 !-------------------------------------------------------------- … … 693 647 ! check that nowhere it has started to melt 694 648 ! zdti_max is a measure of error, it has to be under zdti_bnd 695 zdti_max = 0._wp 696 697 DO jk = 1, nlay_s 698 DO ji = 1, npti 699 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 700 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 649 650 651 DO jk = 1, nlay_s 652 tt_s_1d(jk,ji) = MAX( MIN( tt_s_1d(jk,ji), rt0 ), rt0 - 100._wp ) 653 zdti_max = MAX( zdti_max, ABS( tt_s_1d(jk,ji) - ztsb(jk,ji) ) ) 701 654 END DO 702 END DO703 655 704 DO jk = 1, nlay_i 705 DO ji = 1, npti 706 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 707 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 708 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 656 DO jk = 1, nlay_i 657 ztmelts = -rTmlt * tsz_i_1d(jk,ji) + rt0 658 tt_i_1d(jk,ji) = MAX( MIN( tt_i_1d(jk,ji), ztmelts ), rt0 - 100._wp ) 659 zdti_max = MAX( zdti_max, ABS( tt_i_1d(jk,ji) - ztib(jk,ji) ) ) 709 660 END DO 710 END DO 711 712 ! Compute spatial maximum over all errors 713 ! note that this could be optimized substantially by iterating only the non-converging points 714 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 715 716 ENDIF ! k_jules 717 661 662 ENDIF ! k_jules 663 664 ! IF(zdti_max.le.zdti_bnd) THEN 665 ! liter(ji) = .FALSE. 666 ! itot = itot - 1 667 ! ENDIF 668 669 ENDIF ! liter 670 END DO ! ji 671 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 718 672 END DO ! End of the do while iterative procedure 719 720 721 722 723 673 !need to move this part where all processors are available 674 ! IF( ln_icectl .AND. lwp ) THEN 675 ! WRITE(numout,*) ' zdti_max : ', zdti_max 676 ! WRITE(numout,*) ' iconv : ', iconv 677 ! ENDIF 724 678 725 679 ! … … 729 683 ! 730 684 ! --- calculate conduction fluxes (positive downward) 685 DO jk = 1, nlay_s 686 t_s_1d(1:npti, jk) = tt_s_1d(jk, 1:npti) 687 ENDDO 688 DO jk = 1, nlay_i 689 t_i_1d(1:npti, jk) = tt_i_1d(jk, 1:npti) 690 ENDDO 731 691 732 692 DO ji = 1, npti 733 693 ! ! surface ice conduction flux 734 qcn_ice_top_1d(ji) = - isnow(ji) * zkappa_s( ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) &735 & - ( 1._wp - isnow(ji) ) * zkappa_i( ji,0) * zg1 * ( t_i_1d(ji,1) - t_su_1d(ji) )694 qcn_ice_top_1d(ji) = - isnow(ji) * zkappa_s(0, ji) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) & 695 & - ( 1._wp - isnow(ji) ) * zkappa_i(0, ji) * zg1 * ( t_i_1d(ji,1) - t_su_1d(ji) ) 736 696 ! ! bottom ice conduction flux 737 qcn_ice_bot_1d(ji) = - zkappa_i( ji,nlay_i) * zg1 * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) )697 qcn_ice_bot_1d(ji) = - zkappa_i(nlay_i, ji) * zg1 * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 738 698 END DO 739 699 … … 770 730 771 731 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 772 zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji, nlay_i) - qcn_ice_bot_1d(ji) &732 zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji, nlay_i) - qcn_ice_bot_1d(ji) & 773 733 & + zdq * r1_rdtice ) * a_i_1d(ji) 774 734 ELSE ! case T_su = 0degC 775 zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji, nlay_i) - qcn_ice_bot_1d(ji) &735 zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji, nlay_i) - qcn_ice_bot_1d(ji) & 776 736 & + zdq * r1_rdtice ) * a_i_1d(ji) 777 737 ENDIF … … 779 739 ELSEIF( k_jules == np_jules_ACTIVE ) THEN 780 740 781 zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji, nlay_i) - qcn_ice_bot_1d(ji) &741 zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji, nlay_i) - qcn_ice_bot_1d(ji) & 782 742 & + zdq * r1_rdtice ) * a_i_1d(ji) 783 743 … … 800 760 DO ji = 1, npti 801 761 IF( h_s_1d(ji) > 0.1_wp ) THEN 802 cnd_ice_1d(ji) = 2._wp * zkappa_s( ji,0)762 cnd_ice_1d(ji) = 2._wp * zkappa_s(0, ji) 803 763 ELSE 804 764 IF( h_i_1d(ji) > 0.1_wp ) THEN 805 cnd_ice_1d(ji) = 2._wp * zkappa_i( ji,0)765 cnd_ice_1d(ji) = 2._wp * zkappa_i(0, ji) 806 766 ELSE 807 cnd_ice_1d(ji) = 2._wp * ztcond_i( ji,0) * 10._wp767 cnd_ice_1d(ji) = 2._wp * ztcond_i(0, ji) * 10._wp 808 768 ENDIF 809 769 ENDIF … … 813 773 IF( k_jules == np_jules_EMULE ) THEN 814 774 ! Restore temperatures to their initial values 815 t_s_1d (1:npti,:) = ztsold (1:npti,:) 816 t_i_1d (1:npti,:) = ztiold (1:npti,:) 775 DO jk = 1, nlay_s 776 t_s_1d (1:npti,jk) = ztsold (jk, 1:npti) 777 ENDDO 778 DO jk = 1, nlay_i 779 t_i_1d (1:npti,jk) = ztiold (jk, 1:npti) 780 ENDDO 817 781 qcn_ice_1d(1:npti) = qcn_ice_top_1d(1:npti) 818 782 ENDIF … … 822 786 DO ji = 1, npti 823 787 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 824 zfac = rn_cnd_s * zh_i(ji) + ztcond_i( ji,1) * zh_s(ji)788 zfac = rn_cnd_s * zh_i(ji) + ztcond_i(1, ji) * zh_s(ji) 825 789 IF( h_s_1d(ji) >= zhs_min ) THEN 826 790 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + & 827 & ztcond_i( ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac )791 & ztcond_i(1, ji) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac ) 828 792 ELSE 829 793 t_si_1d(ji) = t_su_1d(ji)
Note: See TracChangeset
for help on using the changeset viewer.
