source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/sshwzv.F90 @ 13167

Last change on this file since 13167 was 13167, checked in by techene, 3 months ago

#2385 : pass sette tests without key_qco, oops should be true this time

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