Changeset 5965 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
- Timestamp:
- 2015-12-01T16:35:30+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r4333 r5965 6 6 !! History : LIM ! 2003-05 (M. Vancoppenolle) Original code in 1D 7 7 !! ! 2005-06 (M. Vancoppenolle) 3D version 8 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm_snw & rdm_ice8 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw & wfx_ice 9 9 !! 3.4 ! 2011-02 (G. Madec) dynamical allocation 10 10 !! 3.5 ! 2012-10 (G. Madec & co) salt flux + bug fixes … … 20 20 USE sbc_oce ! Surface boundary condition: ocean fields 21 21 USE ice ! LIM variables 22 USE par_ice ! LIM parameters23 22 USE thd_ice ! LIM thermodynamics 24 23 USE in_out_manager ! I/O manager … … 26 25 USE wrk_nemo ! work arrays 27 26 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 28 27 29 28 IMPLICIT NONE 30 29 PRIVATE 31 30 32 PUBLIC lim_thd_dh ! called by lim_thd 33 34 REAL(wp) :: epsi20 = 1.e-20 ! constant values 35 REAL(wp) :: epsi10 = 1.e-10 ! 36 REAL(wp) :: epsi13 = 1.e-13 ! 37 REAL(wp) :: zzero = 0._wp ! 38 REAL(wp) :: zone = 1._wp ! 31 PUBLIC lim_thd_dh ! called by lim_thd 32 PUBLIC lim_thd_snwblow ! called in sbcblk/sbcclio/sbccpl and here 33 34 INTERFACE lim_thd_snwblow 35 MODULE PROCEDURE lim_thd_snwblow_1d, lim_thd_snwblow_2d 36 END INTERFACE 39 37 40 38 !!---------------------------------------------------------------------- … … 45 43 CONTAINS 46 44 47 SUBROUTINE lim_thd_dh( kideb, kiut , jl)45 SUBROUTINE lim_thd_dh( kideb, kiut ) 48 46 !!------------------------------------------------------------------ 49 47 !! *** ROUTINE lim_thd_dh *** … … 70 68 !!------------------------------------------------------------------ 71 69 INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied 72 INTEGER , INTENT(in) :: jl ! Thickness cateogry number73 70 !! 74 71 INTEGER :: ji , jk ! dummy loop indices 75 72 INTEGER :: ii, ij ! 2D corresponding indices to ji 76 INTEGER :: isnow ! switch for presence (1) or absence (0) of snow77 INTEGER :: isnowic ! snow ice formation not78 INTEGER :: i_ice_switch ! ice thickness above a certain treshold or not79 73 INTEGER :: iter 80 74 81 REAL(wp) :: zzfmass_i, zihgnew ! local scalar 82 REAL(wp) :: zzfmass_s, zhsnew, ztmelts ! local scalar 83 REAL(wp) :: zhn, zdhcf, zdhbf, zhni, zhnfi, zihg ! 84 REAL(wp) :: zdhnm, zhnnew, zhisn, zihic, zzc ! 75 REAL(wp) :: ztmelts ! local scalar 76 REAL(wp) :: zfdum 85 77 REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment 86 REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads 87 REAL(wp) :: zsm_snowice ! snow-ice salinity 78 REAL(wp) :: zs_snic ! snow-ice salinity 88 79 REAL(wp) :: zswi1 ! switch for computation of bottom salinity 89 80 REAL(wp) :: zswi12 ! switch for computation of bottom salinity 90 81 REAL(wp) :: zswi2 ! switch for computation of bottom salinity 91 82 REAL(wp) :: zgrr ! bottom growth rate 92 REAL(wp) :: ztform ! bottom formation temperature 93 ! 94 REAL(wp), POINTER, DIMENSION(:) :: zh_i ! ice layer thickness 95 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness 96 REAL(wp), POINTER, DIMENSION(:) :: ztfs ! melting point 97 REAL(wp), POINTER, DIMENSION(:) :: zhsold ! old snow thickness 98 REAL(wp), POINTER, DIMENSION(:) :: zqprec ! energy of fallen snow 99 REAL(wp), POINTER, DIMENSION(:) :: zqfont_su ! incoming, remaining surface energy 100 REAL(wp), POINTER, DIMENSION(:) :: zqfont_bo ! incoming, bottom energy 101 REAL(wp), POINTER, DIMENSION(:) :: z_f_surf ! surface heat for ablation 102 REAL(wp), POINTER, DIMENSION(:) :: zhgnew ! new ice thickness 103 REAL(wp), POINTER, DIMENSION(:) :: zfmass_i ! 83 REAL(wp) :: zt_i_new ! bottom formation temperature 84 85 REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean 86 REAL(wp) :: zEi ! specific enthalpy of sea ice (J/kg) 87 REAL(wp) :: zEw ! specific enthalpy of exchanged water (J/kg) 88 REAL(wp) :: zdE ! specific enthalpy difference (J/kg) 89 REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean 90 REAL(wp) :: zsstK ! SST in Kelvin 91 92 REAL(wp), POINTER, DIMENSION(:) :: zqprec ! energy of fallen snow (J.m-3) 93 REAL(wp), POINTER, DIMENSION(:) :: zq_su ! heat for surface ablation (J.m-2) 94 REAL(wp), POINTER, DIMENSION(:) :: zq_bo ! heat for bottom ablation (J.m-2) 95 REAL(wp), POINTER, DIMENSION(:) :: zq_rema ! remaining heat at the end of the routine (J.m-2) 96 REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2) 104 97 105 98 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_mel ! snow melt … … 108 101 109 102 REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah 110 111 ! Pathological cases 112 REAL(wp), POINTER, DIMENSION(:) :: zfdt_init ! total incoming heat for ice melt 113 REAL(wp), POINTER, DIMENSION(:) :: zfdt_final ! total remaing heat for ice melt 114 REAL(wp), POINTER, DIMENSION(:) :: zqt_i ! total ice heat content 115 REAL(wp), POINTER, DIMENSION(:) :: zqt_s ! total snow heat content 116 REAL(wp), POINTER, DIMENSION(:) :: zqt_dummy ! dummy heat content 117 118 REAL(wp), POINTER, DIMENSION(:,:) :: zqt_i_lay ! total ice heat content 119 120 ! mass and salt flux (clem) 121 REAL(wp) :: zdvres, zdvsur, zdvbot 122 REAL(wp), POINTER, DIMENSION(:) :: zviold, zvsold ! old ice volume... 103 REAL(wp), POINTER, DIMENSION(:,:) :: zh_i ! ice layer thickness 104 INTEGER , POINTER, DIMENSION(:,:) :: icount ! number of layers vanished by melting 105 106 REAL(wp), POINTER, DIMENSION(:) :: zqh_i ! total ice heat content (J.m-2) 107 REAL(wp), POINTER, DIMENSION(:) :: zqh_s ! total snow heat content (J.m-2) 108 REAL(wp), POINTER, DIMENSION(:) :: zq_s ! total snow enthalpy (J.m-3) 109 REAL(wp), POINTER, DIMENSION(:) :: zsnw ! distribution of snow after wind blowing 110 111 REAL(wp) :: zswitch_sal 123 112 124 113 ! Heat conservation 125 INTEGER :: num_iter_max, numce_dh 126 REAL(wp) :: meance_dh 127 REAL(wp) :: zinda 128 REAL(wp), POINTER, DIMENSION(:) :: zinnermelt 129 REAL(wp), POINTER, DIMENSION(:) :: zfbase, zdq_i 114 INTEGER :: num_iter_max 115 130 116 !!------------------------------------------------------------------ 131 117 132 CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 133 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 134 CALL wrk_alloc( jpij, zinnermelt, zfbase, zdq_i ) 135 CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay ) 136 137 CALL wrk_alloc( jpij, zviold, zvsold ) ! clem 138 139 ftotal_fin(:) = 0._wp 140 zfdt_init (:) = 0._wp 141 zfdt_final(:) = 0._wp 142 143 dh_i_surf (:) = 0._wp 144 dh_i_bott (:) = 0._wp 145 dh_snowice(:) = 0._wp 146 147 DO ji = kideb, kiut 148 old_ht_i_b(ji) = ht_i_b(ji) 149 old_ht_s_b(ji) = ht_s_b(ji) 150 zviold(ji) = a_i_b(ji) * ht_i_b(ji) ! clem 151 zvsold(ji) = a_i_b(ji) * ht_s_b(ji) ! clem 152 END DO 118 ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values) 119 SELECT CASE( nn_icesal ) ! varying salinity or not 120 CASE( 1, 3, 4 ) ; zswitch_sal = 0 ! prescribed salinity profile 121 CASE( 2 ) ; zswitch_sal = 1 ! varying salinity profile 122 END SELECT 123 124 CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw ) 125 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 126 CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 127 CALL wrk_alloc( jpij, nlay_i, icount ) 128 129 dh_i_surf (:) = 0._wp ; dh_i_bott (:) = 0._wp ; dh_snowice(:) = 0._wp 130 dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp 131 132 zqprec (:) = 0._wp ; zq_su (:) = 0._wp ; zq_bo (:) = 0._wp ; zf_tt(:) = 0._wp 133 zq_rema (:) = 0._wp ; zsnw (:) = 0._wp 134 zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zqh_i(:) = 0._wp 135 zqh_s (:) = 0._wp ; zq_s (:) = 0._wp 136 137 zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp 138 icount (:,:) = 0 139 140 141 ! Initialize enthalpy at nlay_i+1 142 DO ji = kideb, kiut 143 q_i_1d(ji,nlay_i+1) = 0._wp 144 END DO 145 146 ! initialize layer thicknesses and enthalpies 147 h_i_old (:,0:nlay_i+1) = 0._wp 148 qh_i_old(:,0:nlay_i+1) = 0._wp 149 DO jk = 1, nlay_i 150 DO ji = kideb, kiut 151 h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i 152 qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) 153 ENDDO 154 ENDDO 153 155 ! 154 156 !------------------------------------------------------------------------------! 155 ! 1) Calculate available heat for surface a blation!157 ! 1) Calculate available heat for surface and bottom ablation ! 156 158 !------------------------------------------------------------------------------! 157 159 ! 158 160 DO ji = kideb, kiut 159 isnow = INT( 1.0 - MAX( 0.0 , SIGN( 1.0 , - ht_s_b(ji) ) ) ) 160 ztfs (ji) = isnow * rtt + ( 1.0 - isnow ) * rtt 161 z_f_surf (ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 162 z_f_surf (ji) = MAX( zzero , z_f_surf(ji) ) * MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) ) 163 zfdt_init(ji) = ( z_f_surf(ji) + MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) * rdt_ice 164 END DO ! ji 165 166 zqfont_su (:) = 0._wp 167 zqfont_bo (:) = 0._wp 168 dsm_i_se_1d(:) = 0._wp 169 dsm_i_si_1d(:) = 0._wp 161 zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 162 zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 163 164 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 165 zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 166 END DO 167 170 168 ! 171 169 !------------------------------------------------------------------------------! 172 ! 2) Computing layer thicknesses and snow and sea-ice enthalpies. ! 170 ! If snow temperature is above freezing point, then snow melts 171 ! (should not happen but sometimes it does) 173 172 !------------------------------------------------------------------------------! 174 ! 175 DO ji = kideb, kiut ! Layer thickness 176 zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 177 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 178 END DO 179 ! 180 zqt_s(:) = 0._wp ! Total enthalpy of the snow 173 DO ji = kideb, kiut 174 IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting 175 ! Contribution to heat flux to the ocean [W.m-2], < 0 176 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 177 ! Contribution to mass flux 178 wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 179 ! updates 180 ht_s_1d(ji) = 0._wp 181 q_s_1d (ji,1) = 0._wp 182 t_s_1d (ji,1) = rt0 183 END IF 184 END DO 185 186 !------------------------------------------------------------! 187 ! 2) Computing layer thicknesses and enthalpies. ! 188 !------------------------------------------------------------! 189 ! 181 190 DO jk = 1, nlay_s 182 191 DO ji = kideb, kiut 183 zq t_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s )192 zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s 184 193 END DO 185 194 END DO 186 195 ! 187 zqt_i(:) = 0._wp ! Total enthalpy of the ice188 196 DO jk = 1, nlay_i 189 197 DO ji = kideb, kiut 190 zzc = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 191 zqt_i(ji) = zqt_i(ji) + zzc 192 zqt_i_lay(ji,jk) = zzc 198 zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i 199 zqh_i(ji) = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) 193 200 END DO 194 201 END DO … … 212 219 ! Martin Vancoppenolle, December 2006 213 220 214 ! Snow fall 215 DO ji = kideb, kiut 216 zcoeff = ( 1.0 - ( 1.0 - at_i_b(ji) )**betas ) / at_i_b(ji) 217 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn 218 END DO 219 zdh_s_mel(:) = 0._wp 220 221 ! Melt of fallen snow 222 DO ji = kideb, kiut 223 ! tatm_ice is now in K 224 zqprec (ji) = rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus ) 225 zqfont_su(ji) = z_f_surf(ji) * rdt_ice 226 zdeltah (ji,1) = MIN( 0.e0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) ) 227 zqfont_su(ji) = MAX( 0.e0 , - zdh_s_pre(ji) - zdeltah(ji,1) ) * zqprec(ji) 228 zdeltah (ji,1) = MAX( - zdh_s_pre(ji) , zdeltah(ji,1) ) 229 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) 230 ! heat conservation 231 qt_s_in(ji,jl) = qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji) 232 zqt_s (ji) = zqt_s (ji) + zqprec(ji) * zdh_s_pre(ji) 233 zqt_s (ji) = MAX( zqt_s(ji) - zqfont_su(ji) , 0.e0 ) 234 END DO 235 236 237 ! Snow melt due to surface heat imbalance 221 CALL lim_thd_snwblow( 1. - at_i_1d(kideb:kiut), zsnw(kideb:kiut) ) ! snow distribution over ice after wind blowing 222 223 zdeltah(:,:) = 0._wp 224 DO ji = kideb, kiut 225 !----------- 226 ! Snow fall 227 !----------- 228 ! thickness change 229 zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji) 230 ! enthalpy of the precip (>0, J.m-3) 231 zqprec (ji) = - qprec_ice_1d(ji) 232 IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 233 ! heat flux from snow precip (>0, W.m-2) 234 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 235 ! mass flux, <0 236 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice 237 238 !--------------------- 239 ! Melt of falling snow 240 !--------------------- 241 ! thickness change 242 rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 243 zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 244 zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting 245 ! heat used to melt snow (W.m-2, >0) 246 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 247 ! snow melting only = water into the ocean (then without snow precip), >0 248 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 249 ! updates available heat + precipitations after melting 250 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) ) 251 zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1) 252 253 ! update thickness 254 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 255 END DO 256 257 ! If heat still available (zq_su > 0), then melt more snow 258 zdeltah(:,:) = 0._wp 238 259 DO jk = 1, nlay_s 239 260 DO ji = kideb, kiut 240 zdeltah (ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk) 241 zqfont_su(ji) = MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * q_s_b(ji,jk) 242 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) 243 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) ! resulting melt of snow 261 ! thickness change 262 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 263 rswitch = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) ) 264 zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) 265 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting 266 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 267 ! heat used to melt snow(W.m-2, >0) 268 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice 269 ! snow melting only = water into the ocean (then without snow precip) 270 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 271 ! updates available heat + thickness 272 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 273 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 244 274 END DO 245 275 END DO 246 276 247 ! Apply snow melt to snow depth248 DO ji = kideb, kiut249 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji)250 ! Old and new snow depths251 zhsold(ji) = ht_s_b(ji)252 zhsnew = ht_s_b(ji) + dh_s_tot(ji)253 ! If snow is still present zhn = 1, else zhn = 0254 zhn = 1.0 - MAX( zzero , SIGN( zone , - zhsnew ) )255 ht_s_b(ji) = MAX( zzero , zhsnew )256 ! we recompute dh_s_tot (clem)257 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji)258 ! Volume and mass variations of snow259 dvsbq_1d (ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) )260 dvsbq_1d (ji) = MIN( zzero, dvsbq_1d(ji) )261 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji)262 END DO ! ji263 264 !--------------------------265 ! 3.2 Surface ice ablation266 !--------------------------267 DO ji = kideb, kiut268 z_f_surf (ji) = zqfont_su(ji) * r1_rdtice ! heat conservation test269 zdq_i (ji) = 0._wp270 END DO ! ji271 272 DO jk = 1, nlay_i273 DO ji = kideb, kiut274 ! ! melt of layer jk275 zdeltah (ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk)276 ! ! recompute heat available277 zqfont_su(ji ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk)278 ! ! melt of layer jk cannot be higher than its thickness279 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) )280 ! ! update surface melt281 dh_i_surf(ji ) = dh_i_surf(ji) + zdeltah(ji,jk)282 ! ! for energy conservation283 zdq_i (ji ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice284 !285 ! clem286 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) &287 & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice288 END DO289 END DO290 291 ! !-------------------292 IF( con_i .AND. jiindex_1d > 0 ) THEN ! Conservation test293 ! !-------------------294 numce_dh = 0295 meance_dh = 0._wp296 DO ji = kideb, kiut297 IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN298 numce_dh = numce_dh + 1299 meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji)300 ENDIF301 IF( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN!302 WRITE(numout,*) ' ALERTE heat loss for surface melt '303 WRITE(numout,*) ' ii, ij, jl :', ii, ij, jl304 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji)305 WRITE(numout,*) ' z_f_surf : ', z_f_surf(ji)306 WRITE(numout,*) ' zdq_i : ', zdq_i(ji)307 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji)308 WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)309 WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji)310 WRITE(numout,*) ' qlbbq_1d : ', qlbbq_1d(ji)311 WRITE(numout,*) ' s_i_new : ', s_i_new(ji)312 WRITE(numout,*) ' sss_m : ', sss_m(ii,ij)313 ENDIF314 END DO315 !316 IF( numce_dh > 0 ) meance_dh = meance_dh / numce_dh317 WRITE(numout,*) ' Error report - Category : ', jl318 WRITE(numout,*) ' ~~~~~~~~~~~~ '319 WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh320 WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh321 !322 ENDIF323 324 277 !---------------------- 325 ! 3. 3 Snow sublimation278 ! 3.2 Snow sublimation 326 279 !---------------------- 327 328 DO ji = kideb, kiut 329 ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 330 #if defined key_coupled 331 zdh_s_sub(ji) = 0._wp ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) 332 #else 333 ! ! forced mode: snow thickness change due to sublimation 334 zdh_s_sub(ji) = - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice 335 #endif 336 dh_s_tot (ji) = dh_s_tot(ji) + zdh_s_sub(ji) 337 zdhcf = ht_s_b(ji) + zdh_s_sub(ji) 338 ht_s_b (ji) = MAX( zzero , zdhcf ) 339 ! we recompute dh_s_tot 340 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) 341 qt_s_in (ji,jl) = qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1) 342 END DO 343 344 zqt_dummy(:) = 0.e0 280 ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 281 ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean) 282 ! clem comment: ice should also sublimate 283 zdeltah(:,:) = 0._wp 284 ! coupled mode: sublimation is set to 0 (evap_ice = 0) until further notice 285 ! forced mode: snow thickness change due to sublimation 286 DO ji = kideb, kiut 287 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 288 ! Heat flux by sublimation [W.m-2], < 0 289 ! sublimate first snow that had fallen, then pre-existing snow 290 zdeltah(ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 291 hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1) & 292 & ) * a_i_1d(ji) * r1_rdtice 293 ! Mass flux by sublimation 294 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 295 ! new snow thickness 296 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 297 ! update precipitations after sublimation and correct sublimation 298 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 299 zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 300 END DO 301 302 ! --- Update snow diags --- ! 303 DO ji = kideb, kiut 304 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 305 END DO 306 307 !------------------------------------------- 308 ! 3.3 Update temperature, energy 309 !------------------------------------------- 310 ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow) 311 zq_s(:) = 0._wp 345 312 DO jk = 1, nlay_s 346 313 DO ji = kideb,kiut 347 q_s_b (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 348 zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) ! heat conservation 314 rswitch = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) ) 315 q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) * & 316 & ( ( zdh_s_pre(ji) ) * zqprec(ji) + & 317 & ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 318 zq_s(ji) = zq_s(ji) + q_s_1d(ji,jk) 349 319 END DO 350 320 END DO 351 321 352 DO jk = 1, nlay_s 322 !-------------------------- 323 ! 3.4 Surface ice ablation 324 !-------------------------- 325 zdeltah(:,:) = 0._wp ! important 326 DO jk = 1, nlay_i 353 327 DO ji = kideb, kiut 354 ! In case of disparition of the snow, we have to update the snow temperatures 355 zhisn = MAX( zzero , SIGN( zone, - ht_s_b(ji) ) ) 356 t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt 357 q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk) 328 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer k [K] 329 330 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 331 332 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 333 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 334 ! set up at 0 since no energy is needed to melt water...(it is already melted) 335 zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 336 ! this should normally not happen, but sometimes, heat diffusion leads to this 337 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 338 339 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt 340 341 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 342 343 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 344 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 345 346 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 347 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 348 349 ! Contribution to mass flux 350 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 351 352 ELSE !!! Surface melting 353 354 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 355 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] 356 zdE = zEi - zEw ! Specific enthalpy difference < 0 357 358 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 359 360 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Melt of layer jk [m, <0] 361 362 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0] 363 364 zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 365 366 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt 367 368 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 369 370 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 371 372 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 373 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 374 375 ! Contribution to heat flux [W.m-2], < 0 376 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 377 378 ! Total heat flux used in this process [W.m-2], > 0 379 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 380 381 ! Contribution to mass flux 382 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 383 384 END IF 385 ! record which layers have disappeared (for bottom melting) 386 ! => icount=0 : no layer has vanished 387 ! => icount=5 : 5 layers have vanished 388 rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 389 icount(ji,jk) = NINT( rswitch ) 390 zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 391 392 ! update heat content (J.m-2) and layer thickness 393 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 394 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 358 395 END DO 396 END DO 397 ! update ice thickness 398 DO ji = kideb, kiut 399 ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) ) 359 400 END DO 360 401 … … 364 405 !------------------------------------------------------------------------------! 365 406 ! 366 ! Ice basal growth / melt is given by the ratio of heat budget over basal 367 ! ice heat content. Basal heat budget is given by the difference between 368 ! the inner conductive flux (fc_bo_i), from the open water heat flux 369 ! (qlbbqb) and the turbulent ocean flux (fbif). 370 ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice 371 372 !----------------------------------------------------- 373 ! 4.1 Basal growth - (a) salinity not varying in time 374 !----------------------------------------------------- 375 IF( num_sal /= 2 ) THEN ! ice salinity constant in time 407 !------------------ 408 ! 4.1 Basal growth 409 !------------------ 410 ! Basal growth is driven by heat imbalance at the ice-ocean interface, 411 ! between the inner conductive flux (fc_bo_i), from the open water heat flux 412 ! (fhld) and the turbulent ocean flux (fhtur). 413 ! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice 414 415 ! If salinity varies in time, an iterative procedure is required, because 416 ! the involved quantities are inter-dependent. 417 ! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi), 418 ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott) 419 ! -> need for an iterative procedure, which converges quickly 420 421 num_iter_max = 1 422 IF( nn_icesal == 2 ) num_iter_max = 5 423 424 ! Iterative procedure 425 DO iter = 1, num_iter_max 376 426 DO ji = kideb, kiut 377 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0._wp ) THEN 378 s_i_new(ji) = sm_i_b(ji) 379 ! Melting point in K 380 ztmelts = - tmut * s_i_new(ji) + rtt 381 ! New ice heat content (Bitz and Lipscomb, 1999) 382 ztform = t_i_b(ji,nlay_i) ! t_bo_b crashes in the 383 ! Baltic 384 q_i_b(ji,nlay_i+1) = rhoic * ( cpic * ( ztmelts - ztform ) & 385 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( ztform - rtt ) ) & 386 & - rcp * ( ztmelts - rtt ) ) 387 ! Basal growth rate = - F*dt / q 388 dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 389 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 427 IF( zf_tt(ji) < 0._wp ) THEN 428 429 ! New bottom ice salinity (Cox & Weeks, JGR88 ) 430 !--- zswi1 if dh/dt < 2.0e-8 431 !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7 432 !--- zswi2 if dh/dt > 3.6e-7 433 zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) ) 434 zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) ) 435 zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 436 zswi1 = 1. - zswi2 * zswi12 437 zfracs = MIN ( zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) & 438 & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 ) 439 440 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 441 442 s_i_new(ji) = zswitch_sal * zfracs * sss_m(ii,ij) & ! New ice salinity 443 + ( 1. - zswitch_sal ) * sm_i_1d(ji) 444 ! New ice growth 445 ztmelts = - tmut * s_i_new(ji) + rt0 ! New ice melting point (K) 446 447 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 448 449 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) 450 & - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) ) & 451 & + rcp * ( ztmelts-rt0 ) 452 453 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 454 455 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 456 457 dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 458 459 q_i_1d(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0) 460 390 461 ENDIF 391 462 END DO 392 ENDIF 393 394 !------------------------------------------------- 395 ! 4.1 Basal growth - (b) salinity varying in time 396 !------------------------------------------------- 397 IF( num_sal == 2 ) THEN 398 ! the growth rate (dh_i_bott) is function of the new ice heat content (q_i_b(nlay_i+1)). 399 ! q_i_b depends on the new ice salinity (snewice). 400 ! snewice depends on dh_i_bott ; it converges quickly, so, no problem 401 ! See Vancoppenolle et al., OM08 for more info on this 402 403 ! Initial value (tested 1D, can be anything between 1 and 20) 404 num_iter_max = 4 405 s_i_new(:) = 4.0 406 407 ! Iterative procedure 408 DO iter = 1, num_iter_max 409 DO ji = kideb, kiut 410 IF( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0 ) THEN 411 ii = MOD( npb(ji) - 1, jpi ) + 1 412 ij = ( npb(ji) - 1 ) / jpi + 1 413 ! Melting point in K 414 ztmelts = - tmut * s_i_new(ji) + rtt 415 ! New ice heat content (Bitz and Lipscomb, 1999) 416 q_i_b(ji,nlay_i+1) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) & 417 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) ) & 418 & - rcp * ( ztmelts-rtt ) ) 419 ! Bottom growth rate = - F*dt / q 420 dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 421 ! New ice salinity ( Cox and Weeks, JGR, 1988 ) 422 ! zswi2 (1) if dh_i_bott/rdt .GT. 3.6e-7 423 ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8 424 ! zswi1 (1) if dh_i_bott/rdt .LT. 2.0e-8 425 zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi13 ) ) 426 zswi2 = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 427 zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 428 zswi1 = 1. - zswi2 * zswi12 429 zfracs = zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) & 430 & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 431 zfracs = MIN( 0.5 , zfracs ) 432 s_i_new(ji) = zfracs * sss_m(ii,ij) 433 ENDIF ! fc_bo_i 434 END DO ! ji 435 END DO ! iter 436 437 ! Final values 438 DO ji = kideb, kiut 439 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN 440 ! New ice salinity must not exceed 20 psu 441 s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) 442 ! Metling point in K 443 ztmelts = - tmut * s_i_new(ji) + rtt 444 ! New ice heat content (Bitz and Lipscomb, 1999) 445 q_i_b(ji,nlay_i+1) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) & 446 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) ) & 447 & - rcp * ( ztmelts - rtt ) ) 448 ! Basal growth rate = - F*dt / q 449 dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 450 ! Salinity update 451 ! entrapment during bottom growth 452 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 453 ENDIF ! heat budget 454 END DO 455 ENDIF 463 END DO 464 465 ! Contribution to Energy and Salt Fluxes 466 DO ji = kideb, kiut 467 IF( zf_tt(ji) < 0._wp ) THEN 468 ! New ice growth 469 470 zfmdt = - rhoic * dh_i_bott(ji) ! Mass flux x time step (kg/m2, < 0) 471 472 ztmelts = - tmut * s_i_new(ji) + rt0 ! New ice melting point (K) 473 474 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 475 476 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) 477 & - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) ) & 478 & + rcp * ( ztmelts-rt0 ) 479 480 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 481 482 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 483 484 ! Contribution to heat flux to the ocean [W.m-2], >0 485 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 486 487 ! Total heat flux used in this process [W.m-2], <0 488 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 489 490 ! Contribution to salt flux, <0 491 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice 492 493 ! Contribution to mass flux, <0 494 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice 495 496 ! update heat content (J.m-2) and layer thickness 497 qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_1d(ji,nlay_i+1) 498 h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji) 499 ENDIF 500 END DO 456 501 457 502 !---------------- 458 503 ! 4.2 Basal melt 459 504 !---------------- 460 meance_dh = 0._wp 461 numce_dh = 0 462 zinnermelt(:) = 0._wp 463 464 DO ji = kideb, kiut 465 ! heat convergence at the surface > 0 466 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp ) THEN 467 s_i_new(ji) = s_i_b(ji,nlay_i) 468 zqfont_bo(ji) = rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) 469 zfbase(ji) = zqfont_bo(ji) * r1_rdtice ! heat conservation test 470 zdq_i(ji) = 0._wp 471 dh_i_bott(ji) = 0._wp 472 ENDIF 473 END DO 474 505 zdeltah(:,:) = 0._wp ! important 475 506 DO jk = nlay_i, 1, -1 476 507 DO ji = kideb, kiut 477 IF( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) >= 0._wp ) THEN 478 ztmelts = - tmut * s_i_b(ji,jk) + rtt 479 IF( t_i_b(ji,jk) >= ztmelts ) THEN !!gm : a comment is needed 480 zdeltah (ji,jk) = - zh_i(ji) 481 dh_i_bott (ji ) = dh_i_bott(ji) + zdeltah(ji,jk) 482 zinnermelt(ji ) = 1._wp 483 ELSE ! normal ablation 484 zdeltah (ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk) 485 zqfont_bo(ji ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 486 zdeltah (ji,jk) = MAX(zdeltah(ji,jk), - zh_i(ji) ) 487 dh_i_bott(ji ) = dh_i_bott(ji) + zdeltah(ji,jk) 488 zdq_i (ji ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 508 IF( zf_tt(ji) > 0._wp .AND. jk > icount(ji,jk) ) THEN ! do not calculate where layer has already disappeared by surface melting 509 510 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer jk (K) 511 512 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 513 514 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 515 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 516 ! set up at 0 since no energy is needed to melt water...(it is already melted) 517 zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 518 ! this should normally not happen, but sometimes, heat diffusion leads to this 519 520 dh_i_bott (ji) = dh_i_bott(ji) + zdeltah(ji,jk) 521 522 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 523 524 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 525 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 526 527 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 528 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 529 530 ! Contribution to mass flux 531 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 532 533 ! update heat content (J.m-2) and layer thickness 534 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 535 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 536 537 ELSE !!! Basal melting 538 539 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 540 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of meltwater (J/kg, <0) 541 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 542 543 zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0) 544 545 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Gross thickness change 546 547 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 548 549 zq_bo(ji) = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors 550 551 dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) ! Update basal melt 552 553 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 554 555 zQm = zfmdt * zEw ! Heat exchanged with ocean 556 557 ! Contribution to heat flux to the ocean [W.m-2], <0 558 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 559 560 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 561 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 562 563 ! Total heat flux used in this process [W.m-2], >0 564 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 565 566 ! Contribution to mass flux 567 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 568 569 ! update heat content (J.m-2) and layer thickness 570 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 571 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 489 572 ENDIF 490 ! clem: contribution to salt flux 491 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) & 492 & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 493 ENDIF 494 END DO ! ji 495 END DO ! jk 496 497 ! !------------------- 498 IF( con_i .AND. jiindex_1d > 0 ) THEN ! Conservation test 499 ! !------------------- 500 DO ji = kideb, kiut 501 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0 ) THEN 502 IF( ( zfbase(ji) + zdq_i(ji) ) >= 1.e-3 ) THEN 503 numce_dh = numce_dh + 1 504 meance_dh = meance_dh + zfbase(ji) + zdq_i(ji) 505 ENDIF 506 IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN 507 WRITE(numout,*) ' ALERTE heat loss for basal melt : ii, ij, jl :', ii, ij, jl 508 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 509 WRITE(numout,*) ' zfbase : ', zfbase(ji) 510 WRITE(numout,*) ' zdq_i : ', zdq_i(ji) 511 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 512 WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji) 513 WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji) 514 WRITE(numout,*) ' qlbbq_1d : ', qlbbq_1d(ji) 515 WRITE(numout,*) ' s_i_new : ', s_i_new(ji) 516 WRITE(numout,*) ' sss_m : ', sss_m(ii,ij) 517 WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 518 WRITE(numout,*) ' innermelt : ', INT( zinnermelt(ji) ) 519 ENDIF 573 520 574 ENDIF 521 575 END DO 522 IF( numce_dh > 0 ) meance_dh = meance_dh / numce_dh 523 WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh 524 WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh 525 WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d) 526 ! 527 ENDIF 528 529 ! 530 !------------------------------------------------------------------------------! 531 ! 5) Pathological cases ! 532 !------------------------------------------------------------------------------! 533 ! 534 !---------------------------------------------- 535 ! 5.1 Excessive ablation in a 1-category model 536 !---------------------------------------------- 537 538 DO ji = kideb, kiut 539 ! ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15) 540 IF( jpl == 1 ) THEN ; zdhbf = MAX( hmelt , dh_i_bott(ji) ) 541 ELSE ; zdhbf = dh_i_bott(ji) 542 ENDIF 543 zdvres = zdhbf - dh_i_bott(ji) 544 dh_i_bott(ji) = zdhbf 545 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice 546 ! ! excessive energy is sent to lateral ablation 547 zinda = MAX( 0._wp, SIGN( 1._wp , 1.0 - at_i_b(ji) - epsi10 ) ) 548 fsup(ji) = zinda * rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi10 ) * zdvres * r1_rdtice 549 END DO 550 551 !----------------------------------- 552 ! 5.2 More than available ice melts 553 !----------------------------------- 554 ! then heat applied minus heat content at previous time step should equal heat remaining 555 ! 556 DO ji = kideb, kiut 557 ! Adapt the remaining energy if too much ice melts 558 !-------------------------------------------------- 559 zdvres = MAX( 0._wp, - ht_i_b(ji) - dh_i_surf(ji) - dh_i_bott(ji) ) 560 zdvsur = MIN( 0._wp, dh_i_surf(ji) + zdvres ) - dh_i_surf(ji) ! fill the surface first 561 zdvbot = MAX( 0._wp, zdvres - zdvsur ) ! then the bottom 562 dh_i_surf (ji) = dh_i_surf(ji) + zdvsur ! clem 563 dh_i_bott (ji) = dh_i_bott(ji) + zdvbot ! clem 564 565 ! new ice thickness (clem) 566 zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 567 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 568 zhgnew(ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 569 570 ! !since ice volume is only used for outputs, we keep it global for all categories 571 dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 572 573 ! remaining heat 574 zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) + zqfont_bo(ji) ) 575 576 ! If snow remains, energy is used to melt snow 577 zhni = ht_s_b(ji) ! snow depth at previous time step 578 zihg = MAX( zzero , SIGN ( zone , - ht_s_b(ji) ) ) ! =0 if snow 579 580 ! energy of melting of remaining snow 581 zinda = MAX( 0._wp, SIGN( 1._wp , zhni - epsi10 ) ) 582 zqt_s(ji) = ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi10 ) * zinda 583 zdhnm = - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 ) 584 zhnfi = zhni + zdhnm 585 zfdt_final(ji) = MAX( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 ) 586 ht_s_b(ji) = MAX( zzero , zhnfi ) 587 zqt_s(ji) = zqt_s(ji) * ht_s_b(ji) 588 ! we recompute dh_s_tot (clem) 589 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) 590 591 ! Mass variations of ice and snow 592 !--------------------------------- 593 ! ! mass variation of the jl category 594 zzfmass_s = - a_i_b(ji) * ( zhni - ht_s_b(ji) ) * rhosn ! snow 595 zzfmass_i = a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic ! ice 596 ! 597 zfmass_i(ji) = zzfmass_i ! ice variation saved to compute salt flux (see below) 598 ! 599 ! ! mass variation cumulated over category 600 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s ! snow 601 !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i ! ice 602 603 ! Remaining heat to the ocean 604 !--------------------------------- 605 focea(ji) = - zfdt_final(ji) * r1_rdtice ! focea is in W.m-2 * dt 606 607 ! residual salt flux (clem) 608 !-------------------------- 609 ! surface 610 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvsur * rhoic * r1_rdtice 611 ! bottom 612 IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) >= 0._wp ) THEN ! melting 613 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 614 ELSE ! growth 615 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 616 ENDIF 617 ! 618 ! diagnostic 619 ii = MOD( npb(ji) - 1, jpi ) + 1 620 ij = ( npb(ji) - 1 ) / jpi + 1 621 diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 622 diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 623 diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 624 END DO 625 626 ftotal_fin (:) = zfdt_final(:) * r1_rdtice 627 628 !--------------------------- 629 ! heat fluxes 630 !--------------------------- 631 DO ji = kideb, kiut 632 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) ! =1 if ice 633 ! 634 ! Heat flux 635 ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1 636 ! excessive total ablation energy (focea) sent to the ocean 637 qfvbq_1d(ji) = qfvbq_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice 638 639 zihic = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) ) ! equals 0 if ht_i = 0, 1 if ht_i gt 0 640 fscbq_1d(ji) = a_i_b(ji) * fstbif_1d(ji) 641 qldif_1d(ji) = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea (ji) * a_i_b(ji) * rdt_ice & 642 & + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice 643 END DO ! ji 644 645 !------------------------------------------- 646 ! Correct temperature, energy and thickness 647 !------------------------------------------- 648 DO ji = kideb, kiut 649 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 650 t_su_b(ji) = zihgnew * t_su_b(ji) + ( 1.0 - zihgnew ) * rtt 651 END DO ! ji 652 653 DO jk = 1, nlay_i 654 DO ji = kideb, kiut 655 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 656 t_i_b(ji,jk) = zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt 657 q_i_b(ji,jk) = zihgnew * q_i_b(ji,jk) 658 END DO 659 END DO ! ji 660 661 DO ji = kideb, kiut 662 ht_i_b(ji) = zhgnew(ji) 663 END DO ! ji 576 END DO 577 578 !------------------------------------------- 579 ! Update temperature, energy 580 !------------------------------------------- 581 DO ji = kideb, kiut 582 ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) ) 583 END DO 584 585 !------------------------------------------- 586 ! 5. What to do with remaining energy 587 !------------------------------------------- 588 ! If heat still available for melting and snow remains, then melt more snow 589 !------------------------------------------- 590 zdeltah(:,:) = 0._wp ! important 591 DO ji = kideb, kiut 592 zq_rema(ji) = zq_su(ji) + zq_bo(ji) 593 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) ! =1 if snow 594 rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) ) 595 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 ) 596 zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 597 dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) 598 ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1) 599 600 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1) ! update available heat (J.m-2) 601 ! heat used to melt snow 602 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * q_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0) 603 ! Contribution to mass flux 604 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 605 ! 606 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 607 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 608 hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 609 610 IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 611 END DO 612 664 613 ! 665 614 !------------------------------------------------------------------------------| … … 670 619 DO ji = kideb, kiut 671 620 ! 672 dh_snowice(ji) = MAX( zzero , ( rhosn * ht_s_b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic ) ) 673 zhgnew(ji) = MAX( zhgnew(ji) , zhgnew(ji) + dh_snowice(ji) ) 674 zhnnew = MIN( ht_s_b(ji) , ht_s_b(ji) - dh_snowice(ji) ) 675 676 ! Changes in ice volume and ice mass. 677 dvnbq_1d (ji) = a_i_b(ji) * ( zhgnew(ji)-ht_i_b(ji) ) 678 dmgwi_1d (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn 679 680 !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic 681 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew - ht_s_b(ji) ) * rhosn 682 683 ! Equivalent salt flux (1) Snow-ice formation component 684 ! ----------------------------------------------------- 685 ii = MOD( npb(ji) - 1, jpi ) + 1 686 ij = ( npb(ji) - 1 ) / jpi + 1 687 688 IF( num_sal == 2 ) THEN ; zsm_snowice = sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic 689 ELSE ; zsm_snowice = sm_i_b(ji) 690 ENDIF 621 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic ) ) 622 623 ht_i_1d(ji) = ht_i_1d(ji) + dh_snowice(ji) 624 ht_s_1d(ji) = ht_s_1d(ji) - dh_snowice(ji) 625 626 ! Salinity of snow ice 627 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 628 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 629 691 630 ! entrapment during snow ice formation 692 ! clem: new salinity difference stored (to be used in limthd_ent.F90)693 IF ( n um_sal == 2 ) THEN694 i_ice_switch = MAX( 0._wp , SIGN( 1._wp , zhgnew(ji) - epsi10 ) )631 ! new salinity difference stored (to be used in limthd_sal.F90) 632 IF ( nn_icesal == 2 ) THEN 633 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 695 634 ! salinity dif due to snow-ice formation 696 dsm_i_si_1d(ji) = ( zs m_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch635 dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch 697 636 ! salinity dif due to bottom growth 698 IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0._wp ) THEN699 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_ b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch637 IF ( zf_tt(ji) < 0._wp ) THEN 638 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch 700 639 ENDIF 701 640 ENDIF 702 641 703 ! Actualize new snow and ice thickness. 704 ht_s_b(ji) = zhnnew 705 ht_i_b(ji) = zhgnew(ji) 706 707 ! Total ablation ! new lines added to debug 708 IF( ht_i_b(ji) <= 0._wp ) a_i_b(ji) = 0._wp 709 710 ! diagnostic ( snow ice growth ) 711 ii = MOD( npb(ji) - 1, jpi ) + 1 712 ij = ( npb(ji) - 1 ) / jpi + 1 713 diag_sni_gr(ii,ij) = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice 714 ! 715 ! salt flux 716 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 717 !-------------------------------- 718 ! Update mass fluxes (clem) 719 !-------------------------------- 720 rdm_ice_1d(ji) = rdm_ice_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic 721 rdm_snw_1d(ji) = rdm_snw_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn 722 723 END DO !ji 724 ! 725 CALL wrk_dealloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 726 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 727 CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i ) 728 CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay ) 729 ! 730 CALL wrk_dealloc( jpij, zviold, zvsold ) ! clem 642 ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0) 643 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 644 zfmdt = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0._wp ) ! <0 645 zsstK = sst_m(ii,ij) + rt0 646 zEw = rcp * ( zsstK - rt0 ) 647 zQm = zfmdt * zEw 648 649 ! Contribution to heat flux 650 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 651 652 ! Contribution to salt flux 653 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice 654 655 ! Contribution to mass flux 656 ! All snow is thrown in the ocean, and seawater is taken to replace the volume 657 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice 658 wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice 659 660 ! update heat content (J.m-2) and layer thickness 661 qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw 662 h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 663 664 END DO 665 666 ! 667 !------------------------------------------- 668 ! Update temperature, energy 669 !------------------------------------------- 670 DO ji = kideb, kiut 671 rswitch = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 672 t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0 673 END DO 674 675 DO jk = 1, nlay_s 676 DO ji = kideb,kiut 677 ! mask enthalpy 678 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) ) 679 q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk) 680 ! recalculate t_s_1d from q_s_1d 681 t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 682 END DO 683 END DO 684 685 ! --- ensure that a_i = 0 where ht_i = 0 --- 686 WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 687 688 CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw ) 689 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 690 CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 691 CALL wrk_dealloc( jpij, nlay_i, icount ) 692 ! 731 693 ! 732 694 END SUBROUTINE lim_thd_dh 695 696 697 !!-------------------------------------------------------------------------- 698 !! INTERFACE lim_thd_snwblow 699 !! ** Purpose : Compute distribution of precip over the ice 700 !!-------------------------------------------------------------------------- 701 SUBROUTINE lim_thd_snwblow_2d( pin, pout ) 702 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pin ! previous fraction lead ( pfrld or (1. - a_i_b) ) 703 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout 704 pout = ( 1._wp - ( pin )**rn_betas ) 705 END SUBROUTINE lim_thd_snwblow_2d 706 707 SUBROUTINE lim_thd_snwblow_1d( pin, pout ) 708 REAL(wp), DIMENSION(:), INTENT(in ) :: pin 709 REAL(wp), DIMENSION(:), INTENT(inout) :: pout 710 pout = ( 1._wp - ( pin )**rn_betas ) 711 END SUBROUTINE lim_thd_snwblow_1d 712 733 713 734 714 #else
Note: See TracChangeset
for help on using the changeset viewer.
