Changeset 5202 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
- Timestamp:
- 2015-04-07T17:40:16+02:00 (9 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r5167 r5202 91 91 REAL(wp), POINTER, DIMENSION(:) :: zq_rema ! remaining heat at the end of the routine (J.m-2) 92 92 REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2) 93 INTEGER , POINTER, DIMENSION(:) :: icount ! number of layers vanished by melting94 93 95 94 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_mel ! snow melt … … 99 98 REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah 100 99 REAL(wp), POINTER, DIMENSION(:,:) :: zh_i ! ice layer thickness 100 INTEGER , POINTER, DIMENSION(:,:) :: icount ! number of layers vanished by melting 101 101 102 102 REAL(wp), POINTER, DIMENSION(:) :: zqh_i ! total ice heat content (J.m-2) … … 119 119 CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 120 120 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 121 CALL wrk_alloc( jpij, nlay_i +1, zdeltah, zh_i )122 CALL wrk_alloc( jpij, icount )121 CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 122 CALL wrk_alloc( jpij, nlay_i, icount ) 123 123 124 124 dh_i_surf (:) = 0._wp ; dh_i_bott (:) = 0._wp ; dh_snowice(:) = 0._wp … … 136 136 zh_i (:,:) = 0._wp 137 137 zdeltah (:,:) = 0._wp 138 icount (: )= 0138 icount (:,:) = 0 139 139 140 140 ! Initialize enthalpy at nlay_i+1 … … 328 328 DO jk = 1, nlay_i 329 329 DO ji = kideb, kiut 330 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 331 332 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer k [K] 333 334 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] 335 336 zdE = zEi - zEw ! Specific enthalpy difference < 0 337 338 IF( zdE < 0._wp ) THEN 339 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 340 ELSE 341 zfmdt = rhoic * zh_i(ji,jk) !!! internal melting 342 zdE = 0._wp ! All the layer melts if t_i(jk) = Tmelt (i.e. zdE = 0) 330 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer k [K] 331 332 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 333 334 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 335 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 336 ! set up at 0 since no energy is needed to melt water...(it is already melted) 337 zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 338 ! this should normally not happen, but sometimes, heat diffusion leads to this 339 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 340 341 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt 342 343 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 344 345 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 346 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 347 348 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 349 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 350 351 ! Contribution to mass flux 352 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 353 354 ELSE !!! Surface melting 355 356 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 357 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] 358 zdE = zEi - zEw ! Specific enthalpy difference < 0 359 360 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 361 362 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Melt of layer jk [m, <0] 363 364 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] 365 366 zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 367 368 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt 369 370 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 371 372 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 373 374 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 375 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 376 377 ! Contribution to heat flux [W.m-2], < 0 378 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 379 380 ! Total heat flux used in this process [W.m-2], > 0 381 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 382 383 ! Contribution to mass flux 384 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 385 343 386 END IF 344 345 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Melt of layer jk [m, <0]346 347 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]348 349 zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat350 351 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt352 353 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0]354 355 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0]356 357 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)358 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice359 360 ! Contribution to heat flux [W.m-2], < 0361 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice362 363 ! Total heat flux used in this process [W.m-2], > 0364 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice365 366 ! Contribution to mass flux367 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice368 369 387 ! record which layers have disappeared (for bottom melting) 370 388 ! => icount=0 : no layer has vanished 371 389 ! => icount=5 : 5 layers have vanished 372 rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )373 icount(ji ) = icount(ji) +NINT( rswitch )374 zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )390 rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 391 icount(ji,jk) = NINT( rswitch ) 392 zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 375 393 376 394 ! update heat content (J.m-2) and layer thickness … … 490 508 DO jk = nlay_i, 1, -1 491 509 DO ji = kideb, kiut 492 IF( zf_tt(ji) > 0._wp .AND. jk > icount(ji ) ) THEN ! do not calculate where layer has already disappeared by surface melting510 IF( zf_tt(ji) > 0._wp .AND. jk > icount(ji,jk) ) THEN ! do not calculate where layer has already disappeared by surface melting 493 511 494 512 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer jk (K) … … 497 515 498 516 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 499 500 !!zEw = rcp * ( t_i_1d(ji,jk) - rt0 ) ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0)501 502 517 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 503 518 ! set up at 0 since no energy is needed to melt water...(it is already melted) 504 505 zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 506 ! this should normally not happen, but sometimes, heat diffusion leads to this 519 zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 520 ! this should normally not happen, but sometimes, heat diffusion leads to this 507 521 508 522 dh_i_bott (ji) = dh_i_bott(ji) + zdeltah(ji,jk) 509 523 510 zfmdt = - zdeltah(ji,jk) * rhoic 524 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 511 525 512 526 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) … … 525 539 ELSE !!! Basal melting 526 540 527 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 528 529 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of meltwater (J/kg, <0) 530 531 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 532 533 zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0) 534 535 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Gross thickness change 536 537 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 541 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 542 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of meltwater (J/kg, <0) 543 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 544 545 zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0) 546 547 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Gross thickness change 548 549 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 538 550 539 zq_bo(ji) 540 541 dh_i_bott(ji) 542 543 zfmdt 544 545 zQm 551 zq_bo(ji) = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors 552 553 dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) ! Update basal melt 554 555 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 556 557 zQm = zfmdt * zEw ! Heat exchanged with ocean 546 558 547 559 ! Contribution to heat flux to the ocean [W.m-2], <0 548 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice560 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 549 561 550 562 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 551 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice563 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 552 564 553 565 ! Total heat flux used in this process [W.m-2], >0 554 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice566 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 555 567 556 568 ! Contribution to mass flux 557 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice569 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 558 570 559 571 ! update heat content (J.m-2) and layer thickness … … 678 690 CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 679 691 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 680 CALL wrk_dealloc( jpij, nlay_i +1, zdeltah, zh_i )681 CALL wrk_dealloc( jpij, icount )692 CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 693 CALL wrk_dealloc( jpij, nlay_i, icount ) 682 694 ! 683 695 !
Note: See TracChangeset
for help on using the changeset viewer.
