- Timestamp:
- 2017-09-25T21:11:19+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
r8534 r8562 84 84 !!------------------------------------------------------------------- 85 85 INTEGER :: ji, jk ! spatial loop index 86 INTEGER :: numeq! current reference number of equation87 INTEGER :: minnumeqmin, maxnumeqmax86 INTEGER :: jm ! current reference number of equation 87 INTEGER :: jm_mint, jm_maxt 88 88 INTEGER :: iconv ! number of iterations in iterative procedure 89 89 INTEGER :: iconv_max = 50 ! max number of iterations in iterative procedure 90 90 91 INTEGER, DIMENSION(jpij) :: numeqmin! reference number of top equation92 INTEGER, DIMENSION(jpij) :: numeqmax! reference number of bottom equation91 INTEGER, DIMENSION(jpij) :: jm_min ! reference number of top equation 92 INTEGER, DIMENSION(jpij) :: jm_max ! reference number of bottom equation 93 93 94 94 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system … … 113 113 REAL(wp), DIMENSION(jpij) :: zfsw ! solar radiation absorbed at the surface 114 114 REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface 115 REAL(wp), DIMENSION(jpij) :: zf 115 REAL(wp), DIMENSION(jpij) :: zfnet ! surface flux function 116 116 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function 117 117 REAL(wp), DIMENSION(jpij) :: zftrice ! solar radiation transmitted through the ice … … 176 176 ! 2) Radiation 177 177 !------------- 178 !179 178 z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0 180 179 DO ji = 1, nidx … … 224 223 iconv = 0 ! number of iterations 225 224 zdti_max = 1000._wp ! maximal value of error on all points 226 ! ! ----------------------------!225 ! !============================! 227 226 DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max ) ! Iterative procedure begins ! 228 ! !----------------------------! 229 ! 227 ! !============================! 230 228 iconv = iconv + 1 231 229 ! … … 272 270 ! Used in mono-category case only to simulate an ITD implicitly 273 271 ! Fichefet and Morales Maqueda, JGR 1997 274 275 272 zghe(1:nidx) = 1._wp 276 273 ! 277 274 SELECT CASE ( nn_monocat ) 278 275 … … 281 278 zepsilon = 0.1_wp 282 279 DO ji = 1, nidx 283 284 ! Mean sea ice thermal conductivity 285 zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp ) 286 287 ! Effective thickness he (zhe) 288 zhe = ( rn_cnd_s * ht_i_1d(ji) + zcnd_i * ht_s_1d(ji) ) / ( rn_cnd_s + zcnd_i ) 289 290 ! G(he) 280 zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp ) ! Mean sea ice thermal conductivity 281 zhe = ( rn_cnd_s * ht_i_1d(ji) + zcnd_i * ht_s_1d(ji) ) / ( rn_cnd_s + zcnd_i ) ! Effective thickness he (zhe) 291 282 IF( zhe >= zepsilon * 0.5_wp * EXP(1._wp) ) THEN 292 zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ))283 zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) ) ! G(he) 293 284 ENDIF 294 295 285 END DO 296 286 … … 303 293 DO jk = 0, nlay_s-1 304 294 DO ji = 1, nidx 305 zkappa_s(ji,jk) 295 zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 306 296 END DO 307 297 END DO … … 322 312 END DO 323 313 DO ji = 1, nidx ! Snow-ice interface 324 zkappa_i(ji,0) 314 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 325 315 END DO 326 316 ! … … 337 327 DO jk = 1, nlay_s 338 328 DO ji = 1, nidx 339 zeta_s(ji,jk) 329 zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji) 340 330 END DO 341 331 END DO … … 352 342 353 343 DO ji = 1, nidx 354 zf (ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH)344 zfnet(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH) 355 345 END DO 356 346 ! … … 359 349 !---------------------------- 360 350 !!layer denotes the number of the layer in the snow or in the ice 361 !! numeqdenotes the reference number of the equation in the tridiagonal351 !!jm denotes the reference number of the equation in the tridiagonal 362 352 !!system, terms of tridiagonal system are indexed as following : 363 353 !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 364 354 365 355 !!ice interior terms (top equation has the same form as the others) 366 367 DO numeq=1,nlay_i+3 368 DO ji = 1, nidx 369 ztrid(ji,numeq,1) = 0. 370 ztrid(ji,numeq,2) = 0. 371 ztrid(ji,numeq,3) = 0. 372 zindterm(ji,numeq)= 0. 373 zindtbis(ji,numeq)= 0. 374 zdiagbis(ji,numeq)= 0. 375 ENDDO 356 ztrid (1:nidx,:,:) = 0._wp 357 zindterm(1:nidx,:) = 0._wp 358 zindtbis(1:nidx,:) = 0._wp 359 zdiagbis(1:nidx,:) = 0._wp 360 361 DO jm = nlay_s + 2, nlay_s + nlay_i 362 DO ji = 1, nidx 363 jk = jm - nlay_s - 1 364 ztrid(ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 365 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 366 ztrid(ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk) 367 zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 368 END DO 376 369 ENDDO 377 370 378 DO numeq = nlay_s + 2, nlay_s + nlay_i 379 DO ji = 1, nidx 380 jk = numeq - nlay_s - 1 381 ztrid(ji,numeq,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 382 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 383 ztrid(ji,numeq,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk) 384 zindterm(ji,numeq) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 385 END DO 386 ENDDO 387 388 numeq = nlay_s + nlay_i + 1 371 jm = nlay_s + nlay_i + 1 389 372 DO ji = 1, nidx 390 373 !!ice bottom term 391 ztrid(ji, numeq,1) =- zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)392 ztrid(ji, numeq,2) =1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )393 ztrid(ji, numeq,3) =0.0394 zindterm(ji, numeq) =ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * &395 & 374 ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) 375 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 376 ztrid(ji,jm,3) = 0.0 377 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 378 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 396 379 ENDDO 397 380 … … 401 384 IF ( ht_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! 402 385 ! !---------------------! 403 ! 404 !!snow interior terms (bottom equation has the same form as the others) 405 DO numeq = 3, nlay_s + 1 406 jk = numeq - 1 407 ztrid(ji,numeq,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 408 ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 409 ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 410 zindterm(ji,numeq) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 386 ! snow interior terms (bottom equation has the same form as the others) 387 DO jm = 3, nlay_s + 1 388 jk = jm - 1 389 ztrid(ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 390 ztrid(ji,jm,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 391 ztrid(ji,jm,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 392 zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 411 393 END DO 412 394 413 ! !case of only one layer in the ice (ice equation is altered)395 ! case of only one layer in the ice (ice equation is altered) 414 396 IF ( nlay_i == 1 ) THEN 415 ztrid(ji,nlay_s+2,3) =0.0416 zindterm(ji,nlay_s+2) =zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)397 ztrid(ji,nlay_s+2,3) = 0.0 398 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 417 399 ENDIF 418 400 419 401 IF ( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting 420 402 421 numeqmin(ji) =1422 numeqmax(ji) =nlay_i + nlay_s + 1423 424 ! !surface equation403 jm_min(ji) = 1 404 jm_max(ji) = nlay_i + nlay_s + 1 405 406 ! surface equation 425 407 ztrid(ji,1,1) = 0.0 426 408 ztrid(ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 427 409 ztrid(ji,1,3) = zg1s * zkappa_s(ji,0) 428 zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zf (ji)429 430 ! !first layer of snow equation431 ztrid(ji,2,1) = 432 ztrid(ji,2,2) = 433 ztrid(ji,2,3) = 434 zindterm(ji,2) = 410 zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 411 412 ! first layer of snow equation 413 ztrid(ji,2,1) = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1) 414 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 415 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 416 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 435 417 436 418 ELSE !-- case 2 : surface is melting 437 419 ! 438 numeqmin(ji) =2439 numeqmax(ji) =nlay_i + nlay_s + 1440 441 ! !first layer of snow equation442 ztrid(ji,2,1) = 443 ztrid(ji,2,2) = 444 ztrid(ji,2,3) = 420 jm_min(ji) = 2 421 jm_max(ji) = nlay_i + nlay_s + 1 422 423 ! first layer of snow equation 424 ztrid(ji,2,1) = 0.0 425 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 426 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 445 427 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * & 446 & 428 & ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 447 429 ENDIF 448 430 ! !---------------------! … … 452 434 IF ( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting 453 435 ! 454 numeqmin(ji) = nlay_s + 1 455 numeqmax(ji) = nlay_i + nlay_s + 1 456 457 !!surface equation 458 ztrid(ji,numeqmin(ji),1) = 0.0 459 ztrid(ji,numeqmin(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1 460 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 461 zindterm(ji,numeqmin(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zf(ji) 462 463 !!first layer of ice equation 464 ztrid(ji,numeqmin(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 465 ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 466 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1) 467 zindterm(ji,numeqmin(ji)+1)= ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 468 469 !!case of only one layer in the ice (surface & ice equations are altered) 470 436 jm_min(ji) = nlay_s + 1 437 jm_max(ji) = nlay_i + nlay_s + 1 438 439 ! surface equation 440 ztrid(ji,jm_min(ji),1) = 0.0 441 ztrid(ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1 442 ztrid(ji,jm_min(ji),3) = zkappa_i(ji,0)*zg1 443 zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji) 444 445 ! first layer of ice equation 446 ztrid(ji,jm_min(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 447 ztrid(ji,jm_min(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 448 ztrid(ji,jm_min(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1) 449 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 450 451 ! case of only one layer in the ice (surface & ice equations are altered) 471 452 IF ( nlay_i == 1 ) THEN 472 ztrid(ji,numeqmin(ji),1) = 0.0 473 ztrid(ji,numeqmin(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 474 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0) * 2.0 475 ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 476 ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 477 ztrid(ji,numeqmin(ji)+1,3) = 0.0 478 479 zindterm(ji,numeqmin(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * & 480 & ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 453 ztrid(ji,jm_min(ji),1) = 0.0 454 ztrid(ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 455 ztrid(ji,jm_min(ji),3) = zkappa_i(ji,0) * 2.0 456 ztrid(ji,jm_min(ji)+1,1) = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 457 ztrid(ji,jm_min(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 458 ztrid(ji,jm_min(ji)+1,3) = 0.0 459 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * & 460 & ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 481 461 ENDIF 482 462 483 463 ELSE !-- case 2 : surface is melting 484 464 485 numeqmin(ji) = nlay_s + 2486 numeqmax(ji) = nlay_i + nlay_s + 1487 488 ! !first layer of ice equation489 ztrid(ji, numeqmin(ji),1) =0.0490 ztrid(ji, numeqmin(ji),2) =1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )491 ztrid(ji, numeqmin(ji),3) =- zeta_i(ji,1) * zkappa_i(ji,1)492 zindterm(ji, numeqmin(ji)) =ztiold(ji,1) + zeta_i(ji,1) * &493 & 494 495 ! !case of only one layer in the ice (surface & ice equations are altered)465 jm_min(ji) = nlay_s + 2 466 jm_max(ji) = nlay_i + nlay_s + 1 467 468 ! first layer of ice equation 469 ztrid(ji,jm_min(ji),1) = 0.0 470 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 471 ztrid(ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 472 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * & 473 & ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 474 475 ! case of only one layer in the ice (surface & ice equations are altered) 496 476 IF ( nlay_i == 1 ) THEN 497 ztrid(ji, numeqmin(ji),1) =0.0498 ztrid(ji, numeqmin(ji),2) =1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )499 ztrid(ji, numeqmin(ji),3) =0.0500 zindterm(ji, numeqmin(ji)) =ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &501 & 477 ztrid(ji,jm_min(ji),1) = 0.0 478 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 479 ztrid(ji,jm_min(ji),3) = 0.0 480 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) & 481 & + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 502 482 ENDIF 503 483 504 484 ENDIF 505 485 ENDIF 506 486 ! 487 zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) 488 zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 489 ! 507 490 END DO 508 491 ! … … 511 494 !------------------------------ 512 495 ! Solve the tridiagonal system with Gauss elimination method. 513 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984. 514 515 maxnumeqmax = 0 516 minnumeqmin = nlay_i+5 517 518 DO ji = 1, nidx 519 zindtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji)) 520 zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2) 521 minnumeqmin = MIN(numeqmin(ji),minnumeqmin) 522 maxnumeqmax = MAX(numeqmax(ji),maxnumeqmax) 523 END DO 524 525 DO jk = minnumeqmin+1, maxnumeqmax 526 DO ji = 1, nidx 527 numeq = min(max(numeqmin(ji)+1,jk),numeqmax(ji)) 528 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3) / zdiagbis(ji,numeq-1) 529 zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1) 496 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 497 jm_maxt = 0 498 jm_mint = nlay_i+5 499 DO ji = 1, nidx 500 jm_mint = MIN(jm_min(ji),jm_mint) 501 jm_maxt = MAX(jm_max(ji),jm_maxt) 502 END DO 503 504 DO jk = jm_mint+1, jm_maxt 505 DO ji = 1, nidx 506 jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 507 zdiagbis(ji,jm) = ztrid(ji,jm,2) - ztrid(ji,jm,1) * ztrid(ji,jm-1,3) / zdiagbis(ji,jm-1) 508 zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 530 509 END DO 531 510 END DO … … 533 512 DO ji = 1, nidx 534 513 ! ice temperatures 535 t_i_1d(ji,nlay_i) = zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji))536 END DO 537 538 DO numeq= nlay_i + nlay_s, nlay_s + 2, -1539 DO ji = 1, nidx 540 jk = numeq- nlay_s - 1541 t_i_1d(ji,jk) = ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq)514 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 515 END DO 516 517 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 518 DO ji = 1, nidx 519 jk = jm - nlay_s - 1 520 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 542 521 END DO 543 522 END DO … … 546 525 ! snow temperatures 547 526 IF( ht_s_1d(ji) > 0._wp ) THEN 548 t_s_1d(ji,nlay_s) = 549 & 527 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) & 528 & / zdiagbis(ji,nlay_s+1) 550 529 ENDIF 551 530 ! surface temperature 552 531 ztsub(ji) = t_su_1d(ji) 553 532 IF( t_su_1d(ji) < rt0 ) THEN 554 t_su_1d(ji) = ( zindtbis(ji, numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) * &555 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji, numeqmin(ji))533 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 534 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 556 535 ENDIF 557 536 END DO … … 599 578 DO ji = 1, nidx 600 579 ! ! surface ice conduction flux 601 fc_su(ji) 602 & 580 fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & 581 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji)) 603 582 ! ! bottom ice conduction flux 604 fc_bo_i(ji) 583 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 605 584 END DO 606 585 … … 637 616 DO ji = 1, nidx 638 617 !--- Conduction fluxes (positive downwards) 639 diag_fc_bo_1d(ji) 640 diag_fc_su_1d(ji) 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) 641 620 642 621 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 643 622 zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 644 623 IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 645 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + &624 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + & 646 625 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 647 626 ELSE
Note: See TracChangeset
for help on using the changeset viewer.
