Changeset 11822 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/iceitd.F90
- Timestamp:
- 2019-10-29T11:41:36+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/iceitd.F90
r10069 r11822 21 21 USE ice1D ! sea-ice: thermodynamic variables 22 22 USE ice ! sea-ice: variables 23 USE icevar ! sea-ice: operations 23 24 USE icectl ! sea-ice: conservation tests 24 25 USE icetab ! sea-ice: convert 1D<=>2D … … 87 88 88 89 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 90 IF( ln_icediachk ) CALL ice_cons2D (0, 'iceitd_rem', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 89 91 90 92 !----------------------------------------------------------------------------------------------- 91 93 ! 1) Identify grid cells with ice 92 94 !----------------------------------------------------------------------------------------------- 95 at_i(:,:) = SUM( a_i, dim=3 ) 96 ! 93 97 npti = 0 ; nptidx(:) = 0 94 98 DO jj = 1, jpj … … 207 211 CALL itd_glinear( zhb0(1:npti) , zhb1(1:npti) , h_ib_1d(1:npti) , a_i_1d(1:npti) , & ! in 208 212 & g0 (1:npti,1), g1 (1:npti,1), hL (1:npti,1), hR (1:npti,1) ) ! out 209 213 ! 210 214 ! Area lost due to melting of thin ice 211 215 DO ji = 1, npti … … 214 218 ! 215 219 zdh0 = h_i_1d(ji) - h_ib_1d(ji) 216 IF( zdh0 < 0.0 ) THEN ! remove area from category 1220 IF( zdh0 < 0.0 ) THEN ! remove area from category 1 217 221 zdh0 = MIN( -zdh0, hi_max(1) ) 218 222 !Integrate g(1) from 0 to dh0 to estimate area melted … … 222 226 zx1 = zetamax 223 227 zx2 = 0.5 * zetamax * zetamax 224 zda0 = g1(ji,1) * zx2 + g0(ji,1) * zx1 228 zda0 = g1(ji,1) * zx2 + g0(ji,1) * zx1 ! ice area removed 225 229 zdamax = a_i_1d(ji) * (1.0 - h_i_1d(ji) / h_ib_1d(ji) ) ! Constrain new thickness <= h_i 226 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting 227 ! of thin ice (zdamax > 0) 230 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting of thin ice (zdamax > 0) 228 231 ! Remove area, conserving volume 229 232 h_i_1d(ji) = h_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 ) … … 249 252 ! --- g(h) for each thickness category --- ! 250 253 CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti) , a_i_1d(1:npti) , & ! in 251 & g0 (1:npti,jl ), g1 (1:npti,jl), hL (1:npti,jl), hR(1:npti,jl) ) ! out254 & g0 (1:npti,jl ), g1 (1:npti,jl), hL (1:npti,jl), hR (1:npti,jl) ) ! out 252 255 ! 253 256 END DO … … 313 316 ! 314 317 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 318 IF( ln_icediachk ) CALL ice_cons2D (1, 'iceitd_rem', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 315 319 ! 316 320 END SUBROUTINE ice_itd_rem … … 344 348 DO ji = 1, npti 345 349 ! 346 IF( paice(ji) > epsi10 .AND. phice(ji) > 0._wp) THEN350 IF( paice(ji) > epsi10 .AND. phice(ji) > epsi10 ) THEN 347 351 ! 348 352 ! Initialize hL and hR … … 389 393 REAL(wp), DIMENSION(:,:), INTENT(in) :: pdvice ! ice volume transferred across boundary 390 394 ! 391 INTEGER :: ji, j j, jl, jk! dummy loop indices392 INTEGER :: ii, ij, jl2, jl1! local integers395 INTEGER :: ji, jl, jk ! dummy loop indices 396 INTEGER :: jl2, jl1 ! local integers 393 397 REAL(wp) :: ztrans ! ice/snow transferred 394 REAL(wp), DIMENSION(jpij) :: zworka, zworkv ! workspace 395 REAL(wp), DIMENSION(jpij,jpl) :: zaTsfn ! - - 398 REAL(wp), DIMENSION(jpij) :: zworka, zworkv ! workspace 399 REAL(wp), DIMENSION(jpij,jpl) :: zaTsfn ! - - 400 REAL(wp), DIMENSION(jpij,nlay_i,jpl) :: ze_i_2d 401 REAL(wp), DIMENSION(jpij,nlay_s,jpl) :: ze_s_2d 396 402 !!------------------------------------------------------------------ 397 403 … … 405 411 CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 406 412 CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 413 DO jl = 1, jpl 414 DO jk = 1, nlay_s 415 CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 416 END DO 417 DO jk = 1, nlay_i 418 CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 419 END DO 420 END DO 421 ! to correct roundoff errors on a_i 422 CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti), rn_amax_2d ) 407 423 408 424 !---------------------------------------------------------------------------------------------- … … 435 451 ELSE ; zworka(ji) = 0._wp 436 452 ENDIF 437 !438 ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20)439 ! because of truncation error ( i.e. 1. - 1. /= 0 )440 ! I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90)441 453 ! 442 454 a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl) ! Ice areas … … 476 488 ! 477 489 DO jk = 1, nlay_s !--- Snow heat content 478 !479 490 DO ji = 1, npti 480 ii = MOD( nptidx(ji) - 1, jpi ) + 1481 ij = ( nptidx(ji) - 1 ) / jpi + 1482 491 ! 483 492 jl1 = kdonor(ji,jl) … … 487 496 ELSE ; jl2 = jl 488 497 ENDIF 489 ! 490 ztrans = e_s(ii,ij,jk,jl1) * zworkv(ji) 491 e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans 492 e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans 498 ztrans = ze_s_2d(ji,jk,jl1) * zworkv(ji) 499 ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) - ztrans 500 ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ztrans 493 501 ENDIF 494 502 END DO … … 497 505 DO jk = 1, nlay_i !--- Ice heat content 498 506 DO ji = 1, npti 499 ii = MOD( nptidx(ji) - 1, jpi ) + 1500 ij = ( nptidx(ji) - 1 ) / jpi + 1501 507 ! 502 508 jl1 = kdonor(ji,jl) … … 506 512 ELSE ; jl2 = jl 507 513 ENDIF 508 ! 509 ztrans = e_i(ii,ij,jk,jl1) * zworkv(ji) 510 e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans 511 e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans 514 ztrans = ze_i_2d(ji,jk,jl1) * zworkv(ji) 515 ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) - ztrans 516 ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + ztrans 512 517 ENDIF 513 518 END DO … … 515 520 ! 516 521 END DO ! boundaries, 1 to jpl-1 522 523 !------------------- 524 ! 3) roundoff errors 525 !------------------- 526 ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 527 ! because of truncation error ( i.e. 1. - 1. /= 0 ) 528 CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 529 530 ! at_i must be <= rn_amax 531 zworka(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 ) 532 DO jl = 1, jpl 533 WHERE( zworka(1:npti) > rn_amax_1d(1:npti) ) & 534 & a_i_2d(1:npti,jl) = a_i_2d(1:npti,jl) * rn_amax_1d(1:npti) / zworka(1:npti) 535 END DO 517 536 518 537 !------------------------------------------------------------------------------- 519 ! 3) Update ice thickness and temperature538 ! 4) Update ice thickness and temperature 520 539 !------------------------------------------------------------------------------- 521 540 WHERE( a_i_2d(1:npti,:) >= epsi20 ) … … 536 555 CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 537 556 CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 557 DO jl = 1, jpl 558 DO jk = 1, nlay_s 559 CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 560 END DO 561 DO jk = 1, nlay_i 562 CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 563 END DO 564 END DO 538 565 ! 539 566 END SUBROUTINE itd_shiftice … … 558 585 ! 559 586 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' 587 ! 588 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 589 IF( ln_icediachk ) CALL ice_cons2D (0, 'iceitd_reb', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 560 590 ! 561 591 jdonor(:,:) = 0 … … 635 665 END DO 636 666 ! 667 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 668 IF( ln_icediachk ) CALL ice_cons2D (1, 'iceitd_reb', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 669 ! 637 670 END SUBROUTINE ice_itd_reb 638 671 … … 655 688 REWIND( numnam_ice_ref ) ! Namelist namitd in reference namelist : Parameters for ice 656 689 READ ( numnam_ice_ref, namitd, IOSTAT = ios, ERR = 901) 657 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namitd in reference namelist' , lwp)690 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namitd in reference namelist' ) 658 691 REWIND( numnam_ice_cfg ) ! Namelist namitd in configuration namelist : Parameters for ice 659 692 READ ( numnam_ice_cfg, namitd, IOSTAT = ios, ERR = 902 ) 660 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namitd in configuration namelist' , lwp)693 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namitd in configuration namelist' ) 661 694 IF(lwm) WRITE( numoni, namitd ) 662 695 !
Note: See TracChangeset
for help on using the changeset viewer.
