Changeset 8939 for branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3
- Timestamp:
- 2017-12-07T17:53:39+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
r8933 r8939 44 44 INTEGER :: nice_zdf ! Choice of the type of vertical heat diffusion formulation 45 45 ! ! associated indices: 46 INTEGER, PARAMETER :: np_BL99 = 1! Bitz and Lipscomb (1999)47 48 INTEGER, PARAMETER :: np_zdf_jules_OFF 49 INTEGER, PARAMETER :: np_zdf_jules_SND 50 INTEGER, PARAMETER :: np_zdf_jules_RCV 46 INTEGER, PARAMETER :: np_BL99 = 1 ! Bitz and Lipscomb (1999) 47 48 INTEGER, PARAMETER :: np_zdf_jules_OFF = 0 ! compute all temperatures from qsr and qns 49 INTEGER, PARAMETER :: np_zdf_jules_SND = 1 ! compute conductive heat flux and surface temperature from qsr and qns 50 INTEGER, PARAMETER :: np_zdf_jules_RCV = 2 ! compute snow and inner ice temperatures from qcnd 51 51 52 52 !!---------------------------------------------------------------------- … … 143 143 REAL(wp) :: zfac ! dummy factor 144 144 ! 145 REAL(wp), DIMENSION(jpij) :: isnow! switch for presence (1) or absence (0) of snow146 REAL(wp), DIMENSION(jpij) :: ztsub! surface temperature at previous iteration147 REAL(wp), DIMENSION(jpij) 148 REAL(wp), DIMENSION(jpij) 149 REAL(wp), DIMENSION(jpij) :: zqns_ice_b! solar radiation absorbed at the surface150 REAL(wp), DIMENSION(jpij) :: zfnet! surface flux function151 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b! derivative of the surface flux function145 REAL(wp), DIMENSION(jpij) :: isnow ! switch for presence (1) or absence (0) of snow 146 REAL(wp), DIMENSION(jpij) :: ztsub ! surface temperature at previous iteration 147 REAL(wp), DIMENSION(jpij) :: zh_i, z1_h_i ! ice layer thickness 148 REAL(wp), DIMENSION(jpij) :: zh_s, z1_h_s ! snow layer thickness 149 REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface 150 REAL(wp), DIMENSION(jpij) :: zfnet ! surface flux function 151 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function 152 152 ! 153 153 REAL(wp), DIMENSION(jpij ) :: ztsuold ! Old surface temperature in the ice … … 265 265 iconv = 0 ! number of iterations 266 266 zdti_max = 1000._wp ! maximal value of error on all points 267 ! 267 268 ! !============================! 268 269 DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max ) ! Iterative procedure begins ! … … 292 293 ! 293 294 DO ji = 1, npti 294 ztcond_i(ji,0) = rcdic + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) &295 ztcond_i(ji,0) = rcdic + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 295 296 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 296 ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) &297 ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 297 298 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 298 299 END DO 299 300 DO jk = 1, nlay_i-1 300 301 DO ji = 1, npti 301 ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / 302 ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 302 303 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) & 303 304 & - 0.011_wp * ( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) … … 321 322 DO ji = 1, npti 322 323 zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp ) ! Mean sea ice thermal conductivity 323 zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i ) ! Effective thickness he (zhe)324 zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i ) ! Effective thickness he (zhe) 324 325 IF( zhe >= zepsilon * 0.5_wp * EXP(1._wp) ) THEN 325 326 zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) ) ! G(he) … … 394 395 !---------------------------- 395 396 396 DO ji = 1, npti397 397 ! update of the non solar flux according to the update in T_su 398 DO ji = 1, npti 398 399 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 399 400 END DO … … 483 484 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 484 485 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 485 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * & 486 & ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 486 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 487 487 ENDIF 488 488 ! !---------------------! … … 515 515 ztrid(ji,jm_min(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 516 516 ztrid(ji,jm_min(ji)+1,3) = 0.0 517 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * & 518 & ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 517 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 519 518 ENDIF 520 519 … … 561 560 562 561 DO jk = jm_mint+1, jm_maxt 563 DO ji = 1, npti 564 jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 565 zdiagbis(ji,jm) = ztrid(ji,jm,2) - ztrid(ji,jm,1) * ztrid(ji,jm-1,3) / zdiagbis(ji,jm-1) 566 zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 567 END DO 568 END DO 569 570 DO ji = 1, npti 562 DO ji = 1, npti 563 jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 564 zdiagbis(ji,jm) = ztrid(ji,jm,2) - ztrid(ji,jm,1) * ztrid(ji,jm-1,3) / zdiagbis(ji,jm-1) 565 zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 566 END DO 567 END DO 568 571 569 ! ice temperatures 572 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 570 DO ji = 1, npti 571 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 573 572 END DO 574 573 575 574 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 576 DO ji = 1, npti 577 jk = jm - nlay_s - 1 578 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 579 END DO 580 END DO 581 582 DO ji = 1, npti 583 ! snow temperatures 584 IF( h_s_1d(ji) > 0._wp ) THEN 585 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) & 586 & / zdiagbis(ji,nlay_s+1) 587 ENDIF 588 ! surface temperature 589 ztsub(ji) = t_su_1d(ji) 590 IF( t_su_1d(ji) < rt0 ) THEN 591 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 592 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 593 ENDIF 575 DO ji = 1, npti 576 jk = jm - nlay_s - 1 577 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 578 END DO 579 END DO 580 581 DO ji = 1, npti 582 ! snow temperatures 583 IF( h_s_1d(ji) > 0._wp ) THEN 584 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) 585 ENDIF 586 ! surface temperature 587 ztsub(ji) = t_su_1d(ji) 588 IF( t_su_1d(ji) < rt0 ) THEN 589 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 590 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 591 ENDIF 594 592 END DO 595 593 ! … … 601 599 zdti_max = 0._wp 602 600 DO ji = 1, npti 603 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )604 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )605 END DO 606 607 DO jk =1, nlay_s608 DO ji = 1, npti609 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )610 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )611 END DO612 END DO 613 614 DO jk =1, nlay_i615 DO ji = 1, npti616 ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0617 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )618 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )619 END DO601 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 602 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 603 END DO 604 605 DO jk = 1, nlay_s 606 DO ji = 1, npti 607 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 608 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 609 END DO 610 END DO 611 612 DO jk = 1, nlay_i 613 DO ji = 1, npti 614 ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0 615 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) 616 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 617 END DO 620 618 END DO 621 619 … … 652 650 653 651 DO jm = nlay_s + 2, nlay_s + nlay_i 654 DO ji = 1, npti655 jk = jm - nlay_s - 1656 ztrid(ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)657 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )658 ztrid(ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk)659 zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)660 END DO652 DO ji = 1, npti 653 jk = jm - nlay_s - 1 654 ztrid(ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 655 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 656 ztrid(ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk) 657 zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 658 END DO 661 659 ENDDO 662 660 663 661 jm = nlay_s + nlay_i + 1 664 662 DO ji = 1, npti 665 !!ice bottom term666 ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)667 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )668 ztrid(ji,jm,3) = 0.0669 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * &670 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) )663 !!ice bottom term 664 ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) 665 ztrid(ji,jm,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 666 ztrid(ji,jm,3) = 0.0 667 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 668 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 671 669 ENDDO 672 670 673 671 674 672 DO ji = 1, npti 675 ! !---------------------! 676 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! 677 ! !---------------------! 678 ! snow interior terms (bottom equation has the same form as the others) 679 DO jm = 3, nlay_s + 1 680 jk = jm - 1 681 ztrid(ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 682 ztrid(ji,jm,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 683 ztrid(ji,jm,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 684 zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 685 END DO 686 687 ! case of only one layer in the ice (ice equation is altered) 688 IF ( nlay_i == 1 ) THEN 689 ztrid(ji,nlay_s+2,3) = 0.0 690 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 673 ! !---------------------! 674 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! 675 ! !---------------------! 676 ! snow interior terms (bottom equation has the same form as the others) 677 DO jm = 3, nlay_s + 1 678 jk = jm - 1 679 ztrid(ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 680 ztrid(ji,jm,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 681 ztrid(ji,jm,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 682 zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 683 END DO 684 685 ! case of only one layer in the ice (ice equation is altered) 686 IF ( nlay_i == 1 ) THEN 687 ztrid(ji,nlay_s+2,3) = 0.0 688 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 689 ENDIF 690 691 jm_min(ji) = 2 692 jm_max(ji) = nlay_i + nlay_s + 1 693 694 ! first layer of snow equation 695 ztrid(ji,2,1) = 0.0 696 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1) 697 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 698 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 699 700 ! !---------------------! 701 ELSE ! cells without snow ! 702 ! !---------------------! 703 jm_min(ji) = nlay_s + 2 704 jm_max(ji) = nlay_i + nlay_s + 1 705 706 ! first layer of ice equation 707 ztrid(ji,jm_min(ji),1) = 0.0 708 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 709 ztrid(ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 710 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) ) 711 712 ! case of only one layer in the ice (surface & ice equations are altered) 713 IF ( nlay_i == 1 ) THEN 714 ztrid(ji,jm_min(ji),1) = 0.0 715 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 716 ztrid(ji,jm_min(ji),3) = 0.0 717 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) & 718 & + qcn_ice_1d(ji) ) 719 ENDIF 720 691 721 ENDIF 692 693 jm_min(ji) = 2 694 jm_max(ji) = nlay_i + nlay_s + 1 695 696 ! first layer of snow equation 697 ztrid(ji,2,1) = 0.0 698 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1) 699 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 700 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * & 701 & ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 702 703 ! !---------------------! 704 ELSE ! cells without snow ! 705 ! !---------------------! 706 707 jm_min(ji) = nlay_s + 2 708 jm_max(ji) = nlay_i + nlay_s + 1 709 710 ! first layer of ice equation 711 ztrid(ji,jm_min(ji),1) = 0.0 712 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 713 ztrid(ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 714 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * & 715 & ( zradab_i(ji,1) + qcn_ice_1d(ji) ) 716 717 ! case of only one layer in the ice (surface & ice equations are altered) 718 IF ( nlay_i == 1 ) THEN 719 ztrid(ji,jm_min(ji),1) = 0.0 720 ztrid(ji,jm_min(ji),2) = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 721 ztrid(ji,jm_min(ji),3) = 0.0 722 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) & 723 & + qcn_ice_1d(ji) ) 724 725 ENDIF 726 727 ENDIF 728 ! 729 zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) 730 zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 731 ! 722 ! 723 zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) 724 zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 725 ! 732 726 END DO 733 727 ! … … 740 734 jm_mint = nlay_i+5 741 735 DO ji = 1, npti 742 jm_mint = MIN(jm_min(ji),jm_mint)743 jm_maxt = MAX(jm_max(ji),jm_maxt)744 END DO 745 736 jm_mint = MIN(jm_min(ji),jm_mint) 737 jm_maxt = MAX(jm_max(ji),jm_maxt) 738 END DO 739 746 740 DO jk = jm_mint+1, jm_maxt 747 DO ji = 1, npti 748 jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 749 zdiagbis(ji,jm) = ztrid(ji,jm,2) - ztrid(ji,jm,1) * ztrid(ji,jm-1,3) / zdiagbis(ji,jm-1) 750 zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 751 END DO 752 END DO 753 754 DO ji = 1, npti 741 DO ji = 1, npti 742 jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 743 zdiagbis(ji,jm) = ztrid(ji,jm,2) - ztrid(ji,jm,1) * ztrid(ji,jm-1,3) / zdiagbis(ji,jm-1) 744 zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 745 END DO 746 END DO 747 755 748 ! ice temperatures 756 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 749 DO ji = 1, npti 750 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 757 751 END DO 758 752 759 753 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 760 DO ji = 1, npti 761 jk = jm - nlay_s - 1 762 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 763 END DO 764 END DO 765 766 DO ji = 1, npti 754 DO ji = 1, npti 755 jk = jm - nlay_s - 1 756 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 757 END DO 758 END DO 759 767 760 ! snow temperatures 768 IF( h_s_1d(ji) > 0._wp ) THEN769 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) &770 &/ zdiagbis(ji,nlay_s+1)771 ENDIF761 DO ji = 1, npti 762 IF( h_s_1d(ji) > 0._wp ) THEN 763 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) 764 ENDIF 772 765 END DO 773 766 ! … … 779 772 zdti_max = 0._wp 780 773 781 DO jk =1, nlay_s782 DO ji = 1, npti783 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )784 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )785 END DO786 END DO 787 788 DO jk =1, nlay_i789 DO ji = 1, npti790 ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0791 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )792 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )793 END DO774 DO jk = 1, nlay_s 775 DO ji = 1, npti 776 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 777 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 778 END DO 779 END DO 780 781 DO jk = 1, nlay_i 782 DO ji = 1, npti 783 ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0 784 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) 785 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 786 END DO 794 787 END DO 795 788 … … 815 808 ! 816 809 DO ji = 1, npti 817 ! ! surface ice conduction flux810 ! ! surface ice conduction flux 818 811 fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & 819 812 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji)) 820 ! ! bottom ice conduction flux813 ! ! bottom ice conduction flux 821 814 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 822 815 END DO … … 892 885 ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes 893 886 !--------------------------------------------------------------------------------------- 894 ! effective conductivity and 1st layer temperature ( Jules coupling)887 ! effective conductivity and 1st layer temperature (needed by Met Office) 895 888 DO ji = 1, npti 896 IF (h_s_1d(ji) > 0.1) THEN897 889 IF( h_s_1d(ji) > 0.1_wp ) THEN 890 cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0) 898 891 ELSE 899 IF (h_i_1d(ji) > 0.1) THEN900 901 902 cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / 0.1903 892 IF( h_i_1d(ji) > 0.1_wp ) THEN 893 cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0) 894 ELSE 895 cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / 0.1_wp 896 ENDIF 904 897 ENDIF 905 t1_ice_1d (ji) = ( isnow(ji) * t_s_1d (ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d (ji,1))898 t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) 906 899 END DO 907 900 ! 908 901 IF ( k_jules == np_zdf_jules_SND ) THEN ! --- Jules coupling in "SND" mode 909 !910 902 ! Restore temperatures to their initial values 911 903 t_s_1d (1:npti,:) = ztsold(1:npti,:) 912 904 t_i_1d (1:npti,:) = ztiold(1:npti,:) 913 905 qcn_ice_1d(1:npti) = fc_su (1:npti) 914 !915 906 ENDIF 916 907 !
Note: See TracChangeset
for help on using the changeset viewer.
