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 5134 for trunk – NEMO

Changeset 5134 for trunk


Ignore:
Timestamp:
2015-03-09T18:27:34+01:00 (9 years ago)
Author:
clem
Message:

LIM3: changes to constrain ice thickness after advection. Ice volume and concentration are advected while ice thickness is deduced. This could result in overly high thicknesses associated with very low concentrations (in the 5th category)

Location:
trunk/NEMOGCM/NEMO/LIM_SRC_3
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r5128 r5134  
    1919   USE ice              ! LIM variables 
    2020   USE dom_ice          ! LIM domain 
    21    USE limthd_lac       ! LIM 
    2221   USE limvar           ! LIM 
    23    USE in_out_manager   ! I/O manager 
    2422   USE lbclnk           ! lateral boundary condition - MPP exchanges 
    2523   USE lib_mpp          ! MPP library 
    … …  
    2725   USE prtctl           ! Print control 
    2826 
     27   USE in_out_manager   ! I/O manager 
    2928   USE iom              ! I/O manager 
    3029   USE lib_fortran      ! glob_sum 
    … …  
    158157      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 
    159158      !-----------------------------------------------------------------------------! 
    160       Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0                ! proport const for PE 
     159      Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0             ! proport const for PE 
    161160      ! 
    162161      CALL lim_itd_me_ridgeprep                                    ! prepare ridging 
    … …  
    275274         ! 3.4 Compute total area of ice plus open water after ridging. 
    276275         !-----------------------------------------------------------------------------! 
    277  
    278          CALL lim_itd_me_asumr 
     276         ! This is in general not equal to one because of divergence during transport 
     277         asum(:,:) = ato_i(:,:) 
     278         DO jl = 1, jpl 
     279            asum(:,:) = asum(:,:) + a_i(:,:,jl) 
     280         END DO 
    279281 
    280282         ! 3.5 Do we keep on iterating ??? 
    … …  
    341343 
    342344      ! Check if there is a ridging error 
    343       DO jj = 1, jpj 
    344          DO ji = 1, jpi 
    345             IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug 
    346                WRITE(numout,*) ' ' 
    347                WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 
    348                WRITE(numout,*) ' limitd_me ' 
    349                WRITE(numout,*) ' POINT : ', ji, jj 
    350                WRITE(numout,*) ' jpl, a_i, athorn ' 
    351                WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0) 
    352                DO jl = 1, jpl 
    353                   WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl) 
    354                END DO 
    355             ENDIF 
    356          END DO 
    357       END DO 
     345      IF( lwp ) THEN 
     346         DO jj = 1, jpj 
     347            DO ji = 1, jpi 
     348               IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug 
     349                  WRITE(numout,*) ' ' 
     350                  WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 
     351                  WRITE(numout,*) ' limitd_me ' 
     352                  WRITE(numout,*) ' POINT : ', ji, jj 
     353                  WRITE(numout,*) ' jpl, a_i, athorn ' 
     354                  WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0) 
     355                  DO jl = 1, jpl 
     356                     WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl) 
     357                  END DO 
     358               ENDIF 
     359            END DO 
     360         END DO 
     361      END IF 
    358362 
    359363      ! Conservation check 
    … …  
    364368      ENDIF 
    365369 
    366       !-----------------------------------------------------------------------------! 
    367       ! 6) Updating state variables and trend terms (done in limupdate) 
    368       !-----------------------------------------------------------------------------! 
     370      ! updates 
    369371      CALL lim_var_glo2eqv 
    370372      CALL lim_var_zapsmall 
    371373      CALL lim_var_agg( 1 )  
    372374 
    373  
    374       IF(ln_ctl) THEN     ! Control print 
     375      !-----------------------------------------------------------------------------! 
     376      ! control prints 
     377      !-----------------------------------------------------------------------------! 
     378      IF(ln_ctl) THEN  
    375379         CALL prt_ctl_info(' ') 
    376380         CALL prt_ctl_info(' - Cell values : ') 
    … …  
    473477                     ! PE gain from ridging ice 
    474478                     !---------------------------- 
    475                      strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl)/krdg(ji,jj,jl)     & 
     479                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) / krdg(ji,jj,jl)     & 
    476480                        * z1_3 * ( hrmax(ji,jj,jl)**2 + hrmin(ji,jj,jl)**2 +  hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )   
    477481                        !!(a**3-b**3)/(a-b) = a*a+ab+b*b                       
    … …  
    528532            DO ji = 2, jpim1 
    529533               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN ! ice is present 
    530                   zworka(ji,jj) = 4.0 * strength(ji,jj)              & 
    531                      &          + strength(ji-1,jj) * tmask(ji-1,jj,1) &   
    532                      &          + strength(ji+1,jj) * tmask(ji+1,jj,1) &   
    533                      &          + strength(ji,jj-1) * tmask(ji,jj-1,1) &   
    534                      &          + strength(ji,jj+1) * tmask(ji,jj+1,1)     
    535  
    536                   zworka(ji,jj) = zworka(ji,jj) /  & 
    537                      &           ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) ) 
     534                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              & 
     535                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &   
     536                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) & 
     537                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) ) 
    538538               ELSE 
    539539                  zworka(ji,jj) = 0._wp 
    … …  
    625625 
    626626      ! Compute total area of ice plus open water. 
    627       CALL lim_itd_me_asumr 
    628       ! This is in general not equal to one  
    629       ! because of divergence during transport 
     627      ! This is in general not equal to one because of divergence during transport 
     628      asum(:,:) = ato_i(:,:) 
     629      DO jl = 1, jpl 
     630         asum(:,:) = asum(:,:) + a_i(:,:,jl) 
     631      END DO 
    630632 
    631633      ! Compute cumulative thickness distribution function 
    … …  
    678680               DO ji = 1, jpi 
    679681                  IF( Gsum(ji,jj,jl) < rn_gstar) THEN 
    680                      athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * & 
    681                         (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari) 
     682                     athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * & 
     683                        &                        ( 2.0 - (Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari ) 
    682684                  ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN 
    683                      athorn(ji,jj,jl) = Gstari * (rn_gstar-Gsum(ji,jj,jl-1)) *  & 
    684                         (2.0 - (Gsum(ji,jj,jl-1)+rn_gstar)*Gstari) 
     685                     athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) *  & 
     686                        &                        ( 2.0 - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari ) 
    685687                  ELSE 
    686688                     athorn(ji,jj,jl) = 0.0 
    … …  
    693695         !                         
    694696         zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array 
    695  
    696697         DO jl = -1, jpl 
    697698            Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 
    698          END DO !jl 
     699         END DO 
    699700         DO jl = 0, jpl 
    700701             athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 
    701702         END DO 
    702703         ! 
    703       ENDIF ! nn_partfun 
     704      ENDIF 
    704705 
    705706      IF( ln_rafting ) THEN      ! Ridging and rafting ice participation functions 
    … …  
    729730      IF( ln_rafting ) THEN 
    730731 
    731          IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 ) THEN 
     732         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 .AND. lwp ) THEN 
    732733            DO jl = 1, jpl 
    733734               DO jj = 1, jpj 
    … …  
    793794               ENDIF 
    794795 
    795             END DO ! ji 
    796          END DO ! jj 
    797       END DO ! jl 
     796            END DO 
     797         END DO 
     798      END DO 
    798799 
    799800      ! Normalization factor : aksum, ensures mass conservation 
    … …  
    834835      INTEGER ::   icells            ! number of cells with aicen > puny 
    835836      REAL(wp) ::   hL, hR, farea, ztmelts    ! left and right limits of integration 
    836       REAL(wp) ::   zsstK            ! SST in Kelvin 
    837837 
    838838      INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices 
    … …  
    910910               ato_i(ji,jj) = 0._wp 
    911911            ENDIF 
    912          END DO !jj 
    913       END DO !ji 
     912         END DO 
     913      END DO 
    914914 
    915915      ! if negative open water area alert it 
    916       IF( neg_ato_i ) THEN       ! there is a bug 
     916      IF( neg_ato_i .AND. lwp ) THEN       ! there is a bug 
    917917         DO jj = 1, jpj  
    918918            DO ji = 1, jpi 
    … …  
    921921                  WRITE(numout,*) 'Ridging error: ato_i < 0' 
    922922                  WRITE(numout,*) 'ato_i : ', ato_i(ji,jj) 
    923                ENDIF               ! ato_i < -epsi10 
     923               ENDIF 
    924924            END DO 
    925925         END DO 
    … …  
    937937         smv_i_init(:,:,jl) = smv_i(:,:,jl) 
    938938         oa_i_init (:,:,jl) = oa_i (:,:,jl) 
    939       END DO !jl 
     939      END DO 
    940940 
    941941      esnwn_init(:,:,:) = e_s(:,:,1,:) 
    … …  
    968968                  indxi(icells) = ji 
    969969                  indxj(icells) = jj 
    970                ENDIF ! test on a_icen_init  
    971             END DO ! ji 
    972          END DO ! jj 
     970               ENDIF 
     971            END DO 
     972         END DO 
    973973 
    974974         large_afrac = .false. 
    … …  
    10941094            dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) 
    10951095 
    1096          END DO                 ! ij 
     1096         END DO 
    10971097 
    10981098         !-------------------------------------------------------------------- 
    … …  
    11121112               ! enthalpy of the trapped seawater (J/m2, >0) 
    11131113               ! clem: if sst>0, then ersw <0 (is that possible?) 
    1114                zsstK  = sst_m(ji,jj) + rt0 
    1115                ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) * r1_nlay_i 
     1114               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * sst_m(ji,jj) * r1_nlay_i 
    11161115 
    11171116               ! heat flux to the ocean 
    … …  
    11211120               erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 
    11221121 
    1123             END DO ! ij 
    1124          END DO !jk 
     1122            END DO 
     1123         END DO 
    11251124 
    11261125 
    … …  
    11311130                  jj = indxj(ij) 
    11321131                  eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk) 
    1133                END DO ! ij 
    1134             END DO !jk 
     1132               END DO 
     1133            END DO 
    11351134         ENDIF 
    11361135 
    1137          IF( large_afrac ) THEN   ! there is a bug 
     1136         IF( large_afrac .AND. lwp ) THEN   ! there is a bug 
    11381137            DO ij = 1, icells 
    11391138               ji = indxi(ij) 
    … …  
    11461145            END DO 
    11471146         ENDIF 
    1148          IF( large_afrft ) THEN  ! there is a bug 
     1147         IF( large_afrft .AND. lwp ) THEN  ! there is a bug 
    11491148            DO ij = 1, icells 
    11501149               ji = indxi(ij) 
    … …  
    11721171               ! Transfer area, volume, and energy accordingly. 
    11731172 
    1174                IF( hrmin(ji,jj,jl1) >= hi_max(jl2)  .OR.        & 
    1175                    hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN 
     1173               IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN 
    11761174                  hL = 0._wp 
    11771175                  hR = 0._wp 
    … …  
    11991197                  ji = indxi(ij) 
    12001198                  jj = indxj(ij) 
    1201                   e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj)*erdg2(ji,jj,jk) 
     1199                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj) * erdg2(ji,jj,jk) 
    12021200               END DO 
    12031201            END DO 
    … …  
    12131211               ! thickness category jl2, transfer area, volume, and energy accordingly. 
    12141212               ! 
    1215                IF( hraft(ji,jj,jl1) <= hi_max(jl2)  .AND.        & 
    1216                    hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN 
     1213               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN 
    12171214                  a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + arft2 (ji,jj) 
    12181215                  v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + virft (ji,jj) 
    … …  
    12301227                  ji = indxi(ij) 
    12311228                  jj = indxj(ij) 
    1232                   IF(  hraft(ji,jj,jl1)  <=  hi_max(jl2)   .AND.        & 
    1233                        hraft(ji,jj,jl1)  >   hi_max(jl2-1)  ) THEN 
     1229                  IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1)  ) THEN 
    12341230                     e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk) 
    12351231                  ENDIF 
    … …  
    12711267      ! 
    12721268   END SUBROUTINE lim_itd_me_ridgeshift 
    1273  
    1274  
    1275    SUBROUTINE lim_itd_me_asumr 
    1276       !!----------------------------------------------------------------------------- 
    1277       !!                ***  ROUTINE lim_itd_me_asumr *** 
    1278       !! 
    1279       !! ** Purpose :   finds total fractional area 
    1280       !! 
    1281       !! ** Method  :   Find the total area of ice plus open water in each grid cell. 
    1282       !!              This is similar to the aggregate_area subroutine except that the 
    1283       !!              total area can be greater than 1, so the open water area is  
    1284       !!              included in the sum instead of being computed as a residual.  
    1285       !!----------------------------------------------------------------------------- 
    1286       INTEGER ::   jl   ! dummy loop index 
    1287       !!----------------------------------------------------------------------------- 
    1288       ! 
    1289       asum(:,:) = ato_i(:,:)                    ! open water 
    1290       DO jl = 1, jpl                            ! ice categories 
    1291          asum(:,:) = asum(:,:) + a_i(:,:,jl) 
    1292       END DO 
    1293       ! 
    1294    END SUBROUTINE lim_itd_me_asumr 
    1295  
    12961269 
    12971270   SUBROUTINE lim_itd_me_init 
    … …  
    13511324   SUBROUTINE lim_itd_me_icestrength 
    13521325   END SUBROUTINE lim_itd_me_icestrength 
    1353    SUBROUTINE lim_itd_me_sort 
    1354    END SUBROUTINE lim_itd_me_sort 
    13551326   SUBROUTINE lim_itd_me_init 
    13561327   END SUBROUTINE lim_itd_me_init 

  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r5123 r5134  
    2525   USE ice              ! LIM-3 variables 
    2626   USE limvar           ! LIM-3 variables 
    27    USE limcons          ! LIM-3 conservation 
    2827   USE prtctl           ! Print control 
    2928   USE in_out_manager   ! I/O manager 
    … …  
    3837   PUBLIC   lim_itd_th_rem 
    3938   PUBLIC   lim_itd_th_reb 
    40    PUBLIC   lim_itd_fitline 
    41    PUBLIC   lim_itd_shiftice 
    4239 
    4340   !!---------------------------------------------------------------------- 
    … …  
    129126         DO jj = 1, jpj 
    130127            DO ji = 1, jpi 
    131                rswitch           = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) )     !0 if no ice and 1 if yes 
     128               rswitch           = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) )     !0 if no ice and 1 if yes 
    132129               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 
    133                rswitch           = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 
     130               rswitch           = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) )    !0 if no ice and 1 if yes 
    134131               zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 
    135                IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)  
     132!clem               IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)  
     133               zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)  
    136134            END DO 
    137135         END DO 
    … …  
    228226 
    229227            IF( a_i(ji,jj,kubnd) > epsi10 ) THEN 
    230                zhbnew(ji,jj,kubnd) = 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) 
     228               zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) ) 
    231229            ELSE 
    232230               zhbnew(ji,jj,kubnd) = hi_max(kubnd)   
    … …  
    235233            ENDIF 
    236234 
    237             IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 
    238  
    239          END DO !jj 
    240       END DO !jj 
     235         END DO 
     236      END DO 
    241237 
    242238      !----------------------------------------------------------------------------------------------- 
    … …  
    252248         ij = nind_j(ji)  
    253249 
    254          !ji 
    255250         IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 
    256251            zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 
    257             ! ji, a_i > epsi10 
    258252            IF( zdh0 < 0.0 ) THEN !remove area from category 1 
    259                ! ji, a_i > epsi10; zdh0 < 0 
    260253               zdh0 = MIN( -zdh0, hi_max(klbnd) ) 
    261254 
    … …  
    276269                  v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ? 
    277270               ENDIF 
    278                ! ji, a_i > epsi10 
    279  
    280             ELSE ! if ice accretion 
    281                ! ji, a_i > epsi10; zdh0 > 0 
     271 
     272            ELSE ! if ice accretion ! a_i > epsi10; zdh0 > 0 
    282273               zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) )  
    283274               ! zhbnew was 0, and is shifted to the right to account for thin ice 
    … …  
    285276            ENDIF ! zdh0  
    286277 
    287             ! a_i > epsi10 
    288278         ENDIF ! a_i > epsi10 
    289279 
    290       END DO ! ji 
     280      END DO 
    291281 
    292282      !- 7.3 g(h) for each thickness category   
    … …  
    359349         ij = nind_j(ji) 
    360350         IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN 
    361             a_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin  
     351            a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin  
    362352            ht_i(ii,ij,1) = rn_himin 
    363353         ENDIF 
    … …  
    515505      END DO 
    516506 
    517       !---------------------------------------------------------------------------------------------- 
    518       ! 2) Check for daice or dvice out of range, allowing for roundoff error 
    519       !---------------------------------------------------------------------------------------------- 
    520       ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl 
    521       ! has a small area, with h(n) very close to a boundary.  Then 
    522       ! the coefficients of g(h) are large, and the computed daice and 
    523       ! dvice can be in error. If this happens, it is best to transfer 
    524       ! either the entire category or nothing at all, depending on which 
    525       ! side of the boundary hice(n) lies. 
    526       !----------------------------------------------------------------- 
    527       DO jl = klbnd, kubnd-1 
    528  
    529          zdaice_negative = .false. 
    530          zdvice_negative = .false. 
    531          zdaice_greater_aicen = .false. 
    532          zdvice_greater_vicen = .false. 
    533  
    534          DO jj = 1, jpj 
    535             DO ji = 1, jpi 
    536  
    537                IF (zdonor(ji,jj,jl) > 0) THEN 
    538                   jl1 = zdonor(ji,jj,jl) 
    539  
    540                   IF (zdaice(ji,jj,jl) < 0.0) THEN 
    541                      IF (zdaice(ji,jj,jl) > -epsi10) THEN 
    542                         IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  & 
    543                              ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 
    544                            zdaice(ji,jj,jl) = a_i(ji,jj,jl1)  ! shift entire category 
    545                            zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
    546                         ELSE 
    547                            zdaice(ji,jj,jl) = 0.0 ! shift no ice 
    548                            zdvice(ji,jj,jl) = 0.0 
    549                         ENDIF 
    550                      ELSE 
    551                         zdaice_negative = .true. 
    552                      ENDIF 
    553                   ENDIF 
    554  
    555                   IF (zdvice(ji,jj,jl) < 0.0) THEN 
    556                      IF (zdvice(ji,jj,jl) > -epsi10 ) THEN 
    557                         IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  & 
    558                              ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 
    559                            zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 
    560                            zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
    561                         ELSE 
    562                            zdaice(ji,jj,jl) = 0.0    ! shift no ice 
    563                            zdvice(ji,jj,jl) = 0.0 
    564                         ENDIF 
    565                      ELSE 
    566                         zdvice_negative = .true. 
    567                      ENDIF 
    568                   ENDIF 
    569  
    570                   ! If daice is close to aicen, set daice = aicen. 
    571                   IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN 
    572                      IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN 
    573                         zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 
    574                         zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
    575                      ELSE 
    576                         zdaice_greater_aicen = .true. 
    577                      ENDIF 
    578                   ENDIF 
    579  
    580                   IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN 
    581                      IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN 
    582                         zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 
    583                         zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
    584                      ELSE 
    585                         zdvice_greater_vicen = .true. 
    586                      ENDIF 
    587                   ENDIF 
    588  
    589                ENDIF               ! donor > 0 
    590             END DO 
    591          END DO 
    592  
    593       END DO 
    594  
     507!clem: I think the following is wrong (if enabled, it creates negative concentration/volume around -epsi10) 
     508!      !---------------------------------------------------------------------------------------------- 
     509!      ! 2) Check for daice or dvice out of range, allowing for roundoff error 
     510!      !---------------------------------------------------------------------------------------------- 
     511!      ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl 
     512!      ! has a small area, with h(n) very close to a boundary.  Then 
     513!      ! the coefficients of g(h) are large, and the computed daice and 
     514!      ! dvice can be in error. If this happens, it is best to transfer 
     515!      ! either the entire category or nothing at all, depending on which 
     516!      ! side of the boundary hice(n) lies. 
     517!      !----------------------------------------------------------------- 
     518!      DO jl = klbnd, kubnd-1 
     519! 
     520!         zdaice_negative = .false. 
     521!         zdvice_negative = .false. 
     522!         zdaice_greater_aicen = .false. 
     523!         zdvice_greater_vicen = .false. 
     524! 
     525!         DO jj = 1, jpj 
     526!            DO ji = 1, jpi 
     527! 
     528!               IF (zdonor(ji,jj,jl) > 0) THEN 
     529!                  jl1 = zdonor(ji,jj,jl) 
     530! 
     531!                  IF (zdaice(ji,jj,jl) < 0.0) THEN 
     532!                     IF (zdaice(ji,jj,jl) > -epsi10) THEN 
     533!                        IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  & 
     534!                             ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 
     535!                           zdaice(ji,jj,jl) = a_i(ji,jj,jl1)  ! shift entire category 
     536!                           zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
     537!                        ELSE 
     538!                           zdaice(ji,jj,jl) = 0.0 ! shift no ice 
     539!                           zdvice(ji,jj,jl) = 0.0 
     540!                        ENDIF 
     541!                     ELSE 
     542!                        zdaice_negative = .true. 
     543!                     ENDIF 
     544!                  ENDIF 
     545! 
     546!                  IF (zdvice(ji,jj,jl) < 0.0) THEN 
     547!                     IF (zdvice(ji,jj,jl) > -epsi10 ) THEN 
     548!                        IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  & 
     549!                             ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 
     550!                           zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 
     551!                           zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
     552!                        ELSE 
     553!                           zdaice(ji,jj,jl) = 0.0    ! shift no ice 
     554!                           zdvice(ji,jj,jl) = 0.0 
     555!                        ENDIF 
     556!                    ELSE 
     557!                       zdvice_negative = .true. 
     558!                    ENDIF 
     559!                 ENDIF 
     560! 
     561!                  ! If daice is close to aicen, set daice = aicen. 
     562!                  IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN 
     563!                     IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN 
     564!                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 
     565!                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
     566!                    ELSE 
     567!                       zdaice_greater_aicen = .true. 
     568!                    ENDIF 
     569!                  ENDIF 
     570! 
     571!                  IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN 
     572!                     IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN 
     573!                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 
     574!                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
     575!                     ELSE 
     576!                        zdvice_greater_vicen = .true. 
     577!                     ENDIF 
     578!                  ENDIF 
     579! 
     580!               ENDIF               ! donor > 0 
     581!            END DO 
     582!         END DO 
     583! 
     584!      END DO 
     585!clem 
    595586      !------------------------------------------------------------------------------- 
    596587      ! 3) Transfer volume and energy between categories 
    … …  
    614605 
    615606            jl1 = zdonor(ii,ij,jl) 
    616             rswitch       = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) ) 
    617             zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch 
     607            rswitch       = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi20 ) ) 
     608            zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi20 ) * rswitch 
    618609            IF( jl1 == jl) THEN   ;   jl2 = jl1+1 
    619610            ELSE                  ;   jl2 = jl  
    … …  
    710701                  ht_i(ji,jj,jl)  =  v_i   (ji,jj,jl) / a_i(ji,jj,jl)  
    711702                  t_su(ji,jj,jl)  =  zaTsfn(ji,jj,jl) / a_i(ji,jj,jl)  
    712                   rswitch         =  1.0 - MAX( 0.0, SIGN( 1.0, -v_s(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 
    713703               ELSE 
    714704                  ht_i(ji,jj,jl)  = 0._wp 
    … …  
    765755         DO jj = 1, jpj 
    766756            DO ji = 1, jpi  
    767                IF( a_i(ji,jj,jl) > epsi10 ) THEN  
    768                   ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    769                ELSE 
    770                   ht_i(ji,jj,jl) = 0._wp 
    771                ENDIF 
     757               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
     758               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 
    772759            END DO 
    773760         END DO 
    … …  
    775762 
    776763      !------------------------------------------------------------------------------ 
    777       ! 2) Make sure thickness of cat klbnd is at least hi_max(klbnd) 
    778       !------------------------------------------------------------------------------ 
    779       DO jj = 1, jpj  
    780          DO ji = 1, jpi  
    781             IF( a_i(ji,jj,klbnd) > epsi10 ) THEN 
    782                IF( ht_i(ji,jj,klbnd) <= hi_max(0) .AND. hi_max(0) > 0._wp ) THEN 
    783                   a_i(ji,jj,klbnd)  = v_i(ji,jj,klbnd) / hi_max(0)  
    784                   ht_i(ji,jj,klbnd) = hi_max(0) 
    785                ENDIF 
    786             ENDIF 
    787          END DO 
    788       END DO 
    789  
    790       !------------------------------------------------------------------------------ 
    791       ! 3) If a category thickness is not in bounds, shift the 
     764      ! 2) If a category thickness is not in bounds, shift the 
    792765      ! entire area, volume, and energy to the neighboring category 
    793766      !------------------------------------------------------------------------------ 
    … …  
    818791                  zdonor(ji,jj,jl)  = jl  
    819792                  ! begin TECLIM change 
    820                   !zdaice(ji,jj,jl)  = a_i(ji,jj,jl) 
    821                   !zdvice(ji,jj,jl)  = v_i(ji,jj,jl) 
    822793                  !zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * 0.5_wp 
    823794                  !zdvice(ji,jj,jl)  = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1)) * 0.5_wp 
    824795                  ! end TECLIM change  
    825796                  ! clem: how much of a_i you send in cat sup is somewhat arbitrary 
    826                   zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi10 ) / ht_i(ji,jj,jl)   
    827                   zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi10 ) 
     797                  zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl)   
     798                  zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi20 ) 
    828799               ENDIF 
    829800            END DO 
    … …  
    839810         ENDIF 
    840811         ! 
    841       END DO                    ! jl 
     812      END DO 
    842813 
    843814      !---------------------------- 
    … …  
    888859 
    889860      !------------------------------------------------------------------------------ 
    890       ! 4) Conservation check 
     861      ! 3) Conservation check 
    891862      !------------------------------------------------------------------------------ 
    892863 

  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r5128 r5134  
    115115               DO ji = 1, jpi 
    116116                  !0 if no ice and 1 if yes 
    117                   rswitch = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi20 )  ) 
     117                  rswitch = MAX(  0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 )  ) 
    118118                  !Energy of melting q(S,T) [J.m-3] 
    119119                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i ) 
    … …  
    125125               DO ji = 1, jpi 
    126126                  !0 if no ice and 1 if yes 
    127                   rswitch = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi20 )  ) 
     127                  rswitch = MAX(  0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 )  ) 
    128128                  !Energy of melting q(S,T) [J.m-3] 
    129129                  e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s ) 
    … …  
    158158      DO jj = 1, jpj 
    159159         DO ji = 1, jpi 
    160             rswitch  = tmask(ji,jj,1) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice 
     160            rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 
    161161            ! 
    162162            !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 
    … …  
    360360      CALL lim_var_eqv2glo 
    361361 
     362      CALL lim_var_zapsmall 
    362363      !-------------------------------------------- 
    363364      ! Diagnostic thermodynamic growth rates 
    … …  
    413414 
    414415      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    415       IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    416416      !------------------------------------------------------------------------------| 
    417417      !  7) Add frazil ice growing in leads. 
    418418      !------------------------------------------------------------------------------| 
     419      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    419420      CALL lim_thd_lac 
    420421      CALL lim_var_glo2eqv    ! only for info 
    … …  
    513514            a_i_1d(ji)  = MAX( 0._wp, a_i_1d(ji) + zda_mel )  
    514515            ! adjust thickness 
    515             rswitch     = 1._wp - MAX( 0._wp , SIGN( 1._wp , - a_i_1d(ji) + epsi20 ) ) 
     516            rswitch     = MAX( 0._wp , SIGN( 1._wp , a_i_1d(ji) - epsi20 ) ) 
    516517            ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 ) 
    517518            ! retrieve total concentration 

  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r5128 r5134  
    244244         ! thickness change 
    245245         IF( zdh_s_pre(ji) > 0._wp ) THEN 
    246          rswitch        = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) ) 
     246         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 
    247247         zdh_s_mel (ji) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 
    248248         zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting  
    … …  
    266266            ! thickness change 
    267267            rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )  
    268             rswitch          = rswitch * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) + epsi20 ) ) )  
     268            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) )  
    269269            zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) 
    270270            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting 
    … …  
    321321      DO jk = 1, nlay_s 
    322322         DO ji = kideb,kiut 
    323             rswitch       =  MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) + epsi20 )  ) 
    324             q_s_1d(ji,jk) = ( 1._wp - rswitch ) / MAX( ht_s_1d(ji), epsi20 ) *             & 
    325               &            ( (   MAX( 0._wp, dh_s_tot(ji) )              ) * zqprec(ji) +  & 
     323            rswitch       = MAX(  0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 )  ) 
     324            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *                          & 
     325              &            ( (   MAX( 0._wp, dh_s_tot(ji) )               ) * zqprec(ji) +  & 
    326326              &              ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
    327327            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk) 
    … …  
    453453               q_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0) 
    454454                
    455             ENDIF ! fc_bo_i 
    456          END DO ! ji 
    457       END DO ! iter 
     455            ENDIF 
     456         END DO 
     457      END DO 
    458458 
    459459      ! Contribution to Energy and Salt Fluxes 
    … …  
    592592         zq_rema(ji)     = zq_su(ji) + zq_bo(ji)  
    593593         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow 
    594          rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, zq_s(ji) - epsi20 ) ) 
    595          zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( zq_s(ji), epsi20 ) 
     594         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) ) 
     595         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 ) 
    596596         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 
    597597         zdh_s_mel(ji)   = zdh_s_mel(ji) + zdeltah(ji,1)     
    … …  
    599599         ht_s_1d   (ji)  = ht_s_1d(ji)   + zdeltah(ji,1) 
    600600         
    601          zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji)                ! update available heat (J.m-2) 
     601         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1)                ! update available heat (J.m-2) 
    602602         ! heat used to melt snow 
    603          hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0) 
     603         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * q_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0) 
    604604         ! Contribution to mass flux 
    605605         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
    … …  
    630630 
    631631         ! entrapment during snow ice formation 
    632          ! new salinity difference stored (to be used in limthd_ent.F90) 
     632         ! new salinity difference stored (to be used in limthd_sal.F90) 
    633633         IF (  nn_icesal == 2  ) THEN 
    634             rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) ) 
     634            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 
    635635            ! salinity dif due to snow-ice formation 
    636             dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi10 ) * rswitch      
     636            dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch      
    637637            ! salinity dif due to bottom growth  
    638638            IF (  zf_tt(ji)  < 0._wp ) THEN 
    639                dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi10 ) * rswitch 
     639               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch 
    640640            ENDIF 
    641641         ENDIF 
    … …  
    663663         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 
    664664          
    665          ! Total ablation (to debug) 
    666          IF( ht_i_1d(ji) <= 0._wp )   a_i_1d(ji) = 0._wp 
    667  
    668       END DO !ji 
     665      END DO 
    669666 
    670667      ! 
    … …  
    676673         rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )  
    677674         t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0 
    678       END DO  ! ji 
     675      END DO 
    679676 
    680677      DO jk = 1, nlay_s 
    681678         DO ji = kideb,kiut 
    682679            ! mask enthalpy 
    683             rswitch       = MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  ) 
    684             q_s_1d(ji,jk) = ( 1.0 - rswitch ) * q_s_1d(ji,jk) 
     680            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  ) 
     681            q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk) 
    685682            ! recalculate t_s_1d from q_s_1d 
    686             t_s_1d(ji,jk) = rt0 + ( 1._wp - rswitch ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 
    687          END DO 
    688       END DO 
     683            t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 
     684         END DO 
     685      END DO 
     686 
     687      ! --- ensure that a_i = 0 where ht_i = 0 --- 
     688      WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 
    689689       
    690690      CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 

  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90

    r5123 r5134  
    132132      DO jk1 = 1, nlay_i 
    133133         DO ji = kideb, kiut 
    134             rswitch      = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi20 ) )  
     134            rswitch      = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) )  
    135135            qnew(ji,jk1) = rswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 ) 
    136136         ENDDO 

  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r5128 r5134  
    130130               DO ji = 1, jpi 
    131131                  !Energy of melting q(S,T) [J.m-3] 
    132                   rswitch          = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi20 )  )   !0 if no ice 
     132                  rswitch          = MAX(  0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 )  )   !0 if no ice 
    133133                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl), epsi20 ) * REAL( nlay_i, wp ) 
    134134               END DO 
    … …  
    482482         DO jl = 1, jpl 
    483483            DO ji = 1, nbpac 
    484                rswitch          = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) )  ! 0 if no ice and 1 if yes 
     484               rswitch          = MAX( 0._wp , SIGN( 1._wp , za_i_1d(ji,jl) - epsi20 ) )  ! 0 if no ice and 1 if yes 
    485485               zoa_i_1d(ji,jl)  = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * rswitch    
    486486            END DO  

  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r5128 r5134  
    6868      CHARACTER(len=80) ::   cltmp 
    6969      ! 
    70       REAL(wp), POINTER, DIMENSION(:,:)      ::   zsm, zs0at 
    71       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e 
    72       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ow 
    73       REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   zs0e 
     70      REAL(wp), POINTER, DIMENSION(:,:)      ::   zsm 
     71      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi 
     72      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0opw 
     73      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   z0ei 
    7474      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zviold, zvsold, zsmvold  ! old ice volume... 
    7575      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax                   ! old ice thickness 
    … …  
    8080      IF( nn_timing == 1 )  CALL timing_start('limtrp') 
    8181 
    82       CALL wrk_alloc( jpi,jpj,           zsm, zs0at, zatold, zeiold, zesold ) 
    83       CALL wrk_alloc( jpi,jpj,jpl,       zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e ) 
    84       CALL wrk_alloc( jpi,jpj,1,         zs0ow ) 
    85       CALL wrk_alloc( jpi,jpj,nlay_i+1,jpl, zs0e ) 
     82      CALL wrk_alloc( jpi,jpj,           zsm, zatold, zeiold, zesold ) 
     83      CALL wrk_alloc( jpi,jpj,jpl,       z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 
     84      CALL wrk_alloc( jpi,jpj,1,         z0opw ) 
     85      CALL wrk_alloc( jpi,jpj,nlay_i+1,jpl, z0ei ) 
    8686      CALL wrk_alloc( jpi,jpj,jpl,       zhimax, zviold, zvsold, zsmvold ) 
    8787 
    … …  
    118118         ! in case advection creates ice too thick. 
    119119         !--------------------------------------------------------------------- 
    120          zhimax(:,:,:) = ht_i(:,:,:) 
     120         zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) 
    121121         DO jl = 1, jpl 
    122122            DO jj = 2, jpjm1 
    123123               DO ji = 2, jpim1 
    124                   zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) ) 
     124                  zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) ) 
    125125               END DO 
    126126            END DO 
    … …  
    155155         ! transported fields                                         
    156156         !------------------------- 
    157          zs0ow(:,:,1) = ato_i(:,:) * e12t(:,:)              ! Open water area  
    158          DO jl = 1, jpl 
    159             zs0sn (:,:,jl)   = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume 
    160             zs0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume 
    161             zs0a  (:,:,jl)   = a_i  (:,:,jl) * e12t(:,:)    ! Ice area 
    162             zs0sm (:,:,jl)   = smv_i(:,:,jl) * e12t(:,:)    ! Salt content 
    163             zs0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content 
    164             zs0c0 (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content 
     157         z0opw(:,:,1) = ato_i(:,:) * e12t(:,:)             ! Open water area  
     158         DO jl = 1, jpl 
     159            z0snw (:,:,jl)  = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume 
     160            z0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume 
     161            z0ai  (:,:,jl)  = a_i  (:,:,jl) * e12t(:,:)    ! Ice area 
     162            z0smi (:,:,jl)  = smv_i(:,:,jl) * e12t(:,:)    ! Salt content 
     163            z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content 
     164            z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content 
    165165            DO jk = 1, nlay_i 
    166                zs0e  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content 
     166               z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content 
    167167            END DO 
    168168         END DO 
    … …  
    171171         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==! 
    172172            DO jt = 1, initad 
    173                CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   &             !--- ice open water area 
     173               CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
    174174                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    175                CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   & 
     175               CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   & 
    176176                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    177177               DO jl = 1, jpl 
    178                   CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
     178                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    179179                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    180                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   & 
     180                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   & 
    181181                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    182                   CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
     182                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    183183                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    184                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   & 
     184                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   & 
    185185                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    186                   CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
     186                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    187187                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    188                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   & 
     188                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   & 
    189189                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    190                   CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---      
     190                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---      
    191191                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    192                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   & 
     192                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl),   & 
    193193                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    194                   CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations --- 
     194                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations --- 
    195195                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    196                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &  
     196                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &  
    197197                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    198                   CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents --- 
     198                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents --- 
    199199                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    200                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   & 
     200                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl),   & 
    201201                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    202202                  DO jk = 1, nlay_i                                                                !--- ice heat contents --- 
    203                      CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     203                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    204204                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    205205                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    206                      CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     206                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    207207                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    208208                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    … …  
    212212         ELSE 
    213213            DO jt = 1, initad 
    214                CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   &             !--- ice open water area 
     214               CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area 
    215215                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    216                CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   & 
     216               CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   & 
    217217                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    218218               DO jl = 1, jpl 
    219                   CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
     219                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    220220                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    221                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   & 
     221                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   & 
    222222                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    223                   CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
     223                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    224224                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    225                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   & 
     225                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   & 
    226226                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    227                   CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
     227                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    228228                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    229                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   & 
     229                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   & 
    230230                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    231231 
    232                   CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
     232                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
    233233                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    234                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   & 
     234                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl),   & 
    235235                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    236                   CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
     236                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
    237237                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    238                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   & 
     238                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
    239239                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    240                   CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
     240                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
    241241                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    242                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   & 
     242                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl),   & 
    243243                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    244244                  DO jk = 1, nlay_i                                                           !--- ice heat contents --- 
    245                      CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     245                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    246246                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    247247                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    248                      CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     248                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    249249                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    250250                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    … …  
    257257         ! Recover the properties from their contents 
    258258         !------------------------------------------- 
    259          ato_i(:,:) = zs0ow(:,:,1) * r1_e12t(:,:) 
    260          DO jl = 1, jpl 
    261             v_i  (:,:,jl)   = zs0ice(:,:,jl) * r1_e12t(:,:) 
    262             v_s  (:,:,jl)   = zs0sn (:,:,jl) * r1_e12t(:,:) 
    263             smv_i(:,:,jl)   = zs0sm (:,:,jl) * r1_e12t(:,:) 
    264             oa_i (:,:,jl)   = zs0oi (:,:,jl) * r1_e12t(:,:) 
    265             a_i  (:,:,jl)   = zs0a (:,:,jl) * r1_e12t(:,:) 
    266             e_s  (:,:,1,jl) = zs0c0 (:,:,jl) * r1_e12t(:,:) 
     259         ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:) 
     260         DO jl = 1, jpl 
     261            v_i  (:,:,jl)   = z0ice(:,:,jl) * r1_e12t(:,:) 
     262            v_s  (:,:,jl)   = z0snw(:,:,jl) * r1_e12t(:,:) 
     263            smv_i(:,:,jl)   = z0smi(:,:,jl) * r1_e12t(:,:) 
     264            oa_i (:,:,jl)   = z0oi (:,:,jl) * r1_e12t(:,:) 
     265            a_i  (:,:,jl)   = z0ai (:,:,jl) * r1_e12t(:,:) 
     266            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:) 
    267267            DO jk = 1, nlay_i 
    268                e_i  (:,:,jk,jl) = zs0e  (:,:,jk,jl) * r1_e12t(:,:) 
     268               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:) 
    269269            END DO 
    270270         END DO 
    … …  
    293293         END DO 
    294294         ! 
    295          CALL lim_hdf( ato_i (:,:) )   ! Diffusion 
     295         CALL lim_hdf( ato_i (:,:) ) 
    296296 
    297297         !------------------------------------ 
    … …  
    356356 
    357357                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN 
     358 
    358359                     zvi  = v_i  (ji,jj,jl) 
    359360                     zvs  = v_s  (ji,jj,jl) 
    … …  
    361362                     zes  = e_s  (ji,jj,1,jl) 
    362363                     zei  = SUM( e_i(ji,jj,1:nlay_i,jl) ) 
    363                      zdv  = v_i(ji,jj,jl) - zviold(ji,jj,jl)    
    364                       
    365                      rswitch = 1._wp 
    366                      IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. & 
    367                         & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN                                           
    368                         ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 
    369                         rswitch        = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 
    370                         a_i(ji,jj,jl)  = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 
    371                      ELSE 
    372                         ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) ) 
    373                         rswitch        = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) ) 
    374                         a_i(ji,jj,jl)  = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 ) 
     364 
     365                     zdv  = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl)   
     366 
     367                     IF ( ( zdv >  0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. & 
     368                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN                                           
     369 
     370                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) ) 
     371                        a_i(ji,jj,jl)  = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 ) 
     372 
     373                        ! small correction due to *rswitch for a_i 
     374                        v_i  (ji,jj,jl)        = rswitch * v_i  (ji,jj,jl) 
     375                        v_s  (ji,jj,jl)        = rswitch * v_s  (ji,jj,jl) 
     376                        smv_i(ji,jj,jl)        = rswitch * smv_i(ji,jj,jl) 
     377                        e_s(ji,jj,1,jl)        = rswitch * e_s(ji,jj,1,jl) 
     378                        e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl) 
     379 
     380                        ! Update mass fluxes 
     381                        wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 
     382                        wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 
     383                        sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice  
     384                        hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 
     385                        hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0 
     386 
    375387                     ENDIF 
    376388 
    377                      ! small correction due to *rswitch for a_i 
    378                      v_i  (ji,jj,jl) = rswitch * v_i  (ji,jj,jl) 
    379                      v_s  (ji,jj,jl) = rswitch * v_s  (ji,jj,jl) 
    380                      smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl) 
    381                      e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl) 
    382                      e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl) 
    383  
    384                      ! Update mass fluxes 
    385                      wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 
    386                      wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 
    387                      sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice  
    388                      hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 
    389                      hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0 
    390389                  ENDIF 
    391390 
    … …  
    423422         vt_s (:,:) = 0._wp 
    424423         at_i (:,:) = 0._wp 
    425          ! 
    426424         DO jl = 1, jpl 
    427425            DO jj = 1, jpj 
    428426               DO ji = 1, jpi 
    429                   ! 
    430                   vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 
    431                   vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 
    432                   at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 
    433                END DO 
    434             END DO 
    435          END DO 
    436          ! ------------------------------------------------- 
    437  
    438          ! open water 
     427                  vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) 
     428                  vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) 
     429                  at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) 
     430               END DO 
     431            END DO 
     432         END DO 
     433 
     434         ! --- open water = 1 if at_i=0 -------------------------------- 
    439435         DO jj = 1, jpj 
    440436            DO ji = 1, jpi 
    441                ! open water = 1 if at_i=0 
    442437               rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 
    443438               ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj) 
    … …  
    457452      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' ) 
    458453      ! 
    459       CALL wrk_dealloc( jpi,jpj,           zsm, zs0at, zatold, zeiold, zesold ) 
    460       CALL wrk_dealloc( jpi,jpj,jpl,       zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e ) 
    461       CALL wrk_dealloc( jpi,jpj,1,         zs0ow ) 
    462       CALL wrk_dealloc( jpi,jpj,nlay_i+1,jpl, zs0e ) 
     454      CALL wrk_dealloc( jpi,jpj,           zsm, zatold, zeiold, zesold ) 
     455      CALL wrk_dealloc( jpi,jpj,jpl,       z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 
     456      CALL wrk_dealloc( jpi,jpj,1,         z0opw ) 
     457      CALL wrk_dealloc( jpi,jpj,nlay_i+1,jpl, z0ei ) 
    463458      CALL wrk_dealloc( jpi,jpj,jpl,       zviold, zvsold, zhimax, zsmvold ) 
    464459      ! 

  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r5123 r5134  
    6969      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    7070 
    71       !----------------- 
    72       ! zap small values 
    73       !----------------- 
    74       CALL lim_var_zapsmall 
    75  
    7671      CALL lim_var_glo2eqv 
    77       
    7872      !---------------------------------------------------- 
    79       ! Rebin categories with thickness out of bounds 
    80       !---------------------------------------------------- 
    81       IF ( jpl > 1 )   CALL lim_itd_th_reb(1, jpl) 
    82  
     73      ! ice concentration should not exceed amax  
     74      !----------------------------------------------------- 
    8375      at_i(:,:) = 0._wp 
    8476      DO jl = 1, jpl 
    … …  
    8678      END DO 
    8779 
    88       !---------------------------------------------------- 
    89       ! ice concentration should not exceed amax  
    90       !----------------------------------------------------- 
    9180      DO jl  = 1, jpl 
    9281         DO jj = 1, jpj 
    … …  
    9483               IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    9584                  a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
    96                   ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    9785               ENDIF 
    9886            END DO 
    9987         END DO 
    10088      END DO 
    101  
    102       at_i(:,:) = 0._wp 
    103       DO jl = 1, jpl 
    104          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    105       END DO 
    10689     
    107       ! -------------------------------------- 
    108       ! Final thickness distribution rebinning 
    109       ! -------------------------------------- 
     90      !---------------------------------------------------- 
     91      ! Rebin categories with thickness out of bounds 
     92      !---------------------------------------------------- 
    11093      IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) 
    11194 

  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r5128 r5134  
    5656      INTEGER, INTENT(in) ::   kt    ! number of iteration 
    5757      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    58       REAL(wp) ::   zh, zsal 
     58      REAL(wp) ::   zsal 
    5959      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    6060      !!------------------------------------------------------------------- 
    … …  
    6969      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    7070 
    71       !----------------- 
    72       ! zap small values 
    73       !----------------- 
    74       CALL lim_var_agg( 1 ) 
    75       CALL lim_var_zapsmall 
    76       CALL lim_var_glo2eqv 
    77  
    78       !---------------------------------------------------- 
    79       ! Rebin categories with thickness out of bounds 
    80       !---------------------------------------------------- 
    81       IF ( jpl > 1 )   CALL lim_itd_th_reb(1, jpl) 
    82  
    8371      !---------------------------------------------------------------------- 
    8472      ! Constrain the thickness of the smallest category above himin 
    8573      !---------------------------------------------------------------------- 
     74      CALL lim_var_glo2eqv 
    8675      DO jj = 1, jpj  
    8776         DO ji = 1, jpi 
    8877            IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN 
    89                zh             = rn_himin / ht_i(ji,jj,1) 
    90                ht_s(ji,jj,1) = ht_s(ji,jj,1) * zh 
    91                ht_i(ji,jj,1) = ht_i(ji,jj,1) * zh 
    92                a_i (ji,jj,1) = a_i(ji,jj,1)  / zh 
     78               a_i (ji,jj,1) = a_i(ji,jj,1) * ht_i(ji,jj,1) / rn_himin 
    9379            ENDIF 
    9480         END DO 
    … …  
    10894               IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    10995                  a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
    110                   ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    11196               ENDIF 
    11297            END DO 
    … …  
    11499      END DO 
    115100 
    116       at_i(:,:) = 0.0 
    117       DO jl = 1, jpl 
    118          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    119       END DO 
    120  
    121       ! -------------------------------------- 
    122       ! Final thickness distribution rebinning 
    123       ! -------------------------------------- 
     101      !---------------------------------------------------- 
     102      ! Rebin categories with thickness out of bounds 
     103      !---------------------------------------------------- 
    124104      IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl ) 
    125105 
    … …  
    140120                  ! salinity stays in bounds 
    141121                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 
    142                   smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) ) !+ rn_simin * ( 1._wp - rswitch ) * v_i(ji,jj,jl) 
     122                  smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) ) 
    143123                  ! associated salt flux 
    144124                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
    145                END DO ! ji 
    146             END DO ! jj 
    147          END DO !jl 
     125               END DO 
     126            END DO 
     127         END DO 
    148128      ENDIF 
    149129 

  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r5123 r5134  
    161161         DO jj = 1, jpj 
    162162            DO ji = 1, jpi 
    163                rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !0 if no ice and 1 if yes 
     163               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )   !0 if no ice and 1 if yes 
    164164               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 
    165165               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 
    … …  
    173173            DO jj = 1, jpj 
    174174               DO ji = 1, jpi 
    175                   rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !0 if no ice and 1 if yes 
     175                  rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) )   !0 if no ice and 1 if yes 
    176176                  sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * rswitch 
    177177               END DO 
    … …  
    190190               DO ji = 1, jpi 
    191191                  !                                                              ! Energy of melting q(S,T) [J.m-3] 
    192                   rswitch   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes 
     192                  rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes 
    193193                  zq_i    = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp)  
    194194                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rt0              ! Ice layer melt temperature 
    … …  
    215215               DO ji = 1, jpi 
    216216                  !Energy of melting q(S,T) [J.m-3] 
    217                   rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes 
     217                  rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes 
    218218                  zq_s  = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp) 
    219219                  ! 
    … …  
    233233            DO jj = 1, jpj 
    234234               DO ji = 1, jpi 
    235                   rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 
     235                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
    236236                  tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   & 
    237237                     &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  ) 
    … …  
    339339                     !                                      ! weighting the profile 
    340340                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) 
    341                   END DO ! ji 
    342                END DO ! jj 
    343             END DO ! jk 
    344          END DO ! jl 
     341                  END DO 
     342               END DO 
     343            END DO 
     344         END DO 
    345345         ! 
    346346      ENDIF ! nn_icesal 
    … …  
    384384            DO jj = 1, jpj 
    385385               DO ji = 1, jpi 
    386                   rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 
     386                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
    387387                  tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   & 
    388388                     &                      * r1_nlay_i / MAX( vt_i(ji,jj) , epsi10 ) 
    … …  
    572572               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch 
    573573 
    574                ! ice salinity must stay in bounds 
    575                IF(  nn_icesal == 2  ) THEN 
    576                   smv_i(ji,jj,jl) = MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) ) 
    577                ENDIF 
    578574               ! update exchanges with ocean 
    579575               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
Note: See TracChangeset for help on using the changeset viewer.