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/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 – NEMO

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)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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 
Note: See TracChangeset for help on using the changeset viewer.