New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 14021 for NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/ICE/iceitd.F90 – NEMO

Ignore:
Timestamp:
2020-12-02T20:53:00+01:00 (3 years ago)
Author:
laurent
Message:

Caught up with trunk rev 14020...

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/ICE/iceitd.F90

    r13618 r14021  
    1818   !!---------------------------------------------------------------------- 
    1919   USE dom_oce        ! ocean domain 
    20    USE phycst         ! physical constants  
     20   USE phycst         ! physical constants 
    2121   USE ice1D          ! sea-ice: thermodynamic variables 
    2222   USE ice            ! sea-ice: variables 
     
    2929   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
    3030   USE prtctl         ! Print control 
     31   USE timing         ! Timing 
    3132 
    3233   IMPLICIT NONE 
     
    6566      !!                after thermodynamic growth of ice thickness 
    6667      !! 
    67       !! ** Method  :   Linear remapping  
     68      !! ** Method  :   Linear remapping 
    6869      !! 
    6970      !! References :   W.H. Lipscomb, JGR 2001 
    7071      !!------------------------------------------------------------------ 
    71       INTEGER , INTENT (in) ::   kt      ! Ocean time step  
     72      INTEGER , INTENT (in) ::   kt      ! Ocean time step 
    7273      ! 
    7374      INTEGER  ::   ji, jj, jl, jcat     ! dummy loop index 
     
    7576      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars 
    7677      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      - 
    77       REAL(wp) ::   zx3         
     78      REAL(wp) ::   zx3 
    7879      REAL(wp) ::   zslope          ! used to compute local thermodynamic "speeds" 
    7980      ! 
     
    8788      REAL(wp), DIMENSION(jpij,0:jpl) ::   zhbnew          ! new boundaries of ice categories 
    8889      !!------------------------------------------------------------------ 
    89  
    90       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution'  
     90      IF( ln_timing )   CALL timing_start('iceitd_rem') 
     91 
     92      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution' 
    9193 
    9294      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     
    105107         ENDIF 
    106108      END_2D 
    107        
     109 
    108110      !----------------------------------------------------------------------------------------------- 
    109111      !  2) Compute new category boundaries 
     
    141143               ELSEIF( a_ib_2d(ji,jl) <= epsi10 .AND. a_ib_2d(ji,jl+1) >  epsi10 ) THEN   ! a(jl)=0 => Hn* = Hn + fn+1*dt 
    142144                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl+1) 
    143                ELSE                                                                       ! a(jl+1) & a(jl) = 0  
     145               ELSE                                                                       ! a(jl+1) & a(jl) = 0 
    144146                  zhbnew(ji,jl) = hi_max(jl) 
    145147               ENDIF 
    146148               ! 
    147149               ! --- 2 conditions for remapping --- ! 
    148                ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi                
    149                !    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 
    150152               !          in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
    151153# if defined key_single 
     
    157159# endif 
    158160               ! 
    159                ! 2) Hn-1 < Hn* < Hn+1   
     161               ! 2) Hn-1 < Hn* < Hn+1 
    160162               IF( zhbnew(ji,jl) < hi_max(jl-1) )   nptidx(ji) = 0 
    161163               IF( zhbnew(ji,jl) > hi_max(jl+1) )   nptidx(ji) = 0 
     
    169171               zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * h_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) ) 
    170172            ELSE 
    171                zhbnew(ji,jpl) = hi_max(jpl)   
     173               zhbnew(ji,jpl) = hi_max(jpl) 
    172174            ENDIF 
    173175            ! 
    174176            ! --- 1 additional condition for remapping (1st category) --- ! 
    175             ! H0+epsi < h1(t) < H1-epsi  
    176             !    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 
    177179            !    in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
    178180# if defined key_single 
     
    200202         ! 
    201203      ENDIF 
    202     
     204 
    203205      !----------------------------------------------------------------------------------------------- 
    204       !  4) Compute g(h)  
     206      !  4) Compute g(h) 
    205207      !----------------------------------------------------------------------------------------------- 
    206208      IF( npti > 0 ) THEN 
    207209         ! 
    208210         zhb0(:) = hi_max(0)   ;   zhb1(:) = hi_max(1) 
    209          g0(:,:) = 0._wp       ;   g1(:,:) = 0._wp  
    210          hL(:,:) = 0._wp       ;   hR(:,:) = 0._wp  
     211         g0(:,:) = 0._wp       ;   g1(:,:) = 0._wp 
     212         hL(:,:) = 0._wp       ;   hR(:,:) = 0._wp 
    211213         ! 
    212214         DO jl = 1, jpl 
     
    218220            ! 
    219221            IF( jl == 1 ) THEN 
    220                !   
     222               ! 
    221223               ! --- g(h) for category 1 --- ! 
    222224               CALL itd_glinear( zhb0(1:npti)  , zhb1(1:npti)  , h_ib_1d(1:npti)  , a_i_1d(1:npti)  ,  &   ! in 
     
    228230                  IF( a_i_1d(ji) > epsi10 ) THEN 
    229231                     ! 
    230                      zdh0 =  h_i_1d(ji) - h_ib_1d(ji)                 
     232                     zdh0 =  h_i_1d(ji) - h_ib_1d(ji) 
    231233                     IF( zdh0 < 0.0 ) THEN      ! remove area from category 1 
    232234                        zdh0 = MIN( -zdh0, hi_max(1) ) 
     
    236238                        IF( zetamax > 0.0 ) THEN 
    237239                           zx1    = zetamax 
    238                            zx2    = 0.5 * zetamax * zetamax  
     240                           zx2    = 0.5 * zetamax * zetamax 
    239241                           zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                ! ice area removed 
    240                            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 
    241243                           zda0   = MIN( zda0, zdamax )                            ! ice area lost due to melting of thin ice (zdamax > 0) 
    242244                           ! Remove area, conserving volume 
     
    248250                     ELSE ! if ice accretion zdh0 > 0 
    249251                        ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 
    250                         zhbnew(ji,0) = MIN( zdh0, hi_max(1) )  
     252                        zhbnew(ji,0) = MIN( zdh0, hi_max(1) ) 
    251253                     ENDIF 
    252254                     ! 
     
    261263            ENDIF ! jl=1 
    262264            ! 
    263             ! --- g(h) for each thickness category --- !   
     265            ! --- g(h) for each thickness category --- ! 
    264266            CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti)   , a_i_1d(1:npti)   ,  &   ! in 
    265267               &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL    (1:npti,jl), hR    (1:npti,jl)   )   ! out 
    266268            ! 
    267269         END DO 
    268           
     270 
    269271         !----------------------------------------------------------------------------------------------- 
    270272         !  5) Compute area and volume to be shifted across each boundary (Eq. 18) 
     
    276278               ! left and right integration limits in eta space 
    277279               IF (zhbnew(ji,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1 
    278                   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 
    279281                  zetamax = MIN( zhbnew(ji,jl), hR(ji,jl) ) - hL(ji,jl)   ! hR - hL 
    280282                  jdonor(ji,jl) = jl 
     
    299301            END DO 
    300302         END DO 
    301           
     303 
    302304         !---------------------------------------------------------------------------------------------- 
    303305         ! 6) Shift ice between categories 
    304306         !---------------------------------------------------------------------------------------------- 
    305307         CALL itd_shiftice ( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) 
    306           
     308 
    307309         !---------------------------------------------------------------------------------------------- 
    308310         ! 7) Make sure h_i >= minimum ice thickness hi_min 
     
    314316         DO ji = 1, npti 
    315317            IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN 
    316                a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin  
    317                IF( ln_pnd_LEV )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 
     318               a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 
     319               IF( ln_pnd_LEV .OR. ln_pnd_TOPO )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 
    318320               h_i_1d(ji) = rn_himin 
    319321            ENDIF 
     
    328330      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    329331      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_rem',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
     332      IF( ln_timing    )   CALL timing_stop ('iceitd_rem') 
    330333      ! 
    331334   END SUBROUTINE ice_itd_rem 
     
    381384            pg1(ji) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5_wp )   ! Eq. 14 
    382385            ! 
    383          ELSE  ! remap_flag = .false. or a_i < epsi10  
     386         ELSE  ! remap_flag = .false. or a_i < epsi10 
    384387            phL(ji) = 0._wp 
    385388            phR(ji) = 0._wp 
     
    412415      REAL(wp), DIMENSION(jpij,nlay_s,jpl) ::   ze_s_2d 
    413416      !!------------------------------------------------------------------ 
    414           
     417 
    415418      CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i  ) 
    416419      CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i  ) 
     
    442445         END DO 
    443446      END DO 
    444        
     447 
    445448      !------------------------------------------------------------------------------- 
    446449      ! 2) Transfer volume and energy between categories 
     
    454457               ! 
    455458               IF ( jl1 == jl  ) THEN   ;   jl2 = jl1+1 
    456                ELSE                     ;   jl2 = jl  
     459               ELSE                     ;   jl2 = jl 
    457460               ENDIF 
    458461               ! 
     
    472475               ztrans         = v_s_2d(ji,jl1) * zworkv(ji)          ! Snow volumes 
    473476               v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans 
    474                v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans  
     477               v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 
    475478               ! 
    476479               ztrans          = oa_i_2d(ji,jl1) * zworka(ji)        ! Ice age 
     
    485488               zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans 
    486489               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans 
    487                !   
    488                IF ( ln_pnd_LEV ) THEN 
     490               ! 
     491               IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    489492                  ztrans          = a_ip_2d(ji,jl1) * zworka(ji)     ! Pond fraction 
    490493                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans 
    491494                  a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans 
    492                   !                                               
    493                   ztrans          = v_ip_2d(ji,jl1) * zworka(ji)     ! Pond volume (also proportional to da/a) 
     495                  ! 
     496                  ztrans          = v_ip_2d(ji,jl1) * zworkv(ji)     ! Pond volume 
    494497                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 
    495498                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 
    496499                  ! 
    497500                  IF ( ln_pnd_lids ) THEN                            ! Pond lid volume 
    498                      ztrans          = v_il_2d(ji,jl1) * zworka(ji) 
     501                     ztrans          = v_il_2d(ji,jl1) * zworkv(ji) 
    499502                     v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - ztrans 
    500503                     v_il_2d(ji,jl2) = v_il_2d(ji,jl2) + ztrans 
     
    552555            &   a_i_2d(1:npti,jl) = a_i_2d(1:npti,jl) * rn_amax_1d(1:npti) / zworka(1:npti) 
    553556      END DO 
    554        
     557 
    555558      !------------------------------------------------------------------------------- 
    556559      ! 4) Update ice thickness and temperature 
     
    561564      WHERE( a_i_2d(1:npti,:) >= epsi20 ) 
    562565# endif 
    563          h_i_2d (1:npti,:)  =  v_i_2d(1:npti,:) / a_i_2d(1:npti,:)  
    564          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,:) 
    565568      ELSEWHERE 
    566569         h_i_2d (1:npti,:)  = 0._wp 
     
    588591      ! 
    589592   END SUBROUTINE itd_shiftice 
    590     
     593 
    591594 
    592595   SUBROUTINE ice_itd_reb( kt ) 
     
    600603      !!              to the neighboring category 
    601604      !!------------------------------------------------------------------ 
    602       INTEGER , INTENT (in) ::   kt      ! Ocean time step  
     605      INTEGER , INTENT (in) ::   kt      ! Ocean time step 
    603606      INTEGER ::   ji, jj, jl   ! dummy loop indices 
    604607      ! 
     
    606609      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred 
    607610      !!------------------------------------------------------------------ 
    608       ! 
    609       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution'  
     611      IF( ln_timing )   CALL timing_start('iceitd_reb') 
     612      ! 
     613      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' 
    610614      ! 
    611615      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     
    623627            IF( a_i(ji,jj,jl) > 0._wp .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN 
    624628               npti = npti + 1 
    625                nptidx( npti ) = (jj - 1) * jpi + ji                   
     629               nptidx( npti ) = (jj - 1) * jpi + ji 
    626630            ENDIF 
    627631         END_2D 
    628632         ! 
    629          IF( npti > 0 ) THEN             
     633         IF( npti > 0 ) THEN 
    630634            !!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 
    631635            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 
     
    633637            ! 
    634638            DO ji = 1, npti 
    635                jdonor(ji,jl)  = jl  
     639               jdonor(ji,jl)  = jl 
    636640               ! how much of a_i you send in cat sup is somewhat arbitrary 
    637                !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    638                !!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji)   
    639                !!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 
    640                !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    641                !!          zdaice(ji,jl)  = a_i_1d(ji) 
    642                !!          zdvice(ji,jl)  = v_i_1d(ji) 
    643                !!clem: these are from UCL and work ok 
    644                zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp 
    645                zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
     641               ! these are from CICE => transfer everything 
     642               !!zdaice(ji,jl)  = a_i_1d(ji) 
     643               !!zdvice(ji,jl)  = v_i_1d(ji) 
     644               ! these are from LLN => transfer only half of the category 
     645               zdaice(ji,jl)  =                       0.5_wp  * a_i_1d(ji) 
     646               zdvice(ji,jl)  = v_i_1d(ji) - (1._wp - 0.5_wp) * a_i_1d(ji) * hi_mean(jl) 
    646647            END DO 
    647648            ! 
     
    662663            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 
    663664               npti = npti + 1 
    664                nptidx( npti ) = (jj - 1) * jpi + ji                   
     665               nptidx( npti ) = (jj - 1) * jpi + ji 
    665666            ENDIF 
    666667         END_2D 
     
    671672            DO ji = 1, npti 
    672673               jdonor(ji,jl) = jl + 1 
    673                zdaice(ji,jl) = a_i_1d(ji)  
     674               zdaice(ji,jl) = a_i_1d(ji) 
    674675               zdvice(ji,jl) = v_i_1d(ji) 
    675676            END DO 
     
    686687      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    687688      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
     689      IF( ln_timing    )   CALL timing_stop ('iceitd_reb') 
    688690      ! 
    689691   END SUBROUTINE ice_itd_reb 
     
    719721         WRITE(numout,*) '         mean ice thickness in the domain                               rn_himean  = ', rn_himean 
    720722         WRITE(numout,*) '      Ice categories are defined by rn_catbnd                           ln_cat_usr = ', ln_cat_usr 
    721          WRITE(numout,*) '      minimum ice thickness allowed                                     rn_himin   = ', rn_himin  
    722          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 
    723725      ENDIF 
    724726      ! 
     
    727729      !-----------------------------------! 
    728730      !                             !== set the choice of ice categories ==! 
    729       ioptio = 0  
     731      ioptio = 0 
    730732      IF( ln_cat_hfn ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_cathfn    ;   ENDIF 
    731733      IF( ln_cat_usr ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_catusr    ;   ENDIF 
Note: See TracChangeset for help on using the changeset viewer.