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 12377 for NEMO/trunk/src/OCE/DYN/sshwzv.F90 – NEMO

Ignore:
Timestamp:
2020-02-12T15:39:06+01:00 (4 years ago)
Author:
acc
Message:

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk

    • Property svn:externals

      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 

  • NEMO/trunk/src/OCE/DYN/sshwzv.F90

    r11414 r12377  
    1010   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work 
    1111   !!            4.0  !  2018-12  (A. Coward) add mixed implicit/explicit advection 
     12   !!            4.1  !  2019-08  (A. Coward, D. Storkey) Rename ssh_nxt -> ssh_atf. Now only does time filtering. 
    1213   !!---------------------------------------------------------------------- 
    1314 
    1415   !!---------------------------------------------------------------------- 
    1516   !!   ssh_nxt       : after ssh 
    16    !!   ssh_swp       : filter ans swap the ssh arrays 
     17   !!   ssh_atf       : time filter the ssh arrays 
    1718   !!   wzv           : compute now vertical velocity 
    1819   !!---------------------------------------------------------------------- 
    1920   USE oce            ! ocean dynamics and tracers variables 
     21   USE isf_oce        ! ice shelf 
    2022   USE dom_oce        ! ocean space and time domain variables  
    2123   USE sbc_oce        ! surface boundary condition: ocean 
    … …  
    4446   PUBLIC   wzv        ! called by step.F90 
    4547   PUBLIC   wAimp      ! called by step.F90 
    46    PUBLIC   ssh_swp    ! called by step.F90 
     48   PUBLIC   ssh_atf    ! called by step.F90 
    4749 
    4850   !! * Substitutions 
    49 #  include "vectopt_loop_substitute.h90" 
     51#  include "do_loop_substitute.h90" 
    5052   !!---------------------------------------------------------------------- 
    5153   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    5557CONTAINS 
    5658 
    57    SUBROUTINE ssh_nxt( kt ) 
     59   SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa ) 
    5860      !!---------------------------------------------------------------------- 
    5961      !!                ***  ROUTINE ssh_nxt  *** 
    6062      !!                    
    61       !! ** Purpose :   compute the after ssh (ssha) 
     63      !! ** Purpose :   compute the after ssh (ssh(Kaa)) 
    6264      !! 
    6365      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment 
    … …  
    6567      !!      by the time step. 
    6668      !! 
    67       !! ** action  :   ssha, after sea surface height 
     69      !! ** action  :   ssh(:,:,Kaa), after sea surface height 
    6870      !! 
    6971      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    7072      !!---------------------------------------------------------------------- 
    71       INTEGER, INTENT(in) ::   kt   ! time step 
     73      INTEGER                         , INTENT(in   ) ::   kt             ! time step 
     74      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level index 
     75      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height 
    7276      !  
    7377      INTEGER  ::   jk            ! dummy loop indice 
    … …  
    9296      !                                           !------------------------------! 
    9397      IF(ln_wd_il) THEN 
    94          CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 
    95       ENDIF 
    96  
    97       CALL div_hor( kt )                               ! Horizontal divergence 
     98         CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), z2dt, Kmm, uu, vv ) 
     99      ENDIF 
     100 
     101      CALL div_hor( kt, Kbb, Kmm )                     ! Horizontal divergence 
    98102      ! 
    99103      zhdiv(:,:) = 0._wp 
    100104      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports 
    101         zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 
     105        zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk) 
    102106      END DO 
    103107      !                                                ! Sea surface elevation time stepping 
    … …  
    105109      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 
    106110      !  
    107       ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
     111      pssh(:,:,Kaa) = (  pssh(:,:,Kbb) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
    108112      ! 
    109113#if defined key_agrif 
    110       CALL agrif_ssh( kt ) 
     114      Kbb_a = Kbb; Kmm_a = Kmm; Krhs_a = Kaa; CALL agrif_ssh( kt ) 
    111115#endif 
    112116      ! 
    113117      IF ( .NOT.ln_dynspg_ts ) THEN 
    114118         IF( ln_bdy ) THEN 
    115             CALL lbc_lnk( 'sshwzv', ssha, 'T', 1. )    ! Not sure that's necessary 
    116             CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries 
     119            CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1. )    ! Not sure that's necessary 
     120            CALL bdy_ssh( pssh(:,:,Kaa) )             ! Duplicate sea level across open boundaries 
    117121         ENDIF 
    118122      ENDIF 
    … …  
    121125      !                                           !------------------------------! 
    122126      ! 
    123       IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask ) 
     127      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa)  - : ', mask1=tmask ) 
    124128      ! 
    125129      IF( ln_timing )   CALL timing_stop('ssh_nxt') 
    … …  
    128132 
    129133    
    130    SUBROUTINE wzv( kt ) 
     134   SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa ) 
    131135      !!---------------------------------------------------------------------- 
    132136      !!                ***  ROUTINE wzv  *** 
    … …  
    139143      !!        The boundary conditions are w=0 at the bottom (no flux) and. 
    140144      !! 
    141       !! ** action  :   wn      : now vertical velocity 
     145      !! ** action  :   pww      : now vertical velocity 
    142146      !! 
    143147      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    144148      !!---------------------------------------------------------------------- 
    145       INTEGER, INTENT(in) ::   kt   ! time step 
     149      INTEGER                         , INTENT(in)    ::   kt             ! time step 
     150      INTEGER                         , INTENT(in)    ::   Kbb, Kmm, Kaa  ! time level indices 
     151      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pww            ! now vertical velocity 
    146152      ! 
    147153      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    … …  
    157163         IF(lwp) WRITE(numout,*) '~~~~~ ' 
    158164         ! 
    159          wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all) 
     165         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all) 
    160166      ENDIF 
    161167      !                                           !------------------------------! 
    … …  
    171177            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 
    172178            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 
    173             DO jj = 2, jpjm1 
    174                DO ji = fs_2, fs_jpim1   ! vector opt. 
    175                   zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 
    176                END DO 
    177             END DO 
     179            DO_2D_00_00 
     180               zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 
     181            END_2D 
    178182         END DO 
    179183         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions" 
    180184         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells? 
    181          !                             ! Same question holds for hdivn. Perhaps just for security 
     185         !                             ! Same question holds for hdiv. Perhaps just for security 
    182186         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
    183187            ! computation of w 
    184             wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)    & 
    185                &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )     ) * tmask(:,:,jk) 
     188            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk) + zhdiv(:,:,jk)    & 
     189               &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )     ) * tmask(:,:,jk) 
    186190         END DO 
    187          !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 
     191         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0 
    188192         DEALLOCATE( zhdiv )  
    189193      ELSE   ! z_star and linear free surface cases 
    190194         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
    191195            ! computation of w 
    192             wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk)                 & 
    193                &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )  ) * tmask(:,:,jk) 
     196            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                 & 
     197               &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk) 
    194198         END DO 
    195199      ENDIF 
    … …  
    197201      IF( ln_bdy ) THEN 
    198202         DO jk = 1, jpkm1 
    199             wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 
     203            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:) 
    200204         END DO 
    201205      ENDIF 
    … …  
    203207#if defined key_agrif  
    204208      IF( .NOT. AGRIF_Root() ) THEN  
    205          IF ((nbondi ==  1).OR.(nbondi == 2)) wn(nlci-1 , :     ,:) = 0.e0      ! east  
    206          IF ((nbondi == -1).OR.(nbondi == 2)) wn(2      , :     ,:) = 0.e0      ! west  
    207          IF ((nbondj ==  1).OR.(nbondj == 2)) wn(:      ,nlcj-1 ,:) = 0.e0      ! north  
    208          IF ((nbondj == -1).OR.(nbondj == 2)) wn(:      ,2      ,:) = 0.e0      ! south  
     209         IF ((nbondi ==  1).OR.(nbondi == 2)) pww(nlci-1 , :     ,:) = 0.e0      ! east  
     210         IF ((nbondi == -1).OR.(nbondi == 2)) pww(2      , :     ,:) = 0.e0      ! west  
     211         IF ((nbondj ==  1).OR.(nbondj == 2)) pww(:      ,nlcj-1 ,:) = 0.e0      ! north  
     212         IF ((nbondj == -1).OR.(nbondj == 2)) pww(:      ,2      ,:) = 0.e0      ! south  
    209213      ENDIF  
    210214#endif  
    … …  
    215219 
    216220 
    217    SUBROUTINE ssh_swp( kt ) 
    218       !!---------------------------------------------------------------------- 
    219       !!                    ***  ROUTINE ssh_nxt  *** 
    220       !! 
    221       !! ** Purpose :   achieve the sea surface  height time stepping by  
    222       !!              applying Asselin time filter and swapping the arrays 
    223       !!              ssha  already computed in ssh_nxt   
     221   SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh ) 
     222      !!---------------------------------------------------------------------- 
     223      !!                    ***  ROUTINE ssh_atf  *** 
     224      !! 
     225      !! ** Purpose :   Apply Asselin time filter to now SSH. 
    224226      !! 
    225227      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing 
    226228      !!              from the filter, see Leclair and Madec 2010) and swap : 
    227       !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 
     229      !!                pssh(:,:,Kmm) = pssh(:,:,Kaa) + atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) ) 
    228230      !!                            - atfp * rdt * ( emp_b - emp ) / rau0 
    229       !!                sshn = ssha 
    230       !! 
    231       !! ** action  : - sshb, sshn   : before & now sea surface height 
    232       !!                               ready for the next time step 
     231      !! 
     232      !! ** action  : - pssh(:,:,Kmm) time filtered 
    233233      !! 
    234234      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    235235      !!---------------------------------------------------------------------- 
    236       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     236      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index 
     237      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! ocean time level indices 
     238      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! SSH field 
    237239      ! 
    238240      REAL(wp) ::   zcoef   ! local scalar 
    239241      !!---------------------------------------------------------------------- 
    240242      ! 
    241       IF( ln_timing )   CALL timing_start('ssh_swp') 
     243      IF( ln_timing )   CALL timing_start('ssh_atf') 
    242244      ! 
    243245      IF( kt == nit000 ) THEN 
    244246         IF(lwp) WRITE(numout,*) 
    245          IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 
     247         IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height' 
    246248         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    247249      ENDIF 
    248250      !              !==  Euler time-stepping: no filter, just swap  ==! 
    249       IF ( neuler == 0 .AND. kt == nit000 ) THEN 
    250          sshn(:,:) = ssha(:,:)                              ! now    <-- after  (before already = now) 
    251          ! 
    252       ELSE           !==  Leap-Frog time-stepping: Asselin filter + swap  ==! 
    253          !                                                  ! before <-- now filtered 
    254          sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) 
    255          IF( .NOT.ln_linssh ) THEN                          ! before <-- with forcing removed 
     251      IF ( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN   ! Only do time filtering for leapfrog timesteps 
     252         !                                                  ! filtered "now" field 
     253         pssh(:,:,Kmm) = pssh(:,:,Kmm) + atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) ) 
     254         IF( .NOT.ln_linssh ) THEN                          ! "now" <-- with forcing removed 
    256255            zcoef = atfp * rdt * r1_rau0 
    257             sshb(:,:) = sshb(:,:) - zcoef * (     emp_b(:,:) - emp   (:,:)   & 
    258                &                             -    rnf_b(:,:) + rnf   (:,:)   & 
    259                &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:) 
     256            pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * (     emp_b(:,:) - emp   (:,:)   & 
     257               &                             - rnf_b(:,:)        + rnf   (:,:)       & 
     258               &                             + fwfisf_cav_b(:,:) - fwfisf_cav(:,:)   & 
     259               &                             + fwfisf_par_b(:,:) - fwfisf_par(:,:)   ) * ssmask(:,:) 
     260 
     261            ! ice sheet coupling 
     262            IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - atfp * rdt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:) 
     263 
    260264         ENDIF 
    261          sshn(:,:) = ssha(:,:)                              ! now <-- after 
    262       ENDIF 
    263       ! 
    264       IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask ) 
    265       ! 
    266       IF( ln_timing )   CALL timing_stop('ssh_swp') 
    267       ! 
    268    END SUBROUTINE ssh_swp 
    269  
    270    SUBROUTINE wAimp( kt ) 
     265      ENDIF 
     266      ! 
     267      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' pssh(:,:,Kmm)  - : ', mask1=tmask ) 
     268      ! 
     269      IF( ln_timing )   CALL timing_stop('ssh_atf') 
     270      ! 
     271   END SUBROUTINE ssh_atf 
     272 
     273   SUBROUTINE wAimp( kt, Kmm ) 
    271274      !!---------------------------------------------------------------------- 
    272275      !!                ***  ROUTINE wAimp  *** 
    … …  
    277280      !! ** Method  : -  
    278281      !! 
    279       !! ** action  :   wn      : now vertical velocity (to be handled explicitly) 
     282      !! ** action  :   ww      : now vertical velocity (to be handled explicitly) 
    280283      !!            :   wi      : now vertical velocity (for implicit treatment) 
    281284      !! 
    … …  
    285288      !!---------------------------------------------------------------------- 
    286289      INTEGER, INTENT(in) ::   kt   ! time step 
     290      INTEGER, INTENT(in) ::   Kmm  ! time level index 
    287291      ! 
    288292      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    … …  
    305309      ! Calculate Courant numbers 
    306310      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
    307          DO jk = 1, jpkm1 
    308             DO jj = 2, jpjm1 
    309                DO ji = 2, fs_jpim1   ! vector opt. 
    310                   z1_e3t = 1._wp / e3t_n(ji,jj,jk) 
    311                   ! 2*rdt and not r2dt (for restartability) 
    312                   Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )                       &   
    313                      &                             + ( MAX( e2u(ji  ,jj)*e3u_n(ji  ,jj,jk)*un(ji  ,jj,jk) + un_td(ji  ,jj,jk), 0._wp ) -   & 
    314                      &                                 MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk) + un_td(ji-1,jj,jk), 0._wp ) )   & 
    315                      &                               * r1_e1e2t(ji,jj)                                                                     & 
    316                      &                             + ( MAX( e1v(ji,jj  )*e3v_n(ji,jj  ,jk)*vn(ji,jj  ,jk) + vn_td(ji,jj  ,jk), 0._wp ) -   & 
    317                      &                                 MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk) + vn_td(ji,jj-1,jk), 0._wp ) )   & 
    318                      &                               * r1_e1e2t(ji,jj)                                                                     & 
    319                      &                             ) * z1_e3t 
    320                END DO 
    321             END DO 
    322          END DO 
     311         DO_3D_00_00( 1, jpkm1 ) 
     312            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 
     313            ! 2*rdt and not r2dt (for restartability) 
     314            Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )                       &   
     315               &                             + ( MAX( e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm) + un_td(ji  ,jj,jk), 0._wp ) -   & 
     316               &                                 MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) )   & 
     317               &                               * r1_e1e2t(ji,jj)                                                                     & 
     318               &                             + ( MAX( e1v(ji,jj  )*e3v(ji,jj  ,jk,Kmm)*vv(ji,jj  ,jk,Kmm) + vn_td(ji,jj  ,jk), 0._wp ) -   & 
     319               &                                 MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) )   & 
     320               &                               * r1_e1e2t(ji,jj)                                                                     & 
     321               &                             ) * z1_e3t 
     322         END_3D 
    323323      ELSE 
    324          DO jk = 1, jpkm1 
    325             DO jj = 2, jpjm1 
    326                DO ji = 2, fs_jpim1   ! vector opt. 
    327                   z1_e3t = 1._wp / e3t_n(ji,jj,jk) 
    328                   ! 2*rdt and not r2dt (for restartability) 
    329                   Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )   &  
    330                      &                             + ( MAX( e2u(ji  ,jj)*e3u_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   & 
    331                      &                                 MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   & 
    332                      &                               * r1_e1e2t(ji,jj)                                                 & 
    333                      &                             + ( MAX( e1v(ji,jj  )*e3v_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   & 
    334                      &                                 MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   & 
    335                      &                               * r1_e1e2t(ji,jj)                                                 & 
    336                      &                             ) * z1_e3t 
    337                END DO 
    338             END DO 
    339          END DO 
     324         DO_3D_00_00( 1, jpkm1 ) 
     325            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 
     326            ! 2*rdt and not r2dt (for restartability) 
     327            Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )   &  
     328               &                             + ( MAX( e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm), 0._wp ) -   & 
     329               &                                 MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) )   & 
     330               &                               * r1_e1e2t(ji,jj)                                                 & 
     331               &                             + ( MAX( e1v(ji,jj  )*e3v(ji,jj  ,jk,Kmm)*vv(ji,jj  ,jk,Kmm), 0._wp ) -   & 
     332               &                                 MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) )   & 
     333               &                               * r1_e1e2t(ji,jj)                                                 & 
     334               &                             ) * z1_e3t 
     335         END_3D 
    340336      ENDIF 
    341337      CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. ) 
    … …  
    344340      ! 
    345341      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere 
    346          DO jk = jpkm1, 2, -1                           ! or scan Courant criterion and partition 
    347             DO jj = 1, jpj                              ! w where necessary 
    348                DO ji = 1, jpi 
    349                   ! 
    350                   zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) ) 
     342         DO_3DS_11_11( jpkm1, 2, -1 ) 
     343            ! 
     344            zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) ) 
    351345! alt: 
    352 !                  IF ( wn(ji,jj,jk) > 0._wp ) THEN  
     346!                  IF ( ww(ji,jj,jk) > 0._wp ) THEN  
    353347!                     zCu =  Cu_adv(ji,jj,jk)  
    354348!                  ELSE 
    355349!                     zCu =  Cu_adv(ji,jj,jk-1) 
    356350!                  ENDIF  
    357                   ! 
    358                   IF( zCu <= Cu_min ) THEN              !<-- Fully explicit 
    359                      zcff = 0._wp 
    360                   ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit 
    361                      zcff = ( zCu - Cu_min )**2 
    362                      zcff = zcff / ( Fcu + zcff ) 
    363                   ELSE                                  !<-- Mostly implicit 
    364                      zcff = ( zCu - Cu_max )/ zCu 
    365                   ENDIF 
    366                   zcff = MIN(1._wp, zcff) 
    367                   ! 
    368                   wi(ji,jj,jk) =           zcff   * wn(ji,jj,jk) 
    369                   wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk) 
    370                   ! 
    371                   Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl 
    372                END DO 
    373             END DO 
    374          END DO 
     351            ! 
     352            IF( zCu <= Cu_min ) THEN              !<-- Fully explicit 
     353               zcff = 0._wp 
     354            ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit 
     355               zcff = ( zCu - Cu_min )**2 
     356               zcff = zcff / ( Fcu + zcff ) 
     357            ELSE                                  !<-- Mostly implicit 
     358               zcff = ( zCu - Cu_max )/ zCu 
     359            ENDIF 
     360            zcff = MIN(1._wp, zcff) 
     361            ! 
     362            wi(ji,jj,jk) =           zcff   * ww(ji,jj,jk) 
     363            ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk) 
     364            ! 
     365            Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl 
     366         END_3D 
    375367         Cu_adv(:,:,1) = 0._wp  
    376368      ELSE 
    … …  
    381373      CALL iom_put("wimp",wi)  
    382374      CALL iom_put("wi_cff",Cu_adv) 
    383       CALL iom_put("wexp",wn) 
     375      CALL iom_put("wexp",ww) 
    384376      ! 
    385377      IF( ln_timing )   CALL timing_stop('wAimp') 
Note: See TracChangeset for help on using the changeset viewer.