source: NEMO/branches/UKMO/NEMO_4.0_GO8_coupled_iodef/src/OCE/DYN/sshwzv.F90 @ 11355

Last change on this file since 11355 was 11355, checked in by dancopsey, 2 years ago

Add lots of print statements printing out values of files throughout the ocean_ice timestep.. New output will be in files crash_stats_419.out and crash_stats_68.out.

File size: 18.0 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   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   ssh_nxt       : after ssh
15   !!   ssh_swp       : filter ans swap the ssh arrays
16   !!   wzv           : compute now vertical velocity
17   !!----------------------------------------------------------------------
18   USE oce            ! ocean dynamics and tracers variables
19   USE dom_oce        ! ocean space and time domain variables
20   USE sbc_oce        ! surface boundary condition: ocean
21   USE domvvl         ! Variable volume
22   USE divhor         ! horizontal divergence
23   USE phycst         ! physical constants
24   USE bdy_oce , ONLY : ln_bdy, bdytmask   ! Open BounDarY
25   USE bdydyn2d       ! bdy_ssh routine
26#if defined key_agrif
27   USE agrif_oce_interp
28#endif
29   !
30   USE iom 
31   USE in_out_manager ! I/O manager
32   USE restart        ! only for lrst_oce
33   USE prtctl         ! Print control
34   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
35   USE lib_mpp        ! MPP library
36   USE timing         ! Timing
37   USE wet_dry        ! Wetting/Drying flux limiting
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   ssh_nxt    ! called by step.F90
43   PUBLIC   wzv        ! called by step.F90
44   PUBLIC   wAimp      ! called by step.F90
45   PUBLIC   ssh_swp    ! called by step.F90
46
47   !! * Substitutions
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
51   !! $Id$
52   !! Software governed by the CeCILL license (see ./LICENSE)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE ssh_nxt( kt )
57      !!----------------------------------------------------------------------
58      !!                ***  ROUTINE ssh_nxt  ***
59      !!                   
60      !! ** Purpose :   compute the after ssh (ssha)
61      !!
62      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
63      !!      is computed by integrating the horizontal divergence and multiply by
64      !!      by the time step.
65      !!
66      !! ** action  :   ssha, after sea surface height
67      !!
68      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
69      !!----------------------------------------------------------------------
70      INTEGER, INTENT(in) ::   kt   ! time step
71      !
72      INTEGER  ::   jk            ! dummy loop indice
73      REAL(wp) ::   z2dt, zcoef   ! local scalars
74      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace
75      !!----------------------------------------------------------------------
76      !
77      IF( ln_timing )   CALL timing_start('ssh_nxt')
78      !
79      IF( kt == nit000 ) THEN
80         IF(lwp) WRITE(numout,*)
81         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
82         IF(lwp) WRITE(numout,*) '~~~~~~~ '
83      ENDIF
84      !
85      z2dt = 2._wp * rdt                          ! set time step size (Euler/Leapfrog)
86      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
87      zcoef = 0.5_wp * r1_rau0
88
89      !                                           !------------------------------!
90      !                                           !   After Sea Surface Height   !
91      !                                           !------------------------------!
92      IF(ln_wd_il) THEN
93         CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)
94      ENDIF
95
96      CALL div_hor( kt )                               ! Horizontal divergence
97      !
98      zhdiv(:,:) = 0._wp
99      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
100        zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
101      END DO
102      !                                                ! Sea surface elevation time stepping
103      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
104      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
105
106      IF(narea == 419) THEN
107         WRITE(9419,*) 'max ssha ssh_nxt before maths = ',MAXVAL(  ABS( ssha(:,:) )  )
108         WRITE(9419,*) 'max sshb ssh_nxt before maths = ',MAXVAL(  ABS( sshb(:,:) )  )
109         WRITE(9419,*) 'max emp_b ssh_nxt before maths = ',MAXVAL(  ABS( emp_b(:,:) )  )
110         WRITE(9419,*) 'max emp ssh_nxt before maths = ',MAXVAL(  ABS( emp(:,:) )  )
111         WRITE(9419,*) 'max zhdiv ssh_nxt before maths = ',MAXVAL(  ABS( zhdiv(:,:) )  )
112         WRITE(9419,*) 'max ssmask ssh_nxt before maths = ',MAXVAL(  ABS( ssmask(:,:) )  )
113      ENDIF
114
115      !
116      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
117
118      IF(narea == 419) THEN
119         WRITE(9419,*) 'max ssha ssh_nxt after maths = ',MAXVAL(  ABS( ssha(:,:) )  )
120      ENDIF
121      !
122#if defined key_agrif
123      CALL agrif_ssh( kt )
124#endif
125      !
126      IF ( .NOT.ln_dynspg_ts ) THEN
127         IF( ln_bdy ) THEN
128            CALL lbc_lnk( 'sshwzv', ssha, 'T', 1. )    ! Not sure that's necessary
129            CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries
130         ENDIF
131      ENDIF
132
133      IF(narea == 419) THEN
134         WRITE(9419,*) 'max ssha ssh_nxt after bdy_ssh = ',MAXVAL(  ABS( ssha(:,:) )  )
135      ENDIF
136
137      !                                           !------------------------------!
138      !                                           !           outputs            !
139      !                                           !------------------------------!
140      !
141      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask )
142      !
143
144      IF(narea == 419) THEN
145         WRITE(9419,*) 'max ssha ssh_nxt after prt_ctl = ',MAXVAL(  ABS( ssha(:,:) )  )
146      ENDIF
147
148      IF( ln_timing )   CALL timing_stop('ssh_nxt')
149      !
150   END SUBROUTINE ssh_nxt
151
152   
153   SUBROUTINE wzv( kt )
154      !!----------------------------------------------------------------------
155      !!                ***  ROUTINE wzv  ***
156      !!                   
157      !! ** Purpose :   compute the now vertical velocity
158      !!
159      !! ** Method  : - Using the incompressibility hypothesis, the vertical
160      !!      velocity is computed by integrating the horizontal divergence 
161      !!      from the bottom to the surface minus the scale factor evolution.
162      !!        The boundary conditions are w=0 at the bottom (no flux) and.
163      !!
164      !! ** action  :   wn      : now vertical velocity
165      !!
166      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
167      !!----------------------------------------------------------------------
168      INTEGER, INTENT(in) ::   kt   ! time step
169      !
170      INTEGER  ::   ji, jj, jk   ! dummy loop indices
171      REAL(wp) ::   z1_2dt       ! local scalars
172      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv
173      !!----------------------------------------------------------------------
174      !
175      IF( ln_timing )   CALL timing_start('wzv')
176      !
177      IF( kt == nit000 ) THEN
178         IF(lwp) WRITE(numout,*)
179         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
180         IF(lwp) WRITE(numout,*) '~~~~~ '
181         !
182         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
183      ENDIF
184      !                                           !------------------------------!
185      !                                           !     Now Vertical Velocity    !
186      !                                           !------------------------------!
187      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
188      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
189      !
190      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
191         ALLOCATE( zhdiv(jpi,jpj,jpk) ) 
192         !
193         DO jk = 1, jpkm1
194            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
195            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
196            DO jj = 2, jpjm1
197               DO ji = fs_2, fs_jpim1   ! vector opt.
198                  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) )
199               END DO
200            END DO
201         END DO
202         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
203         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
204         !                             ! Same question holds for hdivn. Perhaps just for security
205         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
206            ! computation of w
207            wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)    &
208               &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )     ) * tmask(:,:,jk)
209         END DO
210         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
211         DEALLOCATE( zhdiv ) 
212      ELSE   ! z_star and linear free surface cases
213         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
214            ! computation of w
215            wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk)                 &
216               &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )  ) * tmask(:,:,jk)
217         END DO
218      ENDIF
219
220      IF( ln_bdy ) THEN
221         DO jk = 1, jpkm1
222            wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
223         END DO
224      ENDIF
225      !
226#if defined key_agrif 
227      IF( .NOT. AGRIF_Root() ) THEN
228         IF ((nbondi ==  1).OR.(nbondi == 2)) wn(nlci-1 , :     ,:) = 0.e0      ! east
229         IF ((nbondi == -1).OR.(nbondi == 2)) wn(2      , :     ,:) = 0.e0      ! west
230         IF ((nbondj ==  1).OR.(nbondj == 2)) wn(:      ,nlcj-1 ,:) = 0.e0      ! north
231         IF ((nbondj == -1).OR.(nbondj == 2)) wn(:      ,2      ,:) = 0.e0      ! south
232      ENDIF 
233#endif 
234      !
235      IF( ln_timing )   CALL timing_stop('wzv')
236      !
237   END SUBROUTINE wzv
238
239
240   SUBROUTINE ssh_swp( kt )
241      !!----------------------------------------------------------------------
242      !!                    ***  ROUTINE ssh_nxt  ***
243      !!
244      !! ** Purpose :   achieve the sea surface  height time stepping by
245      !!              applying Asselin time filter and swapping the arrays
246      !!              ssha  already computed in ssh_nxt 
247      !!
248      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
249      !!              from the filter, see Leclair and Madec 2010) and swap :
250      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
251      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
252      !!                sshn = ssha
253      !!
254      !! ** action  : - sshb, sshn   : before & now sea surface height
255      !!                               ready for the next time step
256      !!
257      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
258      !!----------------------------------------------------------------------
259      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
260      !
261      REAL(wp) ::   zcoef   ! local scalar
262      !!----------------------------------------------------------------------
263      !
264      IF( ln_timing )   CALL timing_start('ssh_swp')
265      !
266      IF( kt == nit000 ) THEN
267         IF(lwp) WRITE(numout,*)
268         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
269         IF(lwp) WRITE(numout,*) '~~~~~~~ '
270      ENDIF
271      !              !==  Euler time-stepping: no filter, just swap  ==!
272
273      IF(narea == 419) THEN
274         WRITE(9419,*) 'max sshn in sshwzv = ',MAXVAL(  ABS( sshn(:,:) )  )
275         WRITE(9419,*) 'max ssha in sshwzv = ',MAXVAL(  ABS( ssha(:,:) )  )
276      ENDIF
277
278      IF ( neuler == 0 .AND. kt == nit000 ) THEN
279         sshn(:,:) = ssha(:,:)                              ! now    <-- after  (before already = now)
280         !
281      ELSE           !==  Leap-Frog time-stepping: Asselin filter + swap  ==!
282         !                                                  ! before <-- now filtered
283         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )
284         IF( .NOT.ln_linssh ) THEN                          ! before <-- with forcing removed
285            zcoef = atfp * rdt * r1_rau0
286            sshb(:,:) = sshb(:,:) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
287               &                             -    rnf_b(:,:) + rnf   (:,:)   &
288               &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:)
289         ENDIF
290         sshn(:,:) = ssha(:,:)                              ! now <-- after
291      ENDIF
292      !
293      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask )
294      !
295      IF( ln_timing )   CALL timing_stop('ssh_swp')
296      !
297   END SUBROUTINE ssh_swp
298
299   SUBROUTINE wAimp( kt )
300      !!----------------------------------------------------------------------
301      !!                ***  ROUTINE wAimp  ***
302      !!                   
303      !! ** Purpose :   compute the Courant number and partition vertical velocity
304      !!                if a proportion needs to be treated implicitly
305      !!
306      !! ** Method  : -
307      !!
308      !! ** action  :   wn      : now vertical velocity (to be handled explicitly)
309      !!            :   wi      : now vertical velocity (for implicit treatment)
310      !!
311      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
312      !!----------------------------------------------------------------------
313      INTEGER, INTENT(in) ::   kt   ! time step
314      !
315      INTEGER  ::   ji, jj, jk   ! dummy loop indices
316      REAL(wp)             ::   zCu, zcff, z1_e3w                     ! local scalars
317      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters
318      REAL(wp) , PARAMETER ::   Cu_max = 0.27                         ! local parameters
319      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters
320      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters
321      !!----------------------------------------------------------------------
322      !
323      IF( ln_timing )   CALL timing_start('wAimp')
324      !
325      IF( kt == nit000 ) THEN
326         IF(lwp) WRITE(numout,*)
327         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
328         IF(lwp) WRITE(numout,*) '~~~~~ '
329      ENDIF
330      !
331      !
332      DO jk = 1, jpkm1            ! calculate Courant numbers
333         DO jj = 2, jpjm1
334            DO ji = 2, fs_jpim1   ! vector opt.
335               z1_e3w = 1._wp / e3w_n(ji,jj,jk)
336               Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    &  ! 2*rdt and not r2dt (for restartability)
337                  &                             + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   &
338                  &                                 MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   &
339                  &                               * r1_e1e2t(ji,jj)                                                  &
340                  &                             + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   &
341                  &                                 MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   &
342                  &                               * r1_e1e2t(ji,jj)                                                  &
343                  &                             ) * z1_e3w
344            END DO
345         END DO
346      END DO
347      CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. )
348      !
349      CALL iom_put("Courant",Cu_adv)
350      !
351      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
352         DO jk = 1, jpkm1                               ! or scan Courant criterion and partition
353            DO jj = 1, jpj                              ! w where necessary
354               DO ji = 1, jpi
355                  !
356                  zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) )
357                  !
358                  IF( zCu <= Cu_min ) THEN              !<-- Fully explicit
359                     zcff = 0._wp
360                  ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
361                     zcff = ( zCu - Cu_min )**2
362                     zcff = zcff / ( Fcu + zcff )
363                  ELSE                                  !<-- Mostly implicit
364                     zcff = ( zCu - Cu_max )/ zCu
365                  ENDIF
366                  zcff = MIN(1._wp, zcff)
367                  !
368                  wi(ji,jj,jk) =           zcff   * wn(ji,jj,jk)
369                  wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk)
370                  !
371                  Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient
372               END DO
373            END DO
374         END DO
375      ELSE
376         ! Fully explicit everywhere
377         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient
378         wi    (:,:,:) = 0._wp
379      ENDIF
380      CALL iom_put("wimp",wi) 
381      CALL iom_put("wi_cff",Cu_adv)
382      CALL iom_put("wexp",wn)
383      !
384      IF( ln_timing )   CALL timing_stop('wAimp')
385      !
386   END SUBROUTINE wAimp
387   !!======================================================================
388END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.