- Timestamp:
- 2015-02-02T11:28:50+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r5047 r5048 100 100 INTEGER :: nconv ! number of iterations in iterative procedure 101 101 INTEGER :: minnumeqmin, maxnumeqmax 102 102 103 INTEGER, POINTER, DIMENSION(:) :: numeqmin ! reference number of top equation 103 104 INTEGER, POINTER, DIMENSION(:) :: numeqmax ! reference number of bottom equation 104 105 INTEGER, POINTER, DIMENSION(:) :: isnow ! switch for presence (1) or absence (0) of snow 106 105 107 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system 106 108 REAL(wp) :: zg1 = 2._wp ! … … 113 115 REAL(wp) :: zerritmax ! current maximal error on temperature 114 116 REAL(wp) :: zhsu 115 REAL(wp), POINTER, DIMENSION(:) :: ztfs ! ice melting point 116 REAL(wp), POINTER, DIMENSION(:) :: ztsub ! old surface temperature (before the iterative procedure ) 117 REAL(wp), POINTER, DIMENSION(:) :: ztsubit ! surface temperature at previous iteration 118 REAL(wp), POINTER, DIMENSION(:) :: zh_i ! ice layer thickness 119 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness 120 REAL(wp), POINTER, DIMENSION(:) :: zfsw ! solar radiation absorbed at the surface 121 REAL(wp), POINTER, DIMENSION(:) :: zf ! surface flux function 122 REAL(wp), POINTER, DIMENSION(:) :: dzf ! derivative of the surface flux function 123 REAL(wp), POINTER, DIMENSION(:) :: zerrit ! current error on temperature 124 REAL(wp), POINTER, DIMENSION(:) :: zdifcase ! case of the equation resolution (1->4) 125 REAL(wp), POINTER, DIMENSION(:) :: zftrice ! solar radiation transmitted through the ice 126 REAL(wp), POINTER, DIMENSION(:) :: zihic 127 REAL(wp), POINTER, DIMENSION(:,:) :: ztcond_i ! Ice thermal conductivity 128 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_i ! Radiation transmitted through the ice 129 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_i ! Radiation absorbed in the ice 130 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_i ! Kappa factor in the ice 131 REAL(wp), POINTER, DIMENSION(:,:) :: ztib ! Old temperature in the ice 132 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_i ! Eta factor in the ice 133 REAL(wp), POINTER, DIMENSION(:,:) :: ztitemp ! Temporary temperature in the ice to check the convergence 134 REAL(wp), POINTER, DIMENSION(:,:) :: zspeche_i ! Ice specific heat 135 REAL(wp), POINTER, DIMENSION(:,:) :: z_i ! Vertical cotes of the layers in the ice 136 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_s ! Radiation transmited through the snow 137 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_s ! Radiation absorbed in the snow 138 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_s ! Kappa factor in the snow 139 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_s ! Eta factor in the snow 140 REAL(wp), POINTER, DIMENSION(:,:) :: ztstemp ! Temporary temperature in the snow to check the convergence 141 REAL(wp), POINTER, DIMENSION(:,:) :: ztsb ! Temporary temperature in the snow 142 REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! Vertical cotes of the layers in the snow 143 REAL(wp), POINTER, DIMENSION(:,:) :: zswiterm ! Independent term 144 REAL(wp), POINTER, DIMENSION(:,:) :: zswitbis ! temporary independent term 145 REAL(wp), POINTER, DIMENSION(:,:) :: zdiagbis 146 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! tridiagonal system terms 117 118 REAL(wp), POINTER, DIMENSION(:) :: ztfs ! ice melting point 119 REAL(wp), POINTER, DIMENSION(:) :: ztsub ! old surface temperature (before the iterative procedure ) 120 REAL(wp), POINTER, DIMENSION(:) :: ztsubit ! surface temperature at previous iteration 121 REAL(wp), POINTER, DIMENSION(:) :: zh_i ! ice layer thickness 122 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness 123 REAL(wp), POINTER, DIMENSION(:) :: zfsw ! solar radiation absorbed at the surface 124 REAL(wp), POINTER, DIMENSION(:) :: zf ! surface flux function 125 REAL(wp), POINTER, DIMENSION(:) :: dzf ! derivative of the surface flux function 126 REAL(wp), POINTER, DIMENSION(:) :: zerrit ! current error on temperature 127 REAL(wp), POINTER, DIMENSION(:) :: zdifcase ! case of the equation resolution (1->4) 128 REAL(wp), POINTER, DIMENSION(:) :: zftrice ! solar radiation transmitted through the ice 129 REAL(wp), POINTER, DIMENSION(:) :: zihic 130 131 REAL(wp), POINTER, DIMENSION(:,:) :: ztcond_i ! Ice thermal conductivity 132 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_i ! Radiation transmitted through the ice 133 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_i ! Radiation absorbed in the ice 134 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_i ! Kappa factor in the ice 135 REAL(wp), POINTER, DIMENSION(:,:) :: ztib ! Old temperature in the ice 136 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_i ! Eta factor in the ice 137 REAL(wp), POINTER, DIMENSION(:,:) :: ztitemp ! Temporary temperature in the ice to check the convergence 138 REAL(wp), POINTER, DIMENSION(:,:) :: zspeche_i ! Ice specific heat 139 REAL(wp), POINTER, DIMENSION(:,:) :: z_i ! Vertical cotes of the layers in the ice 140 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_s ! Radiation transmited through the snow 141 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_s ! Radiation absorbed in the snow 142 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_s ! Kappa factor in the snow 143 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_s ! Eta factor in the snow 144 REAL(wp), POINTER, DIMENSION(:,:) :: ztstemp ! Temporary temperature in the snow to check the convergence 145 REAL(wp), POINTER, DIMENSION(:,:) :: ztsb ! Temporary temperature in the snow 146 REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! Vertical cotes of the layers in the snow 147 REAL(wp), POINTER, DIMENSION(:,:) :: zindterm ! 'Ind'ependent term 148 REAL(wp), POINTER, DIMENSION(:,:) :: zindtbis ! Temporary 'ind'ependent term 149 REAL(wp), POINTER, DIMENSION(:,:) :: zdiagbis ! Temporary 'dia'gonal term 150 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! Tridiagonal system terms 151 147 152 ! diag errors on heat 148 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err 153 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err 154 155 ! Mono-category 156 REAL(wp) :: zepsilon ! determines thres. above which computation of G(h) is done 157 REAL(wp) :: zratio_s ! dummy factor 158 REAL(wp) :: zratio_i ! dummy factor 159 REAL(wp) :: zh_thres ! thickness thres. for G(h) computation 160 REAL(wp) :: zhe ! dummy factor 161 REAL(wp) :: zswitch ! dummy switch 162 REAL(wp) :: zkimean ! mean sea ice thermal conductivity 163 REAL(wp) :: zfac ! dummy factor 164 REAL(wp) :: zihe ! dummy factor 165 REAL(wp) :: zheshth ! dummy factor 166 167 REAL(wp), POINTER, DIMENSION(:) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 168 149 169 !!------------------------------------------------------------------ 150 170 ! 151 171 CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 152 172 CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 153 CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic )173 CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 154 174 CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0) 155 175 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0) 156 CALL wrk_alloc( jpij, nlay_i+3, z switerm, zswitbis, zdiagbis )176 CALL wrk_alloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis ) 157 177 CALL wrk_alloc( jpij, nlay_i+3, 3, ztrid ) 158 178 … … 200 220 ! 201 221 !------------------------------------------------------------------------------| 202 ! 2) Radiation s|222 ! 2) Radiation | 203 223 !------------------------------------------------------------------------------| 204 224 ! … … 213 233 ! zftrice = io.qsr_ice is below the surface 214 234 ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 235 ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 215 236 zhsu = 0.1_wp ! threshold for the computation of i0 216 !fr1_i0_1d = i0 for a thin ice surface217 !fr1_i0_2d = i0 for a thick ice surface218 237 DO ji = kideb , kiut 219 238 ! switches … … 339 358 END DO 340 359 ENDIF 341 ! 342 !------------------------------------------------------------------------------| 343 ! 5) kappa factors | 344 !------------------------------------------------------------------------------| 345 ! 360 361 ! 362 !------------------------------------------------------------------------------| 363 ! 6) G(he) - enhancement of thermal conductivity in mono-category case | 364 !------------------------------------------------------------------------------| 365 ! 366 ! Computation of effective thermal conductivity G(h) 367 ! Used in mono-category case only to simulate an ITD implicitly 368 ! Fichefet and Morales Maqueda, JGR 1997 369 370 zghe(:) = 0._wp 371 372 SELECT CASE ( nn_monocat ) 373 374 CASE (0,2,4) 375 376 zghe(kideb:kiut) = 1._wp 377 378 CASE (1,3) ! LIM3 379 380 zepsilon = 0.1 381 zh_thres = EXP( 1._wp ) * zepsilon / 2. 382 383 DO ji = kideb, kiut 384 385 ! Mean sea ice thermal conductivity 386 zkimean = SUM( ztcond_i(ji,0:nlay_i) ) / REAL(nlay_i+1,wp) 387 388 ! Effective thickness he (zhe) 389 zfac = 1._wp / ( rcdsn + zkimean ) 390 zratio_s = rcdsn * zfac 391 zratio_i = zkimean * zfac 392 zhe = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji) 393 394 ! G(he) 395 zswitch = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) ) ! =0 if zhe < zh_thres, if > 396 zghe(ji) = ( 1.0 - zswitch ) + zswitch * ( 0.5 + 0.5 * LOG( 2.*zhe / zepsilon ) ) 397 398 ! Impose G(he) < 2. 399 zghe(ji) = MIN( zghe(ji), 2.0 ) 400 401 END DO 402 403 END SELECT 404 405 ! 406 !------------------------------------------------------------------------------| 407 ! 7) kappa factors | 408 !------------------------------------------------------------------------------| 409 ! 410 !--- Snow 346 411 DO ji = kideb, kiut 347 348 !-- Snow kappa factors 349 zkappa_s(ji,0) = rcdsn / MAX(epsi10,zh_s(ji)) 350 zkappa_s(ji,nlay_s) = rcdsn / MAX(epsi10,zh_s(ji)) 412 zfac = 1. / MAX( epsi10 , zh_s(ji) ) 413 zkappa_s(ji,0) = zghe(ji) * rcdsn * zfac 414 zkappa_s(ji,nlay_s) = zghe(ji) * rcdsn * zfac 351 415 END DO 352 416 353 417 DO jk = 1, nlay_s-1 354 418 DO ji = kideb , kiut 355 zkappa_s(ji,jk) = 2.0 * rcdsn / &356 MAX(epsi10,2.0*zh_s(ji))357 358 END DO 359 419 zkappa_s(ji,jk) = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0*zh_s(ji) ) 420 END DO 421 END DO 422 423 !--- Ice 360 424 DO jk = 1, nlay_i-1 361 425 DO ji = kideb , kiut 362 !-- Ice kappa factors 363 zkappa_i(ji,jk) = 2.0*ztcond_i(ji,jk)/ & 364 MAX(epsi10,2.0*zh_i(ji)) 365 END DO 366 END DO 367 368 DO ji = kideb , kiut 369 zkappa_i(ji,0) = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji)) 370 zkappa_i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji)) 371 !-- Interface 372 zkappa_s(ji,nlay_s) = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, & 373 (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 374 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & 375 + zkappa_i(ji,0)*REAL( 1 - isnow(ji) ) 376 END DO 377 ! 378 !------------------------------------------------------------------------------| 379 ! 6) Sea ice specific heat, eta factors | 426 zkappa_i(ji,jk) = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0*zh_i(ji) ) 427 END DO 428 END DO 429 430 !--- Snow-ice interface 431 DO ji = kideb , kiut 432 zfac = 1./ MAX( epsi10 , zh_i(ji) ) 433 zkappa_i(ji,0) = zghe(ji) * ztcond_i(ji,0) * zfac 434 zkappa_i(ji,nlay_i) = zghe(ji) * ztcond_i(ji,nlay_i) * zfac 435 zkappa_s(ji,nlay_s) = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / & 436 & MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) ) 437 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*REAL( isnow(ji), wp ) + zkappa_i(ji,0)*REAL( 1 - isnow(ji), wp ) 438 END DO 439 440 ! 441 !------------------------------------------------------------------------------| 442 ! 8) Sea ice specific heat, eta factors | 380 443 !------------------------------------------------------------------------------| 381 444 ! … … 383 446 DO ji = kideb , kiut 384 447 ztitemp(ji,jk) = t_i_1d(ji,jk) 385 zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ & 386 MAX((t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt),epsi10) 387 zeta_i(ji,jk) = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), & 388 epsi10) 448 zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ MAX( (t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt) , epsi10 ) 449 zeta_i(ji,jk) = rdt_ice / MAX( rhoic*zspeche_i(ji,jk)*zh_i(ji), epsi10 ) 389 450 END DO 390 451 END DO … … 396 457 END DO 397 458 END DO 398 ! 399 !------------------------------------------------------------------------------| 400 ! 7) surface flux computation | 459 460 ! 461 !------------------------------------------------------------------------------| 462 ! 9) surface flux computation | 401 463 !------------------------------------------------------------------------------| 402 464 ! … … 418 480 ! 419 481 !------------------------------------------------------------------------------| 420 ! 8) tridiagonal system terms|482 ! 10) tridiagonal system terms | 421 483 !------------------------------------------------------------------------------| 422 484 ! … … 433 495 ztrid(ji,numeq,2) = 0. 434 496 ztrid(ji,numeq,3) = 0. 435 z switerm(ji,numeq)= 0.436 z switbis(ji,numeq)= 0.497 zindterm(ji,numeq)= 0. 498 zindtbis(ji,numeq)= 0. 437 499 zdiagbis(ji,numeq)= 0. 438 500 ENDDO … … 446 508 zkappa_i(ji,jk)) 447 509 ztrid(ji,numeq,3) = - zeta_i(ji,jk)*zkappa_i(ji,jk) 448 z switerm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk)* &510 zindterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk)* & 449 511 zradab_i(ji,jk) 450 512 END DO … … 458 520 + zkappa_i(ji,nlay_i-1) ) 459 521 ztrid(ji,numeq,3) = 0.0 460 z switerm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* &522 zindterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 461 523 ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 462 524 * t_bo_1d(ji) ) … … 478 540 zkappa_s(ji,jk) ) 479 541 ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 480 z switerm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk)* &542 zindterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk)* & 481 543 zradab_s(ji,jk) 482 544 END DO … … 485 547 IF ( nlay_i.eq.1 ) THEN 486 548 ztrid(ji,nlay_s+2,3) = 0.0 487 z switerm(ji,nlay_s+2) = zswiterm(ji,nlay_s+2) + zkappa_i(ji,1)* &549 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 488 550 t_bo_1d(ji) 489 551 ENDIF … … 502 564 ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 503 565 ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 504 z switerm(ji,1) = dzf(ji)*t_su_1d(ji) - zf(ji)566 zindterm(ji,1) = dzf(ji)*t_su_1d(ji) - zf(ji) 505 567 506 568 !!first layer of snow equation … … 508 570 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) 509 571 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 510 z switerm(ji,2) = ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)572 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 511 573 512 574 ELSE … … 525 587 zkappa_s(ji,0) * zg1s ) 526 588 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 527 z switerm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * &589 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * & 528 590 ( zradab_s(ji,1) + & 529 591 zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) … … 549 611 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1 550 612 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 551 z switerm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji)613 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji) 552 614 553 615 !!first layer of ice equation … … 556 618 + zkappa_i(ji,0) * zg1 ) 557 619 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1) 558 z switerm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)620 zindterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) 559 621 560 622 !!case of only one layer in the ice (surface & ice equations are altered) … … 569 631 ztrid(ji,numeqmin(ji)+1,3) = 0.0 570 632 571 z switerm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1)* &633 zindterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1)* & 572 634 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) ) 573 635 ENDIF … … 589 651 zg1) 590 652 ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 591 z switerm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &653 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 592 654 zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 593 655 … … 598 660 zkappa_i(ji,1)) 599 661 ztrid(ji,numeqmin(ji),3) = 0.0 600 z switerm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)* &662 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)* & 601 663 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) & 602 664 + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 … … 610 672 ! 611 673 !------------------------------------------------------------------------------| 612 ! 9) tridiagonal system solving|674 ! 11) tridiagonal system solving | 613 675 !------------------------------------------------------------------------------| 614 676 ! … … 622 684 623 685 DO ji = kideb , kiut 624 z switbis(ji,numeqmin(ji)) = zswiterm(ji,numeqmin(ji))686 zindtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji)) 625 687 zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2) 626 688 minnumeqmin = MIN(numeqmin(ji),minnumeqmin) … … 633 695 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 634 696 ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) 635 z switbis(ji,numeq) = zswiterm(ji,numeq) - ztrid(ji,numeq,1)* &636 z switbis(ji,numeq-1)/zdiagbis(ji,numeq-1)697 zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1)* & 698 zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 637 699 END DO 638 700 END DO … … 640 702 DO ji = kideb , kiut 641 703 ! ice temperatures 642 t_i_1d(ji,nlay_i) = z switbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))704 t_i_1d(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 643 705 END DO 644 706 … … 646 708 DO ji = kideb , kiut 647 709 jk = numeq - nlay_s - 1 648 t_i_1d(ji,jk) = (z switbis(ji,numeq) - ztrid(ji,numeq,3)* &710 t_i_1d(ji,jk) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & 649 711 t_i_1d(ji,jk+1))/zdiagbis(ji,numeq) 650 712 END DO … … 654 716 ! snow temperatures 655 717 IF (ht_s_1d(ji).GT.0._wp) & 656 t_s_1d(ji,nlay_s) = (z switbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &718 t_s_1d(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 657 719 * t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & 658 720 * MAX(0.0,SIGN(1.0,ht_s_1d(ji))) … … 662 724 ztsubit(ji) = t_su_1d(ji) 663 725 IF( t_su_1d(ji) < ztfs(ji) ) & 664 t_su_1d(ji) = ( z switbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1) &726 t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1) & 665 727 & + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 666 728 END DO 667 729 ! 668 730 !-------------------------------------------------------------------------- 669 ! 1 0) Has the scheme converged ?, end of the iterative procedure |731 ! 12) Has the scheme converged ?, end of the iterative procedure | 670 732 !-------------------------------------------------------------------------- 671 733 ! … … 709 771 ! 710 772 !-------------------------------------------------------------------------! 711 ! 1 1) Fluxes at the interfaces !773 ! 13) Fluxes at the interfaces ! 712 774 !-------------------------------------------------------------------------! 713 775 DO ji = kideb, kiut … … 771 833 CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 772 834 CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 773 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic )835 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 774 836 CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, & 775 837 & ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 776 838 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 777 CALL wrk_dealloc( jpij, nlay_i+3, z switerm, zswitbis, zdiagbis )839 CALL wrk_dealloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis ) 778 840 CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid ) 779 841 CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err )
Note: See TracChangeset
for help on using the changeset viewer.
