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/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/DYN/sshwzv.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 4 years ago

Add neighborhood collectives calls in the NEMO src - ticket #2496
  • Property svn:keywords set to Id
File size: 22.0 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
31   USE agrif_oce_interp
32#endif
33   !
34   USE iom 
35   USE in_out_manager ! I/O manager
36   USE restart        ! only for lrst_oce
37   USE prtctl         ! Print control
38   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
39   USE lib_mpp        ! MPP library
40   USE timing         ! Timing
41   USE wet_dry        ! Wetting/Drying flux limiting
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC   ssh_nxt    ! called by step.F90
47   PUBLIC   wzv        ! called by step.F90
48   PUBLIC   wAimp      ! called by step.F90
49   PUBLIC   ssh_atf    ! called by step.F90
50
51   !! * Substitutions
52#  include "do_loop_substitute.h90"
53#  include "domzgr_substitute.h90"
54
55   !!----------------------------------------------------------------------
56   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
57   !! $Id$
58   !! Software governed by the CeCILL license (see ./LICENSE)
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa )
63      !!----------------------------------------------------------------------
64      !!                ***  ROUTINE ssh_nxt  ***
65      !!                   
66      !! ** Purpose :   compute the after ssh (ssh(Kaa))
67      !!
68      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
69      !!      is computed by integrating the horizontal divergence and multiply by
70      !!      by the time step.
71      !!
72      !! ** action  :   ssh(:,:,Kaa), after sea surface height
73      !!
74      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
75      !!----------------------------------------------------------------------
76      INTEGER                         , INTENT(in   ) ::   kt             ! time step
77      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level index
78      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height
79      !
80      INTEGER  ::   jk      ! dummy loop index
81      REAL(wp) ::   zcoef   ! local scalar
82      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace
83      !!----------------------------------------------------------------------
84      !
85      IF( ln_timing )   CALL timing_start('ssh_nxt')
86      !
87      IF( kt == nit000 ) THEN
88         IF(lwp) WRITE(numout,*)
89         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
90         IF(lwp) WRITE(numout,*) '~~~~~~~ '
91      ENDIF
92      !
93      zcoef = 0.5_wp * r1_rho0
94
95      !                                           !------------------------------!
96      !                                           !   After Sea Surface Height   !
97      !                                           !------------------------------!
98      IF(ln_wd_il) THEN
99         CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, 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) - rDt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
113      !
114#if defined key_agrif
115      Kbb_a = Kbb   ;   Kmm_a = Kmm   ;   Krhs_a = Kaa
116      CALL agrif_ssh( kt )
117#endif
118      !
119      IF ( .NOT.ln_dynspg_ts ) THEN
120         IF( ln_bdy ) THEN
121#if defined key_mpi3
122            CALL lbc_lnk_nc_multi( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp )    ! Not sure that's necessary
123#else
124            CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp )    ! Not sure that's necessary
125#endif
126            CALL bdy_ssh( pssh(:,:,Kaa) )             ! Duplicate sea level across open boundaries
127         ENDIF
128      ENDIF
129      !                                           !------------------------------!
130      !                                           !           outputs            !
131      !                                           !------------------------------!
132      !
133      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa)  - : ', mask1=tmask )
134      !
135      IF( ln_timing )   CALL timing_stop('ssh_nxt')
136      !
137   END SUBROUTINE ssh_nxt
138
139   
140   SUBROUTINE wzv( kt, Kbb, Kmm, Kaa, pww )
141      !!----------------------------------------------------------------------
142      !!                ***  ROUTINE wzv  ***
143      !!                   
144      !! ** Purpose :   compute the now vertical velocity
145      !!
146      !! ** Method  : - Using the incompressibility hypothesis, the vertical
147      !!      velocity is computed by integrating the horizontal divergence 
148      !!      from the bottom to the surface minus the scale factor evolution.
149      !!        The boundary conditions are w=0 at the bottom (no flux) and.
150      !!
151      !! ** action  :   pww      : now vertical velocity
152      !!
153      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
154      !!----------------------------------------------------------------------
155      INTEGER                         , INTENT(in)    ::   kt             ! time step
156      INTEGER                         , INTENT(in)    ::   Kbb, Kmm, Kaa  ! time level indices
157      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pww            ! vertical velocity at Kmm
158      !
159      INTEGER  ::   ji, jj, jk   ! dummy loop indices
160      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv
161      !!----------------------------------------------------------------------
162      !
163      IF( ln_timing )   CALL timing_start('wzv')
164      !
165      IF( kt == nit000 ) THEN
166         IF(lwp) WRITE(numout,*)
167         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
168         IF(lwp) WRITE(numout,*) '~~~~~ '
169         !
170         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
171      ENDIF
172      !                                           !------------------------------!
173      !                                           !     Now Vertical Velocity    !
174      !                                           !------------------------------!
175      !
176      !                                               !===============================!
177      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      !==  z_tilde and layer cases  ==!
178         !                                            !===============================!
179         ALLOCATE( zhdiv(jpi,jpj,jpk) ) 
180         !
181         DO jk = 1, jpkm1
182            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
183            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
184            DO_2D( 0, 0, 0, 0 )
185               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) )
186            END_2D
187         END DO
188#if defined key_mpi3
189         CALL lbc_lnk_nc_multi('sshwzv', zhdiv, 'T', 1.0_wp)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
190#else
191         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
192#endif
193         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
194         !                             ! Same question holds for hdiv. Perhaps just for security
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               &                            +                  zhdiv(:,:,jk)   &
199               &                            + r1_Dt * (  e3t(:,:,jk,Kaa)       &
200               &                                       - e3t(:,:,jk,Kbb) )   ) * tmask(:,:,jk)
201         END DO
202         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
203         DEALLOCATE( zhdiv ) 
204         !                                            !=================================!
205      ELSEIF( ln_linssh )   THEN                      !==  linear free surface cases  ==!
206         !                                            !=================================!
207         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
208            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)  ) * tmask(:,:,jk)
209         END DO
210         !                                            !==========================================!
211      ELSE                                            !==  Quasi-Eulerian vertical coordinate  ==!   ('key_qco')
212         !                                            !==========================================!
213         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
214            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                 &
215               &                            + r1_Dt * (  e3t(:,:,jk,Kaa)        &
216               &                                       - e3t(:,:,jk,Kbb)  )   ) * tmask(:,:,jk)
217         END DO
218      ENDIF
219
220      IF( ln_bdy ) THEN
221         DO jk = 1, jpkm1
222            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
223         END DO
224      ENDIF
225      !
226#if defined key_agrif
227      IF( .NOT. AGRIF_Root() ) THEN
228         !
229         ! Mask vertical velocity at first/last columns/row
230         ! inside computational domain (cosmetic)
231         DO jk = 1, jpkm1
232            IF( lk_west ) THEN                             ! --- West --- !
233               DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
234                  DO jj = 1, jpj
235                     pww(ji,jj,jk) = 0._wp 
236                  END DO
237               END DO
238            ENDIF
239            IF( lk_east ) THEN                             ! --- East --- !
240               DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
241                  DO jj = 1, jpj
242                     pww(ji,jj,jk) = 0._wp
243                  END DO
244               END DO
245            ENDIF
246            IF( lk_south ) THEN                            ! --- South --- !
247               DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
248                  DO ji = 1, jpi
249                     pww(ji,jj,jk) = 0._wp
250                  END DO
251               END DO
252            ENDIF
253            IF( lk_north ) THEN                            ! --- North --- !
254               DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
255                  DO ji = 1, jpi
256                     pww(ji,jj,jk) = 0._wp
257                  END DO
258               END DO
259            ENDIF
260            !
261         END DO
262         !
263      ENDIF 
264#endif
265      !
266      IF( ln_timing )   CALL timing_stop('wzv')
267      !
268   END SUBROUTINE wzv
269
270
271   SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh, pssh_f )
272      !!----------------------------------------------------------------------
273      !!                    ***  ROUTINE ssh_atf  ***
274      !!
275      !! ** Purpose :   Apply Asselin time filter to now SSH.
276      !!
277      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
278      !!              from the filter, see Leclair and Madec 2010) and swap :
279      !!                pssh(:,:,Kmm) = pssh(:,:,Kaa) + rn_atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
280      !!                            - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0
281      !!
282      !! ** action  : - pssh(:,:,Kmm) time filtered
283      !!
284      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
285      !!----------------------------------------------------------------------
286      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index
287      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! ocean time level indices
288      REAL(wp), DIMENSION(jpi,jpj,jpt)          , TARGET, INTENT(inout) ::   pssh           ! SSH field
289      REAL(wp), DIMENSION(jpi,jpj    ), OPTIONAL, TARGET, INTENT(  out) ::   pssh_f         ! filtered SSH field
290      !
291      REAL(wp) ::   zcoef   ! local scalar
292      REAL(wp), POINTER, DIMENSION(:,:) ::   zssh   ! pointer for filtered SSH
293      !!----------------------------------------------------------------------
294      !
295      IF( ln_timing )   CALL timing_start('ssh_atf')
296      !
297      IF( kt == nit000 ) THEN
298         IF(lwp) WRITE(numout,*)
299         IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
300         IF(lwp) WRITE(numout,*) '~~~~~~~ '
301      ENDIF
302      !              !==  Euler time-stepping: no filter, just swap  ==!
303      IF ( .NOT.( l_1st_euler ) ) THEN   ! Only do time filtering for leapfrog timesteps
304         IF( PRESENT( pssh_f ) ) THEN   ;   zssh => pssh_f
305         ELSE                           ;   zssh => pssh(:,:,Kmm)
306         ENDIF
307         !                                                  ! filtered "now" field
308         pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
309         IF( .NOT.ln_linssh ) THEN                          ! "now" <-- with forcing removed
310            zcoef = rn_atfp * rn_Dt * r1_rho0
311            pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
312               &                             - rnf_b(:,:)        + rnf   (:,:)       &
313               &                             + fwfisf_cav_b(:,:) - fwfisf_cav(:,:)   &
314               &                             + fwfisf_par_b(:,:) - fwfisf_par(:,:)   ) * ssmask(:,:)
315
316            ! ice sheet coupling
317            IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
318
319         ENDIF
320      ENDIF
321      !
322      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' pssh(:,:,Kmm)  - : ', mask1=tmask )
323      !
324      IF( ln_timing )   CALL timing_stop('ssh_atf')
325      !
326   END SUBROUTINE ssh_atf
327
328   
329   SUBROUTINE wAimp( kt, Kmm )
330      !!----------------------------------------------------------------------
331      !!                ***  ROUTINE wAimp  ***
332      !!                   
333      !! ** Purpose :   compute the Courant number and partition vertical velocity
334      !!                if a proportion needs to be treated implicitly
335      !!
336      !! ** Method  : -
337      !!
338      !! ** action  :   ww      : now vertical velocity (to be handled explicitly)
339      !!            :   wi      : now vertical velocity (for implicit treatment)
340      !!
341      !! Reference  : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
342      !!              implicit scheme for vertical advection in oceanic modeling.
343      !!              Ocean Modelling, 91, 38-69.
344      !!----------------------------------------------------------------------
345      INTEGER, INTENT(in) ::   kt   ! time step
346      INTEGER, INTENT(in) ::   Kmm  ! time level index
347      !
348      INTEGER  ::   ji, jj, jk   ! dummy loop indices
349      REAL(wp)             ::   zCu, zcff, z1_e3t, zdt                ! local scalars
350      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters
351      REAL(wp) , PARAMETER ::   Cu_max = 0.30_wp                      ! local parameters
352      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters
353      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters
354      !!----------------------------------------------------------------------
355      !
356      IF( ln_timing )   CALL timing_start('wAimp')
357      !
358      IF( kt == nit000 ) THEN
359         IF(lwp) WRITE(numout,*)
360         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
361         IF(lwp) WRITE(numout,*) '~~~~~ '
362         wi(:,:,:) = 0._wp
363      ENDIF
364      !
365      ! Calculate Courant numbers
366      zdt = 2._wp * rn_Dt                            ! 2*rn_Dt and not rDt (for restartability)
367      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
368         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
369            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
370            Cu_adv(ji,jj,jk) =   zdt *                                                         &
371               &  ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )            &
372               &  + ( MAX( e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)                                  &
373               &                        * uu (ji  ,jj,jk,Kmm) + un_td(ji  ,jj,jk), 0._wp ) -   &
374               &      MIN( e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)                                  &
375               &                        * uu (ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) )   &
376               &                               * r1_e1e2t(ji,jj)                                                                     &
377               &  + ( MAX( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)                                  &
378               &                        * vv (ji,jj  ,jk,Kmm) + vn_td(ji,jj  ,jk), 0._wp ) -   &
379               &      MIN( e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)                                  &
380               &                        * vv (ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) )   &
381               &                               * r1_e1e2t(ji,jj)                                                                     &
382               &                             ) * z1_e3t
383         END_3D
384      ELSE
385         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
386            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
387            Cu_adv(ji,jj,jk) =   zdt *                                                      &
388               &  ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )         &
389               &                             + ( MAX( e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm), 0._wp ) -   &
390               &                                 MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) )   &
391               &                               * r1_e1e2t(ji,jj)                                                 &
392               &                             + ( MAX( e1v(ji,jj  )*e3v(ji,jj  ,jk,Kmm)*vv(ji,jj  ,jk,Kmm), 0._wp ) -   &
393               &                                 MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) )   &
394               &                               * r1_e1e2t(ji,jj)                                                 &
395               &                             ) * z1_e3t
396         END_3D
397      ENDIF
398#if defined key_mpi3
399      CALL lbc_lnk_nc_multi( 'sshwzv', Cu_adv, 'T', 1.0_wp )
400#else
401      CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1.0_wp )
402#endif
403      !
404      CALL iom_put("Courant",Cu_adv)
405      !
406      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
407         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )             ! or scan Courant criterion and partition ! w where necessary
408            !
409            zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
410! alt:
411!                  IF ( ww(ji,jj,jk) > 0._wp ) THEN
412!                     zCu =  Cu_adv(ji,jj,jk)
413!                  ELSE
414!                     zCu =  Cu_adv(ji,jj,jk-1)
415!                  ENDIF
416            !
417            IF( zCu <= Cu_min ) THEN              !<-- Fully explicit
418               zcff = 0._wp
419            ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
420               zcff = ( zCu - Cu_min )**2
421               zcff = zcff / ( Fcu + zcff )
422            ELSE                                  !<-- Mostly implicit
423               zcff = ( zCu - Cu_max )/ zCu
424            ENDIF
425            zcff = MIN(1._wp, zcff)
426            !
427            wi(ji,jj,jk) =           zcff   * ww(ji,jj,jk)
428            ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
429            !
430            Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl
431         END_3D
432         Cu_adv(:,:,1) = 0._wp 
433      ELSE
434         ! Fully explicit everywhere
435         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient below and in stp_ctl
436         wi    (:,:,:) = 0._wp
437      ENDIF
438      CALL iom_put("wimp",wi) 
439      CALL iom_put("wi_cff",Cu_adv)
440      CALL iom_put("wexp",ww)
441      !
442      IF( ln_timing )   CALL timing_stop('wAimp')
443      !
444   END SUBROUTINE wAimp
445   !!======================================================================
446END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.