- Timestamp:
- 2015-06-04T16:12:19+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r5134_UKMO4_CF_compliance/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r5134 r5350 86 86 REAL(wp) :: zsstK ! SST in Kelvin 87 87 88 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness89 88 REAL(wp), POINTER, DIMENSION(:) :: zqprec ! energy of fallen snow (J.m-3) 90 89 REAL(wp), POINTER, DIMENSION(:) :: zq_su ! heat for surface ablation (J.m-2) … … 92 91 REAL(wp), POINTER, DIMENSION(:) :: zq_rema ! remaining heat at the end of the routine (J.m-2) 93 92 REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2) 94 INTEGER , POINTER, DIMENSION(:) :: icount ! number of layers vanished by melting95 93 96 94 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_mel ! snow melt … … 100 98 REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah 101 99 REAL(wp), POINTER, DIMENSION(:,:) :: zh_i ! ice layer thickness 100 INTEGER , POINTER, DIMENSION(:,:) :: icount ! number of layers vanished by melting 102 101 103 102 REAL(wp), POINTER, DIMENSION(:) :: zqh_i ! total ice heat content (J.m-2) … … 118 117 END SELECT 119 118 120 CALL wrk_alloc( jpij, z h_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema )119 CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 121 120 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 122 CALL wrk_alloc( jpij, nlay_i +1, zdeltah, zh_i )123 CALL wrk_alloc( jpij, icount )121 CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 122 CALL wrk_alloc( jpij, nlay_i, icount ) 124 123 125 124 dh_i_surf (:) = 0._wp ; dh_i_bott (:) = 0._wp ; dh_snowice(:) = 0._wp … … 129 128 zq_rema(:) = 0._wp 130 129 131 zh_s (:) = 0._wp132 130 zdh_s_pre(:) = 0._wp 133 131 zdh_s_mel(:) = 0._wp … … 138 136 zh_i (:,:) = 0._wp 139 137 zdeltah (:,:) = 0._wp 140 icount (:) = 0 138 icount (:,:) = 0 139 140 ! Initialize enthalpy at nlay_i+1 141 DO ji = kideb, kiut 142 q_i_1d(ji,nlay_i+1) = 0._wp 143 END DO 141 144 142 145 ! initialize layer thicknesses and enthalpies … … 155 158 ! 156 159 DO ji = kideb, kiut 157 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )158 ztmelts = rswitch * rt0 + ( 1._wp - rswitch ) * rt0159 160 160 zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 161 161 zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 162 162 163 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts) )163 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 164 164 zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 165 165 END DO … … 187 187 !------------------------------------------------------------! 188 188 ! 189 DO ji = kideb, kiut190 zh_s(ji) = ht_s_1d(ji) * r1_nlay_s191 END DO192 !193 189 DO jk = 1, nlay_s 194 190 DO ji = kideb, kiut 195 zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji)191 zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s 196 192 END DO 197 193 END DO … … 222 218 ! Martin Vancoppenolle, December 2006 223 219 220 zdeltah(:,:) = 0._wp 224 221 DO ji = kideb, kiut 225 222 !----------- … … 236 233 ! mass flux, <0 237 234 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice 238 ! update thickness239 ht_s_1d (ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) )240 235 241 236 !--------------------- … … 243 238 !--------------------- 244 239 ! thickness change 245 IF( zdh_s_pre(ji) > 0._wp ) THEN246 240 rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 247 zd h_s_mel (ji) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )248 zd h_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting241 zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 242 zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting 249 243 ! heat used to melt snow (W.m-2, >0) 250 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zd h_s_mel(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice244 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 251 245 ! snow melting only = water into the ocean (then without snow precip), >0 252 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_mel(ji) * r1_rdtice 253 254 ! updates available heat + thickness 255 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) ) 256 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) ) 257 zh_s (ji) = ht_s_1d(ji) * r1_nlay_s 258 259 ENDIF 260 END DO 261 262 ! If heat still available, then melt more snow 263 zdeltah(:,:) = 0._wp ! important 246 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 247 ! updates available heat + precipitations after melting 248 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) ) 249 zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1) 250 251 ! update thickness 252 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 253 END DO 254 255 ! If heat still available (zq_su > 0), then melt more snow 256 zdeltah(:,:) = 0._wp 264 257 DO jk = 1, nlay_s 265 258 DO ji = kideb, kiut … … 268 261 rswitch = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) ) 269 262 zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) 270 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting263 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting 271 264 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 272 265 ! heat used to melt snow(W.m-2, >0) … … 274 267 ! snow melting only = water into the ocean (then without snow precip) 275 268 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 276 277 269 ! updates available heat + thickness 278 270 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 279 271 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 280 281 272 END DO 282 273 END DO … … 286 277 !---------------------- 287 278 ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 288 ! clem comment: not counted in mass exchange in limsbc since this is an exchange with atm. (not ocean)279 ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean) 289 280 ! clem comment: ice should also sublimate 281 zdeltah(:,:) = 0._wp 290 282 IF( lk_cpl ) THEN 291 283 ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) … … 294 286 ! forced mode: snow thickness change due to sublimation 295 287 DO ji = kideb, kiut 296 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - parsub *qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice )288 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 297 289 ! Heat flux by sublimation [W.m-2], < 0 298 290 ! sublimate first snow that had fallen, then pre-existing snow 299 zcoeff = ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) * zqprec(ji) + & 300 & ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_1d(ji,1) ) & 301 & * a_i_1d(ji) * r1_rdtice 302 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff 291 zdeltah(ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 292 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) & 293 & ) * a_i_1d(ji) * r1_rdtice 303 294 ! Mass flux by sublimation 304 295 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 305 296 ! new snow thickness 306 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 297 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 298 ! update precipitations after sublimation and correct sublimation 299 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 300 zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 307 301 END DO 308 302 ENDIF … … 310 304 ! --- Update snow diags --- ! 311 305 DO ji = kideb, kiut 312 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 313 zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 314 END DO ! ji 306 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 307 END DO 315 308 316 309 !------------------------------------------- … … 323 316 rswitch = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) ) 324 317 q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) * & 325 & ( ( MAX( 0._wp, dh_s_tot(ji) )) * zqprec(ji) + &326 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )318 & ( ( zdh_s_pre(ji) ) * zqprec(ji) + & 319 & ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 327 320 zq_s(ji) = zq_s(ji) + q_s_1d(ji,jk) 328 321 END DO … … 334 327 zdeltah(:,:) = 0._wp ! important 335 328 DO jk = 1, nlay_i 336 DO ji = kideb, kiut 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] 340 341 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] 342 343 zdE = zEi - zEw ! Specific enthalpy difference < 0 344 345 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 346 347 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Melt of layer jk [m, <0] 348 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] 350 351 zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 352 353 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt 354 355 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 356 357 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 358 359 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 360 sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 361 362 ! Contribution to heat flux [W.m-2], < 0 363 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 364 365 ! Total heat flux used in this process [W.m-2], > 0 366 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 367 368 ! Contribution to mass flux 369 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 370 329 DO ji = kideb, kiut 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 386 END IF 371 387 ! record which layers have disappeared (for bottom melting) 372 388 ! => icount=0 : no layer has vanished 373 389 ! => icount=5 : 5 layers have vanished 374 rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )375 icount(ji ) = icount(ji) +NINT( rswitch )376 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) ) 377 393 378 394 ! update heat content (J.m-2) and layer thickness … … 405 421 ! -> need for an iterative procedure, which converges quickly 406 422 407 IF ( nn_icesal == 2 ) THEN 408 num_iter_max = 5 409 ELSE 410 num_iter_max = 1 411 ENDIF 412 413 ! Just to be sure that enthalpy at nlay_i+1 is null 414 DO ji = kideb, kiut 415 q_i_1d(ji,nlay_i+1) = 0._wp 416 END DO 423 num_iter_max = 1 424 IF( nn_icesal == 2 ) num_iter_max = 5 417 425 418 426 ! Iterative procedure … … 483 491 484 492 ! Contribution to salt flux, <0 485 sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_1d(ji) * zfmdt* r1_rdtice493 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice 486 494 487 495 ! Contribution to mass flux, <0 … … 500 508 DO jk = nlay_i, 1, -1 501 509 DO ji = kideb, kiut 502 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 503 511 504 512 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer jk (K) … … 507 515 508 516 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)511 512 517 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 513 518 ! set up at 0 since no energy is needed to melt water...(it is already melted) 514 515 zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 516 ! 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 517 521 518 522 dh_i_bott (ji) = dh_i_bott(ji) + zdeltah(ji,jk) 519 523 520 zfmdt = - zdeltah(ji,jk) * rhoic 524 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 521 525 522 526 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) … … 524 528 525 529 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 526 sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic* r1_rdtice530 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 527 531 528 532 ! Contribution to mass flux … … 535 539 ELSE !!! Basal melting 536 540 537 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 546 547 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 548 550 549 zq_bo(ji) 550 551 dh_i_bott(ji) 552 553 zfmdt 554 555 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 556 558 557 559 ! Contribution to heat flux to the ocean [W.m-2], <0 558 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 559 561 560 562 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 561 sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic* r1_rdtice563 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 562 564 563 565 ! Total heat flux used in this process [W.m-2], >0 564 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 565 567 566 568 ! Contribution to mass flux 567 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 568 570 569 571 ! update heat content (J.m-2) and layer thickness … … 595 597 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 ) 596 598 zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 597 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1)598 599 dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) 599 600 ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1) … … 622 623 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic ) ) 623 624 624 ht_i_1d(ji) 625 ht_s_1d(ji) 625 ht_i_1d(ji) = ht_i_1d(ji) + dh_snowice(ji) 626 ht_s_1d(ji) = ht_s_1d(ji) - dh_snowice(ji) 626 627 627 628 ! Salinity of snow ice … … 669 670 ! Update temperature, energy 670 671 !------------------------------------------- 671 !clem bug: we should take snow into account here672 672 DO ji = kideb, kiut 673 673 rswitch = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) … … 688 688 WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 689 689 690 CALL wrk_dealloc( jpij, z h_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema )690 CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 691 691 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 692 CALL wrk_dealloc( jpij, nlay_i +1, zdeltah, zh_i )693 CALL wrk_dealloc( jpij, icount )692 CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 693 CALL wrk_dealloc( jpij, nlay_i, icount ) 694 694 ! 695 695 !
Note: See TracChangeset
for help on using the changeset viewer.
