Changeset 13148
- Timestamp:
- 2020-06-23T22:31:14+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90
r12919 r13148 95 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 REAL(wp) :: zh_min = 1.e-3_wp ! minimum ice/snow thickness for conduction 97 98 REAL(wp) :: ztmelts ! ice melting temperature 98 99 REAL(wp) :: zdti_max ! current maximal error on temperature … … 123 124 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zkappa_s ! Kappa factor in the snow 124 125 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeta_s ! Eta factor in the snow 126 REAL(wp), DIMENSION(jpij) :: zkappa_comb ! Combined snow and ice surface conductivity 125 127 REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindterm ! 'Ind'ependent term 126 128 REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindtbis ! Temporary 'ind'ependent term … … 130 132 REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 131 133 REAL(wp), DIMENSION(jpij) :: za_s_fra ! ice fraction covered by snow 134 REAL(wp), DIMENSION(jpij) :: isnow ! snow presence (1) or not (0) 132 135 ! 133 136 ! Mono-category … … 143 146 END DO 144 147 145 CALL ice_var_snwfra( h_s_1d(1:npti), za_s_fra(1:npti) ) ! calculate ice fraction covered by snow 148 ! calculate ice fraction covered by snow for radiation 149 CALL ice_var_snwfra( h_s_1d(1:npti), za_s_fra(1:npti) ) 146 150 147 151 !------------------ … … 158 162 ENDIF 159 163 ! 160 ! layersthicknesses164 ! thicknesses 161 165 DO ji = 1, npti 162 zh_s(ji) = MAX( zhs_ssl , h_s_1d(ji) ) * r1_nlay_s ! set a minimum snw thickness for conduction 163 zh_i(ji) = MAX( zhi_ssl , h_i_1d(ji) ) * r1_nlay_i ! set a minimum ice thickness for conduction 166 ! ice thickness 167 IF( h_i_1d(ji) > 0._wp ) THEN 168 zh_i (ji) = MAX( zh_min , h_i_1d(ji) ) * r1_nlay_i ! set a minimum thickness for conduction 169 z1_h_i(ji) = 1._wp / zh_i(ji) ! it must be very small 170 ELSE 171 zh_i (ji) = 0._wp 172 z1_h_i(ji) = 0._wp 173 ENDIF 174 ! snow thickness 175 IF( h_s_1d(ji) > 0._wp ) THEN 176 zh_s (ji) = MAX( zh_min , h_s_1d(ji) ) * r1_nlay_s ! set a minimum thickness for conduction 177 z1_h_s(ji) = 1._wp / zh_s(ji) ! it must be very small 178 isnow (ji) = 1._wp 179 ELSE 180 zh_s (ji) = 0._wp 181 z1_h_s(ji) = 0._wp 182 isnow (ji) = 0._wp 183 ENDIF 164 184 END DO 165 ! 166 WHERE( h_i_1d(1:npti) >= zhi_ssl ) ; z1_h_i(1:npti) = 1._wp / zh_i(1:npti) 167 ELSEWHERE ; z1_h_i(1:npti) = 0._wp ! put 0 if ice thick < scattering layer 168 END WHERE 169 ! 170 WHERE( h_s_1d(1:npti) >= zhs_ssl ) ; z1_h_s(1:npti) = 1._wp / zh_s(1:npti) 171 ELSEWHERE ; z1_h_s(1:npti) = 0._wp ! put 0 if snw thick < scattering layer 172 END WHERE 185 ! clem: we should apply correction on snow thickness to take into account snow fraction 186 ! it must be a distribution, so it is a bit complicated 173 187 ! 174 188 ! Store initial temperatures and non solar heat fluxes 175 189 IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN 176 !177 190 ztsub (1:npti) = t_su_1d(1:npti) ! surface temperature at iteration n-1 178 191 ztsuold (1:npti) = t_su_1d(1:npti) ! surface temperature initial value … … 200 213 END DO 201 214 ! 202 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))215 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) ) 203 216 DO jk = 1, nlay_i 204 217 DO ji = 1, npti 205 218 ! ! radiation transmitted below the layer-th ice layer 206 219 zradtr_i(ji,jk) = za_s_fra(ji) * zradtr_s(ji,nlay_s) & ! part covered by snow 207 & * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) &220 & * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zh_min ) ) & 208 221 & + ( 1._wp - za_s_fra(ji) ) * qtr_ice_top_1d(ji) & ! part snow free 209 & * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) ) 210 211 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) ) 222 & * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) ) 212 223 ! ! radiation absorbed by the layer-th ice layer 213 224 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) … … 304 315 END DO 305 316 DO ji = 1, npti ! Snow-ice interface 306 IF ( .NOT. l_T_converged(ji) ) THEN 307 IF( h_s_1d(ji) >= zhs_ssl ) THEN 308 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / & 309 & ( 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) ) 310 ELSE 311 zkappa_s(ji,nlay_s) = 0._wp 312 ENDIF 313 ENDIF 317 IF ( .NOT. l_T_converged(ji) ) & 318 zkappa_s(ji,nlay_s) = isnow(ji) * zghe(ji) * rn_cnd_s * ztcond_i(ji,0) & 319 & / ( 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) ) 314 320 END DO 315 321 … … 324 330 END DO 325 331 DO ji = 1, npti ! Snow-ice interface 326 IF ( .NOT. l_T_converged(ji) ) & 327 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * za_s_fra(ji) + zkappa_i(ji,0) * (1._wp - za_s_fra(ji)) 332 IF ( .NOT. l_T_converged(ji) ) THEN 333 ! Calculate combined surface snow and ice conductivity to pass through the coupler 334 zkappa_comb(ji) = isnow(ji) * zkappa_s(ji,0) + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) 335 ! If there is snow then use the same snow-ice interface conductivity for the top layer of ice 336 IF( h_s_1d(ji) > 0._wp ) zkappa_i(ji,0) = zkappa_s(ji,nlay_s) 337 ENDIF 328 338 END DO 329 339 ! … … 334 344 DO ji = 1, npti 335 345 zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 336 zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi )346 zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / zcpi 337 347 END DO 338 348 END DO … … 558 568 ztsub(ji) = t_su_1d(ji) 559 569 IF( t_su_1d(ji) < rt0 ) THEN 560 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * &561 & ( za_s_fra(ji) * t_s_1d(ji,1) + (1._wp - za_s_fra(ji)) *t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))570 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 571 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 562 572 ENDIF 563 573 ENDIF 564 574 END DO 575 !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1) 565 576 ! 566 577 !-------------------------------------------------------------- … … 575 586 576 587 IF ( .NOT. l_T_converged(ji) ) THEN 588 577 589 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 578 590 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 579 591 580 t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 581 zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 592 IF( h_s_1d(ji) > 0._wp ) THEN 593 DO jk = 1, nlay_s 594 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 595 zdti_max = MAX ( zdti_max , ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 596 END DO 597 ENDIF 582 598 583 599 DO jk = 1, nlay_i … … 586 602 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 587 603 END DO 604 588 605 ! convergence test 589 606 IF( ln_zdf_chkcvg ) THEN … … 745 762 ENDIF 746 763 END DO 764 !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1) 747 765 ! 748 766 !-------------------------------------------------------------- … … 757 775 758 776 IF ( .NOT. l_T_converged(ji) ) THEN 759 ! t_s 760 t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 761 zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 762 ! t_i 777 778 IF( h_s_1d(ji) > 0._wp ) THEN 779 DO jk = 1, nlay_s 780 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 781 zdti_max = MAX ( zdti_max , ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 782 END DO 783 ENDIF 784 763 785 DO jk = 1, nlay_i 764 786 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 … … 766 788 zdti_max = MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 767 789 END DO 790 768 791 ! convergence test 769 792 IF( ln_zdf_chkcvg ) THEN … … 789 812 ! bottom ice conduction flux 790 813 DO ji = 1, npti 791 qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1 814 qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1 * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 792 815 END DO 793 816 ! surface ice conduction flux … … 795 818 ! 796 819 DO ji = 1, npti 797 qcn_ice_top_1d(ji) = - za_s_fra(ji) * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )&798 & - ( 1._wp - za_s_fra(ji) ) * zkappa_i(ji,0) * zg1 * ( t_i_1d(ji,1) - t_su_1d(ji) )820 qcn_ice_top_1d(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) & 821 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * ( t_i_1d(ji,1) - t_su_1d(ji) ) 799 822 END DO 800 823 ! … … 810 833 ! 811 834 DO ji = 1, npti 812 t_su_1d(ji) = ( qcn_ice_top_1d(ji) & ! calculate surface temperature 813 & + za_s_fra(ji) * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) & 814 & + (1._wp - za_s_fra(ji)) * zkappa_i(ji,0) * zg1 * t_i_1d(ji,1) & 815 & ) / MAX( epsi10, za_s_fra(ji) * zkappa_s(ji,0) * zg1s + (1._wp - za_s_fra(ji)) * zkappa_i(ji,0) * zg1 ) 835 t_su_1d(ji) = ( qcn_ice_top_1d(ji) + isnow(ji) * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) + & 836 & ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * t_i_1d(ji,1) ) & 837 & / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 ) 816 838 t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp ) ! cap t_su 817 839 END DO … … 873 895 ! this is a conductivity at mid-layer, hence the factor 2 874 896 DO ji = 1, npti 875 IF( h_i_1d(ji) >= zhi_ssl ) THEN ! zkappa_i=0 when h_i<zhi_ssl (but ztcond_i/=0) 897 IF( h_i_1d(ji) >= zhi_ssl ) THEN 898 !!clem metO cnd_ice_1d(ji) = 2._wp * zkappa_comb(ji) 876 899 cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0) 877 900 ELSE 878 901 cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / zhi_ssl ! cnd_ice is capped by: cond_i/zhi_ssl 879 902 ENDIF 880 t1_ice_1d(ji) = za_s_fra(ji) * t_s_1d(ji,1) + ( 1._wp - za_s_fra(ji) ) * t_i_1d(ji,1)903 t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) 881 904 END DO 882 905 ! … … 893 916 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 894 917 IF( h_s_1d(ji) >= zhs_ssl ) THEN 895 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / & 896 & ( rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) ) 918 t_si_1d(ji) = ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,1) & 919 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) & 920 & / ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i & 921 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s ) 897 922 ELSE 898 923 t_si_1d(ji) = t_su_1d(ji)
Note: See TracChangeset
for help on using the changeset viewer.