Changeset 14072 for NEMO/trunk/src/ICE/iceitd.F90
- Timestamp:
- 2020-12-04T08:48:38+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/iceitd.F90 (modified) (34 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/iceitd.F90
r14005 r14072 18 18 !!---------------------------------------------------------------------- 19 19 USE dom_oce ! ocean domain 20 USE phycst ! physical constants 20 USE phycst ! physical constants 21 21 USE ice1D ! sea-ice: thermodynamic variables 22 22 USE ice ! sea-ice: variables … … 66 66 !! after thermodynamic growth of ice thickness 67 67 !! 68 !! ** Method : Linear remapping 68 !! ** Method : Linear remapping 69 69 !! 70 70 !! References : W.H. Lipscomb, JGR 2001 71 71 !!------------------------------------------------------------------ 72 INTEGER , INTENT (in) :: kt ! Ocean time step 72 INTEGER , INTENT (in) :: kt ! Ocean time step 73 73 ! 74 74 INTEGER :: ji, jj, jl, jcat ! dummy loop index … … 76 76 REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars 77 77 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - 78 REAL(wp) :: zx3 78 REAL(wp) :: zx3 79 79 REAL(wp) :: zslope ! used to compute local thermodynamic "speeds" 80 80 ! … … 90 90 IF( ln_timing ) CALL timing_start('iceitd_rem') 91 91 92 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution' 92 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution' 93 93 94 94 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) … … 107 107 ENDIF 108 108 END_2D 109 109 110 110 !----------------------------------------------------------------------------------------------- 111 111 ! 2) Compute new category boundaries … … 143 143 ELSEIF( a_ib_2d(ji,jl) <= epsi10 .AND. a_ib_2d(ji,jl+1) > epsi10 ) THEN ! a(jl)=0 => Hn* = Hn + fn+1*dt 144 144 zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl+1) 145 ELSE ! a(jl+1) & a(jl) = 0 145 ELSE ! a(jl+1) & a(jl) = 0 146 146 zhbnew(ji,jl) = hi_max(jl) 147 147 ENDIF 148 148 ! 149 149 ! --- 2 conditions for remapping --- ! 150 ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi 151 ! Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 150 ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi 151 ! Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 152 152 ! in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 153 153 # if defined key_single … … 159 159 # endif 160 160 ! 161 ! 2) Hn-1 < Hn* < Hn+1 161 ! 2) Hn-1 < Hn* < Hn+1 162 162 IF( zhbnew(ji,jl) < hi_max(jl-1) ) nptidx(ji) = 0 163 163 IF( zhbnew(ji,jl) > hi_max(jl+1) ) nptidx(ji) = 0 … … 171 171 zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * h_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) ) 172 172 ELSE 173 zhbnew(ji,jpl) = hi_max(jpl) 173 zhbnew(ji,jpl) = hi_max(jpl) 174 174 ENDIF 175 175 ! 176 176 ! --- 1 additional condition for remapping (1st category) --- ! 177 ! H0+epsi < h1(t) < H1-epsi 178 ! h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 177 ! H0+epsi < h1(t) < H1-epsi 178 ! h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 179 179 ! in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 180 180 # if defined key_single … … 202 202 ! 203 203 ENDIF 204 204 205 205 !----------------------------------------------------------------------------------------------- 206 ! 4) Compute g(h) 206 ! 4) Compute g(h) 207 207 !----------------------------------------------------------------------------------------------- 208 208 IF( npti > 0 ) THEN 209 209 ! 210 210 zhb0(:) = hi_max(0) ; zhb1(:) = hi_max(1) 211 g0(:,:) = 0._wp ; g1(:,:) = 0._wp 212 hL(:,:) = 0._wp ; hR(:,:) = 0._wp 211 g0(:,:) = 0._wp ; g1(:,:) = 0._wp 212 hL(:,:) = 0._wp ; hR(:,:) = 0._wp 213 213 ! 214 214 DO jl = 1, jpl … … 220 220 ! 221 221 IF( jl == 1 ) THEN 222 ! 222 ! 223 223 ! --- g(h) for category 1 --- ! 224 224 CALL itd_glinear( zhb0(1:npti) , zhb1(1:npti) , h_ib_1d(1:npti) , a_i_1d(1:npti) , & ! in … … 230 230 IF( a_i_1d(ji) > epsi10 ) THEN 231 231 ! 232 zdh0 = h_i_1d(ji) - h_ib_1d(ji) 232 zdh0 = h_i_1d(ji) - h_ib_1d(ji) 233 233 IF( zdh0 < 0.0 ) THEN ! remove area from category 1 234 234 zdh0 = MIN( -zdh0, hi_max(1) ) … … 238 238 IF( zetamax > 0.0 ) THEN 239 239 zx1 = zetamax 240 zx2 = 0.5 * zetamax * zetamax 240 zx2 = 0.5 * zetamax * zetamax 241 241 zda0 = g1(ji,1) * zx2 + g0(ji,1) * zx1 ! ice area removed 242 zdamax = a_i_1d(ji) * (1.0 - h_i_1d(ji) / h_ib_1d(ji) ) ! Constrain new thickness <= h_i 242 zdamax = a_i_1d(ji) * (1.0 - h_i_1d(ji) / h_ib_1d(ji) ) ! Constrain new thickness <= h_i 243 243 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting of thin ice (zdamax > 0) 244 244 ! Remove area, conserving volume … … 250 250 ELSE ! if ice accretion zdh0 > 0 251 251 ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 252 zhbnew(ji,0) = MIN( zdh0, hi_max(1) ) 252 zhbnew(ji,0) = MIN( zdh0, hi_max(1) ) 253 253 ENDIF 254 254 ! … … 263 263 ENDIF ! jl=1 264 264 ! 265 ! --- g(h) for each thickness category --- ! 265 ! --- g(h) for each thickness category --- ! 266 266 CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti) , a_i_1d(1:npti) , & ! in 267 267 & g0 (1:npti,jl ), g1 (1:npti,jl), hL (1:npti,jl), hR (1:npti,jl) ) ! out 268 268 ! 269 269 END DO 270 270 271 271 !----------------------------------------------------------------------------------------------- 272 272 ! 5) Compute area and volume to be shifted across each boundary (Eq. 18) … … 278 278 ! left and right integration limits in eta space 279 279 IF (zhbnew(ji,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1 280 zetamin = MAX( hi_max(jl) , hL(ji,jl) ) - hL(ji,jl) ! hi_max(jl) - hL 280 zetamin = MAX( hi_max(jl) , hL(ji,jl) ) - hL(ji,jl) ! hi_max(jl) - hL 281 281 zetamax = MIN( zhbnew(ji,jl), hR(ji,jl) ) - hL(ji,jl) ! hR - hL 282 282 jdonor(ji,jl) = jl … … 301 301 END DO 302 302 END DO 303 303 304 304 !---------------------------------------------------------------------------------------------- 305 305 ! 6) Shift ice between categories 306 306 !---------------------------------------------------------------------------------------------- 307 307 CALL itd_shiftice ( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) 308 308 309 309 !---------------------------------------------------------------------------------------------- 310 310 ! 7) Make sure h_i >= minimum ice thickness hi_min … … 316 316 DO ji = 1, npti 317 317 IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN 318 a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 318 a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 319 319 IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 320 320 h_i_1d(ji) = rn_himin … … 384 384 pg1(ji) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5_wp ) ! Eq. 14 385 385 ! 386 ELSE ! remap_flag = .false. or a_i < epsi10 386 ELSE ! remap_flag = .false. or a_i < epsi10 387 387 phL(ji) = 0._wp 388 388 phR(ji) = 0._wp … … 415 415 REAL(wp), DIMENSION(jpij,nlay_s,jpl) :: ze_s_2d 416 416 !!------------------------------------------------------------------ 417 417 418 418 CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i ) 419 419 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) … … 445 445 END DO 446 446 END DO 447 447 448 448 !------------------------------------------------------------------------------- 449 449 ! 2) Transfer volume and energy between categories … … 457 457 ! 458 458 IF ( jl1 == jl ) THEN ; jl2 = jl1+1 459 ELSE ; jl2 = jl 459 ELSE ; jl2 = jl 460 460 ENDIF 461 461 ! … … 475 475 ztrans = v_s_2d(ji,jl1) * zworkv(ji) ! Snow volumes 476 476 v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans 477 v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 477 v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 478 478 ! 479 479 ztrans = oa_i_2d(ji,jl1) * zworka(ji) ! Ice age … … 488 488 zaTsfn(ji,jl1) = zaTsfn(ji,jl1) - ztrans 489 489 zaTsfn(ji,jl2) = zaTsfn(ji,jl2) + ztrans 490 ! 490 ! 491 491 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 492 492 ztrans = a_ip_2d(ji,jl1) * zworka(ji) ! Pond fraction 493 493 a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans 494 494 a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans 495 ! 495 ! 496 496 ztrans = v_ip_2d(ji,jl1) * zworkv(ji) ! Pond volume 497 497 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans … … 555 555 & a_i_2d(1:npti,jl) = a_i_2d(1:npti,jl) * rn_amax_1d(1:npti) / zworka(1:npti) 556 556 END DO 557 557 558 558 !------------------------------------------------------------------------------- 559 559 ! 4) Update ice thickness and temperature … … 564 564 WHERE( a_i_2d(1:npti,:) >= epsi20 ) 565 565 # endif 566 h_i_2d (1:npti,:) = v_i_2d(1:npti,:) / a_i_2d(1:npti,:) 567 t_su_2d(1:npti,:) = zaTsfn(1:npti,:) / a_i_2d(1:npti,:) 566 h_i_2d (1:npti,:) = v_i_2d(1:npti,:) / a_i_2d(1:npti,:) 567 t_su_2d(1:npti,:) = zaTsfn(1:npti,:) / a_i_2d(1:npti,:) 568 568 ELSEWHERE 569 569 h_i_2d (1:npti,:) = 0._wp … … 591 591 ! 592 592 END SUBROUTINE itd_shiftice 593 593 594 594 595 595 SUBROUTINE ice_itd_reb( kt ) … … 603 603 !! to the neighboring category 604 604 !!------------------------------------------------------------------ 605 INTEGER , INTENT (in) :: kt ! Ocean time step 605 INTEGER , INTENT (in) :: kt ! Ocean time step 606 606 INTEGER :: ji, jj, jl ! dummy loop indices 607 607 ! … … 611 611 IF( ln_timing ) CALL timing_start('iceitd_reb') 612 612 ! 613 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' 613 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' 614 614 ! 615 615 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) … … 627 627 IF( a_i(ji,jj,jl) > 0._wp .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN 628 628 npti = npti + 1 629 nptidx( npti ) = (jj - 1) * jpi + ji 629 nptidx( npti ) = (jj - 1) * jpi + ji 630 630 ENDIF 631 631 END_2D 632 632 ! 633 IF( npti > 0 ) THEN 633 IF( npti > 0 ) THEN 634 634 !!clem CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 635 635 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) … … 637 637 ! 638 638 DO ji = 1, npti 639 jdonor(ji,jl) = jl 639 jdonor(ji,jl) = jl 640 640 ! how much of a_i you send in cat sup is somewhat arbitrary 641 641 ! these are from CICE => transfer everything … … 663 663 IF( a_i(ji,jj,jl+1) > 0._wp .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN 664 664 npti = npti + 1 665 nptidx( npti ) = (jj - 1) * jpi + ji 665 nptidx( npti ) = (jj - 1) * jpi + ji 666 666 ENDIF 667 667 END_2D … … 672 672 DO ji = 1, npti 673 673 jdonor(ji,jl) = jl + 1 674 zdaice(ji,jl) = a_i_1d(ji) 674 zdaice(ji,jl) = a_i_1d(ji) 675 675 zdvice(ji,jl) = v_i_1d(ji) 676 676 END DO … … 721 721 WRITE(numout,*) ' mean ice thickness in the domain rn_himean = ', rn_himean 722 722 WRITE(numout,*) ' Ice categories are defined by rn_catbnd ln_cat_usr = ', ln_cat_usr 723 WRITE(numout,*) ' minimum ice thickness allowed rn_himin = ', rn_himin 724 WRITE(numout,*) ' maximum ice thickness allowed rn_himax = ', rn_himax 723 WRITE(numout,*) ' minimum ice thickness allowed rn_himin = ', rn_himin 724 WRITE(numout,*) ' maximum ice thickness allowed rn_himax = ', rn_himax 725 725 ENDIF 726 726 ! … … 729 729 !-----------------------------------! 730 730 ! !== set the choice of ice categories ==! 731 ioptio = 0 731 ioptio = 0 732 732 IF( ln_cat_hfn ) THEN ; ioptio = ioptio + 1 ; nice_catbnd = np_cathfn ; ENDIF 733 733 IF( ln_cat_usr ) THEN ; ioptio = ioptio + 1 ; nice_catbnd = np_catusr ; ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
