Changeset 14072 for NEMO/trunk/src/ICE/icethd_zdf_bl99.F90
- Timestamp:
- 2020-12-04T08:48:38+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icethd_zdf_bl99.F90 (modified) (37 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icethd_zdf_bl99.F90
r14005 r14072 2 2 !!====================================================================== 3 3 !! *** MODULE icethd_zdf_BL99 *** 4 !! sea-ice: vertical heat diffusion in sea ice (computation of temperatures) 4 !! sea-ice: vertical heat diffusion in sea ice (computation of temperatures) 5 5 !!====================================================================== 6 6 !! History : ! 2003-02 (M. Vancoppenolle) original 1D code … … 15 15 !!---------------------------------------------------------------------- 16 16 USE dom_oce ! ocean space and time domain 17 USE phycst ! physical constants (ocean directory) 17 USE phycst ! physical constants (ocean directory) 18 18 USE ice ! sea-ice: variables 19 19 USE ice1D ! sea-ice: thermodynamics variables … … 44 44 !! 45 45 !! ** Method : solves the heat equation diffusion with a Neumann boundary 46 !! condition at the surface and a Dirichlet one at the bottom. 46 !! condition at the surface and a Dirichlet one at the bottom. 47 47 !! Solar radiation is partially absorbed into the ice. 48 !! The specific heat and thermal conductivities depend on ice 49 !! salinity and temperature to take into account brine pocket 48 !! The specific heat and thermal conductivities depend on ice 49 !! salinity and temperature to take into account brine pocket 50 50 !! melting. The numerical scheme is an iterative Crank-Nicolson 51 51 !! on a non-uniform multilayer grid in the ice and snow system. … … 91 91 REAL(wp) :: zbeta = 0.117_wp ! for thermal conductivity (could be 0.13) 92 92 REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity 93 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C 94 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature 95 REAL(wp) :: zhs_ssl = 0.03_wp ! surface scattering layer in the snow 93 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C 94 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature 95 REAL(wp) :: zhs_ssl = 0.03_wp ! surface scattering layer in the snow 96 96 REAL(wp) :: zhi_ssl = 0.10_wp ! surface scattering layer in the ice 97 97 REAL(wp) :: zh_min = 1.e-3_wp ! minimum ice/snow thickness for conduction 98 98 REAL(wp) :: ztmelts ! ice melting temperature 99 REAL(wp) :: zdti_max ! current maximal error on temperature 99 REAL(wp) :: zdti_max ! current maximal error on temperature 100 100 REAL(wp) :: zcpi ! Ice specific heat 101 101 REAL(wp) :: zhfx_err, zdq ! diag errors on heat … … 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(jpij) :: za_s_fra ! ice fraction covered by snow 130 REAL(wp), DIMENSION(jpij) :: isnow ! snow presence (1) or not (0) 131 REAL(wp), DIMENSION(jpij) :: isnow_comb ! snow presence for met-office 129 REAL(wp), DIMENSION(jpij) :: za_s_fra ! ice fraction covered by snow 130 REAL(wp), DIMENSION(jpij) :: isnow ! snow presence (1) or not (0) 131 REAL(wp), DIMENSION(jpij) :: isnow_comb ! snow presence for met-office 132 132 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zindterm ! 'Ind'ependent term 133 133 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zindtbis ! Temporary 'ind'ependent term … … 139 139 REAL(wp) :: zhe ! dummy factor 140 140 REAL(wp) :: zcnd_i ! mean sea ice thermal conductivity 141 !!------------------------------------------------------------------ 141 !!------------------------------------------------------------------ 142 142 143 143 ! --- diag error on heat diffusion - PART 1 --- ! 144 144 DO ji = 1, npti 145 145 zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & 146 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 146 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 147 147 END DO 148 148 149 149 ! calculate ice fraction covered by snow for radiation 150 150 CALL ice_var_snwfra( h_s_1d(1:npti), za_s_fra(1:npti) ) 151 151 152 152 !------------------ 153 153 ! 1) Initialization … … 155 155 ! 156 156 ! extinction radiation in the snow 157 IF ( nn_qtrice == 0 ) THEN ! constant 157 IF ( nn_qtrice == 0 ) THEN ! constant 158 158 zraext_s(1:npti) = rn_kappa_s 159 159 ELSEIF( nn_qtrice == 1 ) THEN ! depends on melting/freezing conditions … … 166 166 DO ji = 1, npti 167 167 ! ice thickness 168 IF( h_i_1d(ji) > 0._wp ) THEN 168 IF( h_i_1d(ji) > 0._wp ) THEN 169 169 zh_i (ji) = MAX( zh_min , h_i_1d(ji) ) * r1_nlay_i ! set a minimum thickness for conduction 170 170 z1_h_i(ji) = 1._wp / zh_i(ji) ! it must be very small … … 198 198 ztsuold (1:npti) = t_su_1d(1:npti) ! surface temperature initial value 199 199 t_su_1d (1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err ) ! required to leave the choice between melting or not 200 zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti) ! derivative of incoming nonsolar flux 200 zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti) ! derivative of incoming nonsolar flux 201 201 zqns_ice_b (1:npti) = qns_ice_1d(1:npti) ! store previous qns_ice_1d value 202 202 ! … … 221 221 ! 222 222 zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * za_s_fra(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - za_s_fra(1:npti) ) 223 DO jk = 1, nlay_i 223 DO jk = 1, nlay_i 224 224 DO ji = 1, npti 225 225 ! ! radiation transmitted below the layer-th ice layer … … 227 227 & * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zh_min ) ) & 228 228 & + ( 1._wp - za_s_fra(ji) ) * qtr_ice_top_1d(ji) & ! part snow free 229 & * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) ) 229 & * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) ) 230 230 ! ! radiation absorbed by the layer-th ice layer 231 231 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) … … 288 288 DO ji = 1, npti 289 289 IF ( .NOT. l_T_converged(ji) ) & 290 ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) ) 290 ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) ) 291 291 END DO 292 292 ! … … 401 401 zdiagbis(1:npti,:) = 0._wp 402 402 403 DO jm = nlay_s + 2, nlay_s + nlay_i 403 DO jm = nlay_s + 2, nlay_s + nlay_i 404 404 DO ji = 1, npti 405 405 jk = jm - nlay_s - 1 … … 414 414 DO ji = 1, npti 415 415 ! ice bottom term 416 ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1) 416 ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1) 417 417 ztrid (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 ) 418 418 ztrid (ji,jm,3) = 0._wp 419 419 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 420 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 420 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 421 421 END DO 422 422 … … 433 433 zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 434 434 END DO 435 435 436 436 ! case of only one layer in the ice (ice equation is altered) 437 437 IF( nlay_i == 1 ) THEN 438 438 ztrid (ji,nlay_s+2,3) = 0._wp 439 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 440 ENDIF 441 439 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 440 ENDIF 441 442 442 IF( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting 443 443 444 444 jm_min(ji) = 1 445 445 jm_max(ji) = nlay_i + nlay_s + 1 446 446 447 447 ! surface equation 448 448 ztrid (ji,1,1) = 0._wp … … 450 450 ztrid (ji,1,3) = zg1s * zkappa_s(ji,0) 451 451 zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 452 452 453 453 ! first layer of snow equation 454 454 ztrid (ji,2,1) = - zeta_s(ji,1) * zkappa_s(ji,0) * zg1s … … 456 456 ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1) 457 457 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 458 458 459 459 ELSE !-- case 2 : surface is melting 460 460 ! 461 461 jm_min(ji) = 2 462 462 jm_max(ji) = nlay_i + nlay_s + 1 463 463 464 464 ! first layer of snow equation 465 465 ztrid (ji,2,1) = 0._wp 466 466 ztrid (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 467 ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1) 468 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 467 ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1) 468 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 469 469 ENDIF 470 470 ! !---------------------! … … 476 476 jm_min(ji) = nlay_s + 1 477 477 jm_max(ji) = nlay_i + nlay_s + 1 478 479 ! surface equation 478 479 ! surface equation 480 480 ztrid (ji,jm_min(ji),1) = 0._wp 481 ztrid (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * zg1 481 ztrid (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * zg1 482 482 ztrid (ji,jm_min(ji),3) = zkappa_i(ji,0) * zg1 483 483 zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 484 484 485 485 ! first layer of ice equation 486 486 ztrid (ji,jm_min(ji)+1,1) = - zeta_i(ji,1) * zkappa_i(ji,0) * zg1 487 487 ztrid (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 488 ztrid (ji,jm_min(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1) 489 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 490 488 ztrid (ji,jm_min(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1) 489 zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 490 491 491 ! case of only one layer in the ice (surface & ice equations are altered) 492 492 IF( nlay_i == 1 ) THEN … … 499 499 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)) 500 500 ENDIF 501 501 502 502 ELSE !-- case 2 : surface is melting 503 503 504 504 jm_min(ji) = nlay_s + 2 505 505 jm_max(ji) = nlay_i + nlay_s + 1 506 506 507 507 ! first layer of ice equation 508 508 ztrid (ji,jm_min(ji),1) = 0._wp 509 ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 509 ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 510 510 ztrid (ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 511 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji)) 512 511 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji)) 512 513 513 ! case of only one layer in the ice (surface & ice equations are altered) 514 514 IF( nlay_i == 1 ) THEN … … 519 519 & + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2._wp 520 520 ENDIF 521 521 522 522 ENDIF 523 523 ENDIF … … 540 540 !!$ END DO 541 541 !!$ !!clem SNWLAY => check why LIM1D does not get this loop. Is nlay_i+5 correct? 542 !!$ 542 !!$ 543 543 !!$ DO jk = jm_mint+1, jm_maxt 544 544 !!$ DO ji = 1, npti … … 574 574 END DO 575 575 576 ! snow temperatures 576 ! snow temperatures 577 577 DO ji = 1, npti 578 578 ! Variables used after iterations … … 589 589 END DO 590 590 END DO 591 591 592 592 ! surface temperature 593 593 DO ji = 1, npti … … 628 628 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 629 629 END DO 630 630 631 631 ! convergence test 632 632 IF( ln_zdf_chkcvg ) THEN … … 665 665 zdiagbis(1:npti,:) = 0._wp 666 666 667 DO jm = nlay_s + 2, nlay_s + nlay_i 667 DO jm = nlay_s + 2, nlay_s + nlay_i 668 668 DO ji = 1, npti 669 669 jk = jm - nlay_s - 1 … … 678 678 DO ji = 1, npti 679 679 ! ice bottom term 680 ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1) 680 ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1) 681 681 ztrid (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 ) 682 682 ztrid (ji,jm,3) = 0._wp 683 683 zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * & 684 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 684 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 685 685 ENDDO 686 686 … … 697 697 zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 698 698 END DO 699 699 700 700 ! case of only one layer in the ice (ice equation is altered) 701 701 IF ( nlay_i == 1 ) THEN 702 702 ztrid (ji,nlay_s+2,3) = 0._wp 703 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 704 ENDIF 705 703 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 704 ENDIF 705 706 706 jm_min(ji) = 2 707 707 jm_max(ji) = nlay_i + nlay_s + 1 708 708 709 709 ! first layer of snow equation 710 710 ztrid (ji,2,1) = 0._wp 711 711 ztrid (ji,2,2) = 1._wp + zeta_s(ji,1) * zkappa_s(ji,1) 712 ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1) 713 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 714 712 ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1) 713 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 714 715 715 ! !---------------------! 716 716 ELSE ! cells without snow ! … … 718 718 jm_min(ji) = nlay_s + 2 719 719 jm_max(ji) = nlay_i + nlay_s + 1 720 720 721 721 ! first layer of ice equation 722 722 ztrid (ji,jm_min(ji),1) = 0._wp 723 ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1) 723 ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1) 724 724 ztrid (ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 725 725 zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) ) 726 726 727 727 ! case of only one layer in the ice (surface & ice equations are altered) 728 728 IF( nlay_i == 1 ) THEN … … 733 733 & ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) + qcn_ice_1d(ji) ) 734 734 ENDIF 735 735 736 736 ENDIF 737 737 ! … … 752 752 !!$ jm_maxt = MAX(jm_max(ji),jm_maxt) 753 753 !!$ END DO 754 !!$ 754 !!$ 755 755 !!$ DO jk = jm_mint+1, jm_maxt 756 756 !!$ DO ji = 1, npti … … 786 786 END DO 787 787 END DO 788 789 ! snow temperatures 788 789 ! snow temperatures 790 790 DO ji = 1, npti 791 791 ! Variables used after iterations … … 823 823 824 824 DO jk = 1, nlay_i 825 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 825 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 826 826 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 827 827 zdti_max = MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) … … 885 885 ! 886 886 DO ji = 1, npti 887 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 887 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 888 888 END DO 889 889 ! … … 893 893 ! 894 894 IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_ON ) THEN 895 896 CALL ice_var_enthalpy 897 895 896 CALL ice_var_enthalpy 897 898 898 ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 899 899 DO ji = 1, npti 900 900 zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & 901 901 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 902 902 903 903 IF( k_cnd == np_cnd_OFF ) THEN 904 904 905 905 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 906 906 zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji) & … … 910 910 & + zdq * r1_Dt_ice ) * a_i_1d(ji) 911 911 ENDIF 912 912 913 913 ELSEIF( k_cnd == np_cnd_ON ) THEN 914 914 915 915 zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji) & 916 916 & + zdq * r1_Dt_ice ) * a_i_1d(ji) 917 917 918 918 ENDIF 919 919 ! … … 921 921 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 922 922 ! 923 ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2 923 ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2 924 924 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_Dt_ice * a_i_1d(ji) 925 925 ! … … 952 952 ! --- SIMIP diagnostics 953 953 ! 954 DO ji = 1, npti 954 DO ji = 1, npti 955 955 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 956 956 IF( h_s_1d(ji) >= zhs_ssl ) THEN
Note: See TracChangeset
for help on using the changeset viewer.
