Changeset 5208 for branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
- Timestamp:
- 2015-04-13T15:08:59+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r4688 r5208 26 26 USE wrk_nemo ! work arrays 27 27 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 28 USE cpl_oasis3, ONLY : lk_cpl29 28 30 29 IMPLICIT NONE … … 32 31 33 32 PUBLIC lim_thd_dh ! called by lim_thd 34 35 REAL(wp) :: epsi20 = 1.e-20 ! constant values36 REAL(wp) :: epsi10 = 1.e-10 !37 33 38 34 !!---------------------------------------------------------------------- … … 112 108 113 109 ! mass and salt flux (clem) 114 REAL(wp) :: zdvres, zswitch_sal , zswitch110 REAL(wp) :: zdvres, zswitch_sal 115 111 116 112 ! Heat conservation 117 113 INTEGER :: num_iter_max 118 REAL(wp) :: zinda, zindq, zindh119 REAL(wp), POINTER, DIMENSION(:) :: zintermelt ! debug120 114 121 115 !!------------------------------------------------------------------ … … 129 123 CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) 130 124 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 ) 125 CALL wrk_alloc( jpij, nlay_i+1, zdeltah, zh_i ) 133 126 CALL wrk_alloc( jpij, icount ) 134 127 … … 148 141 zh_i (:,:) = 0._wp 149 142 zdeltah (:,:) = 0._wp 150 zintermelt(:) = 0._wp151 143 icount (:) = 0 152 144 … … 156 148 DO jk = 1, nlay_i 157 149 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)150 h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 151 qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) 160 152 ENDDO 161 153 ENDDO … … 166 158 ! 167 159 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) * rtt170 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 ) )160 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 161 ztmelts = rswitch * rtt + ( 1._wp - rswitch ) * rtt 162 163 zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 164 zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 165 166 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts ) ) 175 167 zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 176 168 END DO … … 182 174 !------------------------------------------------------------------------------! 183 175 DO ji = kideb, kiut 184 IF( t_s_ b(ji,1) > rtt ) THEN !!! Internal melting176 IF( t_s_1d(ji,1) > rtt ) THEN !!! Internal melting 185 177 ! 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_rdtice178 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 187 179 ! Contribution to mass flux 188 wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_ b(ji) * a_i_b(ji) * r1_rdtice180 wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 189 181 ! updates 190 ht_s_ b(ji) = 0._wp191 q_s_ b(ji,1) = 0._wp192 t_s_ b(ji,1) = rtt182 ht_s_1d(ji) = 0._wp 183 q_s_1d (ji,1) = 0._wp 184 t_s_1d (ji,1) = rtt 193 185 END IF 194 186 END DO … … 199 191 ! 200 192 DO ji = kideb, kiut 201 zh_s(ji) = ht_s_ b(ji) / REAL( nlay_s )193 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 202 194 END DO 203 195 ! 204 196 DO jk = 1, nlay_s 205 197 DO ji = kideb, kiut 206 zqh_s(ji) = zqh_s(ji) + q_s_ b(ji,jk) * zh_s(ji)198 zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji) 207 199 END DO 208 200 END DO … … 210 202 DO jk = 1, nlay_i 211 203 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)204 zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 205 zqh_i(ji) = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) 214 206 END DO 215 207 END DO … … 238 230 !----------- 239 231 ! thickness change 240 zcoeff = ( 1._wp - ( 1._wp - at_i_ b(ji) )**betas ) / at_i_b(ji)232 zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**betas ) / at_i_1d(ji) 241 233 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn 242 234 ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) … … 244 236 IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 245 237 ! 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_rdtice238 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 247 239 ! mass flux, <0 248 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_ b(ji) * zdh_s_pre(ji) * r1_rdtice240 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice 249 241 ! update thickness 250 ht_s_ b (ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_pre(ji) )242 ht_s_1d (ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 251 243 252 244 !--------------------- … … 255 247 ! thickness change 256 248 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 )249 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) ) 250 zdh_s_mel (ji) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 259 251 zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting 260 252 ! heat used to melt snow (W.m-2, >0) 261 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_ b(ji) * zqprec(ji) * r1_rdtice253 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 262 254 ! 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_rdtice255 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_mel(ji) * r1_rdtice 264 256 265 257 ! updates available heat + thickness 266 258 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 )259 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) ) 260 zh_s (ji) = ht_s_1d(ji) / REAL( nlay_s ) 269 261 270 262 ENDIF … … 276 268 DO ji = kideb, kiut 277 269 ! 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 )270 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 271 rswitch = rswitch * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) + epsi20 ) ) ) 272 zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) 281 273 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting 282 274 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 283 275 ! 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_rdtice276 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice 285 277 ! 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_rdtice278 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 287 279 288 280 ! 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) )281 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 282 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 291 283 292 284 END DO … … 305 297 ! forced mode: snow thickness change due to sublimation 306 298 DO ji = kideb, kiut 307 zdh_s_sub(ji) = MAX( - ht_s_ b(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice )299 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 308 300 ! Heat flux by sublimation [W.m-2], < 0 309 301 ! sublimate first snow that had fallen, then pre-existing snow 310 302 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_rdtice303 & ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_1d(ji,1) ) & 304 & * a_i_1d(ji) * r1_rdtice 313 305 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff 314 306 ! Mass flux by sublimation 315 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_ b(ji) * zdh_s_sub(ji) * r1_rdtice307 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 316 308 ! new snow thickness 317 ht_s_ b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_sub(ji) )309 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 318 310 END DO 319 311 ENDIF … … 322 314 DO ji = kideb, kiut 323 315 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 )316 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 325 317 END DO ! ji 326 318 … … 332 324 DO jk = 1, nlay_s 333 325 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 ) * &326 rswitch = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) + epsi20 ) ) 327 q_s_1d(ji,jk) = ( 1._wp - rswitch ) / MAX( ht_s_1d(ji), epsi20 ) * & 336 328 & ( ( 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)329 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) ) 330 zq_s(ji) = zq_s(ji) + q_s_1d(ji,jk) 339 331 END DO 340 332 END DO … … 346 338 DO jk = 1, nlay_i 347 339 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]340 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0] 341 342 ztmelts = - tmut * s_i_1d(ji,jk) + rtt ! Melting point of layer k [K] 351 343 352 344 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] … … 368 360 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 369 361 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_rdtice362 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 363 sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 372 364 373 365 ! Contribution to heat flux [W.m-2], < 0 374 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice366 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 375 367 376 368 ! 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_rdtice369 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 378 370 379 371 ! Contribution to mass flux 380 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice372 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 381 373 382 374 ! record which layers have disappeared (for bottom melting) 383 375 ! => icount=0 : no layer has vanished 384 376 ! => 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) + zindh377 rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 378 icount(ji) = icount(ji) + NINT( rswitch ) 387 379 zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 388 380 389 381 ! 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)382 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 391 383 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 392 384 END DO … … 394 386 ! update ice thickness 395 387 DO ji = kideb, kiut 396 ht_i_ b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_surf(ji) )388 ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) ) 397 389 END DO 398 390 … … 424 416 !clem debug. Just to be sure that enthalpy at nlay_i+1 is null 425 417 DO ji = kideb, kiut 426 q_i_ b(ji,nlay_i+1) = 0._wp418 q_i_1d(ji,nlay_i+1) = 0._wp 427 419 END DO 428 420 … … 446 438 447 439 s_i_new(ji) = zswitch_sal * zfracs * sss_m(ii,ij) & ! New ice salinity 448 + ( 1. - zswitch_sal ) * sm_i_ b(ji)440 + ( 1. - zswitch_sal ) * sm_i_1d(ji) 449 441 ! New ice growth 450 442 ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K) 451 443 452 zt_i_new = zswitch_sal * t_bo_ b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)444 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 453 445 454 446 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) … … 456 448 & + rcp * ( ztmelts-rtt ) 457 449 458 zEw = rcp * ( t_bo_ b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)450 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 459 451 460 452 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) … … 462 454 dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 463 455 464 q_i_ b(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0)456 q_i_1d(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0) 465 457 466 458 ENDIF ! fc_bo_i … … 477 469 ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K) 478 470 479 zt_i_new = zswitch_sal * t_bo_ b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)471 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 480 472 481 473 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) … … 483 475 & + rcp * ( ztmelts-rtt ) 484 476 485 zEw = rcp * ( t_bo_ b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)477 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 486 478 487 479 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 488 480 489 481 ! 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_rdtice482 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 491 483 492 484 ! 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_rdtice485 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 494 486 495 487 ! Contribution to salt flux, <0 496 sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_ b(ji) * zfmdt * r1_rdtice488 sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_1d(ji) * zfmdt * r1_rdtice 497 489 498 490 ! Contribution to mass flux, <0 499 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_ b(ji) * dh_i_bott(ji) * r1_rdtice491 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice 500 492 501 493 ! 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)494 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 495 h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji) 504 496 ENDIF … … 513 505 IF( zf_tt(ji) >= 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared from surface melting 514 506 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) 507 ztmelts = - tmut * s_i_1d(ji,jk) + rtt ! Melting point of layer jk (K) 508 509 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 510 511 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 512 513 !!zEw = rcp * ( t_i_1d(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) 523 514 524 515 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) … … 533 524 534 525 ! 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_rdtice526 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 527 528 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 529 sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 539 530 540 531 ! Contribution to mass flux 541 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice532 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 542 533 543 534 ! 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)535 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 545 536 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 546 537 547 538 ELSE !!! Basal melting 548 539 549 zEi = - q_i_ b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)540 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 550 541 551 542 zEw = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0) … … 568 559 569 560 ! 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_rdtice561 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 562 563 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 564 sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 574 565 575 566 ! 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_rdtice567 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 577 568 578 569 ! Contribution to mass flux 579 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice570 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 580 571 581 572 ! 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)573 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 583 574 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 584 575 ENDIF … … 603 594 ! 604 595 ! ! 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 >=0596 ! rswitch = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_1d(ji) - epsi20 ) ) 597 ! zq_1cat(ji) = rswitch * rhoic * lfus * at_i_1d(ji) / MAX( 1._wp - at_i_1d(ji) , epsi20 ) * zdvres ! J.m-2 >=0 607 598 ! 608 599 ! ! 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 approximation610 ! wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_ b(ji) * zdvres * r1_rdtice600 ! sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation 601 ! wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdvres * r1_rdtice 611 602 ! ENDIF 612 603 ! END DO … … 617 608 !------------------------------------------- 618 609 DO ji = kideb, kiut 619 ht_i_ b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_bott(ji) )610 ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) ) 620 611 END DO 621 612 … … 628 619 DO ji = kideb, kiut 629 620 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 snow621 ! zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) ! =1 if snow 631 622 ! zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zq_s(ji) + epsi20 ) ) 632 623 ! 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 melting624 ! zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 634 625 ! zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) 635 626 ! dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) 636 ! ht_s_ b (ji) = ht_s_b(ji) + zdeltah(ji,1)627 ! ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1) 637 628 ! 638 629 ! zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji) ! update available heat (J.m-2) 639 630 ! ! 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)631 ! hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0) 641 632 ! ! Contribution to mass flux 642 ! wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_ b(ji) * zdeltah(ji,1) * r1_rdtice633 ! wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 643 634 ! 644 635 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 645 636 ! 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_rdtice637 hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 647 638 648 639 IF( ln_nicep .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) … … 657 648 DO ji = kideb, kiut 658 649 ! 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)650 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic ) ) 651 652 ht_i_1d(ji) = ht_i_1d(ji) + dh_snowice(ji) 653 ht_s_1d(ji) = ht_s_1d(ji) - dh_snowice(ji) 663 654 664 655 ! Salinity of snow ice 665 656 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)657 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 667 658 668 659 ! entrapment during snow ice formation 669 660 ! new salinity difference stored (to be used in limthd_ent.F90) 670 661 IF ( num_sal == 2 ) THEN 671 zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_b(ji) - epsi10 ) )662 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) ) 672 663 ! 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 ) * zswitch664 dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi10 ) * rswitch 674 665 ! salinity dif due to bottom growth 675 666 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 ) * zswitch667 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi10 ) * rswitch 677 668 ENDIF 678 669 ENDIF … … 686 677 687 678 ! Contribution to heat flux 688 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice679 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 689 680 690 681 ! Contribution to salt flux 691 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_ b(ji) * zfmdt * r1_rdtice682 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice 692 683 693 684 ! Contribution to mass flux 694 685 ! 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_rdtice686 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice 687 wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice 697 688 698 689 ! 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 * zEw690 qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw 700 691 h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 701 692 702 693 ! Total ablation (to debug) 703 IF( ht_i_ b(ji) <= 0._wp ) a_i_b(ji) = 0._wp694 IF( ht_i_1d(ji) <= 0._wp ) a_i_1d(ji) = 0._wp 704 695 705 696 END DO !ji … … 711 702 !clem bug: we should take snow into account here 712 703 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 ) * rtt704 rswitch = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 705 t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rtt 715 706 END DO ! ji 716 707 … … 718 709 DO ji = kideb,kiut 719 710 ! 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 )711 rswitch = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) ) 712 q_s_1d(ji,jk) = ( 1.0 - rswitch ) * q_s_1d(ji,jk) 713 ! recalculate t_s_1d from q_s_1d 714 t_s_1d(ji,jk) = rtt + ( 1._wp - rswitch ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 724 715 END DO 725 716 END DO … … 727 718 CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) 728 719 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 ) 720 CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i ) 731 721 CALL wrk_dealloc( jpij, icount ) 732 722 !
Note: See TracChangeset
for help on using the changeset viewer.
