- Timestamp:
- 2017-07-12T16:51:20+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r8316 r8325 104 104 INTEGER , POINTER, DIMENSION(:,:) :: icount ! number of layers vanished by melting 105 105 106 REAL(wp), POINTER, DIMENSION(:) :: z qh_i ! total ice heat content (J.m-2)106 REAL(wp), POINTER, DIMENSION(:) :: zeh_i ! total ice heat content (J.m-2) 107 107 REAL(wp), POINTER, DIMENSION(:) :: zsnw ! distribution of snow after wind blowing 108 108 … … 121 121 122 122 CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema ) 123 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, z qh_i )123 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zeh_i ) 124 124 CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 125 125 CALL wrk_alloc( jpij, nlay_i, icount ) … … 127 127 zqprec (:) = 0._wp ; zq_su (:) = 0._wp ; zq_bo (:) = 0._wp ; zf_tt(:) = 0._wp 128 128 zq_rema (:) = 0._wp ; zsnw (:) = 0._wp ; zevap_rema(:) = 0._wp ; 129 zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; z qh_i(:) = 0._wp129 zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zeh_i(:) = 0._wp 130 130 131 131 zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp … … 134 134 ! Initialize enthalpy at nlay_i+1 135 135 DO ji = kideb, kiut 136 q_i_1d(ji,nlay_i+1) = 0._wp136 e_i_1d(ji,nlay_i+1) = 0._wp 137 137 END DO 138 138 139 139 ! initialize layer thicknesses and enthalpies 140 140 h_i_old (:,0:nlay_i+1) = 0._wp 141 qh_i_old(:,0:nlay_i+1) = 0._wp141 eh_i_old(:,0:nlay_i+1) = 0._wp 142 142 DO jk = 1, nlay_i 143 143 DO ji = kideb, kiut 144 144 h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i 145 qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk)145 eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk) 146 146 ENDDO 147 147 ENDDO … … 167 167 IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting 168 168 ! Contribution to heat flux to the ocean [W.m-2], < 0 169 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice169 hfx_res_1d(ji) = hfx_res_1d(ji) + e_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 170 170 ! Contribution to mass flux 171 171 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 172 172 ! updates 173 173 ht_s_1d(ji) = 0._wp 174 q_s_1d (ji,1) = 0._wp174 e_s_1d (ji,1) = 0._wp 175 175 t_s_1d (ji,1) = rt0 176 176 END IF … … 184 184 DO ji = kideb, kiut 185 185 zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i 186 z qh_i(ji) = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk)186 zeh_i(ji) = zeh_i(ji) + e_i_1d(ji,jk) * zh_i(ji,jk) 187 187 END DO 188 188 END DO … … 248 248 ! thickness change 249 249 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 250 rswitch = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) )251 zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 )250 rswitch = rswitch * ( MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) ) 251 zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( e_s_1d(ji,jk), epsi20 ) 252 252 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting 253 253 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 254 254 ! heat used to melt snow(W.m-2, >0) 255 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice255 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d(ji,jk) * r1_rdtice 256 256 ! snow melting only = water into the ocean (then without snow precip) 257 257 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 258 258 ! updates available heat + thickness 259 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )259 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 260 260 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 261 261 END DO … … 274 274 ! Heat flux by sublimation [W.m-2], < 0 (sublimate first snow that had fallen, then pre-existing snow) 275 275 zdeltah(ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 276 hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1) &276 hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) & 277 277 & ) * a_i_1d(ji) * r1_rdtice 278 278 ! Mass flux by sublimation … … 298 298 DO ji = kideb,kiut 299 299 rswitch = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) ) 300 q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) * &300 e_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) * & 301 301 & ( ( zdh_s_pre(ji) ) * zqprec(ji) + & 302 302 & ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) … … 314 314 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 315 315 316 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0]316 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 317 317 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 318 318 ! set up at 0 since no energy is needed to melt water...(it is already melted) … … 336 336 ELSE !!! Surface melting 337 337 338 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0]338 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 339 339 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] 340 340 zdE = zEi - zEw ! Specific enthalpy difference < 0 … … 377 377 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * sm_i_1d(ji) * r1_rdtice 378 378 ! Heat flux [W.m-2], < 0 379 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * q_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice379 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice 380 380 ! Mass flux > 0 381 381 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice … … 392 392 393 393 ! update heat content (J.m-2) and layer thickness 394 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)394 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 395 395 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 396 396 END DO … … 464 464 dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 465 465 466 q_i_1d(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0)466 e_i_1d(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0) 467 467 468 468 ENDIF … … 502 502 503 503 ! update heat content (J.m-2) and layer thickness 504 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)504 eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * e_i_1d(ji,nlay_i+1) 505 505 h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji) 506 506 ENDIF … … 519 519 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 520 520 521 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)521 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 522 522 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 523 523 ! set up at 0 since no energy is needed to melt water...(it is already melted) … … 539 539 540 540 ! update heat content (J.m-2) and layer thickness 541 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)541 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 542 542 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 543 543 544 544 ELSE !!! Basal melting 545 545 546 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)546 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 547 547 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of meltwater (J/kg, <0) 548 548 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) … … 575 575 576 576 ! update heat content (J.m-2) and layer thickness 577 qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)577 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 578 578 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 579 579 ENDIF … … 599 599 zq_rema(ji) = zq_su(ji) + zq_bo(ji) 600 600 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) ! =1 if snow 601 rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) )602 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 )601 rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) ) 602 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 ) 603 603 zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 604 604 dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) 605 605 ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1) 606 606 607 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1) ! update available heat (J.m-2)607 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1) ! update available heat (J.m-2) 608 608 ! heat used to melt snow 609 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * q_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0)609 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0) 610 610 ! Contribution to mass flux 611 611 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice … … 661 661 662 662 ! update heat content (J.m-2) and layer thickness 663 qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw663 eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw 664 664 h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 665 665 … … 679 679 ! mask enthalpy 680 680 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) ) 681 q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk)682 ! recalculate t_s_1d from q_s_1d683 t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )681 e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk) 682 ! recalculate t_s_1d from e_s_1d 683 t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 684 684 END DO 685 685 END DO … … 689 689 690 690 CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema ) 691 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, z qh_i )691 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zeh_i ) 692 692 CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 693 693 CALL wrk_dealloc( jpij, nlay_i, icount )
Note: See TracChangeset
for help on using the changeset viewer.
