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 @ 15620

Last change on this file since 15620 was 15620, checked in by techene, 2 years ago

#2605 update RK3 loops according to the trunk transformation
  • Property svn:keywords set to Id
File size: 29.3 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  ::   ji, jj, 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_3D( 1, nn_hls, 1, nn_hls, 1, jpkm1 )                                 ! Horizontal divergence of barotropic transports
112        zhdiv(ji,jj) = zhdiv(ji,jj) + e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk)
113      END_3D
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      DO_2D_OVR( 1, nn_hls, 1, nn_hls )                ! Loop bounds limited by hdiv definition in div_hor
119         pssh(ji,jj,Kaa) = (  pssh(ji,jj,Kbb) - rDt * ( zcoef * ( emp_b(ji,jj) + emp(ji,jj) ) + zhdiv(ji,jj) )  ) * ssmask(ji,jj)
120      END_2D
121      ! pssh must be defined everywhere (true for dyn_spg_ts, not for dyn_spg_exp)
122      IF ( .NOT. ln_dynspg_ts .AND. nn_hls == 2 ) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp )
123      !
124#if defined key_agrif
125      Kbb_a = Kbb   ;   Kmm_a = Kmm   ;   Krhs_a = Kaa
126      CALL agrif_ssh( kt )
127#endif
128      !
129      IF ( .NOT.ln_dynspg_ts ) THEN
130         IF( ln_bdy ) THEN
131            IF (nn_hls==1) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp )    ! Not sure that's necessary
132            CALL bdy_ssh( pssh(:,:,Kaa) )             ! Duplicate sea level across open boundaries
133         ENDIF
134      ENDIF
135      !                                           !------------------------------!
136      !                                           !           outputs            !
137      !                                           !------------------------------!
138      !
139      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa)  - : ', mask1=tmask )
140      !
141      IF( ln_timing )   CALL timing_stop('ssh_nxt')
142      !
143   END SUBROUTINE ssh_nxt
144
145
146   SUBROUTINE wzv_MLF( kt, Kbb, Kmm, Kaa, pww )
147      !!----------------------------------------------------------------------
148      !!                ***  ROUTINE wzv_MLF  ***
149      !!                   
150      !! ** Purpose :   compute the now vertical velocity
151      !!
152      !! ** Method  : - Using the incompressibility hypothesis, the vertical
153      !!      velocity is computed by integrating the horizontal divergence 
154      !!      from the bottom to the surface minus the scale factor evolution.
155      !!        The boundary conditions are w=0 at the bottom (no flux) and.
156      !!
157      !! ** action  :   pww      : now vertical velocity
158      !!
159      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
160      !!----------------------------------------------------------------------
161      INTEGER                         , INTENT(in)    ::   kt             ! time step
162      INTEGER                         , INTENT(in)    ::   Kbb, Kmm, Kaa  ! time level indices
163      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pww            ! vertical velocity at Kmm
164      !
165      INTEGER  ::   ji, jj, jk   ! dummy loop indices
166      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv
167      !!----------------------------------------------------------------------
168      !
169      IF( ln_timing )   CALL timing_start('wzv_MLF')
170      !
171      IF( kt == nit000 ) THEN
172         IF(lwp) WRITE(numout,*)
173         IF(lwp) WRITE(numout,*) 'wzv_MLF : now vertical velocity '
174         IF(lwp) WRITE(numout,*) '~~~~~~~'
175         !
176         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
177      ENDIF
178      !                                           !------------------------------!
179      !                                           !     Now Vertical Velocity    !
180      !                                           !------------------------------!
181      !
182      !                                               !===============================!
183      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      !==  z_tilde and layer cases  ==!
184         !                                            !===============================!
185         ALLOCATE( zhdiv(jpi,jpj,jpk) ) 
186         !
187         DO jk = 1, jpkm1
188            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
189            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
190            DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
191               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) )
192            END_2D
193         END DO
194         IF( nn_hls == 1)   CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
195         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
196         !                             ! Same question holds for hdiv. Perhaps just for security
197         !                             ! clem: yes it is a problem because ww is used in many other places where we need the halos
198         !
199         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 )     ! integrate from the bottom the hor. divergence
200            ! computation of w
201            pww(ji,jj,jk) = pww(ji,jj,jk+1) - (   e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk)   &
202               &                                  +                  zhdiv(ji,jj,jk)   &
203               &                                  + r1_Dt * (  e3t(ji,jj,jk,Kaa)       &
204               &                                             - e3t(ji,jj,jk,Kbb) )   ) * tmask(ji,jj,jk)
205         END_3D
206         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
207         DEALLOCATE( zhdiv )
208         !                                            !=================================!
209      ELSEIF( ln_linssh )   THEN                      !==  linear free surface cases  ==!
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)  ) * tmask(ji,jj,jk)
213         END_3D
214         !                                            !==========================================!
215      ELSE                                            !==  Quasi-Eulerian vertical coordinate  ==!   ('key_qco')
216         !                                            !==========================================!
217         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 )     ! integrate from the bottom the hor. divergence
218#if defined key_qco
219!!gm slightly faster
220            pww(ji,jj,jk) = pww(ji,jj,jk+1) - (  e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk)    &
221                 &                               + r1_Dt * e3t_0(ji,jj,jk) * ( r3t(ji,jj,Kaa) - r3t(ji,jj,Kbb) )  ) * tmask(ji,jj,jk)
222#else
223            pww(ji,jj,jk) = pww(ji,jj,jk+1) - (  e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk)    &
224               &                                 + r1_Dt * (  e3t(ji,jj,jk,Kaa)        &
225               &                                            - e3t(ji,jj,jk,Kbb)  )   ) * tmask(ji,jj,jk)
226#endif
227         END_3D
228      ENDIF
229
230      IF( ln_bdy ) THEN
231         DO jk = 1, jpkm1
232            DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
233               pww(ji,jj,jk) = pww(ji,jj,jk) * bdytmask(ji,jj)
234            END_2D
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( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) 
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         IF( nn_hls == 1)   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_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 )     ! integrate from the bottom the hor. divergence
339            pww(ji,jj,jk) = pww(ji,jj,jk+1) - (   ze3div(ji,jj,jk) + zhdiv(ji,jj,jk)   &
340                 &                            + r1_Dt * (  e3t(ji,jj,jk,Kaa)       &
341                 &                                       - e3t(ji,jj,jk,Kbb) )   ) * tmask(ji,jj,jk)
342         END_3D
343         !
344         DEALLOCATE( zhdiv ) 
345         !                                            !=================================!
346      ELSEIF( ln_linssh )   THEN                      !==  linear free surface cases  ==!
347         !                                            !=================================!
348         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 )     ! integrate from the bottom the hor. divergence
349            pww(ji,jj,jk) = pww(ji,jj,jk+1) - ze3div(ji,jj,jk) 
350         END_3D
351         !                                            !==========================================!
352      ELSE                                            !==  Quasi-Eulerian vertical coordinate  ==!   ('key_qco')
353         !                                            !==========================================!
354         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 )     ! integrate from the bottom the hor. divergence
355         !                                                              ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]]
356            pww(ji,jj,jk) = pww(ji,jj,jk+1) - (  ze3div(ji,jj,jk)                            &
357               &                               + r1_Dt * e3t_0(ji,jj,jk) * ( r3t(ji,jj,Kaa) - r3t(ji,jj,Kbb) )  ) * tmask(ji,jj,jk)
358         END_3D
359      ENDIF
360
361      IF( ln_bdy ) THEN
362         DO jk = 1, jpkm1
363            DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
364               pww(ji,jj,jk) = pww(ji,jj,jk) * bdytmask(ji,jj)
365            END_2D
366         END DO
367      ENDIF
368      !
369#if defined key_agrif
370      IF( .NOT. AGRIF_Root() ) THEN
371         !
372         ! Mask vertical velocity at first/last columns/row
373         ! inside computational domain (cosmetic)
374         DO jk = 1, jpkm1
375            IF( lk_west ) THEN                             ! --- West --- !
376               DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
377                  DO jj = 1, jpj
378                     pww(ji,jj,jk) = 0._wp 
379                  END DO
380               END DO
381            ENDIF
382            IF( lk_east ) THEN                             ! --- East --- !
383               DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
384                  DO jj = 1, jpj
385                     pww(ji,jj,jk) = 0._wp
386                  END DO
387               END DO
388            ENDIF
389            IF( lk_south ) THEN                            ! --- South --- !
390               DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
391                  DO ji = 1, jpi
392                     pww(ji,jj,jk) = 0._wp
393                  END DO
394               END DO
395            ENDIF
396            IF( lk_north ) THEN                            ! --- North --- !
397               DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
398                  DO ji = 1, jpi
399                     pww(ji,jj,jk) = 0._wp
400                  END DO
401               END DO
402            ENDIF
403            !
404         END DO
405         !
406      ENDIF 
407#endif
408      !
409      IF( ln_timing )   CALL timing_stop('wzv_RK3')
410      !
411   END SUBROUTINE wzv_RK3
412
413
414   SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )
415      !!----------------------------------------------------------------------
416      !!                    ***  ROUTINE ssh_atf  ***
417      !!
418      !! ** Purpose :   Apply Asselin time filter to now SSH.
419      !!
420      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
421      !!              from the filter, see Leclair and Madec 2010) and swap :
422      !!                pssh(:,:,Kmm) = pssh(:,:,Kaa) + rn_atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
423      !!                            - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0
424      !!
425      !! ** action  : - pssh(:,:,Kmm) time filtered
426      !!
427      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
428      !!----------------------------------------------------------------------
429      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index
430      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! ocean time level indices
431      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! SSH field
432      !
433      REAL(wp) ::   zcoef   ! local scalar
434      !!----------------------------------------------------------------------
435      !
436      IF( ln_timing )   CALL timing_start('ssh_atf')
437      !
438      IF( kt == nit000 ) THEN
439         IF(lwp) WRITE(numout,*)
440         IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
441         IF(lwp) WRITE(numout,*) '~~~~~~~ '
442      ENDIF
443      !
444      IF( .NOT.l_1st_euler ) THEN   ! Apply Asselin time filter on Kmm field (not on euler 1st)
445         !
446         IF( ln_linssh ) THEN                ! filtered "now" field
447            pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
448            !
449         ELSE                                ! filtered "now" field with forcing removed
450            zcoef = rn_atfp * rn_Dt * r1_rho0
451            pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )   &
452               &                          - zcoef   * (         emp_b(:,:) -        emp(:,:)   &
453               &                                              - rnf_b(:,:) +        rnf(:,:)   &
454               &                                       - fwfisf_cav_b(:,:) + fwfisf_cav(:,:)   &
455               &                                       - fwfisf_par_b(:,:) + fwfisf_par(:,:)   ) * ssmask(:,:)
456
457            ! ice sheet coupling
458            IF( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1 )   &
459               &   pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0._wp ) * ssmask(:,:)
460
461         ENDIF
462      ENDIF
463      !
464      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' atf  - pssh(:,:,Kmm): ', mask1=tmask )
465      !
466      IF( ln_timing )   CALL timing_stop('ssh_atf')
467      !
468   END SUBROUTINE ssh_atf
469
470   
471   SUBROUTINE wAimp( kt, Kmm )
472      !!----------------------------------------------------------------------
473      !!                ***  ROUTINE wAimp  ***
474      !!                   
475      !! ** Purpose :   compute the Courant number and partition vertical velocity
476      !!                if a proportion needs to be treated implicitly
477      !!
478      !! ** Method  : -
479      !!
480      !! ** action  :   ww      : now vertical velocity (to be handled explicitly)
481      !!            :   wi      : now vertical velocity (for implicit treatment)
482      !!
483      !! Reference  : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
484      !!              implicit scheme for vertical advection in oceanic modeling.
485      !!              Ocean Modelling, 91, 38-69.
486      !!----------------------------------------------------------------------
487      INTEGER, INTENT(in) ::   kt   ! time step
488      INTEGER, INTENT(in) ::   Kmm  ! time level index
489      !
490      INTEGER  ::   ji, jj, jk   ! dummy loop indices
491      REAL(wp)             ::   zCu, zcff, z1_e3t, zdt                ! local scalars
492      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters
493      REAL(wp) , PARAMETER ::   Cu_max = 0.30_wp                      ! local parameters
494      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters
495      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters
496      !!----------------------------------------------------------------------
497      !
498      IF( ln_timing )   CALL timing_start('wAimp')
499      !
500      IF( kt == nit000 ) THEN
501         IF(lwp) WRITE(numout,*)
502         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
503         IF(lwp) WRITE(numout,*) '~~~~~ '
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( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 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( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 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 iom_put("Courant",Cu_adv)
540      !
541      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
542         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 2, -1 )    ! or scan Courant criterion and partition ! w where necessary
543            !
544            zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
545! alt:
546!                  IF ( ww(ji,jj,jk) > 0._wp ) THEN
547!                     zCu =  Cu_adv(ji,jj,jk)
548!                  ELSE
549!                     zCu =  Cu_adv(ji,jj,jk-1)
550!                  ENDIF
551            !
552            IF( zCu <= Cu_min ) THEN              !<-- Fully explicit
553               zcff = 0._wp
554            ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
555               zcff = ( zCu - Cu_min )**2
556               zcff = zcff / ( Fcu + zcff )
557            ELSE                                  !<-- Mostly implicit
558               zcff = ( zCu - Cu_max )/ zCu
559            ENDIF
560            zcff = MIN(1._wp, zcff)
561            !
562            wi(ji,jj,jk) =           zcff   * ww(ji,jj,jk)
563            ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
564            !
565            Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl
566         END_3D
567         Cu_adv(:,:,1) = 0._wp 
568      ELSE
569         ! Fully explicit everywhere
570         Cu_adv(:,:,:) = 0._wp                    ! Reuse array to output coefficient below and in stp_ctl
571         wi    (:,:,:) = 0._wp
572      ENDIF
573      CALL iom_put("wimp",wi) 
574      CALL iom_put("wi_cff",Cu_adv)
575      CALL iom_put("wexp",ww)
576      !
577      IF( ln_timing )   CALL timing_stop('wAimp')
578      !
579   END SUBROUTINE wAimp
580     
581   !!======================================================================
582END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.