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 11084 for NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE – NEMO

Ignore:
Timestamp:
2019-06-06T16:30:49+02:00 (5 years ago)
Author:
davestorkey
Message:

UKMO/NEMO_4.0_old_tidal_mixing : Update to be relative to rev 11081 of NEMO_4.0_mirror.

Location:
NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/icecor.F90

    r10888 r11084  
    6666         WRITE(numout,*) '~~~~~~~' 
    6767      ENDIF 
    68       ! 
    6968      !                             !----------------------------------------------------- 
    70       !                             !  ice thickness must exceed himin (for ice diff)    ! 
     69      !                             !  ice thickness must exceed himin (for temp. diff.) ! 
    7170      !                             !----------------------------------------------------- 
    7271      WHERE( a_i(:,:,:) >= epsi20 )   ;   h_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:) 
    … …  
    9796         END DO 
    9897      ENDIF 
    99  
    10098      !                             !----------------------------------------------------- 
    10199      !                             !  Rebin categories with thickness out of bounds     ! 

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/icectl.F90

    r10888 r11084  
    6969      !! 
    7070      REAL(wp) ::   zv, zs, zt, zfs, zfv, zft 
    71       REAL(wp) ::   zvmin, zamin, zamax  
     71      REAL(wp) ::   zvmin, zamin, zamax, zeimin, zesmin, zsmin 
    7272      REAL(wp) ::   zvtrp, zetrp 
    7373      REAL(wp) ::   zarea, zv_sill, zs_sill, zt_sill 
    … …  
    141141         zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t  ) * zconv 
    142142 
    143          zvmin = glob_min( 'icectl', v_i ) 
    144          zamax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
    145          zamin = glob_min( 'icectl', a_i ) 
     143         zamax  = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
     144         zvmin  = glob_min( 'icectl', v_i ) 
     145         zamin  = glob_min( 'icectl', a_i ) 
     146         zsmin  = glob_min( 'icectl', sv_i ) 
     147         zeimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 
     148         zesmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) 
    146149 
    147150         ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  
    … …  
    152155 
    153156         IF(lwp) THEN 
    154             IF ( ABS( zv   ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv 
    155             IF ( ABS( zs   ) > zs_sill )   WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs 
    156             IF ( ABS( zt   ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt 
    157             IF ( zvmin < -epsi10 )         WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin 
    158             IF ( zamax > MAX( rn_amax_n, rn_amax_s ) + epsi10   & 
    159                & .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' .AND. cd_routine /= 'Hbig' ) & 
    160                &                           WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
    161             IF ( zamin < -epsi10 )         WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
    162 !clem: the following check fails when using UMx advection scheme (see comments in icedyn_adv.F90) 
     157            ! check conservation issues 
     158            IF ( ABS( zv ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv 
     159            IF ( ABS( zs ) > zs_sill )   WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs 
     160            IF ( ABS( zt ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt 
     161            ! check maximum ice concentration 
     162            IF ( zamax > MAX( rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' )  & 
     163               &                         WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
     164            ! check negative values 
     165            IF ( zvmin  < 0. )           WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin 
     166            IF ( zamin  < 0. )           WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
     167            IF ( zsmin  < 0. )           WRITE(numout,*) 'violation s_i<0               (',cd_routine,') = ',zsmin 
     168            IF ( zeimin < 0. )           WRITE(numout,*) 'violation e_i<0               (',cd_routine,') = ',zeimin 
     169            IF ( zesmin < 0. )           WRITE(numout,*) 'violation e_s<0               (',cd_routine,') = ',zesmin 
     170!clem: the following check fails (I think...) 
    163171!            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'icedyn_adv' ) THEN 
    164172!                                           WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp 

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/icedyn.F90

    r10888 r11084  
    7474      INTEGER, INTENT(in) ::   kt     ! ice time step 
    7575      !! 
    76       INTEGER  ::   ji, jj, jl        ! dummy loop indices 
     76      INTEGER  ::   ji, jj        ! dummy loop indices 
    7777      REAL(wp) ::   zcoefu, zcoefv 
    78       REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max 
    79       REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i 
     78      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdivu_i 
    8079      !!-------------------------------------------------------------------- 
    8180      ! 
    … …  
    8988      ENDIF 
    9089      !                       
    91       IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization 
    92          tau_icebfr(:,:) = 0._wp 
    93          DO jl = 1, jpl 
    94             WHERE( h_i_b(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
    95          END DO 
    96       ENDIF 
    97       ! 
    98       !                                !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig) 
    99       DO jl = 1, jpl 
    100          DO jj = 2, jpjm1 
    101             DO ji = fs_2, fs_jpim1 
    102                zhip_max(ji,jj,jl) = MAX( epsi20, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj  ,jl), h_ip_b(ji  ,jj+1,jl), & 
    103                   &                                                h_ip_b(ji-1,jj  ,jl), h_ip_b(ji  ,jj-1,jl), & 
    104                   &                                                h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), & 
    105                   &                                                h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) ) 
    106                zhi_max (ji,jj,jl) = MAX( epsi20, h_i_b (ji,jj,jl), h_i_b (ji+1,jj  ,jl), h_i_b (ji  ,jj+1,jl), & 
    107                   &                                                h_i_b (ji-1,jj  ,jl), h_i_b (ji  ,jj-1,jl), & 
    108                   &                                                h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), & 
    109                   &                                                h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) ) 
    110                zhs_max (ji,jj,jl) = MAX( epsi20, h_s_b (ji,jj,jl), h_s_b (ji+1,jj  ,jl), h_s_b (ji  ,jj+1,jl), & 
    111                   &                                                h_s_b (ji-1,jj  ,jl), h_s_b (ji  ,jj-1,jl), & 
    112                   &                                                h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), & 
    113                   &                                                h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) ) 
    114             END DO 
    115          END DO 
    116       END DO 
    117       CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. ) 
    118       ! 
    119       ! 
    120       SELECT CASE( nice_dyn )           !-- Set which dynamics is running 
     90      ! retrieve thickness from volume for landfast param. and UMx advection scheme 
     91      WHERE( a_i(:,:,:) >= epsi20 ) 
     92         h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:) 
     93         h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:) 
     94      ELSEWHERE 
     95         h_i(:,:,:) = 0._wp 
     96         h_s(:,:,:) = 0._wp 
     97      END WHERE 
     98      ! 
     99      WHERE( a_ip(:,:,:) >= epsi20 ) 
     100         h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) 
     101      ELSEWHERE 
     102         h_ip(:,:,:) = 0._wp 
     103      END WHERE 
     104      ! 
     105      ! 
     106      SELECT CASE( nice_dyn )          !-- Set which dynamics is running 
    121107 
    122108      CASE ( np_dynALL )           !==  all dynamical processes  ==! 
    123          CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology   
    124          CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
    125          CALL ice_dyn_rdgrft( kt )                                                 ! -- ridging/rafting  
    126          CALL ice_cor       ( kt , 1 )                                             ! -- Corrections 
    127  
     109         ! 
     110         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology   
     111         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice 
     112         CALL ice_dyn_rdgrft( kt )                                          ! -- ridging/rafting  
     113         CALL ice_cor       ( kt , 1 )                                      ! -- Corrections 
     114         ! 
    128115      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==! 
    129          CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology   
    130          CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
    131          CALL Hpiling                                                              ! -- simple pile-up (replaces ridging/rafting) 
    132  
     116         ! 
     117         CALL ice_dyn_rhg   ( kt )                                          ! -- rheology   
     118         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice 
     119         CALL Hpiling                                                       ! -- simple pile-up (replaces ridging/rafting) 
     120         CALL ice_var_zapsmall                                              ! -- zap small areas 
     121         ! 
    133122      CASE ( np_dynADV1D )         !==  pure advection ==!   (1D) 
    134          ALLOCATE( zdivu_i(jpi,jpj) ) 
     123         ! 
    135124         ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- ! 
    136125         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length 
    … …  
    145134         END DO 
    146135         ! --- 
    147          CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice 
    148  
    149          ! diagnostics: divergence at T points  
    150          DO jj = 2, jpjm1 
    151             DO ji = 2, jpim1 
    152                zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    153                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
    154             END DO 
    155          END DO 
    156          CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 
    157          IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) ) 
    158  
    159          DEALLOCATE( zdivu_i ) 
    160  
     136         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice 
     137         ! 
    161138      CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities) 
    162          ALLOCATE( zdivu_i(jpi,jpj) ) 
     139         ! 
    163140         u_ice(:,:) = rn_uice * umask(:,:,1) 
    164141         v_ice(:,:) = rn_vice * vmask(:,:,1) 
    … …  
    166143         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 
    167144         ! --- 
    168          CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice 
    169  
    170          ! diagnostics: divergence at T points  
    171          DO jj = 2, jpjm1 
    172             DO ji = 2, jpim1 
    173                zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    174                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     145         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice 
     146 
     147      END SELECT 
     148      ! 
     149      ! 
     150      ! diagnostics: divergence at T points  
     151      IF( iom_use('icediv') ) THEN 
     152         ! 
     153         SELECT CASE( nice_dyn ) 
     154 
     155         CASE ( np_dynADV1D , np_dynADV2D ) 
     156 
     157            ALLOCATE( zdivu_i(jpi,jpj) ) 
     158            DO jj = 2, jpjm1 
     159               DO ji = 2, jpim1 
     160                  zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     161                     &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     162               END DO 
    175163            END DO 
    176          END DO 
    177          CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 
    178          IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) ) 
    179  
    180          DEALLOCATE( zdivu_i ) 
    181  
    182       END SELECT 
     164            CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 
     165            CALL iom_put( "icediv" , zdivu_i(:,:) ) 
     166            DEALLOCATE( zdivu_i ) 
     167 
     168         END SELECT 
     169         ! 
     170      ENDIF 
    183171      ! 
    184172      ! controls 
    … …  
    188176 
    189177 
    190    SUBROUTINE Hbig( phi_max, phs_max, phip_max ) 
    191       !!------------------------------------------------------------------- 
    192       !!                  ***  ROUTINE Hbig  *** 
    193       !! 
    194       !! ** Purpose : Thickness correction in case advection scheme creates 
    195       !!              abnormally tick ice or snow 
    196       !! 
    197       !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points 
    198       !!                 (before advection) and reduce it by adapting ice concentration 
    199       !!              2- check whether snow thickness is larger than the surrounding 9-points 
    200       !!                 (before advection) and reduce it by sending the excess in the ocean 
    201       !!              3- check whether snow load deplets the snow-ice interface below sea level$ 
    202       !!                 and reduce it by sending the excess in the ocean 
    203       !!              4- correct pond fraction to avoid a_ip > a_i 
    204       !! 
    205       !! ** input   : Max thickness of the surrounding 9-points 
    206       !!------------------------------------------------------------------- 
    207       REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    208       ! 
    209       INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    210       REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra 
    211       !!------------------------------------------------------------------- 
    212       ! controls 
    213       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    214       ! 
    215       CALL ice_var_zapsmall                       !-- zap small areas 
    216       ! 
    217       DO jl = 1, jpl 
    218          DO jj = 1, jpj 
    219             DO ji = 1, jpi 
    220                IF ( v_i(ji,jj,jl) > 0._wp ) THEN 
    221                   ! 
    222                   !                               ! -- check h_ip -- ! 
    223                   ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    224                   IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN 
    225                      zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) ) 
    226                      IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN 
    227                         a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
    228                      ENDIF 
    229                   ENDIF 
    230                   ! 
    231                   !                               ! -- check h_i -- ! 
    232                   ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
    233                   zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    234                   IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN 
    235                      a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
    236                   ENDIF 
    237                   ! 
    238                   !                               ! -- check h_s -- ! 
    239                   ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
    240                   zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl) 
    241                   IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN 
    242                      zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 
    243                      ! 
    244                      wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice 
    245                      hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
    246                      ! 
    247                      e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra 
    248                      v_s(ji,jj,jl)          = a_i(ji,jj,jl) * phs_max(ji,jj,jl) 
    249                   ENDIF            
    250                   ! 
    251                   !                               ! -- check snow load -- ! 
    252                   ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 
    253                   !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 
    254                   !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 
    255                   zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
    256                   IF( zvs_excess > 0._wp ) THEN 
    257                      zfra = ( v_s(ji,jj,jl) - zvs_excess ) / MAX( v_s(ji,jj,jl), epsi20 ) 
    258                      wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice 
    259                      hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
    260                      ! 
    261                      e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra 
    262                      v_s(ji,jj,jl)          = v_s(ji,jj,jl) - zvs_excess 
    263                   ENDIF 
    264                    
    265                ENDIF 
    266             END DO 
    267          END DO 
    268       END DO  
    269       !                                           !-- correct pond fraction to avoid a_ip > a_i 
    270       WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:) 
    271       ! 
    272       ! controls 
    273       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    274       ! 
    275    END SUBROUTINE Hbig 
    276  
    277  
    278178   SUBROUTINE Hpiling 
    279179      !!------------------------------------------------------------------- 
    … …  
    290190      ! controls 
    291191      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    292       ! 
    293       CALL ice_var_zapsmall                       !-- zap small areas 
    294192      ! 
    295193      at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    … …  
    337235         WRITE(numout,*) '~~~~~~~~~~~~' 
    338236         WRITE(numout,*) '   Namelist namdyn:' 
    339          WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynALL       = ', ln_dynALL 
    340          WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV    = ', ln_dynRHGADV 
    341          WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)     ln_dynADV1D     = ', ln_dynADV1D 
    342          WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                 ln_dynADV2D     = ', ln_dynADV2D 
    343          WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)    = (', rn_uice,',', rn_vice,')' 
    344          WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat       = ', rn_ishlat 
    345          WRITE(numout,*) '      Landfast: param from Lemieux 2016                       ln_landfast_L16 = ', ln_landfast_L16 
    346          WRITE(numout,*) '      Landfast: param from home made                          ln_landfast_home= ', ln_landfast_home 
    347          WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_depfra       = ', rn_depfra 
    348          WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr       = ', rn_icebfr 
    349          WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax      = ', rn_lfrelax 
    350          WRITE(numout,*) '         isotropic tensile strength                           rn_tensile      = ', rn_tensile 
     237         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr) ln_dynALL       = ', ln_dynALL 
     238         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                     ln_dynRHGADV    = ', ln_dynRHGADV 
     239         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)    ln_dynADV1D     = ', ln_dynADV1D 
     240         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                ln_dynADV2D     = ', ln_dynADV2D 
     241         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)   = (', rn_uice,',',rn_vice,')' 
     242         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat 
     243         WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16 
     244         WRITE(numout,*) '      Landfast: param from home made                         ln_landfast_home= ', ln_landfast_home 
     245         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra 
     246         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr 
     247         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax      = ', rn_lfrelax 
     248         WRITE(numout,*) '         isotropic tensile strength                          rn_tensile      = ', rn_tensile 
    351249         WRITE(numout,*) 
    352250      ENDIF 

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/icedyn_adv.F90

    r10888 r11084  
    6464      !!---------------------------------------------------------------------- 
    6565      INTEGER, INTENT(in) ::   kt   ! number of iteration 
    66       ! 
    67       INTEGER ::   jl   ! dummy loop indice 
    68       REAL(wp), DIMENSION(jpi,jpj) ::   zmask  ! fraction of time step with v_i < 0 
    6966      !!--------------------------------------------------------------------- 
    7067      ! 
    71       IF( ln_timing )   CALL timing_start('icedyn_adv') 
     68      ! controls 
     69      IF( ln_timing    )   CALL timing_start('icedyn_adv')                                                             ! timing 
     70      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    7271      ! 
    7372      IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    7675         WRITE(numout,*) '~~~~~~~~~~~' 
    7776      ENDIF 
    78        
    79       ! conservation test 
    80       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    81                       
    82       !---------- 
    83       ! Advection 
    84       !---------- 
     77      ! 
     78      !---------------! 
     79      !== Advection ==! 
     80      !---------------! 
    8581      SELECT CASE( nice_adv ) 
    8682      !                                !-----------------------! 
    8783      CASE( np_advUMx )                ! ULTIMATE-MACHO scheme ! 
    8884         !                             !-----------------------! 
    89          CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    90       !                                !-----------------------! 
     85         CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, & 
     86            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     87         !                             !-----------------------! 
    9188      CASE( np_advPRA )                ! PRATHER scheme        ! 
    9289         !                             !-----------------------! 
    93          CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     90         CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, & 
     91            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    9492      END SELECT 
    95  
    96       !---------------------------- 
    97       ! Debug the advection schemes 
    98       !---------------------------- 
    99       ! clem: At least one advection scheme above is not strictly positive => UMx 
    100       !       In Prather, I am not sure if the fields are bounded by 0 or not (it seems yes) 
    101       !       In UMx    , advected fields are not perfectly bounded and negative values can appear. 
    102       !                   These values are usually very small but in some occasions they can also be non-negligible 
    103       !                   Therefore one needs to bound the advected fields by 0 (though this is not a clean fix) 
    104       ! 
    105       ! record the negative values resulting from UMx 
    106       zmask(:,:) = 0._wp ! keep the init to 0 here 
    107       DO jl = 1, jpl 
    108          WHERE( v_i(:,:,jl) < 0._wp )   zmask(:,:) = 1._wp 
    109       END DO 
    110       IF( iom_use('iceneg_pres') )   CALL iom_put("iceneg_pres", zmask                                      )  ! fraction of time step with v_i < 0 
    111       IF( iom_use('iceneg_volu') )   CALL iom_put("iceneg_volu", SUM(MIN( v_i, 0. ), dim=3 )                )  ! negative ice volume (only) 
    112       IF( iom_use('iceneg_hfx' ) )   CALL iom_put("iceneg_hfx" , ( SUM(SUM( MIN( e_i(:,:,1:nlay_i,:), 0. )  &  ! negative ice heat content (only) 
    113          &                                                                  , dim=4 ), dim=3 ) )* r1_rdtice )  ! -- eq. heat flux [W/m2] 
    114       ! 
    115       ! ==> conservation is ensured by calling this routine below, 
    116       !     however the global ice volume is then changed by advection (but errors are small)  
    117       CALL ice_var_zapneg( ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    11893 
    11994      !------------ 

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/icedyn_adv_umx.F90

    r10888 r11084  
    1111   !!   'key_si3'                                       SI3 sea-ice model 
    1212   !!---------------------------------------------------------------------- 
    13    !!   ice_dyn_adv_umx   : update the tracer trend with the 3D advection trends using a TVD scheme 
     13   !!   ice_dyn_adv_umx   : update the tracer fields 
    1414   !!   ultimate_x(_y)    : compute a tracer value at velocity points using ULTIMATE scheme at various orders 
    15    !!   macho             : ??? 
    16    !!   nonosc_ice        : compute monotonic tracer fluxes by a non-oscillatory algorithm  
     15   !!   macho             : compute the fluxes 
     16   !!   nonosc_ice        : limit the fluxes using a non-oscillatory algorithm  
    1717   !!---------------------------------------------------------------------- 
    1818   USE phycst         ! physical constant 
    … …  
    3232 
    3333   PUBLIC   ice_dyn_adv_umx   ! called by icedyn_adv.F90 
    34        
    35    REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6 
    36    REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120 
    37     
    38    ! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid) 
    39    INTEGER ::   kn_limiter = 1 
    40  
    41    ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme 
    42    !   clem: if set to true, the 2D test case "diagonal advection" does not work (I do not understand why) 
    43    !         but in realistic cases, it avoids having very negative ice temperature (-50) at low ice concentration  
    44    LOGICAL ::   ll_neg = .TRUE. 
    45     
    46    ! alternate directions for upstream 
    47    LOGICAL ::   ll_upsxy = .TRUE. 
    48  
    49    ! alternate directions for high order 
    50    LOGICAL ::   ll_hoxy = .TRUE. 
    51     
    52    ! prelimiter: use it to avoid overshoot in H 
    53    !   clem: if set to true, the 2D test case "diagnoal advection" does not work (I do not understand why) 
    54    LOGICAL ::   ll_prelimiter_zalesak = .FALSE.  ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D 
    55  
    56  
     34   ! 
     35   INTEGER, PARAMETER ::   np_advS = 1         ! advection for S and T:    dVS/dt = -div(      uVS     ) => np_advS = 1 
     36   !                                                                    or dVS/dt = -div( uA * uHS / u ) => np_advS = 2 
     37   !                                                                    or dVS/dt = -div( uV * uS  / u ) => np_advS = 3 
     38   INTEGER, PARAMETER ::   np_limiter = 1      ! limiter: 1 = nonosc 
     39   !                                                      2 = superbee 
     40   !                                                      3 = h3 
     41   LOGICAL            ::   ll_upsxy  = .TRUE.   ! alternate directions for upstream 
     42   LOGICAL            ::   ll_hoxy   = .TRUE.   ! alternate directions for high order 
     43   LOGICAL            ::   ll_neg    = .TRUE.   ! if T interpolated at u/v points is negative or v_i < 1.e-6 
     44   !                                                 then interpolate T at u/v points using the upstream scheme 
     45   LOGICAL            ::   ll_prelim = .FALSE.  ! prelimiter from: Zalesak(1979) eq. 14 => not well defined in 2D 
     46   ! 
     47   REAL(wp)           ::   z1_6   = 1._wp /   6._wp   ! =1/6 
     48   REAL(wp)           ::   z1_120 = 1._wp / 120._wp   ! =1/120 
     49   ! 
     50   INTEGER, ALLOCATABLE, DIMENSION(:,:,:) ::   imsk_small, jmsk_small 
     51   ! 
    5752   !! * Substitutions 
    5853#  include "vectopt_loop_substitute.h90" 
    … …  
    6459CONTAINS 
    6560 
    66    SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, & 
     61   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, & 
    6762      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    6863      !!---------------------------------------------------------------------- 
    … …  
    7974      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity 
    8075      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity 
     76      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_i       ! ice thickness 
     77      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_s       ! snw thickness 
     78      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_ip      ! ice pond thickness 
    8179      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area 
    8280      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume 
    … …  
    9391      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
    9492      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
    95       REAL(wp) ::   zdt 
    96       REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! send zcflnow and receive zcflprv 
    97       REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box  
     93      REAL(wp) ::   zdt, zvi_cen 
     94      REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! for global communication 
     95      REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box 
    9896      REAL(wp), DIMENSION(jpi,jpj)     ::   zati1, zati2 
    99       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho 
    100       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai, z1_aip 
    101       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhvar 
     97      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zu_cat, zv_cat 
     98      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho, zua_ups, zva_ups 
     99      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai , z1_aip, zhvar 
     100      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max 
     101      ! 
     102      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs  
    102103      !!---------------------------------------------------------------------- 
    103104      ! 
    104105      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 
    105106      ! 
    106       ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! 
    107       !     When needed, the advection split is applied at the next time-step in order to avoid blocking global comm. 
    108       !     ...this should not affect too much the stability... Was ok on the tests we did... 
     107      ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 
     108      DO jl = 1, jpl 
     109         DO jj = 2, jpjm1 
     110            DO ji = fs_2, fs_jpim1 
     111               zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
     112                  &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
     113                  &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
     114                  &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
     115               zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
     116                  &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
     117                  &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
     118                  &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
     119               zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
     120                  &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
     121                  &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
     122                  &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
     123            END DO 
     124         END DO 
     125      END DO 
     126      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. ) 
     127      ! 
     128      ! 
     129      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- ! 
     130      !        Note: the advection split is applied at the next time-step in order to avoid blocking global comm. 
     131      !              this should not affect too much the stability 
    109132      zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 
    110133      zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 
    … …  
    116139      ELSE                         ;   icycle = 1 
    117140      ENDIF 
    118        
    119141      zdt = rdt_ice / REAL(icycle) 
    120142 
    … …  
    122144      zudy(:,:) = pu_ice(:,:) * e2u(:,:) 
    123145      zvdx(:,:) = pv_ice(:,:) * e1v(:,:) 
    124  
     146      ! 
     147      ! setup transport for each ice cat  
     148      DO jl = 1, jpl 
     149         zu_cat(:,:,jl) = zudy(:,:) 
     150         zv_cat(:,:,jl) = zvdx(:,:) 
     151      END DO 
     152      ! 
    125153      ! --- define velocity for advection: u*grad(H) --- ! 
    126154      DO jj = 2, jpjm1 
    … …  
    154182         END WHERE 
    155183         ! 
    156          ! set u*a=u for advection of A only  
    157          DO jl = 1, jpl 
    158             zua_ho(:,:,jl) = zudy(:,:) 
    159             zva_ho(:,:,jl) = zvdx(:,:) 
    160          END DO 
    161           
     184         ! setup a mask where advection will be upstream 
     185         IF( ll_neg ) THEN 
     186            IF( .NOT. ALLOCATED(imsk_small) )   ALLOCATE( imsk_small(jpi,jpj,jpl) )  
     187            IF( .NOT. ALLOCATED(jmsk_small) )   ALLOCATE( jmsk_small(jpi,jpj,jpl) )  
     188            DO jl = 1, jpl 
     189               DO jj = 1, jpjm1 
     190                  DO ji = 1, jpim1 
     191                     zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) ) 
     192                     IF( zvi_cen < epsi06) THEN   ;   imsk_small(ji,jj,jl) = 0 
     193                     ELSE                         ;   imsk_small(ji,jj,jl) = 1   ;   ENDIF 
     194                     zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) ) 
     195                     IF( zvi_cen < epsi06) THEN   ;   jmsk_small(ji,jj,jl) = 0 
     196                     ELSE                         ;   jmsk_small(ji,jj,jl) = 1   ;   ENDIF 
     197                  END DO 
     198               END DO 
     199            END DO 
     200         ENDIF 
     201         ! 
     202         ! ----------------------- ! 
     203         ! ==> start advection <== ! 
     204         ! ----------------------- ! 
     205         ! 
     206         !== Ice area ==! 
    162207         zamsk = 1._wp 
    163          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) !-- Ice area 
    164          zamsk = 0._wp 
    165          ! 
    166          zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
    167          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_i )                !-- Ice volume 
    168          ! 
    169          zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
    170          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_s )                !-- Snw volume 
    171          ! 
    172          zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 
    173          CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, psv_i )               !-- Salt content 
    174          ! 
    175          DO jk = 1, nlay_i 
    176             zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 
    177             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) )   !-- Ice heat content 
    178          END DO 
    179          ! 
    180          DO jk = 1, nlay_s 
    181             zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 
    182             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) )   !-- Snw heat content 
    183          END DO 
    184          ! 
    185          IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN                                                              !-- Ice Age 
    186             ! clem: in theory we should use the formulation below to advect the ice age, but the code is unable to deal with 
    187             !       fields that do not depend on volume (here oa_i depends on concentration). It creates abnormal ages that 
    188             !       spread into the domain. Therefore we cheat and consider that ice age should be advected as ice concentration 
    189             !!zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:) 
    190             !!CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i ) 
    191             ! set u*a=u for advection of ice age 
    192             DO jl = 1, jpl 
    193                zua_ho(:,:,jl) = zudy(:,:) 
    194                zva_ho(:,:,jl) = zvdx(:,:) 
    195             END DO 
     208         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zu_cat , zv_cat , zcu_box, zcv_box, & 
     209            &                                      pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho ) 
     210         ! 
     211         !                             ! --------------------------------- ! 
     212         IF( np_advS == 1 ) THEN       ! -- advection form: -div( uVS ) -- ! 
     213            !                          ! --------------------------------- ! 
     214            zamsk = 0._wp 
     215            !== Ice volume ==! 
     216            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
     217            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     218               &                                      zhvar, pv_i, zua_ups, zva_ups ) 
     219            !== Snw volume ==!          
     220            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
     221            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     222               &                                      zhvar, pv_s, zua_ups, zva_ups ) 
     223            ! 
    196224            zamsk = 1._wp 
    197             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, poa_i, poa_i ) 
     225            !== Salt content ==! 
     226            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 
     227               &                                      psv_i, psv_i ) 
     228            !== Ice heat content ==! 
     229            DO jk = 1, nlay_i 
     230               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 
     231                  &                                      pe_i(:,:,jk,:), pe_i(:,:,jk,:) ) 
     232            END DO 
     233            !== Snw heat content ==! 
     234            DO jk = 1, nlay_s 
     235               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 
     236                  &                                      pe_s(:,:,jk,:), pe_s(:,:,jk,:) ) 
     237            END DO 
     238            ! 
     239            !                          ! ------------------------------------------ ! 
     240         ELSEIF( np_advS == 2 ) THEN   ! -- advection form: -div( uA * uHS / u ) -- ! 
     241            !                          ! ------------------------------------------ ! 
    198242            zamsk = 0._wp 
     243            !== Ice volume ==! 
     244            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
     245            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     246               &                                      zhvar, pv_i, zua_ups, zva_ups ) 
     247            !== Snw volume ==!          
     248            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
     249            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     250               &                                      zhvar, pv_s, zua_ups, zva_ups ) 
     251            !== Salt content ==! 
     252            zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 
     253            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     254               &                                      zhvar, psv_i, zua_ups, zva_ups ) 
     255            !== Ice heat content ==! 
     256            DO jk = 1, nlay_i 
     257               zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 
     258               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, & 
     259                  &                                      zhvar, pe_i(:,:,jk,:), zua_ups, zva_ups ) 
     260            END DO 
     261            !== Snw heat content ==! 
     262            DO jk = 1, nlay_s 
     263               zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 
     264               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, & 
     265                  &                                      zhvar, pe_s(:,:,jk,:), zua_ups, zva_ups ) 
     266            END DO 
     267            ! 
     268            !                          ! ----------------------------------------- ! 
     269         ELSEIF( np_advS == 3 ) THEN   ! -- advection form: -div( uV * uS / u ) -- ! 
     270            !                          ! ----------------------------------------- ! 
     271            zamsk = 0._wp 
     272            ! 
     273            ALLOCATE( zuv_ho (jpi,jpj,jpl), zvv_ho (jpi,jpj,jpl),  & 
     274               &      zuv_ups(jpi,jpj,jpl), zvv_ups(jpi,jpj,jpl), z1_vi(jpi,jpj,jpl), z1_vs(jpi,jpj,jpl) ) 
     275            ! 
     276            ! inverse of Vi 
     277            WHERE( pv_i(:,:,:) >= epsi20 )   ;   z1_vi(:,:,:) = 1._wp / pv_i(:,:,:) 
     278            ELSEWHERE                        ;   z1_vi(:,:,:) = 0. 
     279            END WHERE 
     280            ! inverse of Vs 
     281            WHERE( pv_s(:,:,:) >= epsi20 )   ;   z1_vs(:,:,:) = 1._wp / pv_s(:,:,:) 
     282            ELSEWHERE                        ;   z1_vs(:,:,:) = 0. 
     283            END WHERE 
     284            ! 
     285            ! It is important to first calculate the ice fields and then the snow fields (because we use the same arrays) 
     286            ! 
     287            !== Ice volume ==! 
     288            zuv_ups = zua_ups 
     289            zvv_ups = zva_ups 
     290            zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 
     291            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     292               &                                      zhvar, pv_i, zuv_ups, zvv_ups, zuv_ho , zvv_ho ) 
     293            !== Salt content ==! 
     294            zhvar(:,:,:) = psv_i(:,:,:) * z1_vi(:,:,:) 
     295            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, & 
     296               &                                      zhvar, psv_i, zuv_ups, zvv_ups ) 
     297            !== Ice heat content ==! 
     298            DO jk = 1, nlay_i 
     299               zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_vi(:,:,:) 
     300               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, & 
     301                  &                                      zhvar, pe_i(:,:,jk,:), zuv_ups, zvv_ups ) 
     302            END DO 
     303            !== Snow volume ==!          
     304            zuv_ups = zua_ups 
     305            zvv_ups = zva_ups 
     306            zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 
     307            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 
     308               &                                      zhvar, pv_s, zuv_ups, zvv_ups, zuv_ho , zvv_ho ) 
     309            !== Snw heat content ==! 
     310            DO jk = 1, nlay_s 
     311               zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_vs(:,:,:) 
     312               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, & 
     313                  &                                      zhvar, pe_s(:,:,jk,:), zuv_ups, zvv_ups ) 
     314            END DO 
     315            ! 
     316            DEALLOCATE( zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs ) 
     317            ! 
    199318         ENDIF 
    200319         ! 
    201          IF ( ln_pnd_H12 ) THEN                                                                                               !-- melt ponds 
    202             ! set u*a=u for advection of Ap only  
    203             DO jl = 1, jpl 
    204                zua_ho(:,:,jl) = zudy(:,:) 
    205                zva_ho(:,:,jl) = zvdx(:,:) 
    206             END DO 
    207             ! 
     320         !== Ice age ==! 
     321         IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN 
    208322            zamsk = 1._wp 
    209             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! fraction 
     323            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 
     324               &                                      poa_i, poa_i ) 
     325         ENDIF 
     326         ! 
     327         !== melt ponds ==! 
     328         IF ( ln_pnd_H12 ) THEN 
     329            ! fraction 
     330            zamsk = 1._wp 
     331            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat , zv_cat , zcu_box, zcv_box, & 
     332               &                                      pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho ) 
     333            ! volume 
    210334            zamsk = 0._wp 
    211             ! 
    212335            zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:) 
    213             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_ip )                 ! volume 
     336            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     337               &                                      zhvar, pv_ip, zua_ups, zva_ups ) 
    214338         ENDIF 
    215339         ! 
     340         !== Open water area ==! 
    216341         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
    217342         DO jj = 2, jpjm1 
    218343            DO ji = fs_2, fs_jpim1 
    219                pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                                              !-- Open water area 
     344               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &  
    220345                  &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
    221346            END DO 
    222347         END DO 
    223          CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T',  1. ) 
    224          ! 
     348         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1. ) 
     349         ! 
     350         ! 
     351         ! --- Ensure non-negative fields and in-bound thicknesses --- ! 
     352         ! Remove negative values (conservation is ensured) 
     353         !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
     354         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     355         ! 
     356         ! Make sure ice thickness is not too big 
     357         !    (because ice thickness can be too large where ice concentration is very small) 
     358         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     359 
    225360      END DO 
    226361      ! 
    … …  
    228363 
    229364    
    230    SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho ) 
     365   SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox,  & 
     366      &                                            pt, ptc, pua_ups, pva_ups, pua_ho, pva_ho ) 
    231367      !!---------------------------------------------------------------------- 
    232368      !!                  ***  ROUTINE adv_umx  *** 
    … …  
    235371      !!                 tracers and add it to the general trend of tracer equations 
    236372      !! 
    237       !! **  Method  :   - calculate upstream fluxes and upstream solution for tracer H 
     373      !! **  Method  :   - calculate upstream fluxes and upstream solution for tracers V/A(=H) etc 
    238374      !!                 - calculate tracer H at u and v points (Ultimate) 
    239       !!                 - calculate the high order fluxes using alterning directions (Macho?) 
     375      !!                 - calculate the high order fluxes using alterning directions (Macho) 
    240376      !!                 - apply a limiter on the fluxes (nonosc_ice) 
    241       !!                 - convert this tracer flux to a tracer content flux (uH -> uV) 
    242       !!                 - calculate the high order solution for tracer content V 
     377      !!                 - convert this tracer flux to a "volume" flux (uH -> uV) 
     378      !!                 - apply a limiter a second time on the volumes fluxes (nonosc_ice) 
     379      !!                 - calculate the high order solution for V 
    243380      !! 
    244       !! ** Action : solve 2 equations => a) da/dt = -div(ua) 
    245       !!                                  b) dV/dt = -div(uV) using dH/dt = -u.grad(H) 
    246       !!             in eq. b), - fluxes uH are evaluated (with UMx) and limited (with nonosc_ice). This step is necessary to get a good H. 
    247       !!                        - then we convert this flux to a "volume" flux this way => uH*ua/u 
    248       !!                             where ua is the flux from eq. a) 
    249       !!                        - at last we estimate dV/dt = -div(uH*ua/u) 
     381      !! ** Action : solve 3 equations => a) dA/dt  = -div(uA) 
     382      !!                                  b) dV/dt  = -div(uV)  using dH/dt = -u.grad(H) 
     383      !!                                  c) dVS/dt = -div(uVS) using either dHS/dt = -u.grad(HS) or dS/dt = -u.grad(S) 
    250384      !! 
    251       !! ** Note : - this method can lead to small negative V (since we only limit H) => corrected in icedyn_adv.F90 conserving mass etc. 
    252       !!           - negative tracers at u-v points can also occur from the Ultimate scheme (usually at the ice edge) and the solution for now 
    253       !!             is to apply an upstream scheme when it occurs. A better solution would be to degrade the order of 
    254       !!             the scheme automatically by applying a mask of the ice cover inside Ultimate (not done). 
     385      !!             in eq. b), - fluxes uH are evaluated (with UMx) and limited with nonosc_ice. This step is necessary to get a good H. 
     386      !!                        - then we convert this flux to a "volume" flux this way => uH * uA / u 
     387      !!                             where uA is the flux from eq. a) 
     388      !!                             this "volume" flux is also limited with nonosc_ice (otherwise overshoots can occur) 
     389      !!                        - at last we estimate dV/dt = -div(uH * uA / u) 
     390      !! 
     391      !!             in eq. c), one can solve the equation for  S (ln_advS=T), then dVS/dt = -div(uV * uS  / u) 
     392      !!                                                or for HS (ln_advS=F), then dVS/dt = -div(uA * uHS / u)  
     393      !! 
     394      !! ** Note : - this method can lead to tiny negative V (-1.e-20) => set it to 0 while conserving mass etc. 
     395      !!           - At the ice edge, Ultimate scheme can lead to: 
     396      !!                              1) negative interpolated tracers at u-v points 
     397      !!                              2) non-zero interpolated tracers at u-v points eventhough there is no ice and velocity is outward 
     398      !!                              Solution for 1): apply an upstream scheme when it occurs. A better solution would be to degrade the order of 
     399      !!                                               the scheme automatically by applying a mask of the ice cover inside Ultimate (not done). 
     400      !!                              Solution for 2): we set it to 0 in this case 
    255401      !!           - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good. 
    256402      !!             Large values of H can appear for very small ice concentration, and when it does it messes the things up since we 
    257       !!             work on H (and not V). It probably comes from the prelimiter of zalesak which is coded for 1D and not 2D. 
     403      !!             work on H (and not V). It is partly related to the multi-category approach 
    258404      !!             Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 
    259       !!             concentration is small). 
    260       !! To-do: expand the prelimiter from zalesak to make it work in 2D 
    261       !!---------------------------------------------------------------------- 
    262       REAL(wp)                        , INTENT(in   )           ::   pamsk          ! advection of concentration (1) or other tracers (0) 
    263       INTEGER                         , INTENT(in   )           ::   kn_umx         ! order of the scheme (1-5=UM or 20=CEN2) 
    264       INTEGER                         , INTENT(in   )           ::   jt             ! number of sub-iteration 
    265       INTEGER                         , INTENT(in   )           ::   kt             ! number of iteration 
    266       REAL(wp)                        , INTENT(in   )           ::   pdt            ! tracer time-step 
    267       REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pu   , pv      ! 2 ice velocity components => u*e2 
    268       REAL(wp), DIMENSION(:,:,:)      , INTENT(in   )           ::   puc  , pvc     ! 2 ice velocity components => u*e2 or u*a*e2u 
    269       REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pubox, pvbox   ! upstream velocity 
    270       REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   pt             ! tracer field 
    271       REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   ptc            ! tracer content field 
    272       REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho ! high order u*a fluxes 
     405      !!             concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 
     406      !!             since sv_i and e_i are still good. 
     407      !!---------------------------------------------------------------------- 
     408      REAL(wp)                        , INTENT(in   )           ::   pamsk            ! advection of concentration (1) or other tracers (0) 
     409      INTEGER                         , INTENT(in   )           ::   kn_umx           ! order of the scheme (1-5=UM or 20=CEN2) 
     410      INTEGER                         , INTENT(in   )           ::   jt               ! number of sub-iteration 
     411      INTEGER                         , INTENT(in   )           ::   kt               ! number of iteration 
     412      REAL(wp)                        , INTENT(in   )           ::   pdt              ! tracer time-step 
     413      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pu   , pv        ! 2 ice velocity components => u*e2 
     414      REAL(wp), DIMENSION(:,:,:)      , INTENT(in   )           ::   puc  , pvc       ! 2 ice velocity components => u*e2 or u*a*e2u 
     415      REAL(wp), DIMENSION(:,:  )      , INTENT(in   )           ::   pubox, pvbox     ! upstream velocity 
     416      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   pt               ! tracer field 
     417      REAL(wp), DIMENSION(:,:,:)      , INTENT(inout)           ::   ptc              ! tracer content field 
     418      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout), OPTIONAL ::   pua_ups, pva_ups ! upstream u*a fluxes 
     419      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out), OPTIONAL ::   pua_ho, pva_ho   ! high order u*a fluxes 
    273420      ! 
    274421      INTEGER  ::   ji, jj, jl       ! dummy loop indices   
    … …  
    303450            DO jj = 1, jpjm1 
    304451               DO ji = 1, fs_jpim1 
    305                   IF( ABS( puc(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp ) THEN 
    306                      zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 
    307                      zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * puc(ji,jj,jl) / pu(ji,jj) 
     452                  IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 
     453                     zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj) 
     454                     zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj) 
    308455                  ELSE 
    309456                     zfu_ho (ji,jj,jl) = 0._wp 
    … …  
    311458                  ENDIF 
    312459                  ! 
    313                   IF( ABS( pvc(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp ) THEN 
    314                      zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 
    315                      zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pvc(ji,jj,jl) / pv(ji,jj) 
     460                  IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 
     461                     zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj) 
     462                     zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj) 
    316463                  ELSE 
    317464                     zfv_ho (ji,jj,jl) = 0._wp   
    … …  
    321468            END DO 
    322469         END DO 
     470 
     471         ! the new "volume" fluxes must also be "flux corrected" 
     472         ! thus we calculate the upstream solution and apply a limiter again 
     473         DO jl = 1, jpl 
     474            DO jj = 2, jpjm1 
     475               DO ji = fs_2, fs_jpim1 
     476                  ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) 
     477                  ! 
     478                  zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 
     479               END DO 
     480            END DO 
     481         END DO 
     482         CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T',  1. ) 
     483         ! 
     484         IF    ( np_limiter == 1 ) THEN 
     485            CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho ) 
     486         ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN 
     487            CALL limiter_x( pdt, pu, ptc, zfu_ups, zfu_ho ) 
     488            CALL limiter_y( pdt, pv, ptc, zfv_ups, zfv_ho ) 
     489         ENDIF 
     490         ! 
    323491      ENDIF 
    324       !                                   --ho 
    325       ! in case of advection of A: output u*a 
    326       ! ------------------------------------- 
     492      !                                   --ho    --ups 
     493      ! in case of advection of A: output u*a and u*a 
     494      ! ----------------------------------------------- 
    327495      IF( PRESENT( pua_ho ) ) THEN 
    328496         DO jl = 1, jpl 
    329497            DO jj = 1, jpjm1 
    330498               DO ji = 1, fs_jpim1 
    331                   pua_ho(ji,jj,jl) = zfu_ho(ji,jj,jl) 
    332                   pva_ho(ji,jj,jl) = zfv_ho(ji,jj,jl) 
    333                END DO 
     499                  pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 
     500                  pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 
     501              END DO 
    334502            END DO 
    335503         END DO 
    … …  
    499667         END DO 
    500668         ! 
    501          IF    ( kn_limiter == 1 ) THEN 
     669         IF    ( np_limiter == 1 ) THEN 
    502670            CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    503          ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN 
     671         ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN 
    504672            CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    505673            CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    … …  
    517685               END DO 
    518686            END DO 
    519             IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     687            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    520688 
    521689            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    … …  
    538706               END DO 
    539707            END DO 
    540             IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     708            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    541709 
    542710         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
    … …  
    549717               END DO 
    550718            END DO 
    551             IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     719            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    552720            ! 
    553721            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    … …  
    570738               END DO 
    571739            END DO 
    572             IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     740            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    573741 
    574742         ENDIF 
    575          IF( kn_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     743         IF( np_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    576744          
    577745      ENDIF 
    … …  
    609777         ! 
    610778         !                                                        !--  ultimate interpolation of pt at u-point  --! 
    611          CALL ultimate_x( kn_umx, pdt, pt, pu, zt_u, pfu_ho ) 
     779         CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho ) 
    612780         !                                                        !--  limiter in x --! 
    613          IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     781         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    614782         !                                                        !--  advective form update in zpt  --! 
    615783         DO jl = 1, jpl 
    … …  
    619787                     &                              + pt   (ji,jj,jl) * ( pu  (ji,jj   ) - pu  (ji-1,jj   ) ) * r1_e1e2t(ji,jj) & 
    620788                     &                                                                                        * pamsk           & 
    621                      &                             ) * pdt ) * tmask(ji,jj,1)  
     789                     &                             ) * pdt ) * tmask(ji,jj,1) 
    622790               END DO 
    623791            END DO 
    … …  
    627795         !                                                        !--  ultimate interpolation of pt at v-point  --! 
    628796         IF( ll_hoxy ) THEN 
    629             CALL ultimate_y( kn_umx, pdt, zpt, pv, zt_v, pfv_ho ) 
     797            CALL ultimate_y( pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho ) 
    630798         ELSE 
    631             CALL ultimate_y( kn_umx, pdt, pt , pv, zt_v, pfv_ho ) 
     799            CALL ultimate_y( pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho ) 
    632800         ENDIF 
    633801         !                                                        !--  limiter in y --! 
    634          IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     802         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    635803         !          
    636804         ! 
    … …  
    638806         ! 
    639807         !                                                        !--  ultimate interpolation of pt at v-point  --! 
    640          CALL ultimate_y( kn_umx, pdt, pt, pv, zt_v, pfv_ho ) 
     808         CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho ) 
    641809         !                                                        !--  limiter in y --! 
    642          IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
     810         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    643811         !                                                        !--  advective form update in zpt  --! 
    644812         DO jl = 1, jpl 
    … …  
    656824         !                                                        !--  ultimate interpolation of pt at u-point  --! 
    657825         IF( ll_hoxy ) THEN 
    658             CALL ultimate_x( kn_umx, pdt, zpt, pu, zt_u, pfu_ho ) 
     826            CALL ultimate_x( pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho ) 
    659827         ELSE 
    660             CALL ultimate_x( kn_umx, pdt, pt , pu, zt_u, pfu_ho ) 
     828            CALL ultimate_x( pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho ) 
    661829         ENDIF 
    662830         !                                                        !--  limiter in x --! 
    663          IF( kn_limiter == 2 .OR. kn_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
     831         IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    664832         ! 
    665833      ENDIF 
    666834 
    667       IF( kn_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
     835      IF( np_limiter == 1 )   CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 
    668836      ! 
    669837   END SUBROUTINE macho 
    670838 
    671839 
    672    SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 
     840   SUBROUTINE ultimate_x( pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 
    673841      !!--------------------------------------------------------------------- 
    674842      !!                    ***  ROUTINE ultimate_x  *** 
    … …  
    680848      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    681849      !!---------------------------------------------------------------------- 
     850      REAL(wp)                        , INTENT(in   ) ::   pamsk     ! advection of concentration (1) or other tracers (0) 
    682851      INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
    683852      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step 
    … …  
    806975            DO jj = 1, jpjm1 
    807976               DO ji = 1, fs_jpim1 
    808                   IF( pt_u(ji,jj,jl) < 0._wp ) THEN 
     977                  IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    809978                     pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    810979                        &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
    … …  
    826995    
    827996  
    828    SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 
     997   SUBROUTINE ultimate_y( pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 
    829998      !!--------------------------------------------------------------------- 
    830999      !!                    ***  ROUTINE ultimate_y  *** 
    … …  
    8361005      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74.  
    8371006      !!---------------------------------------------------------------------- 
     1007      REAL(wp)                        , INTENT(in   ) ::   pamsk     ! advection of concentration (1) or other tracers (0) 
    8381008      INTEGER                         , INTENT(in   ) ::   kn_umx    ! order of the scheme (1-5=UM or 20=CEN2) 
    8391009      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step 
    … …  
    9591129            DO jj = 1, jpjm1 
    9601130               DO ji = 1, fs_jpim1 
    961                   IF( pt_v(ji,jj,jl) < 0._wp ) THEN 
     1131                  IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    9621132                     pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
    9631133                        &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
    … …  
    10231193      !                        |      |      |        |    * 
    10241194      !            t_ups :       i-1     i       i+1       i+2    
    1025       IF( ll_prelimiter_zalesak ) THEN 
     1195      IF( ll_prelim ) THEN 
    10261196          
    10271197         DO jl = 1, jpl 
    … …  
    11021272               ! 
    11031273               !                                  ! up & down beta terms 
    1104                IF( zpos > 0._wp ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 
    1105                ELSE                    ; zbetup(ji,jj,jl) = 0._wp ! zbig 
     1274               ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!) 
     1275               IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 
     1276               ELSE                     ; zbetup(ji,jj,jl) = 0._wp ! zbig 
    11061277               ENDIF 
    11071278               ! 
    1108                IF( zneg > 0._wp ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 
    1109                ELSE                    ; zbetdo(ji,jj,jl) = 0._wp ! zbig 
     1279               IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 
     1280               ELSE                     ; zbetdo(ji,jj,jl) = 0._wp ! zbig 
    11101281               ENDIF 
    11111282               ! 
    … …  
    11491320         END DO 
    11501321 
    1151          ! clem test 
    1152 !!         DO jj = 2, jpjm1 
    1153 !!            DO ji = 2, fs_jpim1   ! vector opt. 
    1154 !!               zzt = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
    1155 !!                  &                           - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj)  & 
    1156 !!                  &                     + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
    1157 !!                  &                     + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 
    1158 !!                  &         ) * tmask(ji,jj,1) 
    1159 !!               IF( zzt < -epsi20 ) THEN 
    1160 !!                  WRITE(numout,*) 'T<0 nonosc_ice',zzt 
    1161 !!               ENDIF 
    1162 !!            END DO 
    1163 !!         END DO 
    1164  
    11651322      END DO 
    11661323      ! 
    … …  
    12031360               Rjp = zslpx(ji+1,jj,jl) 
    12041361 
    1205                IF( kn_limiter == 3 ) THEN 
     1362               IF( np_limiter == 3 ) THEN 
    12061363 
    12071364                  IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
    … …  
    12191376                  pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 
    12201377 
    1221                ELSEIF( kn_limiter == 2 ) THEN 
     1378               ELSEIF( np_limiter == 2 ) THEN 
    12221379                  IF( Rj /= 0. ) THEN 
    12231380                     IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj 
    … …  
    12981455               Rjp = zslpy(ji,jj+1,jl) 
    12991456 
    1300                IF( kn_limiter == 3 ) THEN 
     1457               IF( np_limiter == 3 ) THEN 
    13011458 
    13021459                  IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
    … …  
    13141471                  pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 
    13151472 
    1316                ELSEIF( kn_limiter == 2 ) THEN 
     1473               ELSEIF( np_limiter == 2 ) THEN 
    13171474 
    13181475                  IF( Rj /= 0. ) THEN 
    … …  
    13581515   END SUBROUTINE limiter_y 
    13591516 
     1517 
     1518   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     1519      !!------------------------------------------------------------------- 
     1520      !!                  ***  ROUTINE Hbig  *** 
     1521      !! 
     1522      !! ** Purpose : Thickness correction in case advection scheme creates 
     1523      !!              abnormally tick ice or snow 
     1524      !! 
     1525      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points 
     1526      !!                 (before advection) and reduce it by adapting ice concentration 
     1527      !!              2- check whether snow thickness is larger than the surrounding 9-points 
     1528      !!                 (before advection) and reduce it by sending the excess in the ocean 
     1529      !!              3- check whether snow load deplets the snow-ice interface below sea level$ 
     1530      !!                 and reduce it by sending the excess in the ocean 
     1531      !!              4- correct pond fraction to avoid a_ip > a_i 
     1532      !! 
     1533      !! ** input   : Max thickness of the surrounding 9-points 
     1534      !!------------------------------------------------------------------- 
     1535      REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step 
     1536      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
     1537      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip 
     1538      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
     1539      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i 
     1540      ! 
     1541      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
     1542      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zvs_excess, zfra 
     1543      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch 
     1544      !!------------------------------------------------------------------- 
     1545      ! 
     1546      z1_dt = 1._wp / pdt 
     1547      ! 
     1548      DO jl = 1, jpl 
     1549 
     1550         DO jj = 1, jpj 
     1551            DO ji = 1, jpi 
     1552               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     1553                  ! 
     1554                  !                               ! -- check h_ip -- ! 
     1555                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
     1556                  IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     1557                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
     1558                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
     1559                        pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
     1560                     ENDIF 
     1561                  ENDIF 
     1562                  ! 
     1563                  !                               ! -- check h_i -- ! 
     1564                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
     1565                  zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 
     1566                  IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1567                     pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
     1568                  ENDIF 
     1569                  ! 
     1570                  !                               ! -- check h_s -- ! 
     1571                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
     1572                  zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 
     1573                  IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1574                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 
     1575                     ! 
     1576                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt 
     1577                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     1578                     ! 
     1579                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     1580                     pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 
     1581                  ENDIF            
     1582                  ! 
     1583                  !                               ! -- check snow load -- ! 
     1584                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 
     1585                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 
     1586                  !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 
     1587                  zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
     1588                  IF( zvs_excess > 0._wp ) THEN 
     1589                     zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 
     1590                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 
     1591                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     1592                     ! 
     1593                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     1594                     pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess 
     1595                  ENDIF 
     1596                   
     1597               ENDIF 
     1598            END DO 
     1599         END DO 
     1600      END DO  
     1601      !                                           !-- correct pond fraction to avoid a_ip > a_i 
     1602      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:) 
     1603      ! 
     1604      ! 
     1605   END SUBROUTINE Hbig 
     1606    
    13601607#else 
    13611608   !!---------------------------------------------------------------------- 

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/icedyn_rdgrft.F90

    r10888 r11084  
    142142      INTEGER, PARAMETER ::   jp_itermax = 20     
    143143      !!------------------------------------------------------------------- 
    144       ! clem: The redistribution of ice between categories can lead to small negative values (as for the remapping in ice_itd_rem) 
    145       !       likely due to truncation error ( i.e. 1. - 1. /= 0 ) 
    146       !       I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90) 
    147        
    148144      ! controls 
    149145      IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing 
    … …  
    156152      ENDIF       
    157153 
    158       CALL ice_var_zapsmall   ! Zero out categories with very small areas 
    159  
    160154      !-------------------------------- 
    161155      ! 0) Identify grid cells with ice 
    162156      !-------------------------------- 
     157      at_i(:,:) = SUM( a_i, dim=3 ) 
     158      ! 
    163159      npti = 0   ;   nptidx(:) = 0 
    164160      ipti = 0   ;   iptidx(:) = 0 
    165161      DO jj = 1, jpj 
    166162         DO ji = 1, jpi 
    167             IF ( at_i(ji,jj) > 0._wp ) THEN 
     163            IF ( at_i(ji,jj) > epsi10 ) THEN 
    168164               npti           = npti + 1 
    169165               nptidx( npti ) = (jj - 1) * jpi + ji 
    … …  
    178174         
    179175         ! just needed here 
    180          CALL tab_2d_1d( npti, nptidx(1:npti), zdivu(1:npti), divu_i(:,:) ) 
    181          CALL tab_2d_1d( npti, nptidx(1:npti), zdelt(1:npti), delta_i(:,:) ) 
     176         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i ) 
     177         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i ) 
    182178         ! needed here and in the iteration loop 
    183          CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i(:,:,:) ) 
    184          CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i(:,:,:) ) 
    185          CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i(:,:) ) 
     179         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i  ) 
     180         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  ) 
     181         CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i ) 
    186182 
    187183         DO ji = 1, npti 
    … …  
    310306 
    311307      !                       ! Ice thickness needed for rafting 
    312       WHERE( pa_i(1:npti,:) > 0._wp )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 
    313       ELSEWHERE                         ;   zhi(1:npti,:) = 0._wp 
     308      WHERE( pa_i(1:npti,:) > epsi20 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 
     309      ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp 
    314310      END WHERE 
    315311 
    … …  
    329325      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 
    330326      ! 
    331       WHERE( zasum(1:npti) > 0._wp )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
    332       ELSEWHERE                        ;   z1_asum(1:npti) = 0._wp 
     327      WHERE( zasum(1:npti) > epsi20 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
     328      ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp 
    333329      END WHERE 
    334330      ! 
    … …  
    455451      ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.   
    456452      ! NOTE: 0 < aksum <= 1 
    457       WHERE( zaksum(1:npti) > 0._wp )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 
    458       ELSEWHERE                         ;   closing_gross(1:npti) = 0._wp 
     453      WHERE( zaksum(1:npti) > epsi20 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 
     454      ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp 
    459455      END WHERE 
    460456       
    … …  
    466462         DO ji = 1, npti 
    467463            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 
    468             IF( zfac > pa_i(ji,jl) ) THEN 
     464            IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN 
    469465               closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice 
    470466            ENDIF 
    … …  
    510506      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2 
    511507      REAL(wp), DIMENSION(jpij) ::   z1_ai            ! 1 / a 
     508      REAL(wp), DIMENSION(jpij) ::   zvti             ! sum(v_i) 
    512509      ! 
    513510      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrft     ! snow energy of rafting ice 
    … …  
    518515      INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation 
    519516      !!------------------------------------------------------------------- 
    520  
     517      ! 
     518      zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 )   ! total ice volume 
     519      ! 
    521520      ! 1) Change in open water area due to closing and opening 
    522521      !-------------------------------------------------------- 
    … …  
    535534            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging 
    536535 
    537                z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
    538  
     536               IF( a_i_2d(ji,jl1) > epsi20 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
     537               ELSE                                 ;   z1_ai(ji) = 0._wp 
     538               ENDIF 
     539                
    539540               ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2) 
    540541               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice 
    … …  
    549550 
    550551               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges 
    551                vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg 
     552               IF    ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg                                           ! v <= 10m then porosity = rn_porordg 
     553               ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp                                                                         ! v >= 20m then porosity = 0 
     554               ELSE                           ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg * MAX( 0._wp, 2._wp - 0.1_wp * zvti(ji) ) ! v > 10m and v < 20m then porosity = linear transition to 0 
     555               ENDIF 
    552556               ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?) 
    553557 
    554558               ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei) 
    555559               virdg1     = v_i_2d (ji,jl1)   * afrdg 
    556                virdg2(ji) = v_i_2d (ji,jl1)   * afrdg * ( 1. + rn_porordg ) 
     560               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg + vsw 
    557561               vsrdg(ji)  = v_s_2d (ji,jl1)   * afrdg 
    558562               sirdg1     = sv_i_2d(ji,jl1)   * afrdg 
    … …  
    726730      END DO ! jl1 
    727731      ! 
     732      ! roundoff errors 
     733      !---------------- 
    728734      ! In case ridging/rafting lead to very small negative values (sometimes it happens) 
    729       WHERE( a_i_2d(1:npti,:) < 0._wp )   a_i_2d(1:npti,:) = 0._wp 
    730       WHERE( v_i_2d(1:npti,:) < 0._wp )   v_i_2d(1:npti,:) = 0._wp 
     735      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 
    731736      ! 
    732737   END SUBROUTINE rdgrft_shift 
    … …  
    854859         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
    855860         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) ) 
    856  
     861         ! 
    857862         !                 !---------------------! 
    858863      CASE( 2 )            !==  from 1D to 2D  ==! 
    … …  
    945950         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) 
    946951      ENDIF 
    947       !                              ! allocate tke arrays 
     952      ! 
     953      IF( .NOT. ln_icethd ) THEN 
     954         rn_porordg = 0._wp 
     955         rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp 
     956         rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp 
     957         IF( lwp ) THEN 
     958            WRITE(numout,*) '      ==> only ice dynamics is activated, thus some parameters must be changed' 
     959            WRITE(numout,*) '            rn_porordg   = ', rn_porordg 
     960            WRITE(numout,*) '            rn_fsnwrdg   = ', rn_fsnwrdg  
     961            WRITE(numout,*) '            rn_fpndrdg   = ', rn_fpndrdg  
     962            WRITE(numout,*) '            rn_fsnwrft   = ', rn_fsnwrft  
     963            WRITE(numout,*) '            rn_fpndrft   = ', rn_fpndrft  
     964         ENDIF 
     965      ENDIF 
     966      !                              ! allocate arrays 
    948967      IF( ice_dyn_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' ) 
    949968      ! 

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/icedyn_rhg.F90

    r10888 r11084  
    5858      !!-------------------------------------------------------------------- 
    5959      INTEGER, INTENT(in) ::   kt     ! ice time step 
     60      ! 
     61      INTEGER  ::   jl   ! dummy loop indices 
    6062      !!-------------------------------------------------------------------- 
    6163      ! controls 
    … …  
    6870         WRITE(numout,*)'~~~~~~~~~~~' 
    6971      ENDIF 
    70  
    71       ! -------- 
    72       ! Rheology 
    73       ! --------    
     72      ! 
     73      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization 
     74         tau_icebfr(:,:) = 0._wp 
     75         DO jl = 1, jpl 
     76            WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
     77         END DO 
     78      ENDIF 
     79      ! 
     80      !--------------! 
     81      !== Rheology ==! 
     82      !--------------!    
    7483      SELECT CASE( nice_rhg ) 
    7584      !                                !------------------------! 

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/icedyn_rhg_evp.F90

    r10888 r11084  
    112112      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      ! 
    113113      !! 
    114       LOGICAL, PARAMETER ::   ll_bdy_substep = .FALSE. ! temporary option to call bdy at each sub-time step (T) 
     114      LOGICAL, PARAMETER ::   ll_bdy_substep = .TRUE. ! temporary option to call bdy at each sub-time step (T) 
    115115      !                                                                              or only at the main time step (F) 
    116116      INTEGER ::   ji, jj       ! dummy loop indices 
    … …  
    160160 
    161161      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter 
    162       REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity 
     162      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small 
     163      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small 
    163164      !! --- diags 
    164165      REAL(wp), DIMENSION(jpi,jpj) ::   zswi 
    … …  
    319320 
    320321            ! switches 
    321             zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin 
    322             zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin 
     322            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zswitchU(ji,jj) = 0._wp 
     323            ELSE                                                   ;   zswitchU(ji,jj) = 1._wp   ;   ENDIF 
     324            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zswitchV(ji,jj) = 0._wp 
     325            ELSE                                                   ;   zswitchV(ji,jj) = 1._wp   ;   ENDIF 
    323326 
    324327         END DO 
    … …  
    519522                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    520523                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    521                      &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin 
     524                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    522525                     &           ) * zmaskV(ji,jj) 
    523526                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    … …  
    526529                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
    527530                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    528                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
     531                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    529532                     &            ) * zmaskV(ji,jj) 
    530533                  ENDIF 
    … …  
    567570                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    568571                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0 
    569                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin  
     572                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    570573                     &            ) * zmaskU(ji,jj) 
    571574                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    … …  
    574577                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
    575578                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    576                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin  
     579                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    577580                     &            ) * zmaskU(ji,jj) 
    578581                  ENDIF 
    … …  
    617620                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    618621                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0 
    619                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin  
     622                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    620623                     &            ) * zmaskU(ji,jj) 
    621624                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    … …  
    624627                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
    625628                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    626                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin  
     629                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    627630                     &            ) * zmaskU(ji,jj) 
    628631                  ENDIF 
    … …  
    665668                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    666669                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    667                      &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin 
     670                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    668671                     &           ) * zmaskV(ji,jj) 
    669672                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    … …  
    672675                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
    673676                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    674                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
     677                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    675678                     &            ) * zmaskV(ji,jj) 
    676679                  ENDIF 

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/iceitd.F90

    r10888 r11084  
    2121   USE ice1D          ! sea-ice: thermodynamic variables 
    2222   USE ice            ! sea-ice: variables 
     23   USE icevar         ! sea-ice: operations 
    2324   USE icectl         ! sea-ice: conservation tests 
    2425   USE icetab         ! sea-ice: convert 1D<=>2D 
    … …  
    9192      !  1) Identify grid cells with ice 
    9293      !----------------------------------------------------------------------------------------------- 
     94      at_i(:,:) = SUM( a_i, dim=3 ) 
     95      ! 
    9396      npti = 0   ;   nptidx(:) = 0 
    9497      DO jj = 1, jpj 
    … …  
    249252            ! --- g(h) for each thickness category --- !   
    250253            CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti)   , a_i_1d(1:npti)   ,  &   ! in 
    251                &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL     (1:npti,jl), hR   (1:npti,jl)   )   ! out 
     254               &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL    (1:npti,jl), hR    (1:npti,jl)   )   ! out 
    252255            ! 
    253256         END DO 
    … …  
    389392      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdvice   ! ice volume transferred across boundary 
    390393      ! 
    391       INTEGER  ::   ji, jj, jl, jk     ! dummy loop indices 
    392       INTEGER  ::   ii, ij, jl2, jl1   ! local integers 
     394      INTEGER  ::   ji, jl, jk         ! dummy loop indices 
     395      INTEGER  ::   jl2, jl1           ! local integers 
    393396      REAL(wp) ::   ztrans             ! ice/snow transferred 
    394       REAL(wp), DIMENSION(jpij)     ::   zworka, zworkv   ! workspace 
    395       REAL(wp), DIMENSION(jpij,jpl) ::   zaTsfn           !  -    - 
     397      REAL(wp), DIMENSION(jpij)            ::   zworka, zworkv   ! workspace 
     398      REAL(wp), DIMENSION(jpij,jpl)        ::   zaTsfn           !  -    - 
     399      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_2d 
     400      REAL(wp), DIMENSION(jpij,nlay_s,jpl) ::   ze_s_2d 
    396401      !!------------------------------------------------------------------ 
    397402          
    … …  
    405410      CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 
    406411      CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 
     412      DO jl = 1, jpl 
     413         DO jk = 1, nlay_s 
     414            CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 
     415         END DO 
     416         DO jk = 1, nlay_i 
     417            CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 
     418         END DO 
     419      END DO 
     420      ! to correct roundoff errors on a_i 
     421      CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti), rn_amax_2d ) 
    407422 
    408423      !---------------------------------------------------------------------------------------------- 
    … …  
    435450               ELSE                                  ;   zworka(ji) = 0._wp 
    436451               ENDIF 
    437                ! 
    438                ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 
    439                !       because of truncation error ( i.e. 1. - 1. /= 0 ) 
    440                !       I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90) 
    441452               ! 
    442453               a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas 
    … …  
    476487         ! 
    477488         DO jk = 1, nlay_s         !--- Snow heat content 
    478             ! 
    479489            DO ji = 1, npti 
    480                ii = MOD( nptidx(ji) - 1, jpi ) + 1 
    481                ij = ( nptidx(ji) - 1 ) / jpi + 1 
    482490               ! 
    483491               jl1 = kdonor(ji,jl) 
    … …  
    487495                  ELSE                ;  jl2 = jl 
    488496                  ENDIF 
    489                   ! 
    490                   ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji) 
    491                   e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans 
    492                   e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans 
     497                  ztrans             = ze_s_2d(ji,jk,jl1) * zworkv(ji) 
     498                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) - ztrans 
     499                  ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ztrans 
    493500               ENDIF 
    494501            END DO 
    … …  
    497504         DO jk = 1, nlay_i         !--- Ice heat content 
    498505            DO ji = 1, npti 
    499                ii = MOD( nptidx(ji) - 1, jpi ) + 1 
    500                ij = ( nptidx(ji) - 1 ) / jpi + 1 
    501506               ! 
    502507               jl1 = kdonor(ji,jl) 
    … …  
    506511                  ELSE                ;  jl2 = jl 
    507512                  ENDIF 
    508                   ! 
    509                   ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji) 
    510                   e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans 
    511                   e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans 
     513                  ztrans             = ze_i_2d(ji,jk,jl1) * zworkv(ji) 
     514                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) - ztrans 
     515                  ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + ztrans 
    512516               ENDIF 
    513517            END DO 
    … …  
    515519         ! 
    516520      END DO                   ! boundaries, 1 to jpl-1 
     521 
     522      !------------------- 
     523      ! 3) roundoff errors 
     524      !------------------- 
     525      ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 
     526      !       because of truncation error ( i.e. 1. - 1. /= 0 ) 
     527      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 
     528 
     529      ! at_i must be <= rn_amax 
     530      zworka(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 ) 
     531      DO jl  = 1, jpl 
     532         WHERE( zworka(1:npti) > rn_amax_1d(1:npti) )   & 
     533            &   a_i_2d(1:npti,jl) = a_i_2d(1:npti,jl) * rn_amax_1d(1:npti) / zworka(1:npti) 
     534      END DO 
    517535       
    518536      !------------------------------------------------------------------------------- 
    519       ! 3) Update ice thickness and temperature 
     537      ! 4) Update ice thickness and temperature 
    520538      !------------------------------------------------------------------------------- 
    521539      WHERE( a_i_2d(1:npti,:) >= epsi20 ) 
    … …  
    536554      CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 
    537555      CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 
     556      DO jl = 1, jpl 
     557         DO jk = 1, nlay_s 
     558            CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 
     559         END DO 
     560         DO jk = 1, nlay_i 
     561            CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 
     562         END DO 
     563      END DO 
    538564      ! 
    539565   END SUBROUTINE itd_shiftice 
    … …  
    558584      ! 
    559585      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution'  
     586      ! 
     587      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    560588      ! 
    561589      jdonor(:,:) = 0 
    … …  
    635663      END DO 
    636664      ! 
     665      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     666      ! 
    637667   END SUBROUTINE ice_itd_reb 
    638668 

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/icestp.F90

    r10888 r11084  
    189189         IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics       
    190190         ! 
    191          IF( ln_icethd )                CALL ice_cor( kt , 2 )        ! -- Corrections 
     191                                        CALL ice_cor( kt , 2 )        ! -- Corrections 
    192192         ! 
    193193                                        CALL ice_var_glo2eqv          ! necessary calls (at least for coupling) 
    … …  
    323323         WRITE(numout,*) '         maximum ice concentration for SH                              = ', rn_amax_s 
    324324      ENDIF 
     325      !                                        !--- change max ice concentration for roundoff errors 
     326      rn_amax_n = MIN( rn_amax_n, 1._wp - epsi10 ) 
     327      rn_amax_s = MIN( rn_amax_s, 1._wp - epsi10 ) 
    325328      !                                        !--- check consistency 
    326329      IF ( jpl > 1 .AND. ln_virtual_itd ) THEN 
    … …  
    431434      t_si       (:,:,:) = rt0   ! temp at the ice-snow interface 
    432435 
    433       tau_icebfr(:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here) 
    434       cnd_ice   (:,:,:) = 0._wp   ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T) 
    435       qtr_ice_bot(:,:,:) = 0._wp  ! initialization: part of solar radiation transmitted through the ice needed at least for outputs 
     436      tau_icebfr (:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here) 
     437      cnd_ice    (:,:,:) = 0._wp   ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T) 
     438      qcn_ice    (:,:,:) = 0._wp   ! initialisation: conductive flux (ln_cndflx=T & ln_cndemule=T) 
     439      qtr_ice_bot(:,:,:) = 0._wp   ! initialization: part of solar radiation transmitted through the ice needed at least for outputs 
     440      qsb_ice_bot(:,:)   = 0._wp   ! (needed if ln_icethd=F) 
    436441      ! 
    437442      ! for control checks (ln_icediachk) 

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/icethd.F90

    r10888 r11084  
    102102      ENDIF 
    103103       
    104       CALL ice_var_glo2eqv 
    105  
    106104      !---------------------------------------------! 
    107105      ! computation of friction velocity at T points 
    … …  
    162160            qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr ) 
    163161 
    164             ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting  
    165             IF( zqld > 0._wp ) THEN 
     162            ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting  
     163            ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget 
     164            IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN 
    166165               fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 
    167166               qlead(ji,jj) = 0._wp 
    … …  
    178177      ! In case we bypass open-water ice formation 
    179178      IF( .NOT. ln_icedO )  qlead(:,:) = 0._wp 
    180       ! In case we bypass growing/melting from top and bottom: we suppose ice is impermeable => ocean is isolated from atmosphere 
     179      ! In case we bypass growing/melting from top and bottom 
    181180      IF( .NOT. ln_icedH ) THEN 
    182          qt_atm_oi  (:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:) 
    183181         qsb_ice_bot(:,:) = 0._wp 
    184182         fhld       (:,:) = 0._wp 
    … …  
    221219            dh_i_sub  (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp 
    222220            dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp 
    223             ! 
    224             IF( ln_icedH ) THEN                                     ! --- growing/melting --- ! 
    225                               CALL ice_thd_zdf                             ! Ice/Snow Temperature profile 
    226                               CALL ice_thd_dh                              ! Ice/Snow thickness    
    227                               CALL ice_thd_pnd                             ! Melt ponds formation 
    228                               CALL ice_thd_ent( e_i_1d(1:npti,:) )         ! Ice enthalpy remapping 
     221            !                                       
     222                              CALL ice_thd_zdf                      ! --- Ice-Snow temperature --- ! 
     223            ! 
     224            IF( ln_icedH ) THEN                                     ! --- Growing/Melting --- ! 
     225                              CALL ice_thd_dh                           ! Ice-Snow thickness    
     226                              CALL ice_thd_pnd                          ! Melt ponds formation 
     227                              CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping 
    229228            ENDIF 
    230             ! 
    231229                              CALL ice_thd_sal( ln_icedS )          ! --- Ice salinity --- !     
    232230            ! 
    233                               CALL ice_thd_temp                     ! --- temperature update --- ! 
     231                              CALL ice_thd_temp                     ! --- Temperature update --- ! 
    234232            ! 
    235233            IF( ln_icedH .AND. ln_virtual_itd ) & 
    236                &              CALL ice_thd_mono                     ! --- extra lateral melting if virtual_itd --- ! 
    237             ! 
    238             IF( ln_icedA )    CALL ice_thd_da                       ! --- lateral melting --- ! 
     234               &              CALL ice_thd_mono                     ! --- Extra lateral melting if virtual_itd --- ! 
     235            ! 
     236            IF( ln_icedA )    CALL ice_thd_da                       ! --- Lateral melting --- ! 
    239237            ! 
    240238                              CALL ice_thd_1d2d( jl, 2 )            ! --- Change units of e_i, e_s from J/m3 to J/m2 --- ! 
    241239            !                                                       ! --- & Move to 2D arrays --- ! 
    242             ! 
    243240         ENDIF 
    244241         ! 
    245242      END DO 
    246  
     243      ! 
    247244      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    248       ! 
    249                            CALL ice_var_zapsmall           ! --- remove very small ice concentration (<1e-10) --- ! 
    250       !                                                    !     & make sure at_i=SUM(a_i) & ato_i=1 where at_i=0 
    251245      !                    
    252       IF( jpl > 1      )   CALL ice_itd_rem( kt )          ! --- Transport ice between thickness categories --- ! 
    253       ! 
    254       IF( ln_icedO     )   CALL ice_thd_do                 ! --- frazil ice growing in leads --- ! 
     246      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- ! 
     247      ! 
     248      IF( ln_icedO )          CALL ice_thd_do                       ! --- Frazil ice growth in leads --- ! 
    255249      ! 
    256250      ! controls 

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/icethd_do.F90

    r10888 r11084  
    114114      IF( ln_icediachk )   CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft ) 
    115115 
    116       CALL ice_var_agg(1) 
    117       CALL ice_var_glo2eqv 
    118  
     116      at_i(:,:) = SUM( a_i, dim=3 ) 
    119117      !------------------------------------------------------------------------------! 
    120118      ! 1) Collection thickness of ice formed in leads and polynyas 
    … …  
    319317 
    320318         ! --- lateral ice growth --- ! 
    321          ! If lateral ice growth gives an ice concentration gt 1, then 
     319         ! If lateral ice growth gives an ice concentration > amax, then 
    322320         ! we keep the excessive volume in memory and attribute it later to bottom accretion 
    323321         DO ji = 1, npti 
    324             IF ( za_newice(ji) >  ( rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN 
    325                zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - at_i_1d(ji) ) 
     322            IF ( za_newice(ji) >  MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN ! max is for roundoff error 
     323               zda_res(ji)   = za_newice(ji) - MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) 
    326324               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji)  
    327                za_newice(ji) = za_newice(ji) - zda_res  (ji) 
    328                zv_newice(ji) = zv_newice(ji) - zdv_res  (ji) 
     325               za_newice(ji) = MAX( 0._wp, za_newice(ji) - zda_res  (ji) ) 
     326               zv_newice(ji) = MAX( 0._wp, zv_newice(ji) - zdv_res  (ji) ) 
    329327            ELSE 
    330328               zda_res(ji) = 0._wp 

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/icethd_zdf_bl99.F90

    r10888 r11084  
    206206      ! 
    207207      l_T_converged(:) = .FALSE. 
    208       !                                                          !============================! 
    209208      ! Convergence calculated until all sub-domain grid points have converged 
    210209      ! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation) 
    211210      ! but values are not taken into account (results independant of MPI partitioning) 
    212211      ! 
     212      !                                                                            !============================! 
    213213      DO WHILE ( ( .NOT. ALL (l_T_converged(1:npti)) ) .AND. iconv < iconv_max )   ! Iterative procedure begins ! 
    214          !                                                       !============================! 
     214         !                                                                         !============================! 
    215215         iconv = iconv + 1 
    216216         ! 
    … …  
    742742                  zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 
    743743                  ! t_i 
    744                   DO jk = 0, nlay_i 
     744                  DO jk = 1, nlay_i 
    745745                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0  
    746746                     t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 
    … …  
    856856         t_i_1d    (1:npti,:) = ztiold        (1:npti,:) 
    857857         qcn_ice_1d(1:npti)   = qcn_ice_top_1d(1:npti) 
     858 
     859         !!clem 
     860         ! remettre t_su_1d, qns_ice_1d et dqns_ice_1d comme avant puisqu'on devrait faire comme si on avant conduction = input 
     861         !clem 
    858862      ENDIF 
    859863      ! 

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/icevar.F90

    r10888 r11084  
    4444   !!   ice_var_salprof1d : salinity profile in the ice 1D 
    4545   !!   ice_var_zapsmall  : remove very small area and volume 
    46    !!   ice_var_zapneg    : remove negative ice fields (to debug the advection scheme UM3-5) 
     46   !!   ice_var_zapneg    : remove negative ice fields 
     47   !!   ice_var_roundoff  : remove negative values arising from roundoff erros 
    4748   !!   ice_var_itd       : convert 1-cat to jpl-cat 
    4849   !!   ice_var_itd2      : convert N-cat to jpl-cat 
    … …  
    7172   PUBLIC   ice_var_zapsmall 
    7273   PUBLIC   ice_var_zapneg 
     74   PUBLIC   ice_var_roundoff 
    7375   PUBLIC   ice_var_itd 
    7476   PUBLIC   ice_var_itd2 
    … …  
    229231                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area  
    230232                     ! 
    231                      ze_i             =   e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3] 
     233                     ze_i             =   e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i             ! Energy of melting e(S,T) [J.m-3] 
    232234                     ztmelts          = - sz_i(ji,jj,jk,jl) * rTmlt                                 ! Ice layer melt temperature [C] 
    233235                     ! Conversion q(S,T) -> T (second order equation) 
    … …  
    236238                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts 
    237239                     ! 
    238                   ELSE                                !--- no ice 
     240                  ELSE                                   !--- no ice 
    239241                     t_i(ji,jj,jk,jl) = rt0 
    240242                  ENDIF 
    … …  
    537539 
    538540 
    539    SUBROUTINE ice_var_zapneg( pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     541   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    540542      !!------------------------------------------------------------------- 
    541543      !!                   ***  ROUTINE ice_var_zapneg *** 
    … …  
    543545      !! ** Purpose :   Remove negative sea ice fields and correct fluxes 
    544546      !!------------------------------------------------------------------- 
    545       INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices 
    546       ! 
     547      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step 
    547548      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area 
    548549      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume 
    … …  
    555556      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    556557      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
    557       !!------------------------------------------------------------------- 
    558       ! 
     558      ! 
     559      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices 
     560      REAL(wp) ::   z1_dt 
     561      !!------------------------------------------------------------------- 
     562      ! 
     563      z1_dt = 1._wp / pdt 
    559564      ! 
    560565      DO jl = 1, jpl       !==  loop over the categories  ==! 
    561566         ! 
     567         ! make sure a_i=0 where v_i<=0 
     568         WHERE( pv_i(:,:,:) <= 0._wp )   pa_i(:,:,:) = 0._wp 
     569 
    562570         !---------------------------------------- 
    563571         ! zap ice energy and send it to the ocean 
    … …  
    566574            DO jj = 1 , jpj 
    567575               DO ji = 1 , jpi 
    568                   IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
    569                      hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 >0 
     576                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
     577                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0 
    570578                     pe_i(ji,jj,jk,jl) = 0._wp 
    571579                  ENDIF 
    … …  
    577585            DO jj = 1 , jpj 
    578586               DO ji = 1 , jpi 
    579                   IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
    580                      hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0 
     587                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
     588                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0 
    581589                     pe_s(ji,jj,jk,jl) = 0._wp 
    582590                  ENDIF 
    … …  
    590598         DO jj = 1 , jpj 
    591599            DO ji = 1 , jpi 
    592                IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
    593                   wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * r1_rdtice 
     600               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
     601                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt 
    594602                  pv_i   (ji,jj,jl) = 0._wp 
    595603               ENDIF 
    596                IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
    597                   wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * r1_rdtice 
     604               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
     605                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt 
    598606                  pv_s   (ji,jj,jl) = 0._wp 
    599607               ENDIF 
    600                IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN 
    601                   sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * r1_rdtice 
     608               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
     609                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt 
    602610                  psv_i  (ji,jj,jl) = 0._wp 
    603611               ENDIF 
    … …  
    616624   END SUBROUTINE ice_var_zapneg 
    617625 
     626 
     627   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     628      !!------------------------------------------------------------------- 
     629      !!                   ***  ROUTINE ice_var_roundoff *** 
     630      !! 
     631      !! ** Purpose :   Remove negative sea ice values arising from roundoff errors 
     632      !!------------------------------------------------------------------- 
     633      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_i       ! ice concentration 
     634      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_i       ! ice volume 
     635      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_s       ! snw volume 
     636      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   psv_i      ! salt content 
     637      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   poa_i      ! age content 
     638      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
     639      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     640      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
     641      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     642      !!------------------------------------------------------------------- 
     643      ! 
     644      WHERE( pa_i (1:npti,:)   < 0._wp .AND. pa_i (1:npti,:)   > -epsi10 )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0 
     645      WHERE( pv_i (1:npti,:)   < 0._wp .AND. pv_i (1:npti,:)   > -epsi10 )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0 
     646      WHERE( pv_s (1:npti,:)   < 0._wp .AND. pv_s (1:npti,:)   > -epsi10 )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0 
     647      WHERE( psv_i(1:npti,:)   < 0._wp .AND. psv_i(1:npti,:)   > -epsi10 )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0 
     648      WHERE( poa_i(1:npti,:)   < 0._wp .AND. poa_i(1:npti,:)   > -epsi10 )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0 
     649      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0 
     650      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0 
     651      IF ( ln_pnd_H12 ) THEN 
     652         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
     653         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 
     654      ENDIF 
     655      ! 
     656   END SUBROUTINE ice_var_roundoff 
     657    
    618658    
    619659   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i ) 
    … …  
    650690      INTEGER  :: idim, i_fill, jl0   
    651691      REAL(wp) :: zarg, zV, zconv, zdh, zdv 
    652       REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables 
    653       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i ! output ice/snow variables 
     692      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input  ice/snow variables 
     693      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables 
    654694      INTEGER , DIMENSION(4)                  ::   itest 
    655695      !!------------------------------------------------------------------- 
    … …  
    780820      !! 
    781821      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1 
    782       !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
     822       !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
    783823      !!               
    784824      !!               3) Expand the filling to the empty cat between jlmin and jlmax  

  • NEMO/branches/UKMO/NEMO_4.0_old_tidal_mixing/src/ICE/icewri.F90

    r10888 r11084  
    145145 
    146146        ! record presence of fast ice 
    147          WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk00(:,:) == 1._wp ) ; zfast(:,:) = 1._wp 
     147         WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk15(:,:) == 1._wp ) ; zfast(:,:) = 1._wp 
    148148         ELSEWHERE                                                ; zfast(:,:) = 0._wp 
    149149         END WHERE 
Note: See TracChangeset for help on using the changeset viewer.