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 8361 for branches/2017 – NEMO

Changeset 8361 for branches/2017


Ignore:
Timestamp:
2017-07-23T12:51:10+02:00 (7 years ago)
Author:
clem
Message:

first step toward full 1D for limitd_th

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r8355 r8361  
    1515   !!   lim_itd_th_rem   : 
    1616   !!   lim_itd_th_reb   : 
    17    !!   lim_itd_fitline  : 
     17   !!   lim_itd_glinear  : 
    1818   !!   lim_itd_shiftice : 
    1919   !!---------------------------------------------------------------------- 
     
    2424   USE ice              ! LIM-3 variables 
    2525   USE limvar           ! LIM-3 variables 
     26   USE limcons          ! conservation tests 
     27   ! 
    2628   USE prtctl           ! Print control 
    2729   USE in_out_manager   ! I/O manager 
     
    2931   USE wrk_nemo         ! work arrays 
    3032   USE lib_fortran      ! to use key_nosignedzero 
    31    USE limcons          ! conservation tests 
    3233 
    3334   IMPLICIT NONE 
     
    7172      REAL(wp), POINTER, DIMENSION(:,:,:) ::   hL          ! left boundary for the ITD for each thickness 
    7273      REAL(wp), POINTER, DIMENSION(:,:,:) ::   hR          ! left boundary for the ITD for each thickness 
    73       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zht_i_b     ! old ice thickness 
    7474      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice          ! local increment of ice area and volume 
    75       INTEGER , POINTER, DIMENSION(:)     ::   nind_i, nind_j          ! compressed indices for i/j directions 
    76       INTEGER                             ::   nbrem                   ! number of cells with ice to transfer 
     75      INTEGER , DIMENSION(jpij)  ::   nind_i, nind_j          ! compressed indices for i/j directions 
    7776      REAL(wp)                            ::   zslope                  ! used to compute local thermodynamic "speeds" 
    7877      REAL(wp), POINTER, DIMENSION(:,:)   ::   zhb0, zhb1              ! category boundaries for thinnes categories 
     
    8382      CALL wrk_alloc( jpi,jpj, zremap_flag ) 
    8483      CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) 
    85       CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b ) 
     84      CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR ) 
    8685      CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    8786      CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
    88       CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j )  
    8987      CALL wrk_alloc( jpi,jpj, zhb0, zhb1 ) 
    9088 
    91       !!---------------------------------------------------------------------------------------------- 
    92       !! 1) Compute thickness and changes in each ice category 
    93       !!---------------------------------------------------------------------------------------------- 
    9489      IF( kt == nit000 .AND. lwp) THEN 
    9590         WRITE(numout,*) 
     
    9893      ENDIF 
    9994 
    100       zdhice(:,:,:) = 0._wp 
    101       DO jl = 1, jpl 
    102          DO jj = 1, jpj 
    103             DO ji = 1, jpi 
    104                rswitch           = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) )     !0 if no ice and 1 if yes 
    105                ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 
    106                rswitch           = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) ) 
    107                zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 
    108                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?  
    109             END DO 
    110          END DO 
    111       END DO 
    112  
    113 !!      !----------------------------------------------------------------------------------------------- 
    114 !!      !  2) Compute fractional ice area in each grid cell 
    115 !!      !----------------------------------------------------------------------------------------------- 
    116 !!      at_i(:,:) = 0._wp 
    117 !!      DO jl = 1, jpl 
    118 !!         at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    119 !!      END DO 
    120  
    12195      !----------------------------------------------------------------------------------------------- 
    12296      !  3) Identify grid cells with ice 
    12397      !----------------------------------------------------------------------------------------------- 
    124       nbrem = 0 
     98      nidx = 0 ; idxice(:) = 0 
    12599      DO jj = 1, jpj 
    126100         DO ji = 1, jpi 
    127101            IF ( at_i(ji,jj) > epsi10 ) THEN 
    128                nbrem         = nbrem + 1 
    129                nind_i(nbrem) = ji 
    130                nind_j(nbrem) = jj 
     102               nidx         = nidx + 1 
     103               nind_i(nidx) = ji 
     104               nind_j(nidx) = jj 
     105               idxice( nidx ) = (jj - 1) * jpi + ji 
    131106               zremap_flag(ji,jj) = 1 
    132107            ELSE 
     
    139114      !  4) Compute new category boundaries: 1:jpl-1 
    140115      !----------------------------------------------------------------------------------------------- 
     116      zdhice(:,:,:) = 0._wp 
    141117      zhbnew(:,:,:) = 0._wp 
    142  
    143       IF( nbrem > 0 ) THEN 
     118      zhb0(:,:) = hi_max(0) 
     119      zhb1(:,:) = hi_max(1) 
     120 
     121      IF( nidx > 0 ) THEN 
     122 
     123         ! Compute thickness change in each ice category 
     124         DO jl = 1, jpl 
     125            DO ji = 1, nidx 
     126               ii = nind_i(ji) 
     127               ij = nind_j(ji) 
     128               zdhice(ii,ij,jl) = ht_i(ii,ij,jl) - ht_i_b(ii,ij,jl) 
     129            END DO 
     130         END DO 
     131          
    144132         DO jl = 1, jpl - 1 
    145             DO ji = 1, nbrem 
     133            DO ji = 1, nidx 
    146134               ii = nind_i(ji) 
    147135               ij = nind_j(ji) 
    148136               ! 
    149                ! --- New boundary categories Hn* = Hn + Fn*dt --- 
    150                !!    Fn*dt = ( fn + (fn+1 - fn)/(hn+1 - hn) * (Hn - hn) ) * dt = zdhice + zslope * (Hmax - ht_i_b) 
    151                IF    ( a_i_b(ii,ij,jl) >  epsi10 .AND. a_i_b(ii,ij,jl+1) >  epsi10 ) THEN !! a_i_b(jl+1) & a_i_b(jl) /= 0 
    152                   zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 
    153                   zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 
    154 !!               ELSEIF( a_i_b(ii,ij,jl) >  epsi10 .AND. a_i_b(ii,jj,jl+1) <= epsi10 ) THEN !! a_i_b(jl+1)=0 => Hn* = Hn + fn*dt 
    155                ELSEIF( a_i_b(ii,ij,jl) >  epsi10 ) THEN !! a_i_b(jl+1)=0 => Hn* = Hn + fn*dt 
     137               ! --- New boundary categories Hn* = Hn + Fn*dt --- ! 
     138               !     Fn*dt = ( fn + (fn+1 - fn)/(hn+1 - hn) * (Hn - hn) ) * dt = zdhice + zslope * (Hmax - ht_i_b) 
     139               ! 
     140               IF    ( a_i_b(ii,ij,jl) >  epsi10 .AND. a_i_b(ii,ij,jl+1) >  epsi10 ) THEN   ! a_i_b(jl+1) & a_i_b(jl) /= 0 
     141                  zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( ht_i_b(ii,ij,jl+1) - ht_i_b(ii,ij,jl) ) 
     142                  zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - ht_i_b(ii,ij,jl) ) 
     143               ELSEIF( a_i_b(ii,ij,jl) >  epsi10 .AND. a_i_b(ii,ij,jl+1) <= epsi10 ) THEN   ! a_i_b(jl+1)=0 => Hn* = Hn + fn*dt 
    156144                  zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 
    157 !!               ELSEIF( a_i_b(ii,ij,jl) <= epsi10 .AND. a_i_b(ii,ij,jl+1) >  epsi10 ) THEN !! a_i_b(jl)=0 => Hn* = Hn + fn+1*dt 
    158                ELSEIF( a_i_b(ii,ij,jl+1) > epsi10 ) THEN !! a_i_b(jl)=0 => Hn* = Hn + fn+1*dt 
     145               ELSEIF( a_i_b(ii,ij,jl) <= epsi10 .AND. a_i_b(ii,ij,jl+1) >  epsi10 ) THEN   ! a_i_b(jl)=0 => Hn* = Hn + fn+1*dt 
    159146                  zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 
    160                ELSE                                                                       !! a_i_b(jl+1) & a_i_b(jl) = 0  
     147               ELSE                                                                         ! a_i_b(jl+1) & a_i_b(jl) = 0  
    161148                  zhbnew(ii,ij,jl) = hi_max(jl) 
    162149               ENDIF 
     
    165152               ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi                
    166153               !    Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible  
    167                !          in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
     154               !          in lim_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
    168155               IF( a_i(ii,ij,jl  ) > epsi10 .AND. ht_i(ii,ij,jl  ) > ( zhbnew(ii,ij,jl) - epsi10 ) )   zremap_flag(ii,ij) = 0 
    169156               IF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) )   zremap_flag(ii,ij) = 0 
     
    175162            END DO 
    176163         END DO 
     164          
     165         ! --- New boundary for category jpl --- ! 
     166         DO ji = 1, nidx 
     167            ii = nind_i(ji) 
     168            ij = nind_j(ji) 
     169            IF( a_i(ii,ij,jpl) > epsi10 ) THEN 
     170               zhbnew(ii,ij,jpl) = MAX( hi_max(jpl-1), 3._wp * ht_i(ii,ij,jpl) - 2._wp * zhbnew(ii,ij,jpl-1) ) 
     171            ELSE 
     172               zhbnew(ii,ij,jpl) = hi_max(jpl)   
     173            ENDIF 
     174             
     175            ! 1 condition for remapping for the 1st category 
     176            ! H0+epsi < h1(t) < H1-epsi  
     177            !    h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible  
     178            !    in lim_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
     179            IF( ht_i_b(ii,ij,1) < ( hi_max(0) + epsi10 ) )   zremap_flag(ii,ij) = 0 
     180            IF( ht_i_b(ii,ij,1) > ( hi_max(1) - epsi10 ) )   zremap_flag(ii,ij) = 0 
     181         END DO 
     182          
    177183      ENDIF 
    178184    
     
    180186      !  5) Identify cells where ITD is to be remapped 
    181187      !----------------------------------------------------------------------------------------------- 
    182       nbrem = 0 
     188      nidx = 0 
    183189      DO jj = 1, jpj 
    184190         DO ji = 1, jpi 
    185191            IF( zremap_flag(ji,jj) == 1 ) THEN 
    186                nbrem         = nbrem + 1 
    187                nind_i(nbrem) = ji 
    188                nind_j(nbrem) = jj 
     192               nidx         = nidx + 1 
     193               nind_i(nidx) = ji 
     194               nind_j(nidx) = jj 
    189195            ENDIF 
    190196         END DO  
     
    192198 
    193199      !----------------------------------------------------------------------------------------------- 
    194       !  6) Compute new category boundaries: jpl 
    195       !----------------------------------------------------------------------------------------------- 
    196       DO jj = 1, jpj 
    197          DO ji = 1, jpi 
    198             zhb0(ji,jj) = hi_max(0) 
    199             zhb1(ji,jj) = hi_max(1) 
    200  
    201             ! boundary for jpl 
    202             IF( a_i(ji,jj,jpl) > epsi10 ) THEN 
    203                zhbnew(ji,jj,jpl) = MAX( hi_max(jpl-1), 3._wp * ht_i(ji,jj,jpl) - 2._wp * zhbnew(ji,jj,jpl-1) ) 
    204             ELSE 
    205 !clem bug               zhbnew(ji,jj,jpl) = hi_max(jpl)   
    206                zhbnew(ji,jj,jpl) = hi_max(jpl-1) ! not used anyway 
     200      !  7) Compute g(h)  
     201      !----------------------------------------------------------------------------------------------- 
     202      IF( nidx > 0 ) THEN 
     203 
     204         !- 7.1 g(h) for category 1 at start of time step 
     205         CALL lim_itd_glinear( 1, zhb0, zhb1, ht_i_b(:,:,1), g0(:,:,1), g1(:,:,1), hL(:,:,1),   & 
     206            &                  hR(:,:,1), zremap_flag ) 
     207          
     208         !- 7.2 Area lost due to melting of thin ice (first category,  1) 
     209         DO ji = 1, nidx 
     210            ii = nind_i(ji)  
     211            ij = nind_j(ji)  
     212             
     213            IF( a_i(ii,ij,1) > epsi10 ) THEN 
     214                
     215               zdh0 = zdhice(ii,ij,1) !decrease of ice thickness in the lower category 
     216                
     217               IF( zdh0 < 0.0 ) THEN      !remove area from category 1 
     218                  zdh0 = MIN( -zdh0, hi_max(1) ) 
     219                  !Integrate g(1) from 0 to dh0 to estimate area melted 
     220                  zetamax = MIN( zdh0, hR(ii,ij,1) ) - hL(ii,ij,1) 
     221                   
     222                  IF( zetamax > 0.0 ) THEN 
     223                     zx1    = zetamax 
     224                     zx2    = 0.5 * zetamax * zetamax  
     225                     zda0   = g1(ii,ij,1) * zx2 + g0(ii,ij,1) * zx1                        ! ice area removed 
     226                     zdamax = a_i(ii,ij,1) * (1.0 - ht_i(ii,ij,1) / ht_i_b(ii,ij,1) ) ! Constrain new thickness <= ht_i                 
     227                     zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting  
     228                     !     of thin ice (zdamax > 0) 
     229                     ! Remove area, conserving volume 
     230                     ht_i(ii,ij,1) = ht_i(ii,ij,1) * a_i(ii,ij,1) / ( a_i(ii,ij,1) - zda0 ) 
     231                     a_i(ii,ij,1)  = a_i(ii,ij,1) - zda0 
     232                     v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) ! clem-useless ? 
     233                  ENDIF 
     234                   
     235               ELSE ! if ice accretion zdh0 > 0 
     236                  ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 
     237                  zhbnew(ii,ij,0) = MIN( zdh0, hi_max(1) )  
     238               ENDIF 
     239                
    207240            ENDIF 
    208  
    209             ! 1 condition for remapping for the 1st category 
    210             ! H0+epsi < h1(t) < H1-epsi  
    211             !    h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible  
    212             !    in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
    213             IF( zht_i_b(ji,jj,1) < ( zhb0(ji,jj) + epsi10 ) )   zremap_flag(ji,jj) = 0 
    214             IF( zht_i_b(ji,jj,1) > ( zhb1(ji,jj) - epsi10 ) )   zremap_flag(ji,jj) = 0 
    215  
    216          END DO 
    217       END DO 
    218  
    219       !----------------------------------------------------------------------------------------------- 
    220       !  7) Compute g(h)  
    221       !----------------------------------------------------------------------------------------------- 
    222       !- 7.1 g(h) for category 1 at start of time step 
    223       CALL lim_itd_fitline( 1, zhb0, zhb1, zht_i_b(:,:,1), g0(:,:,1), g1(:,:,1), hL(:,:,1),   & 
    224          &                  hR(:,:,1), zremap_flag ) 
    225  
    226       !- 7.2 Area lost due to melting of thin ice (first category,  1) 
    227       DO ji = 1, nbrem 
    228          ii = nind_i(ji)  
    229          ij = nind_j(ji)  
    230  
    231          IF( a_i(ii,ij,1) > epsi10 ) THEN 
    232  
    233             zdh0 = zdhice(ii,ij,1) !decrease of ice thickness in the lower category 
    234  
    235             IF( zdh0 < 0.0 ) THEN      !remove area from category 1 
    236                zdh0 = MIN( -zdh0, hi_max(1) ) 
    237                !Integrate g(1) from 0 to dh0 to estimate area melted 
    238                zetamax = MIN( zdh0, hR(ii,ij,1) ) - hL(ii,ij,1) 
    239  
    240                IF( zetamax > 0.0 ) THEN 
    241                   zx1    = zetamax 
    242                   zx2    = 0.5 * zetamax * zetamax  
    243                   zda0   = g1(ii,ij,1) * zx2 + g0(ii,ij,1) * zx1                        ! ice area removed 
    244                   zdamax = a_i(ii,ij,1) * (1.0 - ht_i(ii,ij,1) / zht_i_b(ii,ij,1) ) ! Constrain new thickness <= ht_i                 
    245                   zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting  
    246                                                                                                 !     of thin ice (zdamax > 0) 
    247                   ! Remove area, conserving volume 
    248                   ht_i(ii,ij,1) = ht_i(ii,ij,1) * a_i(ii,ij,1) / ( a_i(ii,ij,1) - zda0 ) 
    249                   a_i(ii,ij,1)  = a_i(ii,ij,1) - zda0 
    250                   v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) ! clem-useless ? 
    251                ENDIF 
    252  
    253             ELSE ! if ice accretion zdh0 > 0 
    254                ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 
    255                zhbnew(ii,ij,0) = MIN( zdh0, hi_max(1) )  
    256             ENDIF 
    257  
    258          ENDIF 
    259  
    260       END DO 
    261  
    262       !- 7.3 g(h) for each thickness category   
    263       DO jl = 1, jpl 
    264          CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),  & 
    265             &                  g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) 
    266       END DO 
    267  
    268       !----------------------------------------------------------------------------------------------- 
    269       !  8) Compute area and volume to be shifted across each boundary (Eq. 18) 
    270       !----------------------------------------------------------------------------------------------- 
    271  
    272       DO jl = 1, jpl - 1 
    273          DO jj = 1, jpj 
    274             DO ji = 1, jpi 
    275                zdonor(ji,jj,jl) = 0 
    276                zdaice(ji,jj,jl) = 0.0 
    277                zdvice(ji,jj,jl) = 0.0 
    278             END DO 
    279          END DO 
    280  
    281          DO ji = 1, nbrem 
     241             
     242         END DO 
     243          
     244         !- 7.3 g(h) for each thickness category   
     245         DO jl = 1, jpl 
     246            CALL lim_itd_glinear( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),  & 
     247               &                  g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) 
     248         END DO 
     249          
     250         !----------------------------------------------------------------------------------------------- 
     251         !  8) Compute area and volume to be shifted across each boundary (Eq. 18) 
     252         !----------------------------------------------------------------------------------------------- 
     253          
     254         DO jl = 1, jpl - 1 
     255            DO jj = 1, jpj 
     256               DO ji = 1, jpi 
     257                  zdonor(ji,jj,jl) = 0 
     258                  zdaice(ji,jj,jl) = 0.0 
     259                  zdvice(ji,jj,jl) = 0.0 
     260               END DO 
     261            END DO 
     262             
     263            DO ji = 1, nidx 
     264               ii = nind_i(ji) 
     265               ij = nind_j(ji) 
     266                
     267               ! left and right integration limits in eta space 
     268               IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1 
     269                  zetamin = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl)       ! hi_max(jl) - hL  
     270                  zetamax = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) ! hR - hL 
     271                  zdonor(ii,ij,jl) = jl 
     272               ELSE                                    ! Hn* <= Hn => transfer from jl+1 to jl 
     273                  zetamin = 0.0 
     274                  zetamax = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1)   ! hi_max(jl) - hL 
     275                  zdonor(ii,ij,jl) = jl + 1 
     276               ENDIF 
     277               zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin 
     278                
     279               zx1  = zetamax - zetamin 
     280               zwk1 = zetamin * zetamin 
     281               zwk2 = zetamax * zetamax 
     282               zx2  = 0.5 * ( zwk2 - zwk1 ) 
     283               zwk1 = zwk1 * zetamin 
     284               zwk2 = zwk2 * zetamax 
     285               zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 ) 
     286               nd   = zdonor(ii,ij,jl) 
     287               zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 
     288               zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 
     289                
     290            END DO 
     291         END DO 
     292          
     293         !!---------------------------------------------------------------------------------------------- 
     294         !! 9) Shift ice between categories 
     295         !!---------------------------------------------------------------------------------------------- 
     296         CALL lim_itd_shiftice ( zdonor, zdaice, zdvice ) 
     297          
     298         !!---------------------------------------------------------------------------------------------- 
     299         !! 10) Make sure ht_i >= minimum ice thickness hi_min 
     300         !!---------------------------------------------------------------------------------------------- 
     301          
     302         DO ji = 1, nidx 
    282303            ii = nind_i(ji) 
    283304            ij = nind_j(ji) 
    284  
    285             ! left and right integration limits in eta space 
    286             IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1 
    287                zetamin = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl)       ! hi_max(jl) - hL  
    288                zetamax = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) ! hR - hL 
    289                zdonor(ii,ij,jl) = jl 
    290             ELSE                                    ! Hn* <= Hn => transfer from jl+1 to jl 
    291                zetamin = 0.0 
    292                zetamax = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1)   ! hi_max(jl) - hL 
    293                zdonor(ii,ij,jl) = jl + 1 
     305            IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN 
     306               a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin  
     307               ! MV MP 2016 
     308               IF ( nn_pnd_scheme > 0 ) THEN 
     309                  a_ip(ii,ij,1) = a_ip(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 
     310               ENDIF 
     311               ! END MV MP 2016 
     312               ht_i(ii,ij,1) = rn_himin 
    294313            ENDIF 
    295             zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin 
    296  
    297             zx1  = zetamax - zetamin 
    298             zwk1 = zetamin * zetamin 
    299             zwk2 = zetamax * zetamax 
    300             zx2  = 0.5 * ( zwk2 - zwk1 ) 
    301             zwk1 = zwk1 * zetamin 
    302             zwk2 = zwk2 * zetamax 
    303             zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 ) 
    304             nd   = zdonor(ii,ij,jl) 
    305             zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 
    306             zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 
    307  
    308          END DO 
    309       END DO 
    310  
    311       !!---------------------------------------------------------------------------------------------- 
    312       !! 9) Shift ice between categories 
    313       !!---------------------------------------------------------------------------------------------- 
    314       CALL lim_itd_shiftice ( zdonor, zdaice, zdvice ) 
    315  
    316       !!---------------------------------------------------------------------------------------------- 
    317       !! 10) Make sure ht_i >= minimum ice thickness hi_min 
    318       !!---------------------------------------------------------------------------------------------- 
    319  
    320       DO ji = 1, nbrem 
    321          ii = nind_i(ji) 
    322          ij = nind_j(ji) 
    323          IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN 
    324             a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin  
    325             ! MV MP 2016 
    326             IF ( nn_pnd_scheme > 0 ) THEN 
    327                a_ip(ii,ij,1) = a_ip(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 
    328             ENDIF 
    329             ! END MV MP 2016 
    330             ht_i(ii,ij,1) = rn_himin 
    331          ENDIF 
    332       END DO 
     314         END DO 
     315          
     316      ENDIF 
    333317 
    334318      CALL wrk_dealloc( jpi,jpj, zremap_flag ) 
    335319      CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) 
    336       CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b ) 
     320      CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR ) 
    337321      CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    338322      CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
    339       CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j )  
    340323      CALL wrk_dealloc( jpi,jpj, zhb0, zhb1 ) 
    341324 
     
    343326 
    344327 
    345    SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 
    346       !!------------------------------------------------------------------ 
    347       !!                ***  ROUTINE lim_itd_fitline *** 
     328   SUBROUTINE lim_itd_glinear( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 
     329      !!------------------------------------------------------------------ 
     330      !!                ***  ROUTINE lim_itd_glinear *** 
    348331      !! 
    349332      !! ** Purpose :   build g(h) satisfying area and volume constraints (Eq. 6 and 9) 
     
    404387      END DO 
    405388      ! 
    406    END SUBROUTINE lim_itd_fitline 
     389   END SUBROUTINE lim_itd_glinear 
    407390 
    408391 
     
    421404 
    422405      INTEGER ::   ji, jj, jl, jl2, jl1, jk   ! dummy loop indices 
    423       INTEGER ::   ii, ij                     ! indices when changing from 2D-1D is done 
    424  
    425       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zaTsfn 
    426       REAL(wp), POINTER, DIMENSION(:,:)   ::   zworka            ! temporary array used here 
    427  
    428       REAL(wp) ::   zdvsnow, zdesnow   ! snow volume and energy transferred 
    429       REAL(wp) ::   zdeice             ! ice energy transferred 
    430       REAL(wp) ::   zdsm_vice          ! ice salinity times volume transferred 
    431       REAL(wp) ::   zdo_aice           ! ice age times volume transferred 
    432       REAL(wp) ::   zdaTsf             ! aicen*Tsfcn transferred 
    433       ! MV MP 2016  
    434       REAL(wp) ::   zdapnd             ! pond fraction transferred 
    435       REAL(wp) ::   zdvpnd             ! pond volume transferred 
    436       ! END MV MP 2016  
    437  
    438       INTEGER, POINTER, DIMENSION(:) ::   nind_i, nind_j   ! compressed indices for i/j directions 
    439  
    440       INTEGER  ::   nbrem             ! number of cells with ice to transfer 
    441       !!------------------------------------------------------------------ 
    442  
    443       CALL wrk_alloc( jpi,jpj,jpl, zaTsfn ) 
    444       CALL wrk_alloc( jpi,jpj, zworka ) 
    445       CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 
     406 
     407      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zaTsfn 
     408      REAL(wp), DIMENSION(jpi,jpj)     ::   zworka            ! temporary array used here 
     409 
     410      REAL(wp) ::   ztrans             ! ice/snow transferred 
     411      !!------------------------------------------------------------------ 
    446412 
    447413      !---------------------------------------------------------------------------------------------- 
    448414      ! 1) Define a variable equal to a_i*T_su 
    449415      !---------------------------------------------------------------------------------------------- 
    450  
    451416      DO jl = 1, jpl 
    452          zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) 
     417         DO jj = 1, jpj 
     418            DO ji = 1, jpi 
     419               zaTsfn(ji,jj,jl) = a_i(ji,jj,jl) * t_su(ji,jj,jl) 
     420            END DO 
     421         END DO 
    453422      END DO 
    454  
     423       
    455424      !------------------------------------------------------------------------------- 
    456425      ! 2) Transfer volume and energy between categories 
    457426      !------------------------------------------------------------------------------- 
    458  
    459427      DO jl = 1, jpl - 1 
    460          nbrem = 0 
    461428         DO jj = 1, jpj 
    462429            DO ji = 1, jpi 
    463                IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny 
    464                   nbrem = nbrem + 1 
    465                   nind_i(nbrem) = ji 
    466                   nind_j(nbrem) = jj 
    467                ENDIF 
    468             END DO 
    469          END DO 
    470  
    471          DO ji = 1, nbrem  
    472             ii = nind_i(ji) 
    473             ij = nind_j(ji) 
    474  
    475             jl1 = zdonor(ii,ij,jl) 
    476             rswitch       = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi10 ) ) 
    477             zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch 
    478             IF( jl1 == jl) THEN   ;   jl2 = jl1+1 
    479             ELSE                  ;   jl2 = jl  
    480             ENDIF 
    481  
    482             ! Ice areas 
    483             a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 
    484             a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) 
    485  
    486             ! Ice volumes 
    487             v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl)  
    488             v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) 
    489  
    490             ! Snow volumes 
    491             zdvsnow        = v_s(ii,ij,jl1) * zworka(ii,ij) 
    492             v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 
    493             v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow  
    494  
    495             ! Snow heat content   
    496             zdesnow            = e_s(ii,ij,1,jl1) * zworka(ii,ij) 
    497             e_s(ii,ij,1,jl1)   = e_s(ii,ij,1,jl1) - zdesnow 
    498             e_s(ii,ij,1,jl2)   = e_s(ii,ij,1,jl2) + zdesnow 
    499  
    500             ! Ice age  
    501             zdo_aice           = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 
    502             oa_i(ii,ij,jl1)    = oa_i(ii,ij,jl1) - zdo_aice 
    503             oa_i(ii,ij,jl2)    = oa_i(ii,ij,jl2) + zdo_aice 
    504  
    505             ! Ice salinity 
    506             zdsm_vice          = smv_i(ii,ij,jl1) * zworka(ii,ij) 
    507             smv_i(ii,ij,jl1)   = smv_i(ii,ij,jl1) - zdsm_vice 
    508             smv_i(ii,ij,jl2)   = smv_i(ii,ij,jl2) + zdsm_vice 
    509  
    510             ! Surface temperature 
    511             zdaTsf             = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 
    512             zaTsfn(ii,ij,jl1)  = zaTsfn(ii,ij,jl1) - zdaTsf 
    513             zaTsfn(ii,ij,jl2)  = zaTsfn(ii,ij,jl2) + zdaTsf  
    514  
    515             ! MV MP 2016  
    516             IF ( nn_pnd_scheme > 0 ) THEN 
    517             ! Pond fraction 
    518                zdapnd             = a_ip(ii,ij,jl1) * zdaice(ii,ij,jl) 
    519                a_ip(ii,ij,jl1)    = a_ip(ii,ij,jl1) - zdapnd 
    520                a_ip(ii,ij,jl2)    = a_ip(ii,ij,jl2) + zdapnd 
    521  
    522             ! Pond volume 
    523                zdvpnd             = v_ip(ii,ij,jl1) * zdaice(ii,ij,jl) 
    524                v_ip(ii,ij,jl1)    = v_ip(ii,ij,jl1) - zdvpnd 
    525                v_ip(ii,ij,jl2)    = v_ip(ii,ij,jl2) + zdvpnd 
    526  
    527             ENDIF 
    528             ! END MV MP 2016  
    529  
    530          END DO 
    531  
     430 
     431               jl1 = zdonor(ji,jj,jl) 
     432               rswitch       = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl1) - epsi10 ) ) 
     433               zworka(ji,jj) = zdvice(ji,jj,jl) / MAX( v_i(ji,jj,jl1), epsi10 ) * rswitch 
     434 
     435               IF( jl1 == jl) THEN   ;   jl2 = jl1+1 
     436               ELSE                  ;   jl2 = jl  
     437               ENDIF 
     438                
     439               ! Ice areas 
     440               a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - zdaice(ji,jj,jl) 
     441               a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + zdaice(ji,jj,jl) 
     442                
     443               ! Ice volumes 
     444               v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - zdvice(ji,jj,jl)  
     445               v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + zdvice(ji,jj,jl) 
     446                
     447               ! Snow volumes 
     448               ztrans         = v_s(ji,jj,jl1) * zworka(ji,jj) 
     449               v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - ztrans 
     450               v_s(ji,jj,jl2) = v_s(ji,jj,jl2) + ztrans  
     451                
     452               ! Snow heat content   
     453               ztrans             = e_s(ji,jj,1,jl1) * zworka(ji,jj) 
     454               e_s(ji,jj,1,jl1)   = e_s(ji,jj,1,jl1) - ztrans 
     455               e_s(ji,jj,1,jl2)   = e_s(ji,jj,1,jl2) + ztrans 
     456                
     457               ! Ice age  
     458               ztrans             = oa_i(ji,jj,jl1) * zdaice(ji,jj,jl) 
     459               oa_i(ji,jj,jl1)    = oa_i(ji,jj,jl1) - ztrans 
     460               oa_i(ji,jj,jl2)    = oa_i(ji,jj,jl2) + ztrans 
     461                
     462               ! Ice salinity 
     463               ztrans             = smv_i(ji,jj,jl1) * zworka(ji,jj) 
     464               smv_i(ji,jj,jl1)   = smv_i(ji,jj,jl1) - ztrans 
     465               smv_i(ji,jj,jl2)   = smv_i(ji,jj,jl2) + ztrans 
     466                
     467               ! Surface temperature 
     468               ztrans             = t_su(ji,jj,jl1) * zdaice(ji,jj,jl) 
     469               zaTsfn(ji,jj,jl1)  = zaTsfn(ji,jj,jl1) - ztrans 
     470               zaTsfn(ji,jj,jl2)  = zaTsfn(ji,jj,jl2) + ztrans 
     471                
     472               ! MV MP 2016  
     473               IF ( nn_pnd_scheme > 0 ) THEN 
     474                  ! Pond fraction 
     475                  ztrans          = a_ip(ji,jj,jl1) * zdaice(ji,jj,jl) 
     476                  a_ip(ji,jj,jl1) = a_ip(ji,jj,jl1) - ztrans 
     477                  a_ip(ji,jj,jl2) = a_ip(ji,jj,jl2) + ztrans 
     478                   
     479                  ! Pond volume 
     480                  ztrans          = v_ip(ji,jj,jl1) * zdaice(ji,jj,jl) 
     481                  v_ip(ji,jj,jl1) = v_ip(ji,jj,jl1) - ztrans 
     482                  v_ip(ji,jj,jl2) = v_ip(ji,jj,jl2) + ztrans 
     483               ENDIF 
     484               ! END MV MP 2016  
     485                
     486            END DO 
     487         END DO 
     488          
    532489         ! Ice heat content 
    533490         DO jk = 1, nlay_i 
    534             DO ji = 1, nbrem 
    535                ii = nind_i(ji) 
    536                ij = nind_j(ji) 
    537  
    538                jl1 = zdonor(ii,ij,jl) 
    539                IF (jl1 == jl) THEN 
    540                   jl2 = jl+1 
    541                ELSE             ! n1 = n+1 
    542                   jl2 = jl  
    543                ENDIF 
    544  
    545                zdeice = e_i(ii,ij,jk,jl1) * zworka(ii,ij) 
    546                e_i(ii,ij,jk,jl1) =  e_i(ii,ij,jk,jl1) - zdeice 
    547                e_i(ii,ij,jk,jl2) =  e_i(ii,ij,jk,jl2) + zdeice  
    548             END DO 
    549          END DO 
    550  
    551       END DO                   ! boundaries, 1 to ncat-1 
    552  
    553       !----------------------------------------------------------------- 
     491            DO jj = 1, jpj 
     492               DO ji = 1, jpi 
     493                   
     494                  jl1 = zdonor(ji,jj,jl) 
     495                  IF(jl1 == jl) THEN  ;  jl2 = jl+1 
     496                  ELSE                ;  jl2 = jl 
     497                  ENDIF 
     498                   
     499                  ztrans            = e_i(ji,jj,jk,jl1) * zworka(ji,jj) 
     500                  e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - ztrans 
     501                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + ztrans  
     502               END DO 
     503            END DO 
     504         END DO 
     505          
     506      END DO                   ! boundaries, 1 to jpl-1 
     507       
    554508      ! Update ice thickness and temperature 
    555       !----------------------------------------------------------------- 
    556  
    557509      DO jl = 1, jpl 
    558510         DO jj = 1, jpj 
    559             DO ji = 1, jpi  
     511            DO ji = 1, jpi 
    560512               IF ( a_i(ji,jj,jl) > epsi10 ) THEN  
    561513                  ht_i(ji,jj,jl)  =  v_i   (ji,jj,jl) / a_i(ji,jj,jl)  
     
    567519            END DO 
    568520         END DO 
    569       END DO 
    570       ! 
    571       CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn ) 
    572       CALL wrk_dealloc( jpi,jpj, zworka ) 
    573       CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 
     521      END DO       
    574522      ! 
    575523   END SUBROUTINE lim_itd_shiftice 
     
    584532      !! ** Method  : 
    585533      !!------------------------------------------------------------------ 
    586       INTEGER ::   ji,jj, jl   ! dummy loop indices 
     534      INTEGER ::   ji, jj, jl   ! dummy loop indices 
    587535      INTEGER ::   zshiftflag          ! = .true. if ice must be shifted 
    588536 
    589       INTEGER , POINTER, DIMENSION(:,:,:) ::   zdonor           ! donor category index 
    590       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice   ! ice area and volume transferred 
     537      INTEGER , DIMENSION(jpi,jpj,jpl) ::   zdonor           ! donor category index 
     538      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zdaice, zdvice   ! ice area and volume transferred 
    591539      !!------------------------------------------------------------------ 
    592540       
    593       CALL wrk_alloc( jpi,jpj,jpl, zdonor )   ! interger 
    594       CALL wrk_alloc( jpi,jpj,jpl, zdaice, zdvice ) 
    595       !      
    596541      !------------------------------------------------------------------------------ 
    597542      ! 1) Compute ice thickness. 
     
    635580                  zshiftflag        = 1 
    636581                  zdonor(ji,jj,jl)  = jl  
    637                   ! begin TECLIM change 
    638                   !zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * 0.5_wp 
    639                   !zdvice(ji,jj,jl)  = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1)) * 0.5_wp 
    640                   ! end TECLIM change  
    641582                  ! clem: how much of a_i you send in cat sup is somewhat arbitrary 
    642583                  zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl)   
     
    645586            END DO 
    646587         END DO 
    647          IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
     588         IF(lk_mpp)   CALL mpp_max( zshiftflag ) ! clem: for reproducibility ??? 
    648589 
    649590         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
     
    680621         END DO 
    681622 
    682          IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
     623         IF(lk_mpp)   CALL mpp_max( zshiftflag ) ! clem: for reproducibility ??? 
    683624          
    684625         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
     
    692633      END DO 
    693634      ! 
    694       CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) 
    695       CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice ) 
    696  
    697635   END SUBROUTINE lim_itd_th_reb 
    698636 
     
    704642   SUBROUTINE lim_itd_th_rem 
    705643   END SUBROUTINE lim_itd_th_rem 
    706    SUBROUTINE lim_itd_fitline 
    707    END SUBROUTINE lim_itd_fitline 
     644   SUBROUTINE lim_itd_glinear 
     645   END SUBROUTINE lim_itd_glinear 
    708646   SUBROUTINE lim_itd_shiftice 
    709647   END SUBROUTINE lim_itd_shiftice 
Note: See TracChangeset for help on using the changeset viewer.