Ignore:
Timestamp:
2020-07-01T16:09:00+02:00 (3 months ago)
Author:
gsamson
Message:

merge with trunk@r13136 with a more recent svn version; pass all SETTE tests; results identical to trunk@r13136; ticket #2419

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/VORTEX/MY_SRC/domvvl.F90

    r12489 r13197  
    6363   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence 
    6464 
     65   !! * Substitutions 
     66#  include "do_loop_substitute.h90" 
    6567   !!---------------------------------------------------------------------- 
    6668   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    188190      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 
    189191      gdepw(:,:,1,Kbb) = 0.0_wp 
    190       DO jk = 2, jpk                               ! vertical sum 
    191          DO jj = 1,jpj 
    192             DO ji = 1,jpi 
    193                !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    194                !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
    195                !                             ! 0.5 where jk = mikt      
     192      DO_3D_11_11( 2, jpk ) 
     193         !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     194         !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
     195         !                             ! 0.5 where jk = mikt      
    196196!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ?? 
    197                zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
    198                gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
    199                gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
    200                   &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))  
    201                gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    202                gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 
    203                gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  & 
    204                   &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))  
    205             END DO 
    206          END DO 
    207       END DO 
     197         zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
     198         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     199         gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
     200            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))  
     201         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
     202         gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 
     203         gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  & 
     204            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))  
     205      END_3D 
    208206      ! 
    209207      !                    !==  thickness of the water column  !!   (ocean portion only) 
     
    240238         ENDIF 
    241239         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator 
    242             DO jj = 1, jpj 
    243                DO ji = 1, jpi 
     240            DO_2D_11_11 
    244241!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default 
    245                   IF( ABS(gphit(ji,jj)) >= 6.) THEN 
    246                      ! values outside the equatorial band and transition zone (ztilde) 
    247                      frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp ) 
    248                      frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
    249                   ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star 
    250                      ! values inside the equatorial band (ztilde as zstar) 
    251                      frq_rst_e3t(ji,jj) =  0.0_wp 
    252                      frq_rst_hdv(ji,jj) =  1.0_wp / rn_Dt 
    253                   ELSE                                      ! transition band (2.5 to 6 degrees N/S) 
    254                      !                                      ! (linearly transition from z-tilde to z-star) 
    255                      frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   & 
    256                         &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    257                         &                                          * 180._wp / 3.5_wp ) ) 
    258                      frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt)                                & 
    259                         &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp   & 
    260                         &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    261                         &                                          * 180._wp / 3.5_wp ) ) 
    262                   ENDIF 
    263                END DO 
    264             END DO 
     242               IF( ABS(gphit(ji,jj)) >= 6.) THEN 
     243                  ! values outside the equatorial band and transition zone (ztilde) 
     244                  frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp ) 
     245                  frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
     246               ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star 
     247                  ! values inside the equatorial band (ztilde as zstar) 
     248                  frq_rst_e3t(ji,jj) =  0.0_wp 
     249                  frq_rst_hdv(ji,jj) =  1.0_wp / rn_Dt 
     250               ELSE                                      ! transition band (2.5 to 6 degrees N/S) 
     251                  !                                      ! (linearly transition from z-tilde to z-star) 
     252                  frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   & 
     253                     &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
     254                     &                                          * 180._wp / 3.5_wp ) ) 
     255                  frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt)                                & 
     256                     &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp   & 
     257                     &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
     258                     &                                          * 180._wp / 3.5_wp ) ) 
     259               ENDIF 
     260            END_2D 
    265261            IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN 
    266262               IF( nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 
     
    357353      END DO 
    358354      ! 
    359       IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
    360          !                                                            ! ------baroclinic part------ ! 
     355      IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
     356         !                                                               ! ------baroclinic part------ ! 
    361357         ! I - initialization 
    362358         ! ================== 
     
    411407         zwu(:,:) = 0._wp 
    412408         zwv(:,:) = 0._wp 
    413          DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes 
    414             DO jj = 1, jpjm1 
    415                DO ji = 1, jpim1   ! vector opt. 
    416                   un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
    417                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
    418                   vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
    419                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    420                   zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    421                   zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
    422                END DO 
    423             END DO 
    424          END DO 
    425          DO jj = 1, jpj          ! b - correction for last oceanic u-v points 
    426             DO ji = 1, jpi 
    427                un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
    428                vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
    429             END DO 
    430          END DO 
    431          DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes 
    432             DO jj = 2, jpjm1 
    433                DO ji = 2, jpim1   ! vector opt. 
    434                   tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
    435                      &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
    436                      &                                            ) * r1_e1e2t(ji,jj) 
    437                END DO 
    438             END DO 
    439          END DO 
     409         DO_3D_10_10( 1, jpkm1 ) 
     410            un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
     411               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     412            vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
     413               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
     414            zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
     415            zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
     416         END_3D 
     417         DO_2D_11_11 
     418            un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
     419            vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
     420         END_2D 
     421         DO_3D_00_00( 1, jpkm1 ) 
     422            tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
     423               &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
     424               &                                            ) * r1_e1e2t(ji,jj) 
     425         END_3D 
    440426         !                       ! d - thickness diffusion transport: boundary conditions 
    441427         !                             (stored for tracer advction and continuity equation) 
     
    444430         ! 4 - Time stepping of baroclinic scale factors 
    445431         ! --------------------------------------------- 
    446          ! Leapfrog time stepping 
    447          ! ~~~~~~~~~~~~~~~~~~~~~~ 
    448432         CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp ) 
    449433         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 
     
    646630      ! Horizontal scale factor interpolations 
    647631      ! -------------------------------------- 
    648       ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are allready computed in dynnxt 
     632      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 
    649633      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 
    650634       
     
    663647      gdepw(:,:,1,Kmm) = 0.0_wp 
    664648      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 
    665       DO jk = 2, jpk 
    666          DO jj = 1,jpj 
    667             DO ji = 1,jpi 
    668               !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    669                                                                  ! 1 for jk = mikt 
    670                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    671                gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
    672                gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  & 
    673                    &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )  
    674                gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    675             END DO 
    676          END DO 
    677       END DO 
     649      DO_3D_11_11( 2, jpk ) 
     650        !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     651                                                           ! 1 for jk = mikt 
     652         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     653         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     654         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  & 
     655             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )  
     656         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
     657      END_3D 
    678658 
    679659      ! Local depth and Inverse of the local depth of the water 
     
    722702         ! 
    723703      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean 
    724          DO jk = 1, jpk 
    725             DO jj = 1, jpjm1 
    726                DO ji = 1, jpim1   ! vector opt. 
    727                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
    728                      &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    729                      &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
    730                END DO 
    731             END DO 
    732          END DO 
     704         DO_3D_10_10( 1, jpk ) 
     705            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
     706               &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
     707               &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     708         END_3D 
    733709         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 
    734710         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 
    735711         ! 
    736712      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean 
    737          DO jk = 1, jpk 
    738             DO jj = 1, jpjm1 
    739                DO ji = 1, jpim1   ! vector opt. 
    740                   pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
    741                      &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    742                      &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
    743                END DO 
    744             END DO 
    745          END DO 
     713         DO_3D_10_10( 1, jpk ) 
     714            pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
     715               &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
     716               &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
     717         END_3D 
    746718         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 
    747719         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 
    748720         ! 
    749721      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean 
    750          DO jk = 1, jpk 
    751             DO jj = 1, jpjm1 
    752                DO ji = 1, jpim1   ! vector opt. 
    753                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
    754                      &                       *    r1_e1e2f(ji,jj)                                                  & 
    755                      &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    756                      &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
    757                END DO 
    758             END DO 
    759          END DO 
     722         DO_3D_10_10( 1, jpk ) 
     723            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
     724               &                       *    r1_e1e2f(ji,jj)                                                  & 
     725               &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
     726               &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
     727         END_3D 
    760728         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 
    761729         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 
     
    832800            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 
    833801            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 
     802            ! 
    834803            !                             ! --------- ! 
    835804            !                             ! all cases ! 
    836805            !                             ! --------- ! 
     806            ! 
    837807            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
    838808               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
     
    850820               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 
    851821               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 
    852                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
     822               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    853823               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
    854824               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
     
    857827               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 
    858828               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 
    859                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
     829               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    860830               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
    861831               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
     
    864834               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 
    865835               IF(lwp) write(numout,*) 'Compute scale factor from sshn' 
    866                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
     836               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    867837               DO jk = 1, jpk 
    868838                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 
     
    917887                  ssh(:,:,Kbb) = -ssh_ref 
    918888 
    919                   DO jj = 1, jpj 
    920                      DO ji = 1, jpi 
    921                         IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
    922                            ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
    923                            ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
    924                         ENDIF 
    925                      ENDDO 
    926                   ENDDO 
     889                  DO_2D_11_11 
     890                     IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
     891                        ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
     892                        ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
     893                     ENDIF 
     894                  END_2D 
    927895               ENDIF !If test case else 
    928896 
     
    935903               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    936904 
    937                DO ji = 1, jpi 
    938                   DO jj = 1, jpj 
    939                      IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
    940                        CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 
    941                      ENDIF 
    942                   END DO  
    943                END DO  
     905               DO_2D_11_11 
     906                  IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
     907                     CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 
     908                  ENDIF 
     909               END_2D 
    944910               ! 
    945911            ELSE 
Note: See TracChangeset for help on using the changeset viewer.