source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/DYN/sshwzv.F90 @ 13220

Last change on this file since 13220 was 13220, checked in by orioltp, 3 months ago

dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation: updating from trunk r13218

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