- Timestamp:
- 2014-11-27T17:13:38+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r4901 r4902 129 129 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 130 130 CALL wrk_alloc( jpij, zintermelt ) 131 CALL wrk_alloc( jpij, jkmax, zdeltah, zh_i )131 CALL wrk_alloc( jpij, nlay_i+1, zdeltah, zh_i ) 132 132 CALL wrk_alloc( jpij, icount ) 133 133 … … 155 155 DO jk = 1, nlay_i 156 156 DO ji = kideb, kiut 157 h_i_old (ji,jk) = ht_i_ b(ji) / REAL( nlay_i )158 qh_i_old(ji,jk) = q_i_ b(ji,jk) * h_i_old(ji,jk)157 h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 158 qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) 159 159 ENDDO 160 160 ENDDO … … 165 165 ! 166 166 DO ji = kideb, kiut 167 zinda = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )168 ztmelts = zinda * rtt + ( 1._wp - zinda ) * rtt167 zinda = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 168 ztmelts = zinda * rtt + ( 1._wp - zinda ) * rtt 169 169 170 170 zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 171 171 zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 172 172 173 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_ b(ji) - ztmelts ) )173 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts ) ) 174 174 zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 175 175 END DO … … 181 181 !------------------------------------------------------------------------------! 182 182 DO ji = kideb, kiut 183 IF( t_s_ b(ji,1) > rtt ) THEN !!! Internal melting183 IF( t_s_1d(ji,1) > rtt ) THEN !!! Internal melting 184 184 ! Contribution to heat flux to the ocean [W.m-2], < 0 185 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_ b(ji,1) * ht_s_b(ji) * a_i_b(ji) * r1_rdtice185 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 186 186 ! Contribution to mass flux 187 wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_ b(ji) * a_i_b(ji) * r1_rdtice187 wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 188 188 ! updates 189 ht_s_ b(ji) = 0._wp190 q_s_ b(ji,1) = 0._wp191 t_s_ b(ji,1) = rtt189 ht_s_1d(ji) = 0._wp 190 q_s_1d (ji,1) = 0._wp 191 t_s_1d (ji,1) = rtt 192 192 END IF 193 193 END DO … … 198 198 ! 199 199 DO ji = kideb, kiut 200 zh_s(ji) = ht_s_ b(ji) / REAL( nlay_s )200 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 201 201 END DO 202 202 ! 203 203 DO jk = 1, nlay_s 204 204 DO ji = kideb, kiut 205 zqh_s(ji) = zqh_s(ji) + q_s_ b(ji,jk) * zh_s(ji)205 zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji) 206 206 END DO 207 207 END DO … … 209 209 DO jk = 1, nlay_i 210 210 DO ji = kideb, kiut 211 zh_i(ji,jk) = ht_i_ b(ji) / REAL( nlay_i )212 zqh_i(ji) = zqh_i(ji) + q_i_ b(ji,jk) * zh_i(ji,jk)211 zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i ) 212 zqh_i(ji) = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) 213 213 END DO 214 214 END DO … … 237 237 !----------- 238 238 ! thickness change 239 zcoeff = ( 1._wp - ( 1._wp - at_i_ b(ji) )**betas ) / at_i_b(ji)239 zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**betas ) / at_i_1d(ji) 240 240 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn 241 241 ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) … … 243 243 IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 244 244 ! heat flux from snow precip (>0, W.m-2) 245 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_ b(ji) * zqprec(ji) * r1_rdtice245 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 246 246 ! mass flux, <0 247 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_ b(ji) * zdh_s_pre(ji) * r1_rdtice247 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice 248 248 ! update thickness 249 ht_s_ b (ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_pre(ji) )249 ht_s_1d (ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 250 250 251 251 !--------------------- … … 258 258 zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting 259 259 ! heat used to melt snow (W.m-2, >0) 260 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_ b(ji) * zqprec(ji) * r1_rdtice260 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 261 261 ! snow melting only = water into the ocean (then without snow precip), >0 262 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_ b(ji) * zdh_s_mel(ji) * r1_rdtice262 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_mel(ji) * r1_rdtice 263 263 264 264 ! updates available heat + thickness 265 265 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) ) 266 ht_s_ b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_mel(ji) )267 zh_s (ji) = ht_s_ b(ji) / REAL( nlay_s )266 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) ) 267 zh_s (ji) = ht_s_1d(ji) / REAL( nlay_s ) 268 268 269 269 ENDIF … … 275 275 DO ji = kideb, kiut 276 276 ! thickness change 277 zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_ b(ji) ) )278 zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_ b(ji,jk) + epsi20 ) )279 zdeltah (ji,jk) = - zindh * zindq * zq_su(ji) / MAX( q_s_ b(ji,jk), epsi20 )277 zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 278 zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) + epsi20 ) ) 279 zdeltah (ji,jk) = - zindh * zindq * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) 280 280 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting 281 281 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 282 282 ! heat used to melt snow(W.m-2, >0) 283 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_ b(ji) * q_s_b(ji,jk) * r1_rdtice283 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice 284 284 ! snow melting only = water into the ocean (then without snow precip) 285 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice285 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 286 286 287 287 ! updates available heat + thickness 288 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_ b(ji,jk) )289 ht_s_ b(ji) = MAX( 0._wp , ht_s_b(ji) + zdeltah(ji,jk) )288 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 289 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 290 290 291 291 END DO … … 304 304 ! forced mode: snow thickness change due to sublimation 305 305 DO ji = kideb, kiut 306 zdh_s_sub(ji) = MAX( - ht_s_ b(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice )306 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 307 307 ! Heat flux by sublimation [W.m-2], < 0 308 308 ! sublimate first snow that had fallen, then pre-existing snow 309 309 zcoeff = ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) * zqprec(ji) + & 310 & ( 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) ) &311 & * a_i_ b(ji) * r1_rdtice310 & ( 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) ) & 311 & * a_i_1d(ji) * r1_rdtice 312 312 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff 313 313 ! Mass flux by sublimation 314 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_ b(ji) * zdh_s_sub(ji) * r1_rdtice314 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 315 315 ! new snow thickness 316 ht_s_ b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_sub(ji) )316 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 317 317 END DO 318 318 ENDIF … … 321 321 DO ji = kideb, kiut 322 322 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 323 zh_s(ji) = ht_s_ b(ji) / REAL( nlay_s )323 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 324 324 END DO ! ji 325 325 … … 331 331 DO jk = 1, nlay_s 332 332 DO ji = kideb,kiut 333 zindh = MAX( 0._wp , SIGN( 1._wp, - ht_s_ b(ji) + epsi20 ) )334 q_s_ b(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_b(ji), epsi20 ) * &333 zindh = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) + epsi20 ) ) 334 q_s_1d(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_1d(ji), epsi20 ) * & 335 335 & ( ( MAX( 0._wp, dh_s_tot(ji) ) ) * zqprec(ji) + & 336 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_ b(ji) ) * rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) )337 zq_s(ji) = zq_s(ji) + q_s_ b(ji,jk)336 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) ) 337 zq_s(ji) = zq_s(ji) + q_s_1d(ji,jk) 338 338 END DO 339 339 END DO … … 345 345 DO jk = 1, nlay_i 346 346 DO ji = kideb, kiut 347 zEi = - q_i_ b(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0]348 349 ztmelts = - tmut * s_i_ b(ji,jk) + rtt ! Melting point of layer k [K]347 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0] 348 349 ztmelts = - tmut * s_i_1d(ji,jk) + rtt ! Melting point of layer k [K] 350 350 351 351 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] … … 367 367 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 368 368 369 ! Contribution to salt flux (clem: using sm_i_ b and not s_i_b(jk) is ok)370 sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_ b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice369 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 370 sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 371 371 372 372 ! Contribution to heat flux [W.m-2], < 0 373 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice373 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 374 374 375 375 ! Total heat flux used in this process [W.m-2], > 0 376 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_ b(ji) * zdE * r1_rdtice376 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 377 377 378 378 ! Contribution to mass flux 379 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice379 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 380 380 381 381 ! record which layers have disappeared (for bottom melting) … … 387 387 388 388 ! update heat content (J.m-2) and layer thickness 389 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_ b(ji,jk)389 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 390 390 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 391 391 END DO … … 393 393 ! update ice thickness 394 394 DO ji = kideb, kiut 395 ht_i_ b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_surf(ji) )395 ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) ) 396 396 END DO 397 397 … … 423 423 !clem debug. Just to be sure that enthalpy at nlay_i+1 is null 424 424 DO ji = kideb, kiut 425 q_i_ b(ji,nlay_i+1) = 0._wp425 q_i_1d(ji,nlay_i+1) = 0._wp 426 426 END DO 427 427 … … 445 445 446 446 s_i_new(ji) = zswitch_sal * zfracs * sss_m(ii,ij) & ! New ice salinity 447 + ( 1. - zswitch_sal ) * sm_i_ b(ji)447 + ( 1. - zswitch_sal ) * sm_i_1d(ji) 448 448 ! New ice growth 449 449 ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K) 450 450 451 zt_i_new = zswitch_sal * t_bo_ b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)451 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 452 452 453 453 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) … … 455 455 & + rcp * ( ztmelts-rtt ) 456 456 457 zEw = rcp * ( t_bo_ b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)457 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 458 458 459 459 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) … … 461 461 dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 462 462 463 q_i_ b(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0)463 q_i_1d(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0) 464 464 465 465 ENDIF ! fc_bo_i … … 476 476 ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K) 477 477 478 zt_i_new = zswitch_sal * t_bo_ b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)478 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 479 479 480 480 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) … … 482 482 & + rcp * ( ztmelts-rtt ) 483 483 484 zEw = rcp * ( t_bo_ b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)484 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 485 485 486 486 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 487 487 488 488 ! Contribution to heat flux to the ocean [W.m-2], >0 489 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice489 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 490 490 491 491 ! Total heat flux used in this process [W.m-2], <0 492 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_ b(ji) * zdE * r1_rdtice492 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 493 493 494 494 ! Contribution to salt flux, <0 495 sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_ b(ji) * zfmdt * r1_rdtice495 sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_1d(ji) * zfmdt * r1_rdtice 496 496 497 497 ! Contribution to mass flux, <0 498 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_ b(ji) * dh_i_bott(ji) * r1_rdtice498 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice 499 499 500 500 ! update heat content (J.m-2) and layer thickness 501 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)501 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) 502 502 h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji) 503 503 ENDIF … … 512 512 IF( zf_tt(ji) >= 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared from surface melting 513 513 514 ztmelts = - tmut * s_i_ b(ji,jk) + rtt ! Melting point of layer jk (K)515 516 IF( t_i_ b(ji,jk) >= ztmelts ) THEN !!! Internal melting514 ztmelts = - tmut * s_i_1d(ji,jk) + rtt ! Melting point of layer jk (K) 515 516 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 517 517 zintermelt(ji) = 1._wp 518 518 519 zEi = - q_i_ b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)520 521 !!zEw = rcp * ( t_i_ b(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_b(J/kg, <0)519 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 520 521 !!zEw = rcp * ( t_i_1d(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) 522 522 523 523 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) … … 532 532 533 533 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 534 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_ b(ji) * zEi * r1_rdtice535 536 ! Contribution to salt flux (clem: using sm_i_ b and not s_i_b(jk) is ok)537 sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_ b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice534 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 535 536 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 537 sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 538 538 539 539 ! Contribution to mass flux 540 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice540 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 541 541 542 542 ! update heat content (J.m-2) and layer thickness 543 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_ b(ji,jk)543 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 544 544 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 545 545 546 546 ELSE !!! Basal melting 547 547 548 zEi = - q_i_ b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)548 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 549 549 550 550 zEw = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0) … … 567 567 568 568 ! Contribution to heat flux to the ocean [W.m-2], <0 569 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice570 571 ! Contribution to salt flux (clem: using sm_i_ b and not s_i_b(jk) is ok)572 sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_ b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice569 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 570 571 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 572 sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 573 573 574 574 ! Total heat flux used in this process [W.m-2], >0 575 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_ b(ji) * zdE * r1_rdtice575 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 576 576 577 577 ! Contribution to mass flux 578 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_ b(ji) * zdeltah(ji,jk) * r1_rdtice578 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 579 579 580 580 ! update heat content (J.m-2) and layer thickness 581 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_ b(ji,jk)581 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 582 582 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 583 583 ENDIF … … 602 602 ! 603 603 ! ! excessive energy is sent to lateral ablation 604 ! zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_ b(ji) - epsi20 ) )605 ! zq_1cat(ji) = zinda * rhoic * lfus * at_i_ b(ji) / MAX( 1._wp - at_i_b(ji) , epsi20 ) * zdvres ! J.m-2 >=0604 ! zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_1d(ji) - epsi20 ) ) 605 ! zq_1cat(ji) = zinda * rhoic * lfus * at_i_1d(ji) / MAX( 1._wp - at_i_1d(ji) , epsi20 ) * zdvres ! J.m-2 >=0 606 606 ! 607 607 ! ! correct salt and mass fluxes 608 ! 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 approximation609 ! wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_ b(ji) * zdvres * r1_rdtice608 ! 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 609 ! wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdvres * r1_rdtice 610 610 ! ENDIF 611 611 ! END DO … … 616 616 !------------------------------------------- 617 617 DO ji = kideb, kiut 618 ht_i_ b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_bott(ji) )618 ht_i_1d(ji) = MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) ) 619 619 END DO 620 620 … … 627 627 DO ji = kideb, kiut 628 628 zq_rema(ji) = zq_su(ji) + zq_bo(ji) 629 ! zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_ b(ji) ) ) ! =1 if snow629 ! zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) ! =1 if snow 630 630 ! zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zq_s(ji) + epsi20 ) ) 631 631 ! zdeltah (ji,1) = - zindh * zindq * zq_rema(ji) / MAX( zq_s(ji), epsi20 ) 632 ! zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_ b(ji) ) ) ! bound melting632 ! zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 633 633 ! zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) 634 634 ! dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) 635 ! ht_s_ b (ji) = ht_s_b(ji) + zdeltah(ji,1)635 ! ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1) 636 636 ! 637 637 ! zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji) ! update available heat (J.m-2) 638 638 ! ! heat used to melt snow 639 ! hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_ b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0)639 ! hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0) 640 640 ! ! Contribution to mass flux 641 ! wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_ b(ji) * zdeltah(ji,1) * r1_rdtice641 ! wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 642 642 ! 643 643 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 644 644 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 645 hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_ b(ji) ) * r1_rdtice645 hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 646 646 647 647 IF( ln_nicep .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) … … 656 656 DO ji = kideb, kiut 657 657 ! 658 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_ b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic ) )659 660 ht_i_ b(ji) = ht_i_b(ji) + dh_snowice(ji)661 ht_s_ b(ji) = ht_s_b(ji) - dh_snowice(ji)658 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic ) ) 659 660 ht_i_1d(ji) = ht_i_1d(ji) + dh_snowice(ji) 661 ht_s_1d(ji) = ht_s_1d(ji) - dh_snowice(ji) 662 662 663 663 ! Salinity of snow ice 664 664 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 665 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_ b(ji)665 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 666 666 667 667 ! entrapment during snow ice formation 668 668 ! new salinity difference stored (to be used in limthd_ent.F90) 669 669 IF ( num_sal == 2 ) THEN 670 zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_ b(ji) - epsi10 ) )670 zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) ) 671 671 ! salinity dif due to snow-ice formation 672 dsm_i_si_1d(ji) = ( zs_snic - sm_i_ b(ji) ) * dh_snowice(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch672 dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi10 ) * zswitch 673 673 ! salinity dif due to bottom growth 674 674 IF ( zf_tt(ji) < 0._wp ) THEN 675 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_ b(ji) ) * dh_i_bott(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch675 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi10 ) * zswitch 676 676 ENDIF 677 677 ENDIF … … 685 685 686 686 ! Contribution to heat flux 687 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_ b(ji) * zEw * r1_rdtice687 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 688 688 689 689 ! Contribution to salt flux 690 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_ b(ji) * zfmdt * r1_rdtice690 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice 691 691 692 692 ! Contribution to mass flux 693 693 ! All snow is thrown in the ocean, and seawater is taken to replace the volume 694 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_ b(ji) * dh_snowice(ji) * rhoic * r1_rdtice695 wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_ b(ji) * dh_snowice(ji) * rhosn * r1_rdtice694 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice 695 wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice 696 696 697 697 ! update heat content (J.m-2) and layer thickness 698 qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_ b(ji,1) + zfmdt * zEw698 qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw 699 699 h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 700 700 701 701 ! Total ablation (to debug) 702 IF( ht_i_ b(ji) <= 0._wp ) a_i_b(ji) = 0._wp702 IF( ht_i_1d(ji) <= 0._wp ) a_i_1d(ji) = 0._wp 703 703 704 704 END DO !ji … … 710 710 !clem bug: we should take snow into account here 711 711 DO ji = kideb, kiut 712 zindh = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_ b(ji) ) )713 t_su_ b(ji) = zindh * t_su_b(ji) + ( 1.0 - zindh ) * rtt712 zindh = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 713 t_su_1d(ji) = zindh * t_su_1d(ji) + ( 1.0 - zindh ) * rtt 714 714 END DO ! ji 715 715 … … 717 717 DO ji = kideb,kiut 718 718 ! mask enthalpy 719 zinda = MAX( 0._wp , SIGN( 1._wp, - ht_s_ b(ji) ) )720 q_s_ b(ji,jk) = ( 1.0 - zinda ) * q_s_b(ji,jk)721 ! recalculate t_s_ b from q_s_b722 t_s_ b(ji,jk) = rtt + ( 1._wp - zinda ) * ( - q_s_b(ji,jk) / ( rhosn * cpic ) + lfus / cpic )719 zinda = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) ) 720 q_s_1d(ji,jk) = ( 1.0 - zinda ) * q_s_1d(ji,jk) 721 ! recalculate t_s_1d from q_s_1d 722 t_s_1d(ji,jk) = rtt + ( 1._wp - zinda ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 723 723 END DO 724 724 END DO … … 727 727 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 728 728 CALL wrk_dealloc( jpij, zintermelt ) 729 CALL wrk_dealloc( jpij, jkmax, zdeltah, zh_i )729 CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i ) 730 730 CALL wrk_dealloc( jpij, icount ) 731 731 !
Note: See TracChangeset
for help on using the changeset viewer.
