- Timestamp:
- 2018-07-13T09:28:50+02:00 (6 years ago)
- Location:
- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3
- Files:
-
- 1 edited
- 1 copied
-
. (copied) (copied from NEMO/trunk)
-
src/ICE/icethd_dh.F90 (modified) (22 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icethd_dh.F90
r9767 r9939 76 76 REAL(wp) :: zgrr ! bottom growth rate 77 77 REAL(wp) :: zt_i_new ! bottom formation temperature 78 REAL(wp) :: z1_rho ! 1/(rhos n+rau0-rhoic)78 REAL(wp) :: z1_rho ! 1/(rhos+rho0-rhoi) 79 79 80 80 REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean … … 181 181 DO ji = 1, npti 182 182 IF( t_s_1d(ji,jk) > rt0 ) THEN 183 hfx_res_1d (ji) = hfx_res_1d (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_ rdtice ! heat flux to the ocean [W.m-2], < 0184 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos n * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice ! mass flux183 hfx_res_1d (ji) = hfx_res_1d (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0 184 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! mass flux 185 185 ! updates 186 186 dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk) … … 202 202 ! 203 203 ! --- precipitation --- 204 zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhos n / at_i_1d(ji)! thickness change204 zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhos / at_i_1d(ji) ! thickness change 205 205 zqprec (ji) = - qprec_ice_1d(ji) ! enthalpy of the precip (>0, J.m-3) 206 206 ! 207 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_ rdtice ! heat flux from snow precip (>0, W.m-2)208 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos n * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice ! mass flux, <0207 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_Dt_ice ! heat flux from snow precip (>0, W.m-2) 208 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos * a_i_1d(ji) * zdh_s_pre(ji) * r1_Dt_ice ! mass flux, <0 209 209 210 210 ! --- melt of falling snow --- … … 212 212 zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) ! thickness change 213 213 zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting 214 hfx_snw_1d (ji) = hfx_snw_1d (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_ rdtice ! heat used to melt snow (W.m-2, >0)215 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos n * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice ! snow melting only = water into the ocean (then without snow precip), >0214 hfx_snw_1d (ji) = hfx_snw_1d (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_Dt_ice ! heat used to melt snow (W.m-2, >0) 215 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice ! snow melting only = water into the ocean (then without snow precip), >0 216 216 217 217 ! updates available heat + precipitations after melting … … 252 252 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 253 253 254 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_ rdtice ! heat used to melt snow(W.m-2, >0)255 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos n * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! snow melting only = water into the ocean (then without snow precip)254 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_Dt_ice ! heat used to melt snow(W.m-2, >0) 255 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! snow melting only = water into the ocean (then without snow precip) 256 256 257 257 ! updates available heat + thickness … … 273 273 IF( evap_ice_1d(ji) > 0._wp ) THEN 274 274 ! 275 zdh_s_sub (ji) = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos n* rdt_ice )276 zevap_rema(ji) = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhos n! remaining evap in kg.m-2 (used for ice melting later on)275 zdh_s_sub (ji) = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos * rdt_ice ) 276 zevap_rema(ji) = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhos ! remaining evap in kg.m-2 (used for ice melting later on) 277 277 zdeltah (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 278 278 279 279 hfx_sub_1d (ji) = hfx_sub_1d(ji) + & ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow) 280 280 & ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) ) & 281 & * a_i_1d(ji) * r1_ rdtice282 wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos n * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice ! Mass flux by sublimation281 & * a_i_1d(ji) * r1_Dt_ice 282 wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * a_i_1d(ji) * zdh_s_sub(ji) * r1_Dt_ice ! Mass flux by sublimation 283 283 284 284 ! new snow thickness … … 309 309 e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) * & 310 310 & ( ( zdh_s_pre(ji) ) * zqprec(ji) + & 311 & ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos n * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )311 & ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) ) 312 312 END DO 313 313 END DO … … 326 326 IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting 327 327 328 zEi = - e_i_1d(ji,jk) * r1_rhoi c! Specific enthalpy of layer k [J/kg, <0]328 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of layer k [J/kg, <0] 329 329 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 330 330 ! set up at 0 since no energy is needed to melt water...(it is already melted) 331 331 zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 332 332 ! this should normally not happen, but sometimes, heat diffusion leads to this 333 zfmdt = - zdeltah(ji,jk) * rhoi c! Mass flux x time step > 0333 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 334 334 335 335 dh_i_itm(ji) = dh_i_itm(ji) + zdeltah(ji,jk) ! Cumulate internal melting 336 336 337 zfmdt = - rhoi c * zdeltah(ji,jk)! Recompute mass flux [kg/m2, >0]338 339 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_ rdtice ! Heat flux to the ocean [W.m-2], <0337 zfmdt = - rhoi * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 338 339 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 340 340 ! ice enthalpy zEi is "sent" to the ocean 341 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi c * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux341 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux 342 342 ! using s_i_1d and not sz_i_1d(jk) is ok 343 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi c * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux343 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 344 344 345 345 ELSE !-- Surface melting 346 346 347 zEi = - e_i_1d(ji,jk) * r1_rhoi c! Specific enthalpy of layer k [J/kg, <0]347 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of layer k [J/kg, <0] 348 348 zEw = rcp * ztmelts ! Specific enthalpy of resulting meltwater [J/kg, <0] 349 349 zdE = zEi - zEw ! Specific enthalpy difference < 0 … … 351 351 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 352 352 353 zdeltah(ji,jk) = - zfmdt * r1_rhoi c! Melt of layer jk [m, <0]353 zdeltah(ji,jk) = - zfmdt * r1_rhoi ! Melt of layer jk [m, <0] 354 354 355 355 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] 356 356 357 zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoi c* zdE ) ! update available heat357 zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat 358 358 359 359 dh_i_sum(ji) = dh_i_sum(ji) + zdeltah(ji,jk) ! Cumulate surface melt 360 360 361 zfmdt = - rhoi c * zdeltah(ji,jk)! Recompute mass flux [kg/m2, >0]361 zfmdt = - rhoi * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 362 362 363 363 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 364 364 365 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi c * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux >0365 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0 366 366 ! using s_i_1d and not sz_i_1d(jk) is ok) 367 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_ rdtice ! Heat flux [W.m-2], < 0368 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_ rdtice ! Heat flux used in this process [W.m-2], > 0367 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux [W.m-2], < 0 368 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice ! Heat flux used in this process [W.m-2], > 0 369 369 ! 370 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi c * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux370 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 371 371 372 372 END IF … … 374 374 ! Ice sublimation 375 375 ! --------------- 376 zdum = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi c)376 zdum = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi ) 377 377 zdeltah (ji,jk) = zdeltah (ji,jk) + zdum 378 378 dh_i_sub(ji) = dh_i_sub(ji) + zdum 379 379 380 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoi c * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice ! Salt flux >0380 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0 381 381 ! clem: flux is sent to the ocean for simplicity 382 382 ! but salt should remain in the ice except 383 383 ! if all ice is melted. => must be corrected 384 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_ rdtice ! Heat flux [W.m-2], < 0385 386 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi c * a_i_1d(ji) * zdum * r1_rdtice ! Mass flux > 0384 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! Heat flux [W.m-2], < 0 385 386 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_Dt_ice ! Mass flux > 0 387 387 388 388 ! update remaining mass flux 389 zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi c389 zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi 390 390 391 391 ! record which layers have disappeared (for bottom melting) … … 409 409 ! remaining "potential" evap is sent to ocean 410 410 DO ji = 1, npti 411 wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_ rdtice ! <=0 (net evap for the ocean in kg.m-2.s-1)411 wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_Dt_ice ! <=0 (net evap for the ocean in kg.m-2.s-1) 412 412 END DO 413 413 … … 437 437 !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7 438 438 !--- zswi2 if dh/dt > 3.6e-7 439 zgrr = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_ rdtice , epsi10 ) )439 zgrr = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_Dt_ice , epsi10 ) ) 440 440 zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) ) 441 441 zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) … … 450 450 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 451 451 452 zEi = cpic* ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0)453 & - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp* ztmelts452 zEi = rcpi * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0) 453 & - rLfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp * ztmelts 454 454 455 455 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) … … 457 457 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 458 458 459 dh_i_bog(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi c) )459 dh_i_bog(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) ) 460 460 461 461 END DO 462 462 ! Contribution to Energy and Salt Fluxes 463 zfmdt = - rhoi c* dh_i_bog(ji) ! Mass flux x time step (kg/m2, < 0)464 465 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_ rdtice ! Heat flux to the ocean [W.m-2], >0466 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_ rdtice ! Heat flux used in this process [W.m-2], <0467 468 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi c * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_rdtice ! Salt flux, <0469 470 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi c * a_i_1d(ji) * dh_i_bog(ji) * r1_rdtice ! Mass flux, <0463 zfmdt = - rhoi * dh_i_bog(ji) ! Mass flux x time step (kg/m2, < 0) 464 465 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux to the ocean [W.m-2], >0 466 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice ! Heat flux used in this process [W.m-2], <0 467 468 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_Dt_ice ! Salt flux, <0 469 470 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_Dt_ice ! Mass flux, <0 471 471 472 472 ! update heat content (J.m-2) and layer thickness 473 eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi c)473 eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi) 474 474 h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bog(ji) 475 475 … … 489 489 IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting 490 490 491 zEi = - e_i_1d(ji,jk) * r1_rhoi c! Specific enthalpy of melting ice (J/kg, <0)491 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0) 492 492 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 493 493 ! set up at 0 since no energy is needed to melt water...(it is already melted) … … 497 497 dh_i_itm (ji) = dh_i_itm(ji) + zdeltah(ji,jk) 498 498 499 zfmdt = - zdeltah(ji,jk) * rhoi c! Mass flux x time step > 0500 501 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_ rdtice ! Heat flux to the ocean [W.m-2], <0499 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 500 501 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 502 502 ! ice enthalpy zEi is "sent" to the ocean 503 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi c * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux503 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux 504 504 ! using s_i_1d and not sz_i_1d(jk) is ok 505 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi c * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux505 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 506 506 507 507 ! update heat content (J.m-2) and layer thickness … … 511 511 ELSE !-- Basal melting 512 512 513 zEi = - e_i_1d(ji,jk) * r1_rhoi c! Specific enthalpy of melting ice (J/kg, <0)513 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0) 514 514 zEw = rcp * ztmelts ! Specific enthalpy of meltwater (J/kg, <0) 515 515 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) … … 517 517 zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0) 518 518 519 zdeltah(ji,jk) = - zfmdt * r1_rhoi c! Gross thickness change519 zdeltah(ji,jk) = - zfmdt * r1_rhoi ! Gross thickness change 520 520 521 521 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 522 522 523 zq_bo(ji) = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoi c * zdE )! update available heat. MAX is necessary for roundup errors524 525 dh_i_bom(ji) = dh_i_bom(ji) + zdeltah(ji,jk) ! Update basal melt526 527 zfmdt = - zdeltah(ji,jk) * rhoi c! Mass flux x time step > 0523 zq_bo(ji) = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat. MAX is necessary for roundup errors 524 525 dh_i_bom(ji) = dh_i_bom(ji) + zdeltah(ji,jk) ! Update basal melt 526 527 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 528 528 529 529 zQm = zfmdt * zEw ! Heat exchanged with ocean 530 530 531 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_ rdtice ! Heat flux to the ocean [W.m-2], <0532 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_ rdtice ! Heat used in this process [W.m-2], >0533 534 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoi c * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux531 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 532 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice ! Heat used in this process [W.m-2], >0 533 534 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux 535 535 ! using s_i_1d and not sz_i_1d(jk) is ok 536 536 537 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi c * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux537 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 538 538 539 539 ! update heat content (J.m-2) and layer thickness … … 565 565 566 566 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1) ! update available heat (J.m-2) 567 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_ rdtice ! Heat used to melt snow, W.m-2 (>0)568 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos n * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice ! Mass flux567 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_Dt_ice ! Heat used to melt snow, W.m-2 (>0) 568 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice ! Mass flux 569 569 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,1) 570 570 ! 571 571 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 572 hfx_out_1d(ji) = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_ rdtice572 hfx_out_1d(ji) = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice 573 573 574 574 IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) … … 580 580 ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level, 581 581 ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice 582 z1_rho = 1._wp / ( rhos n+rau0-rhoic)582 z1_rho = 1._wp / ( rhos + rho0 - rhoi ) 583 583 DO ji = 1, npti 584 584 ! 585 dh_snowice(ji) = MAX( 0._wp , ( rhos n * h_s_1d(ji) + (rhoic-rau0) * h_i_1d(ji) ) * z1_rho )585 dh_snowice(ji) = MAX( 0._wp , ( rhos * h_s_1d(ji) + (rhoi-rho0) * h_i_1d(ji) ) * z1_rho ) 586 586 587 587 h_i_1d(ji) = h_i_1d(ji) + dh_snowice(ji) … … 589 589 590 590 ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0) 591 zfmdt = ( rhos n - rhoic) * dh_snowice(ji) ! <0591 zfmdt = ( rhos - rhoi ) * dh_snowice(ji) ! <0 592 592 zEw = rcp * sst_1d(ji) 593 593 zQm = zfmdt * zEw 594 594 595 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_ rdtice ! Heat flux596 597 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_ rdtice ! Salt flux595 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux 596 597 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_Dt_ice ! Salt flux 598 598 599 599 ! Case constant salinity in time: virtual salt flux to keep salinity constant 600 600 IF( nn_icesal /= 2 ) THEN 601 sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt * r1_rdtice & ! put back sss_m into the ocean602 & - s_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoi c * r1_rdtice ! and get rn_icesal from the ocean601 sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt * r1_Dt_ice & ! put back sss_m into the ocean 602 & - s_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice ! and get rn_icesal from the ocean 603 603 ENDIF 604 604 605 605 ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume 606 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoi c * r1_rdtice607 wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos n * r1_rdtice606 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice 607 wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_Dt_ice 608 608 609 609 ! update heat content (J.m-2) and layer thickness … … 627 627 e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk) 628 628 ! recalculate t_s_1d from e_s_1d 629 t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos n * r1_cpic + lfus * r1_cpic)629 t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_cpi + rLfus * r1_cpi ) 630 630 END DO 631 631 END DO
Note: See TracChangeset
for help on using the changeset viewer.
