- Timestamp:
- 2018-11-07T18:25:49+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9866_HPC_03_globcom/src/ICE/icethd_zdf_bl99.F90
- Property svn:keywords set to Id
r9656 r10288 31 31 !!---------------------------------------------------------------------- 32 32 !! NEMO/ICE 4.0 , NEMO Consortium (2018) 33 !! $Id : icethd_zdf.F90 8420 2017-08-08 12:18:46Z clem$34 !! Software governed by the CeCILL licen ce (./LICENSE)33 !! $Id$ 34 !! Software governed by the CeCILL license (see ./LICENSE) 35 35 !!---------------------------------------------------------------------- 36 36 CONTAINS … … 93 93 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature 94 94 REAL(wp) :: zhs_min = 0.01_wp ! minimum snow thickness for conductivity calculation 95 REAL(wp) :: ztmelt _i! ice melting temperature95 REAL(wp) :: ztmelts ! ice melting temperature 96 96 REAL(wp) :: zdti_max ! current maximal error on temperature 97 97 REAL(wp) :: zcpi ! Ice specific heat … … 178 178 !------------- 179 179 ! --- Transmission/absorption of solar radiation in the ice --- ! 180 zradtr_s(1:npti,0) = q sr_ice_tr_1d(1:npti)180 zradtr_s(1:npti,0) = qtr_ice_top_1d(1:npti) 181 181 DO jk = 1, nlay_s 182 182 DO ji = 1, npti … … 188 188 END DO 189 189 ! 190 zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + q sr_ice_tr_1d(1:npti) * ( 1._wp - isnow(1:npti) )190 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 191 DO jk = 1, nlay_i 192 192 DO ji = 1, npti … … 198 198 END DO 199 199 ! 200 ftr_ice_1d(1:npti) = zradtr_i(1:npti,nlay_i) ! record radiation transmitted below the ice200 qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti,nlay_i) ! record radiation transmitted below the ice 201 201 ! 202 202 iconv = 0 ! number of iterations … … 217 217 ! 218 218 DO ji = 1, npti 219 ztcond_i(ji,0) = rc dic+ zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )220 ztcond_i(ji,nlay_i) = rc dic+ zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )219 ztcond_i(ji,0) = rcnd_i + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 220 ztcond_i(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) 221 221 END DO 222 222 DO jk = 1, nlay_i-1 223 223 DO ji = 1, npti 224 ztcond_i(ji,jk) = rc dic+ zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &225 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )224 ztcond_i(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 225 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 226 226 END DO 227 227 END DO … … 230 230 ! 231 231 DO ji = 1, npti 232 ztcond_i(ji,0) = rc dic+ 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) &233 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 )234 ztcond_i(ji,nlay_i) = rc dic+ 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) &235 & - 0.011_wp * ( t_bo_1d(ji) - rt0 )232 ztcond_i(ji,0) = rcnd_i + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 233 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 234 ztcond_i(ji,nlay_i) = rcnd_i + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 235 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 236 236 END DO 237 237 DO jk = 1, nlay_i-1 238 238 DO ji = 1, npti 239 ztcond_i(ji,jk) = rc dic+ 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &240 & MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) &241 & - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )239 ztcond_i(ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 240 & MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) & 241 & - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 242 242 END DO 243 243 END DO … … 299 299 DO jk = 1, nlay_i 300 300 DO ji = 1, npti 301 zcpi = cpic+ zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )302 zeta_i(ji,jk) = rdt_ice * r1_rhoi c* z1_h_i(ji) / MAX( epsi10, zcpi )301 zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 302 zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi ) 303 303 END DO 304 304 END DO … … 306 306 DO jk = 1, nlay_s 307 307 DO ji = 1, npti 308 zeta_s(ji,jk) = rdt_ice * r1_rhos n * r1_cpic* z1_h_s(ji)308 zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji) 309 309 END DO 310 310 END DO … … 330 330 331 331 DO ji = 1, npti 332 zfnet(ji) = qsr_ice_1d(ji) - q sr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar332 zfnet(ji) = qsr_ice_1d(ji) - qtr_ice_top_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar 333 333 END DO 334 334 ! … … 539 539 DO jk = 1, nlay_i 540 540 DO ji = 1, npti 541 ztmelt _i = -tmut * sz_i_1d(ji,jk) + rt0542 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt _i), rt0 - 100._wp )541 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 542 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 543 543 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 544 544 END DO … … 704 704 DO jk = 1, nlay_i 705 705 DO ji = 1, npti 706 ztmelt _i = -tmut * sz_i_1d(ji,jk) + rt0707 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt _i), rt0 - 100._wp )706 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 707 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 708 708 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 709 709 END DO … … 728 728 !----------------------------- 729 729 ! 730 ! --- update conduction fluxes731 ! 730 ! --- calculate conduction fluxes (positive downward) 731 732 732 DO ji = 1, npti 733 733 ! ! surface ice conduction flux 734 fc_su(ji) = - isnow(ji) * zkappa_s(ji,0)* zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) &735 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0)* zg1 * ( t_i_1d(ji,1) - t_su_1d(ji) )734 qcn_ice_top_1d(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) & 735 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * ( t_i_1d(ji,1) - t_su_1d(ji) ) 736 736 ! ! bottom ice conduction flux 737 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * zg1 * ( t_bo_1d(ji) - t_i_1d(ji,nlay_i) )737 qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1 * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 738 738 END DO 739 739 … … 750 750 ! 751 751 DO ji = 1, npti 752 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji)- qcn_ice_1d(ji) ) * a_i_1d(ji)752 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qcn_ice_top_1d(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji) 753 753 END DO 754 754 ! … … 770 770 771 771 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 772 zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice )*a_i_1d(ji) 772 zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji) & 773 & + zdq * r1_rdtice ) * a_i_1d(ji) 773 774 ELSE ! case T_su = 0degC 774 zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice )*a_i_1d(ji) 775 zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji) & 776 & + zdq * r1_rdtice ) * a_i_1d(ji) 775 777 ENDIF 776 778 777 779 ELSEIF( k_jules == np_jules_ACTIVE ) THEN 778 780 779 zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 781 zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji) & 782 & + zdq * r1_rdtice ) * a_i_1d(ji) 780 783 781 784 ENDIF … … 787 790 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 788 791 ! 789 END DO790 !791 ! --- SIMIP diagnostics792 !793 DO ji = 1, npti794 !--- Conduction fluxes (positive downwards)795 diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)796 diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su (ji) * a_i_1d(ji) / at_i_1d(ji)797 798 !--- Snow-ice interfacial temperature (diagnostic SIMIP)799 zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)800 IF( h_s_1d(ji) >= zhs_min ) THEN801 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + &802 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac )803 ELSE804 t_si_1d(ji) = t_su_1d(ji)805 ENDIF806 792 END DO 807 793 ! … … 819 805 cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0) 820 806 ELSE 821 cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / 0.1_wp807 cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp 822 808 ENDIF 823 809 ENDIF … … 827 813 IF( k_jules == np_jules_EMULE ) THEN 828 814 ! Restore temperatures to their initial values 829 t_s_1d (1:npti,:) = ztsold (1:npti,:)830 t_i_1d (1:npti,:) = ztiold (1:npti,:)831 qcn_ice_1d(1:npti) = fc_su(1:npti)815 t_s_1d (1:npti,:) = ztsold (1:npti,:) 816 t_i_1d (1:npti,:) = ztiold (1:npti,:) 817 qcn_ice_1d(1:npti) = qcn_ice_top_1d(1:npti) 832 818 ENDIF 833 819 ! 820 ! --- SIMIP diagnostics 821 ! 822 DO ji = 1, npti 823 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 824 zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 825 IF( h_s_1d(ji) >= zhs_min ) THEN 826 t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + & 827 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac ) 828 ELSE 829 t_si_1d(ji) = t_su_1d(ji) 830 ENDIF 831 END DO 832 ! 834 833 END SUBROUTINE ice_thd_zdf_BL99 835 836 834 837 835 #else
Note: See TracChangeset
for help on using the changeset viewer.
