Changeset 5123 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
- Timestamp:
- 2015-03-04T17:06:03+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r4990 r5123 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 … … 70 69 71 70 REAL(wp) :: ztmelts ! local scalar 72 REAL(wp) :: z dh, zfdum !71 REAL(wp) :: zfdum 73 72 REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment 74 73 REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads … … 91 90 REAL(wp), POINTER, DIMENSION(:) :: zq_su ! heat for surface ablation (J.m-2) 92 91 REAL(wp), POINTER, DIMENSION(:) :: zq_bo ! heat for bottom ablation (J.m-2) 93 REAL(wp), POINTER, DIMENSION(:) :: zq_1cat ! corrected heat in case 1-cat and hmelt>15cm (J.m-2)94 92 REAL(wp), POINTER, DIMENSION(:) :: zq_rema ! remaining heat at the end of the routine (J.m-2) 95 REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2)93 REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2) 96 94 INTEGER , POINTER, DIMENSION(:) :: icount ! number of layers vanished by melting 97 95 … … 107 105 REAL(wp), POINTER, DIMENSION(:) :: zq_s ! total snow enthalpy (J.m-3) 108 106 109 ! mass and salt flux (clem) 110 REAL(wp) :: zdvres, zswitch_sal 107 REAL(wp) :: zswitch_sal 111 108 112 109 ! Heat conservation … … 115 112 !!------------------------------------------------------------------ 116 113 117 ! Discriminate between varying salinity (n um_sal=2) and prescribed cases (other values)118 SELECT CASE( n um_sal ) ! varying salinity or not114 ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values) 115 SELECT CASE( nn_icesal ) ! varying salinity or not 119 116 CASE( 1, 3, 4 ) ; zswitch_sal = 0 ! prescribed salinity profile 120 117 CASE( 2 ) ; zswitch_sal = 1 ! varying salinity profile 121 118 END SELECT 122 119 123 CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_ 1cat, zq_rema )120 CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 124 121 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 125 122 CALL wrk_alloc( jpij, nlay_i+1, zdeltah, zh_i ) … … 130 127 131 128 zqprec (:) = 0._wp ; zq_su (:) = 0._wp ; zq_bo (:) = 0._wp ; zf_tt (:) = 0._wp 132 zq_ 1cat(:) = 0._wp ; zq_rema(:) = 0._wp129 zq_rema(:) = 0._wp 133 130 134 131 zh_s (:) = 0._wp … … 148 145 DO jk = 1, nlay_i 149 146 DO ji = kideb, kiut 150 h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i )147 h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i 151 148 qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) 152 149 ENDDO … … 159 156 DO ji = kideb, kiut 160 157 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 161 ztmelts = rswitch * rt t + ( 1._wp - rswitch ) * rtt158 ztmelts = rswitch * rt0 + ( 1._wp - rswitch ) * rt0 162 159 163 160 zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) … … 174 171 !------------------------------------------------------------------------------! 175 172 DO ji = kideb, kiut 176 IF( t_s_1d(ji,1) > rt t) THEN !!! Internal melting173 IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting 177 174 ! Contribution to heat flux to the ocean [W.m-2], < 0 178 175 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice … … 182 179 ht_s_1d(ji) = 0._wp 183 180 q_s_1d (ji,1) = 0._wp 184 t_s_1d (ji,1) = rt t181 t_s_1d (ji,1) = rt0 185 182 END IF 186 183 END DO … … 191 188 ! 192 189 DO ji = kideb, kiut 193 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s )190 zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 194 191 END DO 195 192 ! … … 202 199 DO jk = 1, nlay_i 203 200 DO ji = kideb, kiut 204 zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i )201 zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i 205 202 zqh_i(ji) = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) 206 203 END DO … … 230 227 !----------- 231 228 ! thickness change 232 zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )** betas ) / at_i_1d(ji)233 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice /rhosn229 zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**rn_betas ) / at_i_1d(ji) 230 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice * r1_rhosn 234 231 ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) 235 zqprec (ji) = rhosn * ( cpic * ( rt t- MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )232 zqprec (ji) = rhosn * ( cpic * ( rt0 - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus ) 236 233 IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 237 234 ! heat flux from snow precip (>0, W.m-2) … … 258 255 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) ) 259 256 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 )257 zh_s (ji) = ht_s_1d(ji) * r1_nlay_s 261 258 262 259 ENDIF … … 279 276 280 277 ! updates available heat + thickness 281 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )278 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 282 279 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 283 280 … … 314 311 DO ji = kideb, kiut 315 312 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 316 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s )313 zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 317 314 END DO ! ji 318 315 … … 327 324 q_s_1d(ji,jk) = ( 1._wp - rswitch ) / MAX( ht_s_1d(ji), epsi20 ) * & 328 325 & ( ( MAX( 0._wp, dh_s_tot(ji) ) ) * zqprec(ji) + & 329 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rt t- t_s_1d(ji,jk) ) + lfus ) )326 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 330 327 zq_s(ji) = zq_s(ji) + q_s_1d(ji,jk) 331 328 END DO … … 338 335 DO jk = 1, nlay_i 339 336 DO ji = kideb, kiut 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) + rt t! Melting point of layer k [K]337 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 338 339 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer k [K] 343 340 344 341 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] … … 348 345 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 349 346 350 zdeltah(ji,jk) = - zfmdt / rhoic! Melt of layer jk [m, <0]347 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Melt of layer jk [m, <0] 351 348 352 349 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] … … 408 405 ! -> need for an iterative procedure, which converges quickly 409 406 410 IF ( n um_sal == 2 ) THEN407 IF ( nn_icesal == 2 ) THEN 411 408 num_iter_max = 5 412 409 ELSE … … 414 411 ENDIF 415 412 416 ! clem debug.Just to be sure that enthalpy at nlay_i+1 is null413 ! Just to be sure that enthalpy at nlay_i+1 is null 417 414 DO ji = kideb, kiut 418 415 q_i_1d(ji,nlay_i+1) = 0._wp … … 440 437 + ( 1. - zswitch_sal ) * sm_i_1d(ji) 441 438 ! New ice growth 442 ztmelts = - tmut * s_i_new(ji) + rt t! New ice melting point (K)439 ztmelts = - tmut * s_i_new(ji) + rt0 ! New ice melting point (K) 443 440 444 441 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 445 442 446 443 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) 447 & - lfus * ( 1.0 - ( ztmelts - rt t ) / ( zt_i_new - rtt) ) &448 & + rcp * ( ztmelts-rt t)444 & - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) ) & 445 & + rcp * ( ztmelts-rt0 ) 449 446 450 447 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) … … 467 464 zfmdt = - rhoic * dh_i_bott(ji) ! Mass flux x time step (kg/m2, < 0) 468 465 469 ztmelts = - tmut * s_i_new(ji) + rt t! New ice melting point (K)466 ztmelts = - tmut * s_i_new(ji) + rt0 ! New ice melting point (K) 470 467 471 468 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 472 469 473 470 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) 474 & - lfus * ( 1.0 - ( ztmelts - rt t ) / ( zt_i_new - rtt) ) &475 & + rcp * ( ztmelts-rt t)471 & - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) ) & 472 & + rcp * ( ztmelts-rt0 ) 476 473 477 474 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) … … 503 500 DO jk = nlay_i, 1, -1 504 501 DO ji = kideb, kiut 505 IF( zf_tt(ji) >= 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared fromsurface melting506 507 ztmelts = - tmut * s_i_1d(ji,jk) + rt t! Melting point of layer jk (K)502 IF( zf_tt(ji) >= 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared by surface melting 503 504 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer jk (K) 508 505 509 506 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 510 507 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) - rt t) ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0)508 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 509 510 !!zEw = rcp * ( t_i_1d(ji,jk) - rt0 ) ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) 514 511 515 512 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) … … 538 535 ELSE !!! Basal melting 539 536 540 zEi = - q_i_1d(ji,jk) /rhoic ! Specific enthalpy of melting ice (J/kg, <0)541 542 zEw = rcp * ( ztmelts - rt t )! Specific enthalpy of meltwater (J/kg, <0)543 544 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)545 546 zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0)547 548 zdeltah(ji,jk) = - zfmdt / rhoic! Gross thickness change537 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 538 539 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of meltwater (J/kg, <0) 540 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 549 546 550 547 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change … … 576 573 577 574 ENDIF 578 END DO ! ji 579 END DO ! jk 580 581 !------------------------------------------------------------------------------! 582 ! Excessive ablation in a 1-category model 583 ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15) 584 !------------------------------------------------------------------------------! 585 ! ??? keep ??? 586 ! clem bug: I think this should be included above, so we would not have to 587 ! track heat/salt/mass fluxes backwards 588 ! IF( jpl == 1 ) THEN 589 ! DO ji = kideb, kiut 590 ! IF( zf_tt(ji) >= 0._wp ) THEN 591 ! zdh = MAX( hmelt , dh_i_bott(ji) ) 592 ! zdvres = zdh - dh_i_bott(ji) ! >=0 593 ! dh_i_bott(ji) = zdh 594 ! 595 ! ! excessive energy is sent to lateral ablation 596 ! 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 598 ! 599 ! ! correct salt and mass fluxes 600 ! 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 602 ! ENDIF 603 ! END DO 604 ! ENDIF 575 END DO 576 END DO 605 577 606 578 !------------------------------------------- … … 635 607 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 636 608 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 637 hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_ 1cat(ji) + zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice609 hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 638 610 639 611 IF( ln_nicep .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) … … 655 627 ! Salinity of snow ice 656 628 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 657 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) /rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji)629 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 658 630 659 631 ! entrapment during snow ice formation 660 632 ! new salinity difference stored (to be used in limthd_ent.F90) 661 IF ( n um_sal == 2 ) THEN633 IF ( nn_icesal == 2 ) THEN 662 634 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) ) 663 635 ! salinity dif due to snow-ice formation … … 703 675 DO ji = kideb, kiut 704 676 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 ) * rt t677 t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0 706 678 END DO ! ji 707 679 … … 712 684 q_s_1d(ji,jk) = ( 1.0 - rswitch ) * q_s_1d(ji,jk) 713 685 ! recalculate t_s_1d from q_s_1d 714 t_s_1d(ji,jk) = rt t+ ( 1._wp - rswitch ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )686 t_s_1d(ji,jk) = rt0 + ( 1._wp - rswitch ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 715 687 END DO 716 688 END DO 717 718 CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_ 1cat, zq_rema )689 690 CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 719 691 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 720 692 CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i )
Note: See TracChangeset
for help on using the changeset viewer.