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/trunk/src/OCE/DYN – NEMO

source: NEMO/trunk/src/OCE/DYN/sshwzv.F90

Last change on this file was 15150, checked in by smasson, 3 years ago

trunk: remove some of the changes done in [15145]

  • Property svn:keywords set to Id
File size: 22.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   !!             -   !  2020-08  (S. Techene, G. Madec)  add here ssh initiatlisation
14   !!----------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------
17   !!   ssh_nxt       : after ssh
18   !!   ssh_atf       : time filter the ssh arrays
19   !!   wzv           : compute now vertical velocity
20   !!----------------------------------------------------------------------
21   USE oce            ! ocean dynamics and tracers variables
22   USE isf_oce        ! ice shelf
23   USE dom_oce        ! ocean space and time domain variables
24   USE sbc_oce        ! surface boundary condition: ocean
25   USE domvvl         ! Variable volume
26   USE divhor         ! horizontal divergence
27   USE phycst         ! physical constants
28   USE bdy_oce , ONLY : ln_bdy, bdytmask   ! Open BounDarY
29   USE bdydyn2d       ! bdy_ssh routine
30   USE wet_dry        ! Wetting/Drying flux limiting
31#if defined key_agrif
32   USE agrif_oce
33   USE agrif_oce_interp
34#endif
35   !
36   USE iom 
37   USE in_out_manager ! I/O manager
38   USE restart        ! only for lrst_oce
39   USE prtctl         ! Print control
40   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
41   USE lib_mpp        ! MPP library
42   USE timing         ! Timing
43   
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC   ssh_nxt        ! called by step.F90
48   PUBLIC   wzv            ! called by step.F90
49   PUBLIC   wAimp          ! called by step.F90
50   PUBLIC   ssh_atf        ! called by step.F90
51
52   !! * Substitutions
53#  include "do_loop_substitute.h90"
54#  include "domzgr_substitute.h90"
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  ::   ji, jj, 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_3D( 1, nn_hls, 1, nn_hls, 1, jpkm1 )                                 ! Horizontal divergence of barotropic transports
106        zhdiv(ji,jj) = zhdiv(ji,jj) + e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk)
107      END_3D
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      DO_2D_OVR( 1, nn_hls, 1, nn_hls )                ! Loop bounds limited by hdiv definition in div_hor
113         pssh(ji,jj,Kaa) = (  pssh(ji,jj,Kbb) - rDt * ( zcoef * ( emp_b(ji,jj) + emp(ji,jj) ) + zhdiv(ji,jj) )  ) * ssmask(ji,jj)
114      END_2D
115      ! pssh must be defined everywhere (true for dyn_spg_ts, not for dyn_spg_exp)
116      IF ( .NOT. ln_dynspg_ts .AND. nn_hls == 2 ) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp )
117      !
118#if defined key_agrif
119      Kbb_a = Kbb   ;   Kmm_a = Kmm   ;   Krhs_a = Kaa
120      CALL agrif_ssh( kt )
121#endif
122      !
123      IF ( .NOT.ln_dynspg_ts ) THEN
124         IF( ln_bdy ) THEN
125            IF (nn_hls==1) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp )    ! Not sure that's necessary
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( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
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( nn_hls == 1)   CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
189         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
190         !                             ! Same question holds for hdiv. Perhaps just for security
191         !                             ! clem: yes it is a problem because ww is used in many other places where we need the halos
192         !
193         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 )     ! integrate from the bottom the hor. divergence
194            ! computation of w
195            pww(ji,jj,jk) = pww(ji,jj,jk+1) - (   e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk)   &
196               &                                  +                  zhdiv(ji,jj,jk)   &
197               &                                  + r1_Dt * (  e3t(ji,jj,jk,Kaa)       &
198               &                                             - e3t(ji,jj,jk,Kbb) )   ) * tmask(ji,jj,jk)
199         END_3D
200         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
201         DEALLOCATE( zhdiv ) 
202         !                                            !=================================!
203      ELSEIF( ln_linssh )   THEN                      !==  linear free surface cases  ==!
204         !                                            !=================================!
205         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 )     ! integrate from the bottom the hor. divergence
206            pww(ji,jj,jk) = pww(ji,jj,jk+1) - (  e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk)  ) * tmask(ji,jj,jk)
207         END_3D
208         !                                            !==========================================!
209      ELSE                                            !==  Quasi-Eulerian vertical coordinate  ==!   ('key_qco')
210         !                                            !==========================================!
211         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 )     ! integrate from the bottom the hor. divergence
212            pww(ji,jj,jk) = pww(ji,jj,jk+1) - (  e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk)    &
213               &                                 + r1_Dt * (  e3t(ji,jj,jk,Kaa)        &
214               &                                            - e3t(ji,jj,jk,Kbb)  )   ) * tmask(ji,jj,jk)
215         END_3D
216      ENDIF
217
218      IF( ln_bdy ) THEN
219         DO jk = 1, jpkm1
220            DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
221               pww(ji,jj,jk) = pww(ji,jj,jk) * bdytmask(ji,jj)
222            END_2D
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 )
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), INTENT(inout) ::   pssh           ! SSH field
289      !
290      REAL(wp) ::   zcoef   ! local scalar
291      !!----------------------------------------------------------------------
292      !
293      IF( ln_timing )   CALL timing_start('ssh_atf')
294      !
295      IF( kt == nit000 ) THEN
296         IF(lwp) WRITE(numout,*)
297         IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
298         IF(lwp) WRITE(numout,*) '~~~~~~~ '
299      ENDIF
300      !
301      IF( .NOT.l_1st_euler ) THEN   ! Apply Asselin time filter on Kmm field (not on euler 1st)
302         !
303         IF( ln_linssh ) THEN                ! filtered "now" field
304            pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
305            !
306         ELSE                                ! filtered "now" field with forcing removed
307            zcoef = rn_atfp * rn_Dt * r1_rho0
308            pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )   &
309               &                          - zcoef   * (         emp_b(:,:) -        emp(:,:)   &
310               &                                              - rnf_b(:,:) +        rnf(:,:)   &
311               &                                       - fwfisf_cav_b(:,:) + fwfisf_cav(:,:)   &
312               &                                       - fwfisf_par_b(:,:) + fwfisf_par(:,:)   ) * ssmask(:,:)
313
314            ! ice sheet coupling
315            IF( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1 )   &
316               &   pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0._wp ) * ssmask(:,:)
317
318         ENDIF
319      ENDIF
320      !
321      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' atf  - pssh(:,:,Kmm): ', mask1=tmask )
322      !
323      IF( ln_timing )   CALL timing_stop('ssh_atf')
324      !
325   END SUBROUTINE ssh_atf
326
327   
328   SUBROUTINE wAimp( kt, Kmm )
329      !!----------------------------------------------------------------------
330      !!                ***  ROUTINE wAimp  ***
331      !!                   
332      !! ** Purpose :   compute the Courant number and partition vertical velocity
333      !!                if a proportion needs to be treated implicitly
334      !!
335      !! ** Method  : -
336      !!
337      !! ** action  :   ww      : now vertical velocity (to be handled explicitly)
338      !!            :   wi      : now vertical velocity (for implicit treatment)
339      !!
340      !! Reference  : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
341      !!              implicit scheme for vertical advection in oceanic modeling.
342      !!              Ocean Modelling, 91, 38-69.
343      !!----------------------------------------------------------------------
344      INTEGER, INTENT(in) ::   kt   ! time step
345      INTEGER, INTENT(in) ::   Kmm  ! time level index
346      !
347      INTEGER  ::   ji, jj, jk   ! dummy loop indices
348      REAL(wp)             ::   zCu, zcff, z1_e3t, zdt                ! local scalars
349      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters
350      REAL(wp) , PARAMETER ::   Cu_max = 0.30_wp                      ! local parameters
351      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters
352      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters
353      !!----------------------------------------------------------------------
354      !
355      IF( ln_timing )   CALL timing_start('wAimp')
356      !
357      IF( kt == nit000 ) THEN
358         IF(lwp) WRITE(numout,*)
359         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
360         IF(lwp) WRITE(numout,*) '~~~~~ '
361      ENDIF
362      !
363      ! Calculate Courant numbers
364      zdt = 2._wp * rn_Dt                            ! 2*rn_Dt and not rDt (for restartability)
365      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
366         DO_3D( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )
367            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
368            Cu_adv(ji,jj,jk) =   zdt *                                                         &
369               &  ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )            &
370               &  + ( MAX( e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)                                  &
371               &                        * uu (ji  ,jj,jk,Kmm) + un_td(ji  ,jj,jk), 0._wp ) -   &
372               &      MIN( e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)                                  &
373               &                        * uu (ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) )   &
374               &                               * r1_e1e2t(ji,jj)                                                                     &
375               &  + ( MAX( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)                                  &
376               &                        * vv (ji,jj  ,jk,Kmm) + vn_td(ji,jj  ,jk), 0._wp ) -   &
377               &      MIN( e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)                                  &
378               &                        * vv (ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) )   &
379               &                               * r1_e1e2t(ji,jj)                                                                     &
380               &                             ) * z1_e3t
381         END_3D
382      ELSE
383         DO_3D( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )
384            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
385            Cu_adv(ji,jj,jk) =   zdt *                                                      &
386               &  ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )         &
387               &                             + ( MAX( e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm), 0._wp ) -   &
388               &                                 MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) )   &
389               &                               * r1_e1e2t(ji,jj)                                                 &
390               &                             + ( MAX( e1v(ji,jj  )*e3v(ji,jj  ,jk,Kmm)*vv(ji,jj  ,jk,Kmm), 0._wp ) -   &
391               &                                 MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) )   &
392               &                               * r1_e1e2t(ji,jj)                                                 &
393               &                             ) * z1_e3t
394         END_3D
395      ENDIF
396      CALL iom_put("Courant",Cu_adv)
397      !
398      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
399         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 2, -1 )    ! or scan Courant criterion and partition ! w where necessary
400            !
401            zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
402! alt:
403!                  IF ( ww(ji,jj,jk) > 0._wp ) THEN
404!                     zCu =  Cu_adv(ji,jj,jk)
405!                  ELSE
406!                     zCu =  Cu_adv(ji,jj,jk-1)
407!                  ENDIF
408            !
409            IF( zCu <= Cu_min ) THEN              !<-- Fully explicit
410               zcff = 0._wp
411            ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
412               zcff = ( zCu - Cu_min )**2
413               zcff = zcff / ( Fcu + zcff )
414            ELSE                                  !<-- Mostly implicit
415               zcff = ( zCu - Cu_max )/ zCu
416            ENDIF
417            zcff = MIN(1._wp, zcff)
418            !
419            wi(ji,jj,jk) =           zcff   * ww(ji,jj,jk)
420            ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
421            !
422            Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl
423         END_3D
424         Cu_adv(:,:,1) = 0._wp 
425      ELSE
426         ! Fully explicit everywhere
427         Cu_adv(:,:,:) = 0._wp                    ! Reuse array to output coefficient below and in stp_ctl
428         wi    (:,:,:) = 0._wp
429      ENDIF
430      CALL iom_put("wimp",wi) 
431      CALL iom_put("wi_cff",Cu_adv)
432      CALL iom_put("wexp",ww)
433      !
434      IF( ln_timing )   CALL timing_stop('wAimp')
435      !
436   END SUBROUTINE wAimp
437     
438   !!======================================================================
439END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.