Changeset 4990 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
- Timestamp:
- 2014-12-15T17:42:49+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r4873 r4990 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 125 CALL wrk_alloc( jpij, nlay_i+1, zdeltah, zh_i ) 133 126 CALL wrk_alloc( jpij, icount ) … … 148 141 zh_i (:,:) = 0._wp 149 142 zdeltah (:,:) = 0._wp 150 zintermelt(:) = 0._wp151 143 icount (:) = 0 152 144 … … 166 158 ! 167 159 DO ji = kideb, kiut 168 zinda= 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(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)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) 173 165 174 166 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts ) ) … … 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) … … 276 268 DO ji = kideb, kiut 277 269 ! thickness change 278 zindh= 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )279 zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) + epsi20) )280 zdeltah (ji,jk) = - zindh * zindq* zq_su(ji) / MAX( q_s_1d(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) … … 332 324 DO jk = 1, nlay_s 333 325 DO ji = kideb,kiut 334 zindh= MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) + epsi20 ) )335 q_s_1d(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_1d(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 329 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) ) … … 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 … … 516 508 517 509 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 518 zintermelt(ji) = 1._wp519 510 520 511 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) … … 603 594 ! 604 595 ! ! excessive energy is sent to lateral ablation 605 ! zinda= MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_1d(ji) - epsi20 ) )606 ! zq_1cat(ji) = zinda* rhoic * lfus * at_i_1d(ji) / MAX( 1._wp - at_i_1d(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 … … 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_1d(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_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(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_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(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 … … 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_1d(ji) ) )714 t_su_1d(ji) = zindh * t_su_1d(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_1d(ji) ) )721 q_s_1d(ji,jk) = ( 1.0 - zinda) * q_s_1d(ji,jk)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) 722 713 ! recalculate t_s_1d from q_s_1d 723 t_s_1d(ji,jk) = rtt + ( 1._wp - zinda) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )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 720 CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i ) 731 721 CALL wrk_dealloc( jpij, icount )
Note: See TracChangeset
for help on using the changeset viewer.