Changeset 4990 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
- Timestamp:
- 2014-12-15T17:42:49+01:00 (10 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 (modified) (30 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r4873 r4990 25 25 USE wrk_nemo ! work arrays 26 26 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 27 USE cpl_oasis3, ONLY : lk_cpl27 USE sbc_oce, ONLY : lk_cpl 28 28 29 29 IMPLICIT NONE … … 32 32 PUBLIC lim_thd_dif ! called by lim_thd 33 33 34 REAL(wp) :: epsi10 = 1.e-10_wp !35 34 !!---------------------------------------------------------------------- 36 35 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 108 107 REAL(wp) :: zgamma = 18009._wp ! for specific heat 109 108 REAL(wp) :: zbeta = 0.117_wp ! for thermal conductivity (could be 0.13) 110 REAL(wp) :: zraext_s = 1 .e+8_wp! extinction coefficient of radiation in the snow109 REAL(wp) :: zraext_s = 10._wp ! extinction coefficient of radiation in the snow 111 110 REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity 112 111 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered as 0°C … … 141 140 REAL(wp), POINTER, DIMENSION(:,:) :: ztsb ! Temporary temperature in the snow 142 141 REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! Vertical cotes of the layers in the snow 143 REAL(wp), POINTER, DIMENSION(:,:) :: z indterm ! Independent term144 REAL(wp), POINTER, DIMENSION(:,:) :: z indtbis ! temporary independent term142 REAL(wp), POINTER, DIMENSION(:,:) :: zswiterm ! Independent term 143 REAL(wp), POINTER, DIMENSION(:,:) :: zswitbis ! temporary independent term 145 144 REAL(wp), POINTER, DIMENSION(:,:) :: zdiagbis 146 145 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! tridiagonal system terms 147 146 ! diag errors on heat 148 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini 149 REAL(wp) :: zhfx_err 147 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err 150 148 !!------------------------------------------------------------------ 151 149 ! … … 155 153 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) 156 154 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0) 157 CALL wrk_alloc( jpij, nlay_i+3, z indterm, zindtbis, zdiagbis )155 CALL wrk_alloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis ) 158 156 CALL wrk_alloc( jpij, nlay_i+3, 3, ztrid ) 159 157 160 CALL wrk_alloc( jpij, zdq, zq_ini )158 CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err ) 161 159 162 160 ! --- diag error on heat diffusion - PART 1 --- ! … … 272 270 273 271 DO ji = kideb, kiut ! Radiation transmitted below the ice 274 !!!ftr_ice_1d(ji) = ftr_ice_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_1d(ji) / at_i_1d(ji) ! clem modif275 272 ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 276 273 END DO … … 408 405 !------------------------------------------------------------------------------| 409 406 ! 410 DO ji = kideb , kiut 411 ! update of the non solar flux according to the update in T_su 412 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) ) 413 407 IF( .NOT. lk_cpl ) THEN !--- forced atmosphere case 408 DO ji = kideb , kiut 409 ! update of the non solar flux according to the update in T_su 410 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) ) 411 END DO 412 ENDIF 413 414 ! Update incoming flux 415 DO ji = kideb , kiut 414 416 ! update incoming flux 415 417 zf(ji) = zfsw(ji) & ! net absorbed solar radiation 416 + qns_ice_1d(ji) ! non solar total flux418 + qns_ice_1d(ji) ! non solar total flux 417 419 ! (LWup, LWdw, SH, LH) 418 420 END DO … … 435 437 ztrid(ji,numeq,2) = 0. 436 438 ztrid(ji,numeq,3) = 0. 437 z indterm(ji,numeq)= 0.438 z indtbis(ji,numeq)= 0.439 zswiterm(ji,numeq)= 0. 440 zswitbis(ji,numeq)= 0. 439 441 zdiagbis(ji,numeq)= 0. 440 442 ENDDO … … 448 450 zkappa_i(ji,jk)) 449 451 ztrid(ji,numeq,3) = - zeta_i(ji,jk)*zkappa_i(ji,jk) 450 z indterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk)* &452 zswiterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk)* & 451 453 zradab_i(ji,jk) 452 454 END DO … … 460 462 + zkappa_i(ji,nlay_i-1) ) 461 463 ztrid(ji,numeq,3) = 0.0 462 z indterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* &464 zswiterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 463 465 ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 464 466 * t_bo_1d(ji) ) … … 480 482 zkappa_s(ji,jk) ) 481 483 ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 482 z indterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk)* &484 zswiterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk)* & 483 485 zradab_s(ji,jk) 484 486 END DO … … 487 489 IF ( nlay_i.eq.1 ) THEN 488 490 ztrid(ji,nlay_s+2,3) = 0.0 489 z indterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* &491 zswiterm(ji,nlay_s+2) = zswiterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 490 492 t_bo_1d(ji) 491 493 ENDIF … … 504 506 ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 505 507 ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 506 z indterm(ji,1) = dzf(ji)*t_su_1d(ji) - zf(ji)508 zswiterm(ji,1) = dzf(ji)*t_su_1d(ji) - zf(ji) 507 509 508 510 !!first layer of snow equation … … 510 512 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) 511 513 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 512 z indterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)514 zswiterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 513 515 514 516 ELSE … … 527 529 zkappa_s(ji,0) * zg1s ) 528 530 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 529 z indterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * &531 zswiterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * & 530 532 ( zradab_s(ji,1) + & 531 533 zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) … … 551 553 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1 552 554 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 553 z indterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji)555 zswiterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji) 554 556 555 557 !!first layer of ice equation … … 558 560 + zkappa_i(ji,0) * zg1 ) 559 561 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1) 560 z indterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)562 zswiterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) 561 563 562 564 !!case of only one layer in the ice (surface & ice equations are altered) … … 571 573 ztrid(ji,numeqmin(ji)+1,3) = 0.0 572 574 573 z indterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1)* &575 zswiterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1)* & 574 576 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) ) 575 577 ENDIF … … 591 593 zg1) 592 594 ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 593 z indterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &595 zswiterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 594 596 zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 595 597 … … 600 602 zkappa_i(ji,1)) 601 603 ztrid(ji,numeqmin(ji),3) = 0.0 602 z indterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)* &604 zswiterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)* & 603 605 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) & 604 606 + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 … … 624 626 625 627 DO ji = kideb , kiut 626 z indtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji))628 zswitbis(ji,numeqmin(ji)) = zswiterm(ji,numeqmin(ji)) 627 629 zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2) 628 630 minnumeqmin = MIN(numeqmin(ji),minnumeqmin) … … 635 637 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 636 638 ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) 637 z indtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1)* &638 z indtbis(ji,numeq-1)/zdiagbis(ji,numeq-1)639 zswitbis(ji,numeq) = zswiterm(ji,numeq) - ztrid(ji,numeq,1)* & 640 zswitbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 639 641 END DO 640 642 END DO … … 642 644 DO ji = kideb , kiut 643 645 ! ice temperatures 644 t_i_1d(ji,nlay_i) = z indtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))646 t_i_1d(ji,nlay_i) = zswitbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 645 647 END DO 646 648 … … 648 650 DO ji = kideb , kiut 649 651 jk = numeq - nlay_s - 1 650 t_i_1d(ji,jk) = (z indtbis(ji,numeq) - ztrid(ji,numeq,3)* &652 t_i_1d(ji,jk) = (zswitbis(ji,numeq) - ztrid(ji,numeq,3)* & 651 653 t_i_1d(ji,jk+1))/zdiagbis(ji,numeq) 652 654 END DO … … 656 658 ! snow temperatures 657 659 IF (ht_s_1d(ji).GT.0._wp) & 658 t_s_1d(ji,nlay_s) = (z indtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &660 t_s_1d(ji,nlay_s) = (zswitbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 659 661 * t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & 660 662 * MAX(0.0,SIGN(1.0,ht_s_1d(ji))) … … 664 666 ztsubit(ji) = t_su_1d(ji) 665 667 IF( t_su_1d(ji) < ztfs(ji) ) & 666 t_su_1d(ji) = ( z indtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1) &668 t_su_1d(ji) = ( zswitbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1) & 667 669 & + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 668 670 END DO … … 740 742 CALL lim_thd_enmelt( kideb, kiut ) 741 743 742 ! --- diag erroron heat diffusion - PART 2 --- !744 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 743 745 DO ji = kideb, kiut 744 746 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) + & 745 747 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 746 zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 747 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err * a_i_1d(ji) 748 ! --- correction of qns_ice and surface conduction flux --- ! 749 qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err 750 fc_su (ji) = fc_su (ji) - zhfx_err 751 ! --- Heat flux at the ice surface in W.m-2 --- ! 748 zhfx_err(ji) = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 749 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 750 END DO 751 752 ! diagnose external surface (forced case) or bottom (forced case) from heat conservation 753 IF( .NOT. lk_cpl ) THEN ! --- forced case: qns_ice and fc_su are diagnosed 754 ! 755 DO ji = kideb, kiut 756 qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err(ji) 757 fc_su (ji) = fc_su(ji) - zhfx_err(ji) 758 END DO 759 ! 760 ELSE ! --- coupled case: ocean turbulent heat flux is diagnosed 761 ! 762 DO ji = kideb, kiut 763 fhtur_1d (ji) = fhtur_1d(ji) - zhfx_err(ji) 764 END DO 765 ! 766 ENDIF 767 768 ! --- compute diagnostic net heat flux at the surface of the snow-ice system (W.m2) 769 DO ji = kideb, kiut 752 770 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 753 771 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_1d(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) … … 761 779 & ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 762 780 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 763 CALL wrk_dealloc( jpij, nlay_i+3, z indterm, zindtbis, zdiagbis )781 CALL wrk_dealloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis ) 764 782 CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid ) 765 CALL wrk_dealloc( jpij, zdq, zq_ini )783 CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) 766 784 767 785 END SUBROUTINE lim_thd_dif … … 778 796 ! 779 797 INTEGER :: ji, jk ! dummy loop indices 780 REAL(wp) :: ztmelts , zindb! local scalar798 REAL(wp) :: ztmelts ! local scalar 781 799 !!------------------------------------------------------------------- 782 800 ! … … 784 802 DO ji = kideb, kiut 785 803 ztmelts = - tmut * s_i_1d(ji,jk) + rtt 786 zindb= MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi10 ) )804 rswitch = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi10 ) ) 787 805 q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) & 788 & + lfus * ( 1.0 - zindb* ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk)-rtt, -epsi10 ) ) &806 & + lfus * ( 1.0 - rswitch * ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk)-rtt, -epsi10 ) ) & 789 807 & - rcp * ( ztmelts-rtt ) ) 790 808 END DO
Note: See TracChangeset
for help on using the changeset viewer.
