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/CANAL/MY_SRC/domvvl.F90

    r12489 r13197  
    3737 
    3838   PUBLIC  dom_vvl_init       ! called by domain.F90 
     39   PUBLIC  dom_vvl_zgr        ! called by isfcpl.F90 
    3940   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90 
    4041   PUBLIC  dom_vvl_sf_update  ! called by step.F90 
     
    6263   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence 
    6364 
     65   !! * Substitutions 
     66#  include "do_loop_substitute.h90" 
    6467   !!---------------------------------------------------------------------- 
    6568   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    116119      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
    117120      ! 
     121      IF(lwp) WRITE(numout,*) 
     122      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 
     123      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     124      ! 
     125      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer) 
     126      ! 
     127      !                    ! Allocate module arrays 
     128      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 
     129      ! 
     130      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
     131      CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 
     132      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
     133      ! 
     134      CALL dom_vvl_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column 
     135      ! 
     136   END SUBROUTINE dom_vvl_init 
     137   ! 
     138   SUBROUTINE dom_vvl_zgr(Kbb, Kmm, Kaa) 
     139      !!---------------------------------------------------------------------- 
     140      !!                ***  ROUTINE dom_vvl_init  *** 
     141      !!                    
     142      !! ** Purpose :  Interpolation of all scale factors,  
     143      !!               depths and water column heights 
     144      !! 
     145      !! ** Method  :  - interpolate scale factors 
     146      !! 
     147      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b) 
     148      !!              - Regrid: e3(u/v)_n 
     149      !!                        e3(u/v)_b        
     150      !!                        e3w_n            
     151      !!                        e3(u/v)w_b       
     152      !!                        e3(u/v)w_n       
     153      !!                        gdept_n, gdepw_n and gde3w_n 
     154      !!              - h(t/u/v)_0 
     155      !!              - frq_rst_e3t and frq_rst_hdv 
     156      !! 
     157      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling. 
     158      !!---------------------------------------------------------------------- 
     159      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
     160      !!---------------------------------------------------------------------- 
    118161      INTEGER ::   ji, jj, jk 
    119162      INTEGER ::   ii0, ii1, ij0, ij1 
    120163      REAL(wp)::   zcoef 
    121164      !!---------------------------------------------------------------------- 
    122       ! 
    123       IF(lwp) WRITE(numout,*) 
    124       IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 
    125       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    126       ! 
    127       CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer) 
    128       ! 
    129       !                    ! Allocate module arrays 
    130       IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 
    131       ! 
    132       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
    133       CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 
    134       e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
    135165      ! 
    136166      !                    !== Set of all other vertical scale factors  ==!  (now and before) 
     
    160190      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 
    161191      gdepw(:,:,1,Kbb) = 0.0_wp 
    162       DO jk = 2, jpk                               ! vertical sum 
    163          DO jj = 1,jpj 
    164             DO ji = 1,jpi 
    165                !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    166                !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
    167                !                             ! 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      
    168196!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ?? 
    169                zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
    170                gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
    171                gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
    172                   &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))  
    173                gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    174                gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 
    175                gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  & 
    176                   &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))  
    177             END DO 
    178          END DO 
    179       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 
    180206      ! 
    181207      !                    !==  thickness of the water column  !!   (ocean portion only) 
     
    212238         ENDIF 
    213239         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator 
    214             DO jj = 1, jpj 
    215                DO ji = 1, jpi 
     240            DO_2D_11_11 
    216241!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default 
    217                   IF( ABS(gphit(ji,jj)) >= 6.) THEN 
    218                      ! values outside the equatorial band and transition zone (ztilde) 
    219                      frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp ) 
    220                      frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
    221                   ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star 
    222                      ! values inside the equatorial band (ztilde as zstar) 
    223                      frq_rst_e3t(ji,jj) =  0.0_wp 
    224                      frq_rst_hdv(ji,jj) =  1.0_wp / rn_Dt 
    225                   ELSE                                      ! transition band (2.5 to 6 degrees N/S) 
    226                      !                                      ! (linearly transition from z-tilde to z-star) 
    227                      frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   & 
    228                         &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    229                         &                                          * 180._wp / 3.5_wp ) ) 
    230                      frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt)                                & 
    231                         &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp   & 
    232                         &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    233                         &                                          * 180._wp / 3.5_wp ) ) 
    234                   ENDIF 
    235                END DO 
    236             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 
    237261            IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN 
    238262               IF( nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 
     
    264288      ENDIF 
    265289      ! 
    266    END SUBROUTINE dom_vvl_init 
     290   END SUBROUTINE dom_vvl_zgr 
    267291 
    268292 
     
    329353      END DO 
    330354      ! 
    331       IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
    332          !                                                            ! ------baroclinic part------ ! 
     355      IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
     356         !                                                               ! ------baroclinic part------ ! 
    333357         ! I - initialization 
    334358         ! ================== 
     
    383407         zwu(:,:) = 0._wp 
    384408         zwv(:,:) = 0._wp 
    385          DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes 
    386             DO jj = 1, jpjm1 
    387                DO ji = 1, fs_jpim1   ! vector opt. 
    388                   un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
    389                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
    390                   vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
    391                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    392                   zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    393                   zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
    394                END DO 
    395             END DO 
    396          END DO 
    397          DO jj = 1, jpj          ! b - correction for last oceanic u-v points 
    398             DO ji = 1, jpi 
    399                un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
    400                vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
    401             END DO 
    402          END DO 
    403          DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes 
    404             DO jj = 2, jpjm1 
    405                DO ji = fs_2, fs_jpim1   ! vector opt. 
    406                   tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
    407                      &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
    408                      &                                            ) * r1_e1e2t(ji,jj) 
    409                END DO 
    410             END DO 
    411          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 
    412426         !                       ! d - thickness diffusion transport: boundary conditions 
    413427         !                             (stored for tracer advction and continuity equation) 
     
    416430         ! 4 - Time stepping of baroclinic scale factors 
    417431         ! --------------------------------------------- 
    418          ! Leapfrog time stepping 
    419          ! ~~~~~~~~~~~~~~~~~~~~~~ 
    420432         CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp ) 
    421433         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 
     
    613625         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 
    614626      ENDIF 
    615       gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 
    616       gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 
    617  
    618       e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa) 
    619       e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa) 
    620       e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa) 
    621627 
    622628      ! Compute all missing vertical scale factor and depths 
     
    641647      gdepw(:,:,1,Kmm) = 0.0_wp 
    642648      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 
    643       DO jk = 2, jpk 
    644          DO jj = 1,jpj 
    645             DO ji = 1,jpi 
    646               !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    647                                                                  ! 1 for jk = mikt 
    648                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    649                gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
    650                gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  & 
    651                    &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )  
    652                gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    653             END DO 
    654          END DO 
    655       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 
    656658 
    657659      ! Local depth and Inverse of the local depth of the water 
     
    700702         ! 
    701703      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean 
    702          DO jk = 1, jpk 
    703             DO jj = 1, jpjm1 
    704                DO ji = 1, fs_jpim1   ! vector opt. 
    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 DO 
    709             END DO 
    710          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 
    711709         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 
    712710         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 
    713711         ! 
    714712      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean 
    715          DO jk = 1, jpk 
    716             DO jj = 1, jpjm1 
    717                DO ji = 1, fs_jpim1   ! vector opt. 
    718                   pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
    719                      &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    720                      &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
    721                END DO 
    722             END DO 
    723          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 
    724718         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 
    725719         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 
    726720         ! 
    727721      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean 
    728          DO jk = 1, jpk 
    729             DO jj = 1, jpjm1 
    730                DO ji = 1, fs_jpim1   ! vector opt. 
    731                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
    732                      &                       *    r1_e1e2f(ji,jj)                                                  & 
    733                      &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    734                      &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
    735                END DO 
    736             END DO 
    737          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 
    738728         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 
    739729         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 
     
    810800            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 
    811801            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 
     802            ! 
    812803            !                             ! --------- ! 
    813804            !                             ! all cases ! 
    814805            !                             ! --------- ! 
     806            ! 
    815807            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
    816808               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
     
    828820               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 
    829821               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 
    830                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
     822               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    831823               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
    832824               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
     
    835827               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 
    836828               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 
    837                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
     829               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    838830               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
    839831               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
     
    842834               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 
    843835               IF(lwp) write(numout,*) 'Compute scale factor from sshn' 
    844                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
     836               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    845837               DO jk = 1, jpk 
    846838                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 
     
    895887                  ssh(:,:,Kbb) = -ssh_ref 
    896888 
    897                   DO jj = 1, jpj 
    898                      DO ji = 1, jpi 
    899                         IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
    900                            ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
    901                            ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
    902                         ENDIF 
    903                      ENDDO 
    904                   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 
    905895               ENDIF !If test case else 
    906896 
     
    913903               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    914904 
    915                DO ji = 1, jpi 
    916                   DO jj = 1, jpj 
    917                      IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
    918                        CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 
    919                      ENDIF 
    920                   END DO  
    921                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 
    922910               ! 
    923911            ELSE 
    924912               ! 
    925                ! usr_def_istate called here only to get sshb, that is needed to initialize e3t(Kbb) and e3t(Kmm) 
    926                CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
    927                ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn) 
     913               ! usr_def_istate called here only to get ssh(Kbb) needed to initialize e3t(Kbb) and e3t(Kmm) 
     914               ! 
     915               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )   
     916               ! 
     917               ! usr_def_istate will be called again in istate_init to initialize ts, ssh, u and v 
    928918               ! 
    929919               DO jk=1,jpk 
    930                   e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 
    931                     &                              / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    932                     &              + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t_b != 0 on land points 
     920                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 
     921                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
     922                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t(:,:,:,Kbb) != 0 on land points 
    933923               END DO 
    934924               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
    935                ssh(:,:  ,Kmm) = ssh(:,:  ,Kbb)   ! needed later for gde3w 
     925               ssh(:,:,Kmm) = ssh(:,:,Kbb)                                     ! needed later for gde3w 
    936926               ! 
    937927            END IF           ! end of ll_wd edits 
     
    10251015      ! 
    10261016      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 
    1027       IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 
    10281017      ! 
    10291018      IF(lwp) THEN                   ! Print the choice 
Note: See TracChangeset for help on using the changeset viewer.