Changeset 5407 for trunk/NEMOGCM/NEMO/LIM_SRC_3
- Timestamp:
- 2015-06-11T21:13:22+02:00 (9 years ago)
- Location:
- trunk/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 8 edited
-
limistate.F90 (modified) (2 diffs)
-
limitd_th.F90 (modified) (28 diffs)
-
limsbc.F90 (modified) (10 diffs)
-
limthd.F90 (modified) (11 diffs)
-
limthd_dh.F90 (modified) (9 diffs)
-
limthd_dif.F90 (modified) (2 diffs)
-
limwri.F90 (modified) (2 diffs)
-
thd_ice.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r5202 r5407 117 117 118 118 ! basal temperature (considered at freezing point) 119 t_bo(:,:) = ( eos_fzp( tsn(:,:,1,jp_sal) ) + rt0 ) * tmask(:,:,1)119 t_bo(:,:) = ( eos_fzp( sss_m(:,:) ) + rt0 ) * tmask(:,:,1) 120 120 121 121 IF( ln_iceini ) THEN … … 127 127 DO jj = 1, jpj ! ice if sst <= t-freez + ttest 128 128 DO ji = 1, jpi 129 IF( ( tsn(ji,jj,1,jp_tem) - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN129 IF( ( sst_m(ji,jj) - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN 130 130 zswitch(ji,jj) = 0._wp * tmask(ji,jj,1) ! no ice 131 131 ELSE -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r5202 r5407 91 91 !!------------------------------------------------------------------ 92 92 93 CALL wrk_alloc( jpi,jpj, zremap_flag ) ! integer94 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) ! integer93 CALL wrk_alloc( jpi,jpj, zremap_flag ) 94 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) 95 95 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 96 96 CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice ) 97 97 CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 98 98 CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax ) 99 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer99 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 100 100 CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 101 101 … … 128 128 rswitch = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) ) !0 if no ice and 1 if yes 129 129 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 130 rswitch = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) ) !0 if no ice and 1 if yes130 rswitch = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) ) 131 131 zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 132 zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)132 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) ! clem: useless IF statement? 133 133 END DO 134 134 END DO … … 172 172 ! 173 173 zhbnew(ii,ij,jl) = hi_max(jl) 174 IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN174 IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 175 175 !interpolate between adjacent category growth rates 176 176 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 177 177 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 178 ELSEIF 178 ELSEIF( a_i_b(ii,ij,jl) > epsi10) THEN 179 179 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 180 ELSEIF 180 ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN 181 181 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 182 182 ENDIF … … 187 187 ii = nind_i(ji) 188 188 ij = nind_j(ji) 189 IF( a_i(ii,ij,jl) > epsi10 .AND. ht_i(ii,ij,jl) >= zhbnew(ii,ij,jl) ) THEN 189 190 ! clem: we do not want ht_i to be too close to either HR or HL otherwise a division by nearly 0 is possible 191 ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 192 IF ( a_i(ii,ij,jl ) > epsi10 .AND. ht_i(ii,ij,jl ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) THEN 190 193 zremap_flag(ii,ij) = 0 191 ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < = zhbnew(ii,ij,jl) ) THEN194 ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) THEN 192 195 zremap_flag(ii,ij) = 0 193 196 ENDIF 194 197 195 198 !- 4.3 Check that each zhbnew does not exceed maximal values hi_max 199 IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 196 200 IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 197 IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 201 ! clem bug: why is not the following instead? 202 !!IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 203 !!IF( zhbnew(ii,ij,jl) > hi_max(jl ) ) zremap_flag(ii,ij) = 0 204 198 205 END DO 199 206 … … 219 226 DO jj = 1, jpj 220 227 DO ji = 1, jpi 221 zhb0(ji,jj) = hi_max(0) ! 0eme 222 zhb1(ji,jj) = hi_max(1) ! 1er 223 224 zhbnew(ji,jj,klbnd-1) = 0._wp 228 zhb0(ji,jj) = hi_max(0) 229 zhb1(ji,jj) = hi_max(1) 225 230 226 231 IF( a_i(ji,jj,kubnd) > epsi10 ) THEN 227 232 zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) ) 228 233 ELSE 229 zhbnew(ji,jj,kubnd) = hi_max(kubnd) 230 !!? clem bug: since hi_max(jpl)=99, this limit is very high 231 !!? but I think it is erased in fitline subroutine 234 !clem bug zhbnew(ji,jj,kubnd) = hi_max(kubnd) 235 zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) ! not used anyway 236 ENDIF 237 238 ! clem: we do not want ht_i_b to be too close to either HR or HL otherwise a division by nearly 0 is possible 239 ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 240 IF ( zht_i_b(ji,jj,klbnd) < ( zhb0(ji,jj) + epsi10 ) ) THEN 241 zremap_flag(ji,jj) = 0 242 ELSEIF( zht_i_b(ji,jj,klbnd) > ( zhb1(ji,jj) - epsi10 ) ) THEN 243 zremap_flag(ji,jj) = 0 232 244 ENDIF 233 245 … … 248 260 249 261 IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 262 250 263 zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 251 IF( zdh0 < 0.0 ) THEN !remove area from category 1 264 265 IF( zdh0 < 0.0 ) THEN !remove area from category 1 252 266 zdh0 = MIN( -zdh0, hi_max(klbnd) ) 253 254 267 !Integrate g(1) from 0 to dh0 to estimate area melted 255 268 zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd) 269 256 270 IF( zetamax > 0.0 ) THEN 257 zx1 = zetamax 258 zx2 = 0.5 * zetamax * zetamax 259 zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 260 ! Constrain new thickness <= ht_i 261 zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! zdamax > 0 262 !ice area lost due to melting of thin ice 263 zda0 = MIN( zda0, zdamax ) 264 271 zx1 = zetamax 272 zx2 = 0.5 * zetamax * zetamax 273 zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 ! ice area removed 274 zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! Constrain new thickness <= ht_i 275 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting 276 ! of thin ice (zdamax > 0) 265 277 ! Remove area, conserving volume 266 278 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) … … 269 281 ENDIF 270 282 271 ELSE ! if ice accretion ! a_i > epsi10; zdh0 > 0 283 ELSE ! if ice accretion zdh0 > 0 284 ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 272 285 zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) ) 273 ! zhbnew was 0, and is shifted to the right to account for thin ice 274 ! growth in openwater (F0 = f1) 275 ENDIF ! zdh0 276 277 ENDIF ! a_i > epsi10 286 ENDIF 287 288 ENDIF 278 289 279 290 END DO … … 303 314 304 315 IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1 305 306 316 ! left and right integration limits in eta space 307 317 zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl) 308 zvetamax(ji) = MIN (zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl)318 zvetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) 309 319 zdonor(ii,ij,jl) = jl 310 320 311 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 312 321 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 313 322 ! left and right integration limits in eta space 314 323 zvetamin(ji) = 0.0 … … 316 325 zdonor(ii,ij,jl) = jl + 1 317 326 318 ENDIF ! zhbnew(jl) > hi_max(jl)327 ENDIF 319 328 320 329 zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin … … 333 342 334 343 END DO 335 END DO ! jl klbnd -> kubnd - 1344 END DO 336 345 337 346 !!---------------------------------------------------------------------------------------------- … … 375 384 ENDIF 376 385 377 CALL wrk_dealloc( jpi,jpj, zremap_flag ) ! integer378 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) ! integer386 CALL wrk_dealloc( jpi,jpj, zremap_flag ) 387 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) 379 388 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 380 389 CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice ) 381 390 CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 382 391 CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax ) 383 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer392 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 384 393 CALL wrk_dealloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 385 394 … … 406 415 INTEGER , DIMENSION(jpi,jpj), INTENT(in ) :: zremap_flag ! 407 416 ! 408 INTEGER :: ji,jj! horizontal indices417 INTEGER :: ji,jj ! horizontal indices 409 418 REAL(wp) :: zh13 ! HbL + 1/3 * (HbR - HbL) 410 419 REAL(wp) :: zh23 ! HbL + 2/3 * (HbR - HbL) … … 413 422 !!------------------------------------------------------------------ 414 423 ! 415 !416 424 DO jj = 1, jpj 417 425 DO ji = 1, jpi 418 426 ! 419 427 IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10 & 420 & .AND. hice(ji,jj) > 0._wp )THEN428 & .AND. hice(ji,jj) > 0._wp ) THEN 421 429 422 430 ! Initialize hL and hR 423 424 431 hL(ji,jj) = HbL(ji,jj) 425 432 hR(ji,jj) = HbR(ji,jj) 426 433 427 434 ! Change hL or hR if hice falls outside central third of range 428 429 435 zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) ) 430 436 zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) ) … … 435 441 436 442 ! Compute coefficients of g(eta) = g0 + g1*eta 437 438 443 zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj)) 439 444 zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr … … 442 447 g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) 443 448 ! 444 ELSE 449 ELSE ! remap_flag = .false. or a_i < epsi10 445 450 hL(ji,jj) = 0._wp 446 451 hR(ji,jj) = 0._wp 447 452 g0(ji,jj) = 0._wp 448 453 g1(ji,jj) = 0._wp 449 ENDIF ! a_i > epsi10454 ENDIF 450 455 ! 451 456 END DO … … 471 476 472 477 INTEGER :: ji, jj, jl, jl2, jl1, jk ! dummy loop indices 473 INTEGER :: ii, ij ! indices when changing from 2D-1D is done478 INTEGER :: ii, ij ! indices when changing from 2D-1D is done 474 479 475 480 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaTsfn … … 484 489 INTEGER, POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions 485 490 486 INTEGER :: nbrem ! number of cells with ice to transfer 487 488 LOGICAL :: zdaice_negative ! true if daice < -puny 489 LOGICAL :: zdvice_negative ! true if dvice < -puny 490 LOGICAL :: zdaice_greater_aicen ! true if daice > aicen 491 LOGICAL :: zdvice_greater_vicen ! true if dvice > vicen 491 INTEGER :: nbrem ! number of cells with ice to transfer 492 492 !!------------------------------------------------------------------ 493 493 494 494 CALL wrk_alloc( jpi,jpj,jpl, zaTsfn ) 495 495 CALL wrk_alloc( jpi,jpj, zworka ) 496 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer496 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 497 497 498 498 !---------------------------------------------------------------------------------------------- … … 504 504 END DO 505 505 506 !clem: I think the following is wrong (if enabled, it creates negative concentration/volume around -epsi10)507 ! !----------------------------------------------------------------------------------------------508 ! ! 2) Check for daice or dvice out of range, allowing for roundoff error509 ! !----------------------------------------------------------------------------------------------510 ! ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl511 ! ! has a small area, with h(n) very close to a boundary. Then512 ! ! the coefficients of g(h) are large, and the computed daice and513 ! ! dvice can be in error. If this happens, it is best to transfer514 ! ! either the entire category or nothing at all, depending on which515 ! ! side of the boundary hice(n) lies.516 ! !-----------------------------------------------------------------517 ! DO jl = klbnd, kubnd-1518 !519 ! zdaice_negative = .false.520 ! zdvice_negative = .false.521 ! zdaice_greater_aicen = .false.522 ! zdvice_greater_vicen = .false.523 !524 ! DO jj = 1, jpj525 ! DO ji = 1, jpi526 !527 ! IF (zdonor(ji,jj,jl) > 0) THEN528 ! jl1 = zdonor(ji,jj,jl)529 !530 ! IF (zdaice(ji,jj,jl) < 0.0) THEN531 ! IF (zdaice(ji,jj,jl) > -epsi10) THEN532 ! IF ( ( jl1 == jl .AND. ht_i(ji,jj,jl1) > hi_max(jl) ) .OR. &533 ! ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN534 ! zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category535 ! zdvice(ji,jj,jl) = v_i(ji,jj,jl1)536 ! ELSE537 ! zdaice(ji,jj,jl) = 0.0 ! shift no ice538 ! zdvice(ji,jj,jl) = 0.0539 ! ENDIF540 ! ELSE541 ! zdaice_negative = .true.542 ! ENDIF543 ! ENDIF544 !545 ! IF (zdvice(ji,jj,jl) < 0.0) THEN546 ! IF (zdvice(ji,jj,jl) > -epsi10 ) THEN547 ! IF ( ( jl1 == jl .AND. ht_i(ji,jj,jl1) > hi_max(jl) ) .OR. &548 ! ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN549 ! zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category550 ! zdvice(ji,jj,jl) = v_i(ji,jj,jl1)551 ! ELSE552 ! zdaice(ji,jj,jl) = 0.0 ! shift no ice553 ! zdvice(ji,jj,jl) = 0.0554 ! ENDIF555 ! ELSE556 ! zdvice_negative = .true.557 ! ENDIF558 ! ENDIF559 !560 ! ! If daice is close to aicen, set daice = aicen.561 ! IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN562 ! IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN563 ! zdaice(ji,jj,jl) = a_i(ji,jj,jl1)564 ! zdvice(ji,jj,jl) = v_i(ji,jj,jl1)565 ! ELSE566 ! zdaice_greater_aicen = .true.567 ! ENDIF568 ! ENDIF569 !570 ! IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN571 ! IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN572 ! zdaice(ji,jj,jl) = a_i(ji,jj,jl1)573 ! zdvice(ji,jj,jl) = v_i(ji,jj,jl1)574 ! ELSE575 ! zdvice_greater_vicen = .true.576 ! ENDIF577 ! ENDIF578 !579 ! ENDIF ! donor > 0580 ! END DO581 ! END DO582 !583 ! END DO584 !clem585 506 !------------------------------------------------------------------------------- 586 ! 3) Transfer volume and energy between categories507 ! 2) Transfer volume and energy between categories 587 508 !------------------------------------------------------------------------------- 588 509 … … 604 525 605 526 jl1 = zdonor(ii,ij,jl) 606 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi 20 ) )607 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi 20 ) * rswitch527 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi10 ) ) 528 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch 608 529 IF( jl1 == jl) THEN ; jl2 = jl1+1 609 530 ELSE ; jl2 = jl … … 613 534 ! Ice areas 614 535 !-------------- 615 616 536 a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 617 537 a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) … … 620 540 ! Ice volumes 621 541 !-------------- 622 623 542 v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl) 624 543 v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) … … 627 546 ! Snow volumes 628 547 !-------------- 629 630 548 zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij) 631 549 v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow … … 635 553 ! Snow heat content 636 554 !-------------------- 637 638 555 zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij) 639 556 e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow … … 643 560 ! Ice age 644 561 !-------------- 645 646 562 zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 647 563 oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice … … 651 567 ! Ice salinity 652 568 !-------------- 653 654 569 zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij) 655 570 smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice … … 659 574 ! Surface temperature 660 575 !--------------------- 661 662 576 zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 663 577 zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf … … 710 624 CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn ) 711 625 CALL wrk_dealloc( jpi,jpj, zworka ) 712 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer626 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 713 627 ! 714 628 END SUBROUTINE lim_itd_shiftice … … 859 773 ENDIF 860 774 ! 861 CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) ! interger775 CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) 862 776 CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice ) 863 777 CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final ) -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r5187 r5407 30 30 USE sbc_oce ! Surface boundary condition: ocean fields 31 31 USE sbccpl 32 USE oce , ONLY : fraqsr_1lev,sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass32 USE oce , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass 33 33 USE albedo ! albedo parameters 34 34 USE lbclnk ! ocean lateral boundary condition - MPP exchanges … … 94 94 !! - fr_i : ice fraction 95 95 !! - tn_ice : sea-ice surface temperature 96 !! - alb_ice : sea-ice albedo ( lk_cpl=T)96 !! - alb_ice : sea-ice albedo (only useful in coupled mode) 97 97 !! 98 98 !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90. … … 101 101 !! The ref should be Rousset et al., 2015 102 102 !!--------------------------------------------------------------------- 103 INTEGER, INTENT(in) :: kt ! number of iteration 104 INTEGER :: ji, jj, jl, jk ! dummy loop indices 105 REAL(wp) :: zemp ! local scalars 106 REAL(wp) :: zf_mass ! Heat flux associated with mass exchange ice->ocean (W.m-2) 107 REAL(wp) :: zfcm1 ! New solar flux received by the ocean 103 INTEGER, INTENT(in) :: kt ! number of iteration 104 INTEGER :: ji, jj, jl, jk ! dummy loop indices 105 REAL(wp) :: zqmass ! Heat flux associated with mass exchange ice->ocean (W.m-2) 106 REAL(wp) :: zqsr ! New solar flux received by the ocean 108 107 ! 109 108 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb_cs, zalb_os ! 2D/3D workspace … … 111 110 112 111 ! make calls for heat fluxes before it is modified 113 IF( iom_use('qsr_oce') ) CALL iom_put( "qsr_oce" , qsr(:,:) * pfrld(:,:) ) ! solar flux at ocean surface 114 IF( iom_use('qns_oce') ) CALL iom_put( "qns_oce" , qns(:,:) * pfrld(:,:) ) ! non-solar flux at ocean surface 115 IF( iom_use('qsr_ice') ) CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux at ice surface 116 IF( iom_use('qns_ice') ) CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! non-solar flux at ice surface 117 IF( iom_use('qtr_ice') ) CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux transmitted thru ice 118 IF( iom_use('qt_oce' ) ) CALL iom_put( "qt_oce" , ( qsr(:,:) + qns(:,:) ) * pfrld(:,:) ) 119 IF( iom_use('qt_ice' ) ) CALL iom_put( "qt_ice" , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) * a_i_b(:,:,:), dim=3 ) ) 112 IF( iom_use('qsr_oce') ) CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) ) ! solar flux at ocean surface 113 IF( iom_use('qns_oce') ) CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) ) ! non-solar flux at ocean surface 114 IF( iom_use('qsr_ice') ) CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux at ice surface 115 IF( iom_use('qns_ice') ) CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface 116 IF( iom_use('qtr_ice') ) CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux transmitted thru ice 117 IF( iom_use('qt_oce' ) ) CALL iom_put( "qt_oce" , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) ) 118 IF( iom_use('qt_ice' ) ) CALL iom_put( "qt_ice" , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) & 119 & * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) 120 IF( iom_use('qemp_oce' ) ) CALL iom_put( "qemp_oce" , qemp_oce(:,:) ) 121 IF( iom_use('qemp_ice' ) ) CALL iom_put( "qemp_ice" , qemp_ice(:,:) ) 120 122 121 123 ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) … … 126 128 ! heat flux at the ocean surface ! 127 129 !------------------------------------------! 128 ! Solar heat flux reaching the ocean = z fcm1(W.m-2)130 ! Solar heat flux reaching the ocean = zqsr (W.m-2) 129 131 !--------------------------------------------------- 130 IF( lk_cpl ) THEN 131 !!! LIM2 version zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) ) 132 zfcm1 = qsr_tot(ji,jj) 133 DO jl = 1, jpl 134 zfcm1 = zfcm1 + ( ftr_ice(ji,jj,jl) - qsr_ice(ji,jj,jl) ) * a_i_b(ji,jj,jl) 135 END DO 136 ELSE 137 !!! LIM2 version zqsr = pfrld(ji,jj) * qsr(ji,jj) + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj) 138 zfcm1 = pfrld(ji,jj) * qsr(ji,jj) 139 DO jl = 1, jpl 140 zfcm1 = zfcm1 + a_i_b(ji,jj,jl) * ftr_ice(ji,jj,jl) 141 END DO 142 ENDIF 132 zqsr = qsr_tot(ji,jj) 133 DO jl = 1, jpl 134 zqsr = zqsr - a_i_b(ji,jj,jl) * ( qsr_ice(ji,jj,jl) - ftr_ice(ji,jj,jl) ) 135 END DO 143 136 144 137 ! Total heat flux reaching the ocean = hfx_out (W.m-2) 145 138 !--------------------------------------------------- 146 z f_mass= hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC)147 hfx_out(ji,jj) = hfx_out(ji,jj) + z f_mass + zfcm1139 zqmass = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC) 140 hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr 148 141 149 142 ! Add the residual from heat diffusion equation (W.m-2) … … 153 146 ! New qsr and qns used to compute the oceanic heat flux at the next time step 154 147 !--------------------------------------------------- 155 qsr(ji,jj) = z fcm1156 qns(ji,jj) = hfx_out(ji,jj) - z fcm1148 qsr(ji,jj) = zqsr 149 qns(ji,jj) = hfx_out(ji,jj) - zqsr 157 150 158 151 !------------------------------------------! … … 167 160 ! Even if i see Ice melting as a FW and SALT flux 168 161 ! 169 ! computing freshwater exchanges at the ice/ocean interface170 IF( lk_cpl ) THEN171 zemp = emp_tot(ji,jj) & ! net mass flux over grid cell172 & - emp_ice(ji,jj) * ( 1._wp - pfrld(ji,jj) ) & ! minus the mass flux intercepted by sea ice173 & + sprecip(ji,jj) * ( pfrld(ji,jj) - pfrld(ji,jj)**rn_betas ) !174 ELSE175 zemp = emp(ji,jj) * pfrld(ji,jj) & ! evaporation over oceanic fraction176 & - tprecip(ji,jj) * ( 1._wp - pfrld(ji,jj) ) & ! all precipitation reach the ocean177 & + sprecip(ji,jj) * ( 1._wp - pfrld(ji,jj)**rn_betas ) ! except solid precip intercepted by sea-ice178 ENDIF179 180 162 ! mass flux from ice/ocean 181 163 wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) & … … 184 166 ! mass flux at the ocean/ice interface 185 167 fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) ) * r1_rdtice ! F/M mass flux save at least for biogeochemical model 186 emp(ji,jj) = zemp - wfx_ice(ji,jj) - wfx_snw(ji,jj)! mass flux + F/M mass flux (always ice/ocean mass exchange)168 emp(ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) ! mass flux + F/M mass flux (always ice/ocean mass exchange) 187 169 188 170 END DO … … 213 195 tn_ice(:,:,:) = t_su(:,:,:) ! Ice surface temperature 214 196 215 !------------------------------------------------! 216 ! Snow/ice albedo (only if sent to coupler) ! 217 !------------------------------------------------! 218 IF( lk_cpl ) THEN ! coupled case 219 220 CALL wrk_alloc( jpi, jpj, jpl, zalb_cs, zalb_os ) 221 222 CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 223 224 alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 225 226 CALL wrk_dealloc( jpi, jpj, jpl, zalb_cs, zalb_os ) 227 228 ENDIF 197 !------------------------------------------------------------------------! 198 ! Snow/ice albedo (only if sent to coupler, useless in forced mode) ! 199 !------------------------------------------------------------------------! 200 CALL wrk_alloc( jpi, jpj, jpl, zalb_cs, zalb_os ) 201 CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 202 alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 203 CALL wrk_dealloc( jpi, jpj, jpl, zalb_cs, zalb_os ) 229 204 230 205 ! conservation test … … 346 321 sice_0(:,:) = 2._wp 347 322 END WHERE 348 ENDIF349 350 IF( .NOT. ln_rstart ) THEN351 fraqsr_1lev(:,:) = 1._wp352 323 ENDIF 353 324 ! -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r5385 r5407 22 22 USE phycst ! physical constants 23 23 USE dom_oce ! ocean space and time domain variables 24 USE oce , ONLY : fraqsr_1lev25 24 USE ice ! LIM: sea-ice variables 26 25 USE sbc_oce ! Surface boundary condition: ocean fields … … 28 27 USE thd_ice ! LIM thermodynamic sea-ice variables 29 28 USE dom_ice ! LIM sea-ice domain 30 USE domvvl ! domain: variable volume level31 29 USE limthd_dif ! LIM: thermodynamics, vertical diffusion 32 30 USE limthd_dh ! LIM: thermodynamics, ice and snow thickness variation … … 50 48 PRIVATE 51 49 52 PUBLIC lim_thd ! called by limstp module53 PUBLIC lim_thd_init ! called by sbc_lim_init50 PUBLIC lim_thd ! called by limstp module 51 PUBLIC lim_thd_init ! called by sbc_lim_init 54 52 55 53 !! * Substitutions … … 92 90 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient 93 91 ! 94 REAL(wp), POINTER, DIMENSION(:,:) :: zqsr, zqns95 92 !!------------------------------------------------------------------- 96 CALL wrk_alloc( jpi,jpj, zqsr, zqns )97 93 98 94 IF( nn_timing == 1 ) CALL timing_start('limthd') … … 136 132 ! 2) Partial computation of forcing for the thermodynamic sea ice model. ! 137 133 !-----------------------------------------------------------------------------! 138 139 !--- Ocean solar and non solar fluxes to be used in zqld140 IF ( .NOT. lk_cpl ) THEN ! --- forced case, fluxes to the lead are the same as over the ocean141 !142 zqsr(:,:) = qsr(:,:) ; zqns(:,:) = qns(:,:)143 !144 ELSE ! --- coupled case, fluxes to the lead are total - intercepted145 !146 zqsr(:,:) = qsr_tot(:,:) ; zqns(:,:) = qns_tot(:,:)147 !148 DO jl = 1, jpl149 DO jj = 1, jpj150 DO ji = 1, jpi151 zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * a_i_b(ji,jj,jl)152 zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * a_i_b(ji,jj,jl)153 END DO154 END DO155 END DO156 !157 ENDIF158 159 134 DO jj = 1, jpj 160 135 DO ji = 1, jpi … … 167 142 ! ! temperature and turbulent mixing (McPhee, 1992) 168 143 ! 169 170 144 ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- ! 171 ! REMARK valid at least in forced mode from clem 172 ! precip is included in qns but not in qns_ice 173 IF ( lk_cpl ) THEN 174 zqld = tmask(ji,jj,1) * rdt_ice * & 175 & ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) & ! pfrld already included in coupled mode 176 & + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj) * & ! heat content of precip 177 & ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus ) & 178 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) ) 179 ELSE 180 zqld = tmask(ji,jj,1) * rdt_ice * & 181 & ( pfrld(ji,jj) * ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) ) & 182 & + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj) * & ! heat content of precip 183 & ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus ) & 184 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) ) 185 ENDIF 145 zqld = tmask(ji,jj,1) * rdt_ice * & 146 & ( pfrld(ji,jj) * qsr_oce(ji,jj) * frq_m(ji,jj) + pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 186 147 187 148 ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! … … 210 171 ! Net heat flux on top of ice-ocean [W.m-2] 211 172 ! ----------------------------------------- 212 ! heat flux at the ocean surface + precip 213 ! + heat flux at the ice surface 214 hfx_in(ji,jj) = hfx_in(ji,jj) & 215 ! heat flux above the ocean 216 & + pfrld(ji,jj) * ( zqns(ji,jj) + zqsr(ji,jj) ) & 217 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 218 & + ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus ) & 219 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) & 220 ! heat flux above the ice 221 & + SUM( a_i_b(ji,jj,:) * ( qns_ice(ji,jj,:) + qsr_ice(ji,jj,:) ) ) 173 hfx_in(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 222 174 223 175 ! ----------------------------------------------------------------------------- 224 ! Net heat flux that is retroceded to the ocean or taken from the ocean[W.m-2]176 ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2] 225 177 ! ----------------------------------------------------------------------------- 226 178 ! First step here : non solar + precip - qlead - qturb 227 179 ! Second step in limthd_dh : heat remaining if total melt (zq_rema) 228 180 ! Third step in limsbc : heat from ice-ocean mass exchange (zf_mass) + solar 229 hfx_out(ji,jj) = hfx_out(ji,jj) & 230 ! Non solar heat flux received by the ocean 231 & + pfrld(ji,jj) * zqns(ji,jj) & 232 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 233 & + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj) & 234 & * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus ) & 235 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) & 236 ! heat flux taken from the ocean where there is open water ice formation 237 & - qlead(ji,jj) * r1_rdtice & 238 ! heat flux taken from the ocean during bottom growth/melt (fhld should be 0 while bott growth) 239 & - at_i(ji,jj) * fhtur(ji,jj) & 240 & - at_i(ji,jj) * fhld(ji,jj) 241 181 hfx_out(ji,jj) = pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) & ! Non solar heat flux received by the ocean 182 & - qlead(ji,jj) * r1_rdtice & ! heat flux taken from the ocean where there is open water ice formation 183 & - at_i(ji,jj) * fhtur(ji,jj) & ! heat flux taken by turbulence 184 & - at_i(ji,jj) * fhld(ji,jj) ! heat flux taken during bottom growth/melt 185 ! (fhld should be 0 while bott growth) 242 186 END DO 243 187 END DO … … 412 356 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 413 357 414 CALL wrk_dealloc( jpi,jpj, zqsr, zqns )415 416 358 !------------------------------------------------------------------------------| 417 359 ! 6) Transport of ice between thickness categories. | … … 472 414 END SUBROUTINE lim_thd 473 415 416 474 417 SUBROUTINE lim_thd_temp( kideb, kiut ) 475 418 !!----------------------------------------------------------------------- … … 570 513 END DO 571 514 572 CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:), jpi, jpj, npb(1:nbpb) )515 CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) ) 573 516 CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 574 517 CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb), fr1_i0 , jpi, jpj, npb(1:nbpb) ) … … 576 519 CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 577 520 CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 578 IF( .NOT. lk_cpl ) THEN 579 CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 580 CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 581 ENDIF 521 CALL tab_2d_1d( nbpb, evap_ice_1d (1:nbpb), evap_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 582 522 CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 583 523 CALL tab_2d_1d( nbpb, t_bo_1d (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) ) -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r5202 r5407 29 29 PRIVATE 30 30 31 PUBLIC lim_thd_dh ! called by lim_thd 31 PUBLIC lim_thd_dh ! called by lim_thd 32 PUBLIC lim_thd_snwblow ! called in sbcblk/sbccpl and here 33 34 INTERFACE lim_thd_snwblow 35 MODULE PROCEDURE lim_thd_snwblow_1d, lim_thd_snwblow_2d 36 END INTERFACE 32 37 33 38 !!---------------------------------------------------------------------- … … 71 76 REAL(wp) :: zfdum 72 77 REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment 73 REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads 74 REAL(wp) :: zs_snic ! snow-ice salinity 78 REAL(wp) :: zs_snic ! snow-ice salinity 75 79 REAL(wp) :: zswi1 ! switch for computation of bottom salinity 76 80 REAL(wp) :: zswi12 ! switch for computation of bottom salinity … … 103 107 REAL(wp), POINTER, DIMENSION(:) :: zqh_s ! total snow heat content (J.m-2) 104 108 REAL(wp), POINTER, DIMENSION(:) :: zq_s ! total snow enthalpy (J.m-3) 109 REAL(wp), POINTER, DIMENSION(:) :: zsnw ! distribution of snow after wind blowing 105 110 106 111 REAL(wp) :: zswitch_sal … … 118 123 119 124 CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 120 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )125 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s, zsnw ) 121 126 CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 122 127 CALL wrk_alloc( jpij, nlay_i, icount ) … … 219 224 220 225 zdeltah(:,:) = 0._wp 226 CALL lim_thd_snwblow( 1. - at_i_1d, zsnw ) ! snow distribution over ice after wind blowing 221 227 DO ji = kideb, kiut 222 228 !----------- … … 224 230 !----------- 225 231 ! thickness change 226 zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**rn_betas ) / at_i_1d(ji) 227 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice * r1_rhosn 228 ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) 229 zqprec (ji) = rhosn * ( cpic * ( rt0 - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus ) 232 zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji) 233 ! enthalpy of the precip (>0, J.m-3) 234 zqprec (ji) = - qprec_ice_1d(ji) 230 235 IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 231 236 ! heat flux from snow precip (>0, W.m-2) … … 280 285 ! clem comment: ice should also sublimate 281 286 zdeltah(:,:) = 0._wp 282 IF( lk_cpl ) THEN 283 ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) 284 zdh_s_sub(:) = 0._wp 285 ELSE 286 ! forced mode: snow thickness change due to sublimation 287 DO ji = kideb, kiut 288 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 289 ! Heat flux by sublimation [W.m-2], < 0 290 ! sublimate first snow that had fallen, then pre-existing snow 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 294 ! Mass flux by sublimation 295 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 296 ! new snow thickness 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) 301 END DO 302 ENDIF 303 287 ! coupled mode: sublimation is set to 0 (evap_ice = 0) until further notice 288 ! forced mode: snow thickness change due to sublimation 289 DO ji = kideb, kiut 290 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 291 ! Heat flux by sublimation [W.m-2], < 0 292 ! sublimate first snow that had fallen, then pre-existing snow 293 zdeltah(ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 294 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) & 295 & ) * a_i_1d(ji) * r1_rdtice 296 ! Mass flux by sublimation 297 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 298 ! new snow thickness 299 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 300 ! update precipitations after sublimation and correct sublimation 301 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 302 zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 303 END DO 304 304 305 ! --- Update snow diags --- ! 305 306 DO ji = kideb, kiut … … 689 690 690 691 CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 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, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s, zsnw ) 692 693 CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 693 694 CALL wrk_dealloc( jpij, nlay_i, icount ) … … 695 696 ! 696 697 END SUBROUTINE lim_thd_dh 698 699 700 !!-------------------------------------------------------------------------- 701 !! INTERFACE lim_thd_snwblow 702 !! ** Purpose : Compute distribution of precip over the ice 703 !!-------------------------------------------------------------------------- 704 SUBROUTINE lim_thd_snwblow_2d( pin, pout ) 705 REAL(wp), DIMENSION(:,:), INTENT(in) :: pin ! previous fraction lead ( pfrld or (1. - a_i_b) ) 706 REAL(wp), DIMENSION(:,:), INTENT(out) :: pout 707 pout = ( 1._wp - ( pin )**rn_betas ) 708 END SUBROUTINE lim_thd_snwblow_2d 709 710 SUBROUTINE lim_thd_snwblow_1d( pin, pout ) 711 REAL(wp), DIMENSION(:), INTENT(in) :: pin 712 REAL(wp), DIMENSION(:), INTENT(out) :: pout 713 pout = ( 1._wp - ( pin )**rn_betas ) 714 END SUBROUTINE lim_thd_snwblow_1d 715 697 716 698 717 #else -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r5385 r5407 24 24 USE wrk_nemo ! work arrays 25 25 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 26 USE sbc_oce, ONLY : lk_cpl27 26 28 27 IMPLICIT NONE … … 745 744 !-------------------------------------------------------------------------! 746 745 DO ji = kideb, kiut 747 ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)748 IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) )749 746 ! ! surface ice conduction flux 750 747 isnow(ji) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) -
trunk/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r5202 r5407 176 176 CALL iom_put( "utau_ice" , utau_ice ) ! wind stress over ice along i-axis at I-point 177 177 CALL iom_put( "vtau_ice" , vtau_ice ) ! wind stress over ice along j-axis at I-point 178 CALL iom_put( "snowpre" , sprecip 178 CALL iom_put( "snowpre" , sprecip * 86400. ) ! snow precipitation 179 179 CALL iom_put( "micesalt" , smt_i ) ! mean ice salinity 180 180 … … 232 232 CALL iom_put ('hfxdif' , hfx_dif(:,:) ) ! 233 233 CALL iom_put ('hfxopw' , hfx_opw(:,:) ) ! 234 CALL iom_put ('hfxtur' , fhtur(:,:) * at_i(:,:) ) ! turbulent heat flux at ice base234 CALL iom_put ('hfxtur' , fhtur(:,:) * SUM(a_i_b(:,:,:), dim=3) ) ! turbulent heat flux at ice base 235 235 CALL iom_put ('hfxdhc' , diag_heat(:,:) ) ! Heat content variation in snow and ice 236 236 CALL iom_put ('hfxspr' , hfx_spr(:,:) ) ! Heat content of snow precip -
trunk/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
r5167 r5407 89 89 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fhld_1d !: <==> the 2D fhld 90 90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dqns_ice_1d !: <==> the 2D dqns_ice 91 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qla_ice_1d !: <==> the 2D qla_ice 92 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: dqla_ice_1d !: <==> the 2D dqla_ice 93 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: tatm_ice_1d !: <==> the 2D tatm_ice 91 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: evap_ice_1d !: <==> the 2D evap_ice 92 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qprec_ice_1d !: <==> the 2D qprec_ice 94 93 ! ! to reintegrate longwave flux inside the ice thermodynamics 95 94 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: i0 !: fraction of radiation transmitted to the ice … … 153 152 & fhld_1d (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d (jpij) , wfx_bom_1d(jpij) , & 154 153 & wfx_sum_1d(jpij) , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d(jpij) , & 155 & dqns_ice_1d(jpij) , qla_ice_1d (jpij) , dqla_ice_1d(jpij) ,&156 & tatm_ice_1d(jpij), i0 (jpij) , &154 & dqns_ice_1d(jpij) , evap_ice_1d (jpij), & 155 & qprec_ice_1d(jpij), i0 (jpij) , & 157 156 & sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij), & 158 157 & sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , &
Note: See TracChangeset
for help on using the changeset viewer.
