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

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/sshwzv.F90 @ 14618

Last change on this file since 14618 was 14618, checked in by techene, 3 years ago

#2506 RK3 work in progess : for 2 configurag
  • Property svn:keywords set to Id
File size: 28.9 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           : generic interface of vertical velocity calculation
20   !!   wzv_MLF       : MLF: compute NOW vertical velocity
21   !!   wzv_RK3       : RK3: Compute a vertical velocity
22   !!----------------------------------------------------------------------
23   USE oce            ! ocean dynamics and tracers variables
24   USE isf_oce        ! ice shelf
25   USE dom_oce        ! ocean space and time domain variables
26   USE sbc_oce        ! surface boundary condition: ocean
27   USE domvvl         ! Variable volume
28   USE divhor         ! horizontal divergence
29   USE phycst         ! physical constants
30   USE bdy_oce , ONLY : ln_bdy, bdytmask   ! Open BounDarY
31   USE bdydyn2d       ! bdy_ssh routine
32   USE wet_dry        ! Wetting/Drying flux limiting
33#if defined key_agrif
34   USE agrif_oce
35   USE agrif_oce_interp
36#endif
37   !
38   USE iom 
39   USE in_out_manager ! I/O manager
40   USE restart        ! only for lrst_oce
41   USE prtctl         ! Print control
42   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
43   USE lib_mpp        ! MPP library
44   USE timing         ! Timing
45   
46   IMPLICIT NONE
47   PRIVATE
48   !                  !! * Interface
49   INTERFACE wzv
50      MODULE PROCEDURE wzv_MLF, wzv_RK3
51   END INTERFACE
52
53   PUBLIC   ssh_nxt        ! called by step.F90
54   PUBLIC   wzv            ! called by step.F90
55   PUBLIC   wAimp          ! called by step.F90
56   PUBLIC   ssh_atf        ! called by step.F90
57
58   !! * Substitutions
59#  include "do_loop_substitute.h90"
60#  include "domzgr_substitute.h90"
61   !!----------------------------------------------------------------------
62   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
63   !! $Id$
64   !! Software governed by the CeCILL license (see ./LICENSE)
65   !!----------------------------------------------------------------------
66CONTAINS
67
68   SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa )
69      !!----------------------------------------------------------------------
70      !!                ***  ROUTINE ssh_nxt  ***
71      !!                   
72      !! ** Purpose :   compute the after ssh (ssh(Kaa))
73      !!
74      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
75      !!      is computed by integrating the horizontal divergence and multiply by
76      !!      by the time step.
77      !!
78      !! ** action  :   ssh(:,:,Kaa), after sea surface height
79      !!
80      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
81      !!----------------------------------------------------------------------
82      INTEGER                         , INTENT(in   ) ::   kt             ! time step
83      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level index
84      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height
85      !
86      INTEGER  ::   jk      ! dummy loop index
87      REAL(wp) ::   zcoef   ! local scalar
88      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace
89      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace
90      !!----------------------------------------------------------------------
91      !
92      IF( ln_timing )   CALL timing_start('ssh_nxt')
93      !
94      IF( kt == nit000 ) THEN
95         IF(lwp) WRITE(numout,*)
96         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
97         IF(lwp) WRITE(numout,*) '~~~~~~~ '
98      ENDIF
99      !
100      zcoef = 0.5_wp * r1_rho0
101      !                                           !------------------------------!
102      !                                           !   After Sea Surface Height   !
103      !                                           !------------------------------!
104      IF(ln_wd_il) THEN
105         CALL wad_lmt( pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )
106      ENDIF
107
108      CALL div_hor( kt, Kbb, Kmm )                     ! Horizontal divergence
109      !
110      zhdiv(:,:) = 0._wp
111      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
112        zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
113      END DO
114      !                                                ! Sea surface elevation time stepping
115      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
116      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
117      !
118      pssh(:,:,Kaa) = (  pssh(:,:,Kbb) - rDt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
119      !
120#if defined key_agrif
121      Kbb_a = Kbb   ;   Kmm_a = Kmm   ;   Krhs_a = Kaa
122      CALL agrif_ssh( kt )
123#endif
124      !
125      IF ( .NOT.ln_dynspg_ts ) THEN
126         IF( ln_bdy ) THEN
127            CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp )    ! Not sure that's necessary
128            CALL bdy_ssh( pssh(:,:,Kaa) )             ! Duplicate sea level across open boundaries
129         ENDIF
130      ENDIF
131      !                                           !------------------------------!
132      !                                           !           outputs            !
133      !                                           !------------------------------!
134      !
135      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa)  - : ', mask1=tmask )
136      !
137      IF( ln_timing )   CALL timing_stop('ssh_nxt')
138      !
139   END SUBROUTINE ssh_nxt
140
141
142   SUBROUTINE wzv_MLF( kt, Kbb, Kmm, Kaa, pww )
143      !!----------------------------------------------------------------------
144      !!                ***  ROUTINE wzv_MLF  ***
145      !!                   
146      !! ** Purpose :   compute the now vertical velocity
147      !!
148      !! ** Method  : - Using the incompressibility hypothesis, the vertical
149      !!      velocity is computed by integrating the horizontal divergence 
150      !!      from the bottom to the surface minus the scale factor evolution.
151      !!        The boundary conditions are w=0 at the bottom (no flux) and.
152      !!
153      !! ** action  :   pww      : now vertical velocity
154      !!
155      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
156      !!----------------------------------------------------------------------
157      INTEGER                         , INTENT(in)    ::   kt             ! time step
158      INTEGER                         , INTENT(in)    ::   Kbb, Kmm, Kaa  ! time level indices
159      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pww            ! vertical velocity at Kmm
160      !
161      INTEGER  ::   ji, jj, jk   ! dummy loop indices
162      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv
163      !!----------------------------------------------------------------------
164      !
165      IF( ln_timing )   CALL timing_start('wzv_MLF')
166      !
167      IF( kt == nit000 ) THEN
168         IF(lwp) WRITE(numout,*)
169         IF(lwp) WRITE(numout,*) 'wzv_MLF : now vertical velocity '
170         IF(lwp) WRITE(numout,*) '~~~~~~~'
171         !
172         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
173      ENDIF
174      !                                           !------------------------------!
175      !                                           !     Now Vertical Velocity    !
176      !                                           !------------------------------!
177      !
178      !                                               !===============================!
179      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      !==  z_tilde and layer cases  ==!
180         !                                            !===============================!
181         ALLOCATE( zhdiv(jpi,jpj,jpk) ) 
182         !
183         DO jk = 1, jpkm1
184            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
185            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
186            DO_2D( 0, 0, 0, 0 )
187               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) )
188            END_2D
189         END DO
190         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
191         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
192         !                             ! Same question holds for hdiv. Perhaps just for security
193         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
194            ! computation of w
195            pww(:,:,jk) = pww(:,:,jk+1) - (   e3t(:,:,jk,Kmm) * hdiv(:,:,jk)   &
196               &                            +                  zhdiv(:,:,jk)   &
197               &                            + r1_Dt * (  e3t(:,:,jk,Kaa)       &
198               &                                       - e3t(:,:,jk,Kbb) )   ) * tmask(:,:,jk)
199         END DO
200         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
201         DEALLOCATE( zhdiv )
202         !                                            !=================================!
203      ELSEIF( ln_linssh )   THEN                      !==  linear free surface cases  ==!
204         !                                            !=================================!
205         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
206            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)  ) * tmask(:,:,jk)
207         END DO
208         !                                            !==========================================!
209      ELSE                                            !==  Quasi-Eulerian vertical coordinate  ==!   ('key_qco')
210         !                                            !==========================================!
211         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
212            !                                               ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]]
213!!gm old
214!            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)   &
215!               &                            + r1_Dt * (  e3t(:,:,jk,Kaa)      &
216!               &                                       - e3t(:,:,jk,Kbb)  )   ) * tmask(:,:,jk)
217!!gm slightly faster :
218#if defined key_qco
219            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                            &
220               &                            + r1_Dt * e3t_0(:,:,jk) * ( r3t(:,:,Kaa) - r3t(:,:,Kbb) )  ) * tmask(:,:,jk)
221#else
222            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)   &
223               &                            + r1_Dt * (  e3t(:,:,jk,Kaa)      &
224               &                                       - e3t(:,:,jk,Kbb)  )   ) * tmask(:,:,jk)
225#endif
226!!gm end !!st need to be there when not using qco
227
228
229         END DO
230      ENDIF
231
232      IF( ln_bdy ) THEN
233         DO jk = 1, jpkm1
234            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
235         END DO
236      ENDIF
237      !
238#if defined key_agrif
239      IF( .NOT. AGRIF_Root() ) THEN
240         !
241         ! Mask vertical velocity at first/last columns/row
242         ! inside computational domain (cosmetic)
243         DO jk = 1, jpkm1
244            IF( lk_west ) THEN                             ! --- West --- !
245               DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
246                  DO jj = 1, jpj
247                     pww(ji,jj,jk) = 0._wp 
248                  END DO
249               END DO
250            ENDIF
251            IF( lk_east ) THEN                             ! --- East --- !
252               DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
253                  DO jj = 1, jpj
254                     pww(ji,jj,jk) = 0._wp
255                  END DO
256               END DO
257            ENDIF
258            IF( lk_south ) THEN                            ! --- South --- !
259               DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
260                  DO ji = 1, jpi
261                     pww(ji,jj,jk) = 0._wp
262                  END DO
263               END DO
264            ENDIF
265            IF( lk_north ) THEN                            ! --- North --- !
266               DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
267                  DO ji = 1, jpi
268                     pww(ji,jj,jk) = 0._wp
269                  END DO
270               END DO
271            ENDIF
272            !
273         END DO
274         !
275      ENDIF 
276#endif
277      !
278      IF( ln_timing )   CALL timing_stop('wzv_MLF')
279      !
280   END SUBROUTINE wzv_MLF
281
282
283   SUBROUTINE wzv_RK3( kt, Kbb, Kmm, Kaa, puu, pvv, pww )
284      !!----------------------------------------------------------------------
285      !!                ***  ROUTINE wzv_RK3  ***
286      !!                   
287      !! ** Purpose :   compute the now vertical velocity
288      !!
289      !! ** Method  : - Using the incompressibility hypothesis, the vertical
290      !!      velocity is computed by integrating the horizontal divergence 
291      !!      from the bottom to the surface minus the scale factor evolution.
292      !!        The boundary conditions are w=0 at the bottom (no flux) and.
293      !!
294      !! ** action  :   pww      : now vertical velocity
295      !!
296      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
297      !!----------------------------------------------------------------------
298      INTEGER                         , INTENT(in   ) ::   kt             ! time step
299      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level indices
300      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   puu, pvv       ! horizontal velocity at Kmm
301      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::             pww  ! vertical velocity at Kmm
302      !
303      INTEGER  ::   ji, jj, jk   ! dummy loop indices
304      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv
305      REAL(wp)             , DIMENSION(jpi,jpj,jpk) ::   ze3div
306      !!----------------------------------------------------------------------
307      !
308      IF( ln_timing )   CALL timing_start('wzv_RK3')
309      !
310      IF( kt == nit000 ) THEN
311         IF(lwp) WRITE(numout,*)
312         IF(lwp) WRITE(numout,*) 'wzv_RK3 : now vertical velocity '
313         IF(lwp) WRITE(numout,*) '~~~~~ '
314         !
315         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
316      ENDIF
317      !
318      CALL div_hor( kt, Kbb, Kmm, puu, pvv, ze3div )
319      !                                           !------------------------------!
320      !                                           !     Now Vertical Velocity    !
321      !                                           !------------------------------!
322      !
323      !                                               !===============================!
324      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      !==  z_tilde and layer cases  ==!
325         !                                            !===============================!
326         ALLOCATE( zhdiv(jpi,jpj,jpk) ) 
327         !
328         DO jk = 1, jpkm1
329            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
330            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
331            DO_2D( 0, 0, 0, 0 )
332               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) )
333            END_2D
334         END DO
335         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
336         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
337         !                             ! Same question holds for hdiv. Perhaps just for security
338         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
339            ! computation of w
340            pww(:,:,jk) = pww(:,:,jk+1) - (   ze3div(:,:,jk) + zhdiv(:,:,jk)   &
341               &                            + r1_Dt * (  e3t(:,:,jk,Kaa)       &
342               &                                       - e3t(:,:,jk,Kbb) )   ) * tmask(:,:,jk)
343         END DO
344         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
345         DEALLOCATE( zhdiv ) 
346         !                                            !=================================!
347      ELSEIF( ln_linssh )   THEN                      !==  linear free surface cases  ==!
348         !                                            !=================================!
349         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
350            pww(:,:,jk) = pww(:,:,jk+1) - ze3div(:,:,jk)
351         END DO
352         !                                            !==========================================!
353      ELSE                                            !==  Quasi-Eulerian vertical coordinate  ==!   ('key_qco')
354         !                                            !==========================================!
355         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
356            !                                               ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]]
357            pww(:,:,jk) = pww(:,:,jk+1) - (  ze3div(:,:,jk)                            &
358               &                            + r1_Dt * e3t_0(:,:,jk) * ( r3t(:,:,Kaa) - r3t(:,:,Kbb) )  ) * tmask(:,:,jk)
359         END DO
360      ENDIF
361
362      IF( ln_bdy ) THEN
363         DO jk = 1, jpkm1
364            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
365         END DO
366      ENDIF
367      !
368#if defined key_agrif
369      IF( .NOT. AGRIF_Root() ) THEN
370         !
371         ! Mask vertical velocity at first/last columns/row
372         ! inside computational domain (cosmetic)
373         DO jk = 1, jpkm1
374            IF( lk_west ) THEN                             ! --- West --- !
375               DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
376                  DO jj = 1, jpj
377                     pww(ji,jj,jk) = 0._wp 
378                  END DO
379               END DO
380            ENDIF
381            IF( lk_east ) THEN                             ! --- East --- !
382               DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
383                  DO jj = 1, jpj
384                     pww(ji,jj,jk) = 0._wp
385                  END DO
386               END DO
387            ENDIF
388            IF( lk_south ) THEN                            ! --- South --- !
389               DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
390                  DO ji = 1, jpi
391                     pww(ji,jj,jk) = 0._wp
392                  END DO
393               END DO
394            ENDIF
395            IF( lk_north ) THEN                            ! --- North --- !
396               DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
397                  DO ji = 1, jpi
398                     pww(ji,jj,jk) = 0._wp
399                  END DO
400               END DO
401            ENDIF
402            !
403         END DO
404         !
405      ENDIF 
406#endif
407      !
408      IF( ln_timing )   CALL timing_stop('wzv_RK3')
409      !
410   END SUBROUTINE wzv_RK3
411
412
413   SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )
414      !!----------------------------------------------------------------------
415      !!                    ***  ROUTINE ssh_atf  ***
416      !!
417      !! ** Purpose :   Apply Asselin time filter to now SSH.
418      !!
419      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
420      !!              from the filter, see Leclair and Madec 2010) and swap :
421      !!                pssh(:,:,Kmm) = pssh(:,:,Kaa) + rn_atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
422      !!                            - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0
423      !!
424      !! ** action  : - pssh(:,:,Kmm) time filtered
425      !!
426      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
427      !!----------------------------------------------------------------------
428      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index
429      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! ocean time level indices
430      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! SSH field
431      !
432      REAL(wp) ::   zcoef   ! local scalar
433      !!----------------------------------------------------------------------
434      !
435      IF( ln_timing )   CALL timing_start('ssh_atf')
436      !
437      IF( kt == nit000 ) THEN
438         IF(lwp) WRITE(numout,*)
439         IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
440         IF(lwp) WRITE(numout,*) '~~~~~~~ '
441      ENDIF
442      !
443      IF( .NOT.l_1st_euler ) THEN   ! Apply Asselin time filter on Kmm field (not on euler 1st)
444         !
445         IF( ln_linssh ) THEN                ! filtered "now" field
446            pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
447            !
448         ELSE                                ! filtered "now" field with forcing removed
449            zcoef = rn_atfp * rn_Dt * r1_rho0
450            pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )   &
451               &                          - zcoef   * (         emp_b(:,:) -        emp(:,:)   &
452               &                                              - rnf_b(:,:) +        rnf(:,:)   &
453               &                                       + fwfisf_cav_b(:,:) - fwfisf_cav(:,:)   &
454               &                                       + fwfisf_par_b(:,:) - fwfisf_par(:,:)   ) * ssmask(:,:)
455
456            ! ice sheet coupling
457            IF( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1 )   &
458               &   pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0._wp ) * ssmask(:,:)
459
460         ENDIF
461      ENDIF
462      !
463      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' atf  - pssh(:,:,Kmm): ', mask1=tmask )
464      !
465      IF( ln_timing )   CALL timing_stop('ssh_atf')
466      !
467   END SUBROUTINE ssh_atf
468
469   
470   SUBROUTINE wAimp( kt, Kmm )
471      !!----------------------------------------------------------------------
472      !!                ***  ROUTINE wAimp  ***
473      !!                   
474      !! ** Purpose :   compute the Courant number and partition vertical velocity
475      !!                if a proportion needs to be treated implicitly
476      !!
477      !! ** Method  : -
478      !!
479      !! ** action  :   ww      : now vertical velocity (to be handled explicitly)
480      !!            :   wi      : now vertical velocity (for implicit treatment)
481      !!
482      !! Reference  : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
483      !!              implicit scheme for vertical advection in oceanic modeling.
484      !!              Ocean Modelling, 91, 38-69.
485      !!----------------------------------------------------------------------
486      INTEGER, INTENT(in) ::   kt   ! time step
487      INTEGER, INTENT(in) ::   Kmm  ! time level index
488      !
489      INTEGER  ::   ji, jj, jk   ! dummy loop indices
490      REAL(wp)             ::   zCu, zcff, z1_e3t, zdt                ! local scalars
491      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters
492      REAL(wp) , PARAMETER ::   Cu_max = 0.30_wp                      ! local parameters
493      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters
494      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters
495      !!----------------------------------------------------------------------
496      !
497      IF( ln_timing )   CALL timing_start('wAimp')
498      !
499      IF( kt == nit000 ) THEN
500         IF(lwp) WRITE(numout,*)
501         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
502         IF(lwp) WRITE(numout,*) '~~~~~ '
503         wi(:,:,:) = 0._wp
504      ENDIF
505      !
506      ! Calculate Courant numbers
507      zdt = 2._wp * rn_Dt                            ! 2*rn_Dt and not rDt (for restartability)
508      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
509         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
510            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
511            Cu_adv(ji,jj,jk) =   zdt *                                                         &
512               &  ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )            &
513               &  + ( MAX( e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)                                  &
514               &                        * uu (ji  ,jj,jk,Kmm) + un_td(ji  ,jj,jk), 0._wp ) -   &
515               &      MIN( e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)                                  &
516               &                        * uu (ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) )   &
517               &                               * r1_e1e2t(ji,jj)                                                                     &
518               &  + ( MAX( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)                                  &
519               &                        * vv (ji,jj  ,jk,Kmm) + vn_td(ji,jj  ,jk), 0._wp ) -   &
520               &      MIN( e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)                                  &
521               &                        * vv (ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) )   &
522               &                               * r1_e1e2t(ji,jj)                                                                     &
523               &                             ) * z1_e3t
524         END_3D
525      ELSE
526         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
527            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
528            Cu_adv(ji,jj,jk) =   zdt *                                                      &
529               &  ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )         &
530               &                             + ( MAX( e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm), 0._wp ) -   &
531               &                                 MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) )   &
532               &                               * r1_e1e2t(ji,jj)                                                 &
533               &                             + ( MAX( e1v(ji,jj  )*e3v(ji,jj  ,jk,Kmm)*vv(ji,jj  ,jk,Kmm), 0._wp ) -   &
534               &                                 MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) )   &
535               &                               * r1_e1e2t(ji,jj)                                                 &
536               &                             ) * z1_e3t
537         END_3D
538      ENDIF
539      CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1.0_wp )
540      !
541      CALL iom_put("Courant",Cu_adv)
542      !
543      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
544         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )             ! or scan Courant criterion and partition ! w where necessary
545            !
546            zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
547! alt:
548!                  IF ( ww(ji,jj,jk) > 0._wp ) THEN
549!                     zCu =  Cu_adv(ji,jj,jk)
550!                  ELSE
551!                     zCu =  Cu_adv(ji,jj,jk-1)
552!                  ENDIF
553            !
554            IF( zCu <= Cu_min ) THEN              !<-- Fully explicit
555               zcff = 0._wp
556            ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
557               zcff = ( zCu - Cu_min )**2
558               zcff = zcff / ( Fcu + zcff )
559            ELSE                                  !<-- Mostly implicit
560               zcff = ( zCu - Cu_max )/ zCu
561            ENDIF
562            zcff = MIN(1._wp, zcff)
563            !
564            wi(ji,jj,jk) =           zcff   * ww(ji,jj,jk)
565            ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
566            !
567            Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl
568         END_3D
569         Cu_adv(:,:,1) = 0._wp 
570      ELSE
571         ! Fully explicit everywhere
572         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient below and in stp_ctl
573         wi    (:,:,:) = 0._wp
574      ENDIF
575      CALL iom_put("wimp",wi) 
576      CALL iom_put("wi_cff",Cu_adv)
577      CALL iom_put("wexp",ww)
578      !
579      IF( ln_timing )   CALL timing_stop('wAimp')
580      !
581   END SUBROUTINE wAimp
582     
583   !!======================================================================
584END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.