- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r4688 r6225 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 USE cpl_oasis3, ONLY : lk_cpl29 27 30 28 IMPLICIT NONE 31 29 PRIVATE 32 30 33 PUBLIC lim_thd_dh ! called by lim_thd 34 35 REAL(wp) :: epsi20 = 1.e-20 ! constant values 36 REAL(wp) :: epsi10 = 1.e-10 ! 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 37 37 38 38 !!---------------------------------------------------------------------- … … 74 74 75 75 REAL(wp) :: ztmelts ! local scalar 76 REAL(wp) :: z dh, zfdum !76 REAL(wp) :: zfdum 77 77 REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment 78 REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads 79 REAL(wp) :: zs_snic ! snow-ice salinity 78 REAL(wp) :: zs_snic ! snow-ice salinity 80 79 REAL(wp) :: zswi1 ! switch for computation of bottom salinity 81 80 REAL(wp) :: zswi12 ! switch for computation of bottom salinity … … 91 90 REAL(wp) :: zsstK ! SST in Kelvin 92 91 93 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness94 92 REAL(wp), POINTER, DIMENSION(:) :: zqprec ! energy of fallen snow (J.m-3) 95 93 REAL(wp), POINTER, DIMENSION(:) :: zq_su ! heat for surface ablation (J.m-2) 96 94 REAL(wp), POINTER, DIMENSION(:) :: zq_bo ! heat for bottom ablation (J.m-2) 97 REAL(wp), POINTER, DIMENSION(:) :: zq_1cat ! corrected heat in case 1-cat and hmelt>15cm (J.m-2)98 95 REAL(wp), POINTER, DIMENSION(:) :: zq_rema ! remaining heat at the end of the routine (J.m-2) 99 REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2) 100 INTEGER , POINTER, DIMENSION(:) :: icount ! number of layers vanished by melting 96 REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2) 101 97 102 98 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_mel ! snow melt … … 106 102 REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah 107 103 REAL(wp), POINTER, DIMENSION(:,:) :: zh_i ! ice layer thickness 104 INTEGER , POINTER, DIMENSION(:,:) :: icount ! number of layers vanished by melting 108 105 109 106 REAL(wp), POINTER, DIMENSION(:) :: zqh_i ! total ice heat content (J.m-2) 110 107 REAL(wp), POINTER, DIMENSION(:) :: zqh_s ! total snow heat content (J.m-2) 111 108 REAL(wp), POINTER, DIMENSION(:) :: zq_s ! total snow enthalpy (J.m-3) 112 113 ! mass and salt flux (clem) 114 REAL(wp) :: z dvres, zswitch_sal, zswitch109 REAL(wp), POINTER, DIMENSION(:) :: zsnw ! distribution of snow after wind blowing 110 111 REAL(wp) :: zswitch_sal 115 112 116 113 ! Heat conservation 117 114 INTEGER :: num_iter_max 118 REAL(wp) :: zinda, zindq, zindh119 REAL(wp), POINTER, DIMENSION(:) :: zintermelt ! debug120 115 121 116 !!------------------------------------------------------------------ 122 117 123 ! Discriminate between varying salinity (n um_sal=2) and prescribed cases (other values)124 SELECT CASE( n um_sal ) ! varying salinity or not118 ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values) 119 SELECT CASE( nn_icesal ) ! varying salinity or not 125 120 CASE( 1, 3, 4 ) ; zswitch_sal = 0 ! prescribed salinity profile 126 121 CASE( 2 ) ; zswitch_sal = 1 ! varying salinity profile 127 122 END SELECT 128 123 129 CALL wrk_alloc( jpij, z h_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema)124 CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw ) 130 125 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 131 CALL wrk_alloc( jpij, zintermelt ) 132 CALL wrk_alloc( jpij, jkmax, zdeltah, zh_i ) 133 CALL wrk_alloc( jpij, icount ) 134 126 CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 127 CALL wrk_alloc( jpij, nlay_i, icount ) 128 135 129 dh_i_surf (:) = 0._wp ; dh_i_bott (:) = 0._wp ; dh_snowice(:) = 0._wp 136 130 dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp 137 138 zqprec (:) = 0._wp ; zq_su (:) = 0._wp ; zq_bo (:) = 0._wp ; zf_tt (:) = 0._wp 139 zq_1cat(:) = 0._wp ; zq_rema(:) = 0._wp 140 141 zh_s (:) = 0._wp 142 zdh_s_pre(:) = 0._wp 143 zdh_s_mel(:) = 0._wp 144 zdh_s_sub(:) = 0._wp 145 zqh_s (:) = 0._wp 146 zqh_i (:) = 0._wp 147 148 zh_i (:,:) = 0._wp 149 zdeltah (:,:) = 0._wp 150 zintermelt(:) = 0._wp 151 icount (:) = 0 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 152 145 153 146 ! initialize layer thicknesses and enthalpies … … 156 149 DO jk = 1, nlay_i 157 150 DO ji = kideb, kiut 158 h_i_old (ji,jk) = ht_i_ b(ji) / REAL( nlay_i )159 qh_i_old(ji,jk) = q_i_ b(ji,jk) * h_i_old(ji,jk)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) 160 153 ENDDO 161 154 ENDDO … … 166 159 ! 167 160 DO ji = kideb, kiut 168 zinda = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) ) 169 ztmelts = zinda * rtt + ( 1._wp - zinda ) * rtt 170 171 zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 172 zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 173 174 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_b(ji) - ztmelts ) ) 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 ) ) 175 165 zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 176 166 END DO … … 182 172 !------------------------------------------------------------------------------! 183 173 DO ji = kideb, kiut 184 IF( t_s_ b(ji,1) > rtt) THEN !!! Internal melting174 IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting 185 175 ! Contribution to heat flux to the ocean [W.m-2], < 0 186 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_ b(ji,1) * ht_s_b(ji) * a_i_b(ji) * r1_rdtice176 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 187 177 ! Contribution to mass flux 188 wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_ b(ji) * a_i_b(ji) * r1_rdtice178 wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 189 179 ! updates 190 ht_s_ b(ji) = 0._wp191 q_s_ b(ji,1) = 0._wp192 t_s_ b (ji,1) = rtt180 ht_s_1d(ji) = 0._wp 181 q_s_1d (ji,1) = 0._wp 182 t_s_1d (ji,1) = rt0 193 183 END IF 194 184 END DO … … 198 188 !------------------------------------------------------------! 199 189 ! 200 DO ji = kideb, kiut201 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )202 END DO203 !204 190 DO jk = 1, nlay_s 205 191 DO ji = kideb, kiut 206 zqh_s(ji) = zqh_s(ji) + q_s_ b(ji,jk) * zh_s(ji)192 zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s 207 193 END DO 208 194 END DO … … 210 196 DO jk = 1, nlay_i 211 197 DO ji = kideb, kiut 212 zh_i(ji,jk) = ht_i_ b(ji) / REAL( nlay_i )213 zqh_i(ji) = zqh_i(ji) + q_i_ b(ji,jk) * zh_i(ji,jk)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) 214 200 END DO 215 201 END DO … … 233 219 ! Martin Vancoppenolle, December 2006 234 220 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 235 224 DO ji = kideb, kiut 236 225 !----------- … … 238 227 !----------- 239 228 ! thickness change 240 zcoeff = ( 1._wp - ( 1._wp - at_i_b(ji) )**betas ) / at_i_b(ji) 241 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn 242 ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) 243 zqprec (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus ) 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) 244 232 IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 245 233 ! heat flux from snow precip (>0, W.m-2) 246 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_ b(ji) * zqprec(ji) * r1_rdtice234 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 247 235 ! mass flux, <0 248 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_b(ji) * zdh_s_pre(ji) * r1_rdtice 249 ! update thickness 250 ht_s_b (ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_pre(ji) ) 236 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice 251 237 252 238 !--------------------- … … 254 240 !--------------------- 255 241 ! thickness change 256 IF( zdh_s_pre(ji) > 0._wp ) THEN 257 zindq = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) ) 258 zdh_s_mel (ji) = - zindq * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 259 zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting 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 260 245 ! heat used to melt snow (W.m-2, >0) 261 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zd h_s_mel(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice246 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 262 247 ! snow melting only = water into the ocean (then without snow precip), >0 263 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdh_s_mel(ji) * r1_rdtice 264 265 ! updates available heat + thickness 266 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) ) 267 ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_mel(ji) ) 268 zh_s (ji) = ht_s_b(ji) / REAL( nlay_s ) 269 270 ENDIF 271 END DO 272 273 ! If heat still available, then melt more snow 274 zdeltah(:,:) = 0._wp ! important 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 275 259 DO jk = 1, nlay_s 276 260 DO ji = kideb, kiut 277 261 ! thickness change 278 zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_b(ji) ) )279 zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_b(ji,jk) + epsi20) )280 zdeltah (ji,jk) = - zindh * zindq * zq_su(ji) / MAX( q_s_b(ji,jk), epsi20 )281 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting262 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 282 266 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 283 267 ! heat used to melt snow(W.m-2, >0) 284 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_ b(ji) * q_s_b(ji,jk) * r1_rdtice268 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice 285 269 ! snow melting only = water into the ocean (then without snow precip) 286 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice 287 270 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 288 271 ! updates available heat + thickness 289 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_b(ji,jk) ) 290 ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdeltah(ji,jk) ) 291 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) ) 292 274 END DO 293 275 END DO … … 297 279 !---------------------- 298 280 ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 299 ! clem comment: not counted in mass exchange in limsbc since this is an exchange with atm. (not ocean)281 ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean) 300 282 ! clem comment: ice should also sublimate 301 IF( lk_cpl ) THEN 302 ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) 303 zdh_s_sub(:) = 0._wp 304 ELSE 305 ! forced mode: snow thickness change due to sublimation 306 DO ji = kideb, kiut 307 zdh_s_sub(ji) = MAX( - ht_s_b(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 308 ! Heat flux by sublimation [W.m-2], < 0 309 ! sublimate first snow that had fallen, then pre-existing snow 310 zcoeff = ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) * zqprec(ji) + & 311 & ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_b(ji,1) ) & 312 & * a_i_b(ji) * r1_rdtice 313 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff 314 ! Mass flux by sublimation 315 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_b(ji) * zdh_s_sub(ji) * r1_rdtice 316 ! new snow thickness 317 ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_sub(ji) ) 318 END DO 319 ENDIF 320 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 321 302 ! --- Update snow diags --- ! 322 303 DO ji = kideb, kiut 323 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 324 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 325 END DO ! ji 304 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 305 END DO 326 306 327 307 !------------------------------------------- … … 332 312 DO jk = 1, nlay_s 333 313 DO ji = kideb,kiut 334 zindh = MAX( 0._wp , SIGN( 1._wp, - ht_s_b(ji) +epsi20 ) )335 q_s_ b(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_b(ji), epsi20 ) *&336 & ( ( MAX( 0._wp, dh_s_tot(ji) )) * zqprec(ji) + &337 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_b(ji) ) * rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) )338 zq_s(ji) = zq_s(ji) + q_s_ b(ji,jk)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) 339 319 END DO 340 320 END DO … … 345 325 zdeltah(:,:) = 0._wp ! important 346 326 DO jk = 1, nlay_i 347 DO ji = kideb, kiut 348 zEi = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0] 349 350 ztmelts = - tmut * s_i_b(ji,jk) + rtt ! Melting point of layer k [K] 351 352 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] 353 354 zdE = zEi - zEw ! Specific enthalpy difference < 0 355 356 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 357 358 zdeltah(ji,jk) = - zfmdt / rhoic ! Melt of layer jk [m, <0] 359 360 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] 361 362 zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 363 364 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt 365 366 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 367 368 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 369 370 ! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok) 371 sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 372 373 ! Contribution to heat flux [W.m-2], < 0 374 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice 375 376 ! Total heat flux used in this process [W.m-2], > 0 377 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice 378 379 ! Contribution to mass flux 380 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice 381 327 DO ji = kideb, kiut 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 382 385 ! record which layers have disappeared (for bottom melting) 383 386 ! => icount=0 : no layer has vanished 384 387 ! => icount=5 : 5 layers have vanished 385 zindh = NINT( MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk)) ) ) )386 icount(ji ) = icount(ji) + zindh387 zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )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) ) 388 391 389 392 ! update heat content (J.m-2) and layer thickness 390 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_ b(ji,jk)393 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 391 394 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 392 395 END DO … … 394 397 ! update ice thickness 395 398 DO ji = kideb, kiut 396 ht_i_ b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_surf(ji) )399 ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) ) 397 400 END DO 398 401 … … 416 419 ! -> need for an iterative procedure, which converges quickly 417 420 418 IF ( num_sal == 2 ) THEN 419 num_iter_max = 5 420 ELSE 421 num_iter_max = 1 422 ENDIF 423 424 !clem debug. Just to be sure that enthalpy at nlay_i+1 is null 425 DO ji = kideb, kiut 426 q_i_b(ji,nlay_i+1) = 0._wp 427 END DO 421 num_iter_max = 1 422 IF( nn_icesal == 2 ) num_iter_max = 5 428 423 429 424 ! Iterative procedure … … 446 441 447 442 s_i_new(ji) = zswitch_sal * zfracs * sss_m(ii,ij) & ! New ice salinity 448 + ( 1. - zswitch_sal ) * sm_i_ b(ji)443 + ( 1. - zswitch_sal ) * sm_i_1d(ji) 449 444 ! New ice growth 450 ztmelts = - tmut * s_i_new(ji) + rt t! New ice melting point (K)451 452 zt_i_new = zswitch_sal * t_bo_ b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)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) 453 448 454 449 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) 455 & - lfus * ( 1.0 - ( ztmelts - rt t ) / ( zt_i_new - rtt) ) &456 & + rcp * ( ztmelts-rt t)457 458 zEw = rcp * ( t_bo_ b(ji) - rt0 ) ! Specific enthalpy of seawater (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) 459 454 460 455 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) … … 462 457 dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 463 458 464 q_i_ b(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0)465 466 ENDIF ! fc_bo_i467 END DO ! ji468 END DO ! iter459 q_i_1d(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0) 460 461 ENDIF 462 END DO 463 END DO 469 464 470 465 ! Contribution to Energy and Salt Fluxes … … 475 470 zfmdt = - rhoic * dh_i_bott(ji) ! Mass flux x time step (kg/m2, < 0) 476 471 477 ztmelts = - tmut * s_i_new(ji) + rt t! New ice melting point (K)472 ztmelts = - tmut * s_i_new(ji) + rt0 ! New ice melting point (K) 478 473 479 zt_i_new = zswitch_sal * t_bo_ b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)474 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 480 475 481 476 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) 482 & - lfus * ( 1.0 - ( ztmelts - rt t ) / ( zt_i_new - rtt) ) &483 & + rcp * ( ztmelts-rt t)477 & - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) ) & 478 & + rcp * ( ztmelts-rt0 ) 484 479 485 zEw = rcp * ( t_bo_ b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)480 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 486 481 487 482 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 488 483 489 484 ! Contribution to heat flux to the ocean [W.m-2], >0 490 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice485 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 491 486 492 487 ! Total heat flux used in this process [W.m-2], <0 493 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_ b(ji) * zdE * r1_rdtice488 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 494 489 495 490 ! Contribution to salt flux, <0 496 sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_b(ji) * zfmdt* r1_rdtice491 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice 497 492 498 493 ! Contribution to mass flux, <0 499 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_ b(ji) * dh_i_bott(ji) * r1_rdtice494 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice 500 495 501 496 ! update heat content (J.m-2) and layer thickness 502 qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_ b(ji,nlay_i+1)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) 503 498 h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji) 504 499 ENDIF … … 511 506 DO jk = nlay_i, 1, -1 512 507 DO ji = kideb, kiut 513 IF( zf_tt(ji) >= 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared from surface melting 514 515 ztmelts = - tmut * s_i_b(ji,jk) + rtt ! Melting point of layer jk (K) 516 517 IF( t_i_b(ji,jk) >= ztmelts ) THEN !!! Internal melting 518 zintermelt(ji) = 1._wp 519 520 zEi = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 521 522 !!zEw = rcp * ( t_i_b(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_b (J/kg, <0) 523 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) 524 515 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 525 516 ! set up at 0 since no energy is needed to melt water...(it is already melted) 526 527 zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 528 ! this should normally not happen, but sometimes, heat diffusion leads to this 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 529 519 530 520 dh_i_bott (ji) = dh_i_bott(ji) + zdeltah(ji,jk) 531 521 532 zfmdt = - zdeltah(ji,jk) * rhoic 522 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 533 523 534 524 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 535 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_ b(ji) * zEi * r1_rdtice536 537 ! Contribution to salt flux (clem: using sm_i_ b and not s_i_b(jk) is ok)538 sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic* r1_rdtice525 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 539 529 540 530 ! Contribution to mass flux 541 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice531 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 542 532 543 533 ! update heat content (J.m-2) and layer thickness 544 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_ b(ji,jk)534 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 545 535 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 546 536 547 537 ELSE !!! Basal melting 548 538 549 zEi = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 550 551 zEw = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0) 552 553 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 554 555 zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0) 556 557 zdeltah(ji,jk) = - zfmdt / rhoic ! Gross thickness change 558 559 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 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 560 548 561 zq_bo(ji) 562 563 dh_i_bott(ji) 564 565 zfmdt 566 567 zQm 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 568 556 569 557 ! Contribution to heat flux to the ocean [W.m-2], <0 570 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice571 572 ! Contribution to salt flux (clem: using sm_i_ b and not s_i_b(jk) is ok)573 sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic* r1_rdtice558 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 574 562 575 563 ! Total heat flux used in this process [W.m-2], >0 576 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice564 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 577 565 578 566 ! Contribution to mass flux 579 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice567 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 580 568 581 569 ! update heat content (J.m-2) and layer thickness 582 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_ b(ji,jk)570 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 583 571 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 584 572 ENDIF 585 573 586 574 ENDIF 587 END DO ! ji 588 END DO ! jk 589 590 !------------------------------------------------------------------------------! 591 ! Excessive ablation in a 1-category model 592 ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15) 593 !------------------------------------------------------------------------------! 594 ! ??? keep ??? 595 ! clem bug: I think this should be included above, so we would not have to 596 ! track heat/salt/mass fluxes backwards 597 ! IF( jpl == 1 ) THEN 598 ! DO ji = kideb, kiut 599 ! IF( zf_tt(ji) >= 0._wp ) THEN 600 ! zdh = MAX( hmelt , dh_i_bott(ji) ) 601 ! zdvres = zdh - dh_i_bott(ji) ! >=0 602 ! dh_i_bott(ji) = zdh 603 ! 604 ! ! excessive energy is sent to lateral ablation 605 ! zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_b(ji) - epsi20 ) ) 606 ! zq_1cat(ji) = zinda * rhoic * lfus * at_i_b(ji) / MAX( 1._wp - at_i_b(ji) , epsi20 ) * zdvres ! J.m-2 >=0 607 ! 608 ! ! correct salt and mass fluxes 609 ! sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation 610 ! wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_b(ji) * zdvres * r1_rdtice 611 ! ENDIF 612 ! END DO 613 ! ENDIF 575 END DO 576 END DO 614 577 615 578 !------------------------------------------- … … 617 580 !------------------------------------------- 618 581 DO ji = kideb, kiut 619 ht_i_ b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_bott(ji) )582 ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) ) 620 583 END DO 621 584 … … 628 591 DO ji = kideb, kiut 629 592 zq_rema(ji) = zq_su(ji) + zq_bo(ji) 630 ! zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_b(ji) ) ) ! =1 if snow 631 ! zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zq_s(ji) + epsi20 ) ) 632 ! zdeltah (ji,1) = - zindh * zindq * zq_rema(ji) / MAX( zq_s(ji), epsi20 ) 633 ! zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_b(ji) ) ) ! bound melting 634 ! zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) 635 ! dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) 636 ! ht_s_b (ji) = ht_s_b(ji) + zdeltah(ji,1) 637 ! 638 ! zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji) ! update available heat (J.m-2) 639 ! ! heat used to melt snow 640 ! hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0) 641 ! ! Contribution to mass flux 642 ! wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdeltah(ji,1) * r1_rdtice 643 ! 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 ! 644 606 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 645 607 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 646 hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_ 1cat(ji) + zq_rema(ji) * a_i_b(ji) ) * r1_rdtice647 648 IF( ln_ nicep.AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)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) 649 611 END DO 650 612 … … 657 619 DO ji = kideb, kiut 658 620 ! 659 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_ b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic ) )660 661 ht_i_ b(ji) = ht_i_b(ji) + dh_snowice(ji)662 ht_s_ b(ji) = ht_s_b(ji) - dh_snowice(ji)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) 663 625 664 626 ! Salinity of snow ice 665 627 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 666 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_b(ji)628 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 667 629 668 630 ! entrapment during snow ice formation 669 ! new salinity difference stored (to be used in limthd_ ent.F90)670 IF ( n um_sal == 2 ) THEN671 zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_b(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 ) ) 672 634 ! salinity dif due to snow-ice formation 673 dsm_i_si_1d(ji) = ( zs_snic - sm_i_ b(ji) ) * dh_snowice(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch635 dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch 674 636 ! salinity dif due to bottom growth 675 637 IF ( zf_tt(ji) < 0._wp ) THEN 676 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_ b(ji) ) * dh_i_bott(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch638 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch 677 639 ENDIF 678 640 ENDIF … … 686 648 687 649 ! Contribution to heat flux 688 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice650 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 689 651 690 652 ! Contribution to salt flux 691 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_ b(ji) * zfmdt * r1_rdtice653 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice 692 654 693 655 ! Contribution to mass flux 694 656 ! All snow is thrown in the ocean, and seawater is taken to replace the volume 695 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_ b(ji) * dh_snowice(ji) * rhoic * r1_rdtice696 wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_ b(ji) * dh_snowice(ji) * rhosn * r1_rdtice657 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 697 659 698 660 ! update heat content (J.m-2) and layer thickness 699 qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_ b(ji,1) + zfmdt * zEw661 qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw 700 662 h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 701 663 702 ! Total ablation (to debug) 703 IF( ht_i_b(ji) <= 0._wp ) a_i_b(ji) = 0._wp 704 705 END DO !ji 664 END DO 706 665 707 666 ! … … 709 668 ! Update temperature, energy 710 669 !------------------------------------------- 711 !clem bug: we should take snow into account here 712 DO ji = kideb, kiut 713 zindh = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) ) 714 t_su_b(ji) = zindh * t_su_b(ji) + ( 1.0 - zindh ) * rtt 715 END DO ! ji 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 716 674 717 675 DO jk = 1, nlay_s 718 676 DO ji = kideb,kiut 719 677 ! mask enthalpy 720 zinda = MAX( 0._wp , SIGN( 1._wp, - ht_s_b(ji) ) )721 q_s_ b(ji,jk) = ( 1.0 - zinda ) * q_s_b(ji,jk)722 ! recalculate t_s_ b from q_s_b723 t_s_ b(ji,jk) = rtt + ( 1._wp - zinda ) * ( - q_s_b(ji,jk) / ( rhosn * cpic ) + lfus / cpic )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 ) 724 682 END DO 725 683 END DO 726 684 727 CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) 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 ) 728 689 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 729 CALL wrk_dealloc( jpij, zintermelt ) 730 CALL wrk_dealloc( jpij, jkmax, zdeltah, zh_i ) 731 CALL wrk_dealloc( jpij, icount ) 690 CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 691 CALL wrk_dealloc( jpij, nlay_i, icount ) 732 692 ! 733 693 ! 734 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 735 713 736 714 #else
Note: See TracChangeset
for help on using the changeset viewer.
