- Timestamp:
- 2020-05-07T16:11:23+02:00 (4 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
r12854 r12890 94 94 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C 95 95 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature 96 REAL(wp) :: zhs_min = 0.1_wp ! minimum snow thickness for conductivity calculation 96 REAL(wp) :: zhs_ssl = 0.03_wp ! surface scattering layer in the snow 97 REAL(wp) :: zhi_ssl = 0.10_wp ! surface scattering layer in the ice 97 98 REAL(wp) :: ztmelts ! ice melting temperature 98 99 REAL(wp) :: zdti_max ! current maximal error on temperature 99 100 REAL(wp) :: zcpi ! Ice specific heat 100 101 REAL(wp) :: zhfx_err, zdq ! diag errors on heat 101 REAL(wp) :: zfac ! dummy factor 102 ! 103 REAL(wp), DIMENSION(jpij) :: isnow ! fraction of sea ice that is snow covered (for thermodynamic use only) 102 ! 104 103 REAL(wp), DIMENSION(jpij) :: ztsub ! surface temperature at previous iteration 105 104 REAL(wp), DIMENSION(jpij) :: zh_i, z1_h_i ! ice layer thickness … … 150 149 !------------------ 151 150 DO ji = 1, npti 152 153 ! If the snow thickness drops below zhs_min then reduce the snow fraction instead 154 !!IF( h_s_1d(ji) < zhs_min ) THEN 155 !! isnow(ji) = h_s_1d(ji) / zhs_min 156 !! zh_s(ji) = zhs_min * r1_nlay_s 157 !!ELSE 158 !! isnow(ji) = 1.0_wp 159 !! zh_s(ji) = h_s_1d(ji) * r1_nlay_s 160 !!END IF 161 isnow(ji) = za_s_fra(ji) !!clem: 2 variables for the same thing 162 zh_s(ji) = h_s_1d(ji) * r1_nlay_s 163 164 165 ! layer thickness 166 zh_i(ji) = h_i_1d(ji) * r1_nlay_i 167 151 zh_s(ji) = MAX( zhs_ssl , h_s_1d(ji) ) * r1_nlay_s ! set a minimum snw thickness for conduction 152 zh_i(ji) = MAX( zhi_ssl , h_i_1d(ji) ) * r1_nlay_i ! set a minimum ice thickness for conduction 168 153 END DO 169 154 ! 170 WHERE( zh_i(1:npti) >= epsi10) ; z1_h_i(1:npti) = 1._wp / zh_i(1:npti)171 ELSEWHERE ; z1_h_i(1:npti) = 0._wp155 WHERE( h_i_1d(1:npti) >= zhi_ssl ) ; z1_h_i(1:npti) = 1._wp / zh_i(1:npti) 156 ELSEWHERE ; z1_h_i(1:npti) = 0._wp ! put 0 if ice thick < scattering layer 172 157 END WHERE 173 158 ! 174 WHERE( zh_s(1:npti) >= epsi10) ; z1_h_s(1:npti) = 1._wp / zh_s(1:npti)175 ELSEWHERE ; z1_h_s(1:npti) = 0._wp159 WHERE( h_s_1d(1:npti) >= zhs_ssl ) ; z1_h_s(1:npti) = 1._wp / zh_s(1:npti) 160 ELSEWHERE ; z1_h_s(1:npti) = 0._wp ! put 0 if snw thick < scattering layer 176 161 END WHERE 177 162 ! … … 198 183 DO ji = 1, npti 199 184 ! ! radiation transmitted below the layer-th snow layer 200 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * h_s_1d(ji) * r1_nlay_s * REAL(jk) )185 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * MAX( 0._wp, zh_s(ji) * REAL(jk) - zhs_ssl ) ) 201 186 ! ! radiation absorbed by the layer-th snow layer 202 187 zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) … … 204 189 END DO 205 190 ! 206 zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - isnow(1:npti))191 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)) 207 192 DO jk = 1, nlay_i 208 193 DO ji = 1, npti 209 194 ! ! radiation transmitted below the layer-th ice layer 210 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )195 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) ) 211 196 ! ! radiation absorbed by the layer-th ice layer 212 197 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) … … 216 201 qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti,nlay_i) ! record radiation transmitted below the ice 217 202 ! 218 iconv 203 iconv = 0 ! number of iterations 219 204 ! 220 205 l_T_converged(:) = .FALSE. … … 243 228 DO ji = 1, npti 244 229 ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 245 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )230 & MIN( -epsi10, 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 ) 246 231 END DO 247 232 END DO … … 251 236 DO ji = 1, npti 252 237 ztcond_i_cp(ji,0) = rcnd_i + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 253 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 )238 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 254 239 ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 255 & - 0.011_wp * ( t_bo_1d(ji) - rt0 )240 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 256 241 END DO 257 242 DO jk = 1, nlay_i-1 258 243 DO ji = 1, npti 259 ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / 260 & MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )&261 & - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d(ji,jk+1) ) - rt0 )244 ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 245 & MIN( -epsi10, 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 ) & 246 & - 0.011_wp * ( 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 ) 262 247 END DO 263 248 END DO … … 304 289 DO ji = 1, npti ! Snow-ice interface 305 290 IF ( .NOT. l_T_converged(ji) ) THEN 306 zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )307 IF( zfac > epsi10 ) THEN308 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac291 IF( h_s_1d(ji) >= zhs_ssl ) THEN 292 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / & 293 & ( 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) ) 309 294 ELSE 310 295 zkappa_s(ji,nlay_s) = 0._wp … … 324 309 DO ji = 1, npti ! Snow-ice interface 325 310 IF ( .NOT. l_T_converged(ji) ) & 326 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji))311 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * za_s_fra(ji) + zkappa_i(ji,0) * (1._wp - za_s_fra(ji)) 327 312 END DO 328 313 ! … … 558 543 IF( t_su_1d(ji) < rt0 ) THEN 559 544 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 560 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji)) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))545 & ( 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)) 561 546 ENDIF 562 547 ENDIF … … 790 775 ! 791 776 DO ji = 1, npti 792 qcn_ice_top_1d(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) &793 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * ( t_i_1d(ji,1) - t_su_1d(ji) )777 qcn_ice_top_1d(ji) = - za_s_fra(ji) * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) & 778 & - ( 1._wp - za_s_fra(ji) ) * zkappa_i(ji,0) * zg1 * ( t_i_1d(ji,1) - t_su_1d(ji) ) 794 779 END DO 795 780 ! … … 806 791 DO ji = 1, npti 807 792 t_su_1d(ji) = ( qcn_ice_top_1d(ji) & ! calculate surface temperature 808 & + isnow(ji)* zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) &809 & + ( 1._wp - isnow(ji)) * zkappa_i(ji,0) * zg1 * t_i_1d(ji,1) &810 & ) / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji)) * zkappa_i(ji,0) * zg1 )793 & + za_s_fra(ji) * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) & 794 & + (1._wp - za_s_fra(ji)) * zkappa_i(ji,0) * zg1 * t_i_1d(ji,1) & 795 & ) / MAX( epsi10, za_s_fra(ji) * zkappa_s(ji,0) * zg1s + (1._wp - za_s_fra(ji)) * zkappa_i(ji,0) * zg1 ) 811 796 t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp ) ! cap t_su 812 797 END DO … … 866 851 !-------------------------------------------------------------------- 867 852 ! effective conductivity and 1st layer temperature (needed by Met Office) 853 ! this is a conductivity at mid-layer, hence the factor 2 868 854 DO ji = 1, npti 869 IF( h_i_1d(ji) > 0.1_wp ) THEN855 IF( h_i_1d(ji) >= zhi_ssl ) THEN ! zkappa_i=0 when h_i<zhi_ssl (but ztcond_i/=0) 870 856 cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0) 871 857 ELSE 872 cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp ! cnd_ice is capped by: cond_i/0.1m858 cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / zhi_ssl ! cnd_ice is capped by: cond_i/zhi_ssl 873 859 ENDIF 874 t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)860 t1_ice_1d(ji) = za_s_fra(ji) * t_s_1d(ji,1) + ( 1._wp - za_s_fra(ji) ) * t_i_1d(ji,1) 875 861 END DO 876 862 ! … … 886 872 DO ji = 1, npti 887 873 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 888 zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 889 IF( h_s_1d(ji) >= zhs_min ) THEN !!clem change that 890 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + & 891 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac ) 874 IF( h_s_1d(ji) >= zhs_ssl ) THEN 875 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) ) / & 876 & ( rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) ) 892 877 ELSE 893 878 t_si_1d(ji) = t_su_1d(ji)
Note: See TracChangeset
for help on using the changeset viewer.
