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 10413 for NEMO/trunk/src/ICE/iceistate.F90 – NEMO

Ignore:
Timestamp:
2018-12-18T18:59:59+01:00 (5 years ago)
Author:
clem
Message:

merge dev_r9947_SI3_advection with the trunk. All sette tests passed. There is probably a conservation issue with the new advection scheme that I should solve but the routine icedyn_adv_umx.F90 is going to change anyway in the next couple of weeks. At worst, the old routine can be plugged without harm

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/iceistate.F90

    r10069 r10413  
    6464   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    6565   !! $Id$ 
    66    !! Software governed by the CeCILL license (see ./LICENSE) 
     66   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    6767   !!---------------------------------------------------------------------- 
    6868CONTAINS 
     
    174174         ! then we check whether the distribution fullfills 
    175175         ! volume and area conservation, positivity and ice categories bounds 
    176          zh_i_ini(:,:,:) = 0._wp  
    177          za_i_ini(:,:,:) = 0._wp 
    178          ! 
    179          DO jj = 1, jpj 
    180             DO ji = 1, jpi 
    181                ! 
    182                IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN 
    183  
    184                   ! find which category (jl0) the input ice thickness falls into 
    185                   jl0 = jpl 
    186                   DO jl = 1, jpl 
    187                      IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 
    188                         jl0 = jl 
    189                         CYCLE 
    190                      ENDIF 
    191                   END DO 
     176 
     177         IF( jpl == 1 ) THEN 
     178            ! 
     179            zh_i_ini(:,:,1) = zht_i_ini(:,:) 
     180            za_i_ini(:,:,1) = zat_i_ini(:,:)             
     181            ! 
     182         ELSE 
     183            zh_i_ini(:,:,:) = 0._wp  
     184            za_i_ini(:,:,:) = 0._wp 
     185            ! 
     186            DO jj = 1, jpj 
     187               DO ji = 1, jpi 
    192188                  ! 
    193                   itest(:) = 0 
    194                   i_fill   = jpl + 1                                            !------------------------------------ 
    195                   DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories 
    196                      !                                                          !------------------------------------ 
    197                      i_fill = i_fill - 1 
     189                  IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN 
     190 
     191                     ! find which category (jl0) the input ice thickness falls into 
     192                     jl0 = jpl 
     193                     DO jl = 1, jpl 
     194                        IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 
     195                           jl0 = jl 
     196                           CYCLE 
     197                        ENDIF 
     198                     END DO 
    198199                     ! 
    199                      zh_i_ini(ji,jj,:) = 0._wp  
    200                      za_i_ini(ji,jj,:) = 0._wp 
    201200                     itest(:) = 0 
    202                      ! 
    203                      IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1 
    204                         zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj) 
    205                         za_i_ini(ji,jj,1) = zat_i_ini(ji,jj) 
    206                      ELSE                         !-- case ice is thicker: fill categories >1 
    207                         ! thickness 
    208                         DO jl = 1, i_fill-1 
    209                            zh_i_ini(ji,jj,jl) = hi_mean(jl) 
    210                         END DO 
     201                     i_fill   = jpl + 1                                            !------------------------------------ 
     202                     DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories 
     203                        !                                                          !------------------------------------ 
     204                        i_fill = i_fill - 1 
    211205                        ! 
    212                         ! concentration 
    213                         za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 
    214                         DO jl = 1, i_fill - 1 
    215                            IF( jl /= jl0 )THEN 
    216                               zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) ) 
    217                               za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 
     206                        zh_i_ini(ji,jj,:) = 0._wp  
     207                        za_i_ini(ji,jj,:) = 0._wp 
     208                        itest(:) = 0 
     209                        ! 
     210                        IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1 
     211                           zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj) 
     212                           za_i_ini(ji,jj,1) = zat_i_ini(ji,jj) 
     213                        ELSE                         !-- case ice is thicker: fill categories >1 
     214                           ! thickness 
     215                           DO jl = 1, i_fill-1 
     216                              zh_i_ini(ji,jj,jl) = hi_mean(jl) 
     217                           END DO 
     218                           ! 
     219                           ! concentration 
     220                           za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 
     221                           DO jl = 1, i_fill - 1 
     222                              IF( jl /= jl0 )THEN 
     223                                 zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) ) 
     224                                 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 
     225                              ENDIF 
     226                           END DO 
     227 
     228                           ! last category 
     229                           za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) ) 
     230                           zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) ) 
     231                           zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 )  
     232 
     233                           ! correction if concentration of upper cat is greater than lower cat 
     234                           !   (it should be a gaussian around jl0 but sometimes it is not) 
     235                           IF ( jl0 /= jpl ) THEN 
     236                              DO jl = jpl, jl0+1, -1 
     237                                 IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN 
     238                                    zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl) 
     239                                    zh_i_ini(ji,jj,jl    ) = 0._wp 
     240                                    za_i_ini(ji,jj,jl    ) = 0._wp 
     241                                    za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  & 
     242                                       &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 ) 
     243                                 END IF 
     244                              ENDDO 
    218245                           ENDIF 
    219                         END DO 
    220  
    221                         ! last category 
    222                         za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) ) 
    223                         zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) ) 
    224                         zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 )  
    225  
    226                         ! correction if concentration of upper cat is greater than lower cat 
    227                         !   (it should be a gaussian around jl0 but sometimes it is not) 
    228                         IF ( jl0 /= jpl ) THEN 
    229                            DO jl = jpl, jl0+1, -1 
    230                               IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN 
    231                                  zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl) 
    232                                  zh_i_ini(ji,jj,jl    ) = 0._wp 
    233                                  za_i_ini(ji,jj,jl    ) = 0._wp 
    234                                  za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  & 
    235                                     &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 ) 
    236                               END IF 
    237                            ENDDO 
     246                           ! 
    238247                        ENDIF 
    239248                        ! 
     249                        ! Compatibility tests 
     250                        zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )           ! Test 1: area conservation 
     251                        IF ( zconv < epsi06 ) itest(1) = 1 
     252                        ! 
     253                        zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &         ! Test 2: volume conservation 
     254                           &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) ) 
     255                        IF ( zconv < epsi06 ) itest(2) = 1 
     256                        ! 
     257                        IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1           ! Test 3: thickness of the last category is in-bounds ? 
     258                        ! 
     259                        itest(4) = 1 
     260                        DO jl = 1, i_fill 
     261                           IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0                        ! Test 4: positivity of ice concentrations 
     262                        END DO 
     263                        !                                                          !---------------------------- 
     264                     END DO                                                        ! end iteration on categories 
     265                     !                                                             !---------------------------- 
     266                     IF( lwp .AND. SUM(itest) /= 4 ) THEN  
     267                        WRITE(numout,*) 
     268                        WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! ' 
     269                        WRITE(numout,*) ' !!!! Something is wrong in the SI3 initialization procedure ' 
     270                        WRITE(numout,*) 
     271                        WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:) 
     272                        WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 
     273                        WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 
    240274                     ENDIF 
    241275                     ! 
    242                      ! Compatibility tests 
    243                      zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )           ! Test 1: area conservation 
    244                      IF ( zconv < epsi06 ) itest(1) = 1 
    245                      ! 
    246                      zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &         ! Test 2: volume conservation 
    247                         &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) ) 
    248                      IF ( zconv < epsi06 ) itest(2) = 1 
    249                      ! 
    250                      IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1           ! Test 3: thickness of the last category is in-bounds ? 
    251                      ! 
    252                      itest(4) = 1 
    253                      DO jl = 1, i_fill 
    254                         IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0                        ! Test 4: positivity of ice concentrations 
    255                      END DO 
    256                      !                                                          !---------------------------- 
    257                   END DO                                                        ! end iteration on categories 
    258                   !                                                             !---------------------------- 
    259                   IF( lwp .AND. SUM(itest) /= 4 ) THEN  
    260                      WRITE(numout,*) 
    261                      WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! ' 
    262                      WRITE(numout,*) ' !!!! Something is wrong in the SI3 initialization procedure ' 
    263                      WRITE(numout,*) 
    264                      WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:) 
    265                      WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 
    266                      WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 
    267276                  ENDIF 
    268277                  ! 
    269                ENDIF 
    270                ! 
    271             END DO    
    272          END DO    
    273  
     278               END DO 
     279            END DO 
     280         ENDIF 
     281          
    274282         !--------------------------------------------------------------------- 
    275283         ! 4) Fill in sea ice arrays 
Note: See TracChangeset for help on using the changeset viewer.