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.
sshwzv.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/sshwzv.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:


Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

  • Property svn:keywords set to Id
File size: 19.2 KB
Line 
1MODULE sshwzv   
2   !!==============================================================================
3   !!                       ***  MODULE  sshwzv  ***
4   !! Ocean dynamics : sea surface height and vertical velocity
5   !!==============================================================================
6   !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code
7   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
8   !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
9   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
10   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work
11   !!            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.
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   ssh_nxt       : after ssh
17   !!   ssh_atf       : time filter the ssh arrays
18   !!   wzv           : compute now vertical velocity
19   !!----------------------------------------------------------------------
20   USE oce            ! ocean dynamics and tracers variables
21   USE isf_oce        ! ice shelf
22   USE dom_oce        ! ocean space and time domain variables
23   USE sbc_oce        ! surface boundary condition: ocean
24   USE domvvl         ! Variable volume
25   USE divhor         ! horizontal divergence
26   USE phycst         ! physical constants
27   USE bdy_oce , ONLY : ln_bdy, bdytmask   ! Open BounDarY
28   USE bdydyn2d       ! bdy_ssh routine
29#if defined key_agrif
30   USE agrif_oce_interp
31#endif
32   !
33   USE iom 
34   USE in_out_manager ! I/O manager
35   USE restart        ! only for lrst_oce
36   USE prtctl         ! Print control
37   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
38   USE lib_mpp        ! MPP library
39   USE timing         ! Timing
40   USE wet_dry        ! Wetting/Drying flux limiting
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   ssh_nxt    ! called by step.F90
46   PUBLIC   wzv        ! called by step.F90
47   PUBLIC   wAimp      ! called by step.F90
48   PUBLIC   ssh_atf    ! called by step.F90
49
50   !! * Substitutions
51#  include "vectopt_loop_substitute.h90"
52#  include "do_loop_substitute.h90"
53   !!----------------------------------------------------------------------
54   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
55   !! $Id$
56   !! Software governed by the CeCILL license (see ./LICENSE)
57   !!----------------------------------------------------------------------
58CONTAINS
59
60   SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa )
61      !!----------------------------------------------------------------------
62      !!                ***  ROUTINE ssh_nxt  ***
63      !!                   
64      !! ** Purpose :   compute the after ssh (ssh(Kaa))
65      !!
66      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
67      !!      is computed by integrating the horizontal divergence and multiply by
68      !!      by the time step.
69      !!
70      !! ** action  :   ssh(:,:,Kaa), after sea surface height
71      !!
72      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
73      !!----------------------------------------------------------------------
74      INTEGER                         , INTENT(in   ) ::   kt             ! time step
75      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level index
76      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height
77      !
78      INTEGER  ::   jk            ! dummy loop indice
79      REAL(wp) ::   z2dt, zcoef   ! local scalars
80      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace
81      !!----------------------------------------------------------------------
82      !
83      IF( ln_timing )   CALL timing_start('ssh_nxt')
84      !
85      IF( kt == nit000 ) THEN
86         IF(lwp) WRITE(numout,*)
87         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
88         IF(lwp) WRITE(numout,*) '~~~~~~~ '
89      ENDIF
90      !
91      z2dt = 2._wp * rdt                          ! set time step size (Euler/Leapfrog)
92      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
93      zcoef = 0.5_wp * r1_rau0
94
95      !                                           !------------------------------!
96      !                                           !   After Sea Surface Height   !
97      !                                           !------------------------------!
98      IF(ln_wd_il) THEN
99         CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), z2dt, Kmm, uu, vv )
100      ENDIF
101
102      CALL div_hor( kt, Kbb, Kmm )                     ! Horizontal divergence
103      !
104      zhdiv(:,:) = 0._wp
105      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
106        zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
107      END DO
108      !                                                ! Sea surface elevation time stepping
109      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
110      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
111      !
112      pssh(:,:,Kaa) = (  pssh(:,:,Kbb) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
113      !
114#if defined key_agrif
115      Kbb_a = Kbb; Kmm_a = Kmm; Krhs_a = Kaa; CALL agrif_ssh( kt )
116#endif
117      !
118      IF ( .NOT.ln_dynspg_ts ) THEN
119         IF( ln_bdy ) THEN
120            CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1. )    ! Not sure that's necessary
121            CALL bdy_ssh( pssh(:,:,Kaa) )             ! Duplicate sea level across open boundaries
122         ENDIF
123      ENDIF
124      !                                           !------------------------------!
125      !                                           !           outputs            !
126      !                                           !------------------------------!
127      !
128      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa)  - : ', mask1=tmask )
129      !
130      IF( ln_timing )   CALL timing_stop('ssh_nxt')
131      !
132   END SUBROUTINE ssh_nxt
133
134   
135   SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa )
136      !!----------------------------------------------------------------------
137      !!                ***  ROUTINE wzv  ***
138      !!                   
139      !! ** Purpose :   compute the now vertical velocity
140      !!
141      !! ** Method  : - Using the incompressibility hypothesis, the vertical
142      !!      velocity is computed by integrating the horizontal divergence 
143      !!      from the bottom to the surface minus the scale factor evolution.
144      !!        The boundary conditions are w=0 at the bottom (no flux) and.
145      !!
146      !! ** action  :   pww      : now vertical velocity
147      !!
148      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
149      !!----------------------------------------------------------------------
150      INTEGER                         , INTENT(in)    ::   kt             ! time step
151      INTEGER                         , INTENT(in)    ::   Kbb, Kmm, Kaa  ! time level indices
152      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pww            ! now vertical velocity
153      !
154      INTEGER  ::   ji, jj, jk   ! dummy loop indices
155      REAL(wp) ::   z1_2dt       ! local scalars
156      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv
157      !!----------------------------------------------------------------------
158      !
159      IF( ln_timing )   CALL timing_start('wzv')
160      !
161      IF( kt == nit000 ) THEN
162         IF(lwp) WRITE(numout,*)
163         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
164         IF(lwp) WRITE(numout,*) '~~~~~ '
165         !
166         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
167      ENDIF
168      !                                           !------------------------------!
169      !                                           !     Now Vertical Velocity    !
170      !                                           !------------------------------!
171      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
172      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
173      !
174      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
175         ALLOCATE( zhdiv(jpi,jpj,jpk) ) 
176         !
177         DO jk = 1, jpkm1
178            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
179            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
180            DO_2D_00_00
181               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) )
182            END_2D
183         END DO
184         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
185         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
186         !                             ! Same question holds for hdiv. Perhaps just for security
187         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
188            ! computation of w
189            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk) + zhdiv(:,:,jk)    &
190               &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )     ) * tmask(:,:,jk)
191         END DO
192         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
193         DEALLOCATE( zhdiv ) 
194      ELSE   ! z_star and linear free surface cases
195         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
196            ! computation of w
197            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                 &
198               &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk)
199         END DO
200      ENDIF
201
202      IF( ln_bdy ) THEN
203         DO jk = 1, jpkm1
204            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
205         END DO
206      ENDIF
207      !
208#if defined key_agrif 
209      IF( .NOT. AGRIF_Root() ) THEN
210         IF ((nbondi ==  1).OR.(nbondi == 2)) pww(nlci-1 , :     ,:) = 0.e0      ! east
211         IF ((nbondi == -1).OR.(nbondi == 2)) pww(2      , :     ,:) = 0.e0      ! west
212         IF ((nbondj ==  1).OR.(nbondj == 2)) pww(:      ,nlcj-1 ,:) = 0.e0      ! north
213         IF ((nbondj == -1).OR.(nbondj == 2)) pww(:      ,2      ,:) = 0.e0      ! south
214      ENDIF 
215#endif 
216      !
217      IF( ln_timing )   CALL timing_stop('wzv')
218      !
219   END SUBROUTINE wzv
220
221
222   SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )
223      !!----------------------------------------------------------------------
224      !!                    ***  ROUTINE ssh_atf  ***
225      !!
226      !! ** Purpose :   Apply Asselin time filter to now SSH.
227      !!
228      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
229      !!              from the filter, see Leclair and Madec 2010) and swap :
230      !!                pssh(:,:,Kmm) = pssh(:,:,Kaa) + atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
231      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
232      !!
233      !! ** action  : - pssh(:,:,Kmm) time filtered
234      !!
235      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
236      !!----------------------------------------------------------------------
237      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index
238      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! ocean time level indices
239      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! SSH field
240      !
241      REAL(wp) ::   zcoef   ! local scalar
242      !!----------------------------------------------------------------------
243      !
244      IF( ln_timing )   CALL timing_start('ssh_atf')
245      !
246      IF( kt == nit000 ) THEN
247         IF(lwp) WRITE(numout,*)
248         IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
249         IF(lwp) WRITE(numout,*) '~~~~~~~ '
250      ENDIF
251      !              !==  Euler time-stepping: no filter, just swap  ==!
252      IF ( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN   ! Only do time filtering for leapfrog timesteps
253         !                                                  ! filtered "now" field
254         pssh(:,:,Kmm) = pssh(:,:,Kmm) + atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
255         IF( .NOT.ln_linssh ) THEN                          ! "now" <-- with forcing removed
256            zcoef = atfp * rdt * r1_rau0
257            pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
258               &                             - rnf_b(:,:)        + rnf   (:,:)       &
259               &                             + fwfisf_cav_b(:,:) - fwfisf_cav(:,:)   &
260               &                             + fwfisf_par_b(:,:) - fwfisf_par(:,:)   ) * ssmask(:,:)
261
262            ! ice sheet coupling
263            IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - atfp * rdt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
264
265         ENDIF
266      ENDIF
267      !
268      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' pssh(:,:,Kmm)  - : ', mask1=tmask )
269      !
270      IF( ln_timing )   CALL timing_stop('ssh_atf')
271      !
272   END SUBROUTINE ssh_atf
273
274   SUBROUTINE wAimp( kt, Kmm )
275      !!----------------------------------------------------------------------
276      !!                ***  ROUTINE wAimp  ***
277      !!                   
278      !! ** Purpose :   compute the Courant number and partition vertical velocity
279      !!                if a proportion needs to be treated implicitly
280      !!
281      !! ** Method  : -
282      !!
283      !! ** action  :   ww      : now vertical velocity (to be handled explicitly)
284      !!            :   wi      : now vertical velocity (for implicit treatment)
285      !!
286      !! Reference  : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
287      !!              implicit scheme for vertical advection in oceanic modeling.
288      !!              Ocean Modelling, 91, 38-69.
289      !!----------------------------------------------------------------------
290      INTEGER, INTENT(in) ::   kt   ! time step
291      INTEGER, INTENT(in) ::   Kmm  ! time level index
292      !
293      INTEGER  ::   ji, jj, jk   ! dummy loop indices
294      REAL(wp)             ::   zCu, zcff, z1_e3t                     ! local scalars
295      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters
296      REAL(wp) , PARAMETER ::   Cu_max = 0.30_wp                      ! local parameters
297      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters
298      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters
299      !!----------------------------------------------------------------------
300      !
301      IF( ln_timing )   CALL timing_start('wAimp')
302      !
303      IF( kt == nit000 ) THEN
304         IF(lwp) WRITE(numout,*)
305         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
306         IF(lwp) WRITE(numout,*) '~~~~~ '
307         wi(:,:,:) = 0._wp
308      ENDIF
309      !
310      ! Calculate Courant numbers
311      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
312         DO_3D_00_00( 1, jpkm1 )
313            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
314            ! 2*rdt and not r2dt (for restartability)
315            Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )                       & 
316               &                             + ( MAX( e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm) + un_td(ji  ,jj,jk), 0._wp ) -   &
317               &                                 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 ) )   &
318               &                               * r1_e1e2t(ji,jj)                                                                     &
319               &                             + ( MAX( e1v(ji,jj  )*e3v(ji,jj  ,jk,Kmm)*vv(ji,jj  ,jk,Kmm) + vn_td(ji,jj  ,jk), 0._wp ) -   &
320               &                                 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 ) )   &
321               &                               * r1_e1e2t(ji,jj)                                                                     &
322               &                             ) * z1_e3t
323         END_3D
324      ELSE
325         DO_3D_00_00( 1, jpkm1 )
326            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
327            ! 2*rdt and not r2dt (for restartability)
328            Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )   & 
329               &                             + ( MAX( e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm), 0._wp ) -   &
330               &                                 MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) )   &
331               &                               * r1_e1e2t(ji,jj)                                                 &
332               &                             + ( MAX( e1v(ji,jj  )*e3v(ji,jj  ,jk,Kmm)*vv(ji,jj  ,jk,Kmm), 0._wp ) -   &
333               &                                 MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) )   &
334               &                               * r1_e1e2t(ji,jj)                                                 &
335               &                             ) * z1_e3t
336         END_3D
337      ENDIF
338      CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. )
339      !
340      CALL iom_put("Courant",Cu_adv)
341      !
342      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
343         DO_3DS_11_11( jpkm1, 2, -1 )
344            !
345            zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
346! alt:
347!                  IF ( ww(ji,jj,jk) > 0._wp ) THEN
348!                     zCu =  Cu_adv(ji,jj,jk)
349!                  ELSE
350!                     zCu =  Cu_adv(ji,jj,jk-1)
351!                  ENDIF
352            !
353            IF( zCu <= Cu_min ) THEN              !<-- Fully explicit
354               zcff = 0._wp
355            ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
356               zcff = ( zCu - Cu_min )**2
357               zcff = zcff / ( Fcu + zcff )
358            ELSE                                  !<-- Mostly implicit
359               zcff = ( zCu - Cu_max )/ zCu
360            ENDIF
361            zcff = MIN(1._wp, zcff)
362            !
363            wi(ji,jj,jk) =           zcff   * ww(ji,jj,jk)
364            ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
365            !
366            Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl
367         END_3D
368         Cu_adv(:,:,1) = 0._wp 
369      ELSE
370         ! Fully explicit everywhere
371         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient below and in stp_ctl
372         wi    (:,:,:) = 0._wp
373      ENDIF
374      CALL iom_put("wimp",wi) 
375      CALL iom_put("wi_cff",Cu_adv)
376      CALL iom_put("wexp",ww)
377      !
378      IF( ln_timing )   CALL timing_stop('wAimp')
379      !
380   END SUBROUTINE wAimp
381   !!======================================================================
382END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.