source: branches/UKMO/dev_r5518_GO6_package_inc_asm/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 8190

Last change on this file since 8190 was 8190, checked in by jwhile, 3 years ago

Update due to Tim's review

File size: 13.2 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 divcur          ! hor. divergence and curl      (div & cur routines)
23   USE restart         ! only for lrst_oce
24   USE in_out_manager  ! I/O manager
25   USE prtctl          ! Print control
26   USE phycst
27   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
28   USE lib_mpp         ! MPP library
29   USE bdy_oce
30   USE bdy_par         
31   USE bdydyn2d        ! bdy_ssh routine
32#if defined key_agrif
33   USE agrif_opa_interp
34#endif
35#if defined key_asminc   
36   USE asminc          ! Assimilation increment
37#endif
38   USE wrk_nemo        ! Memory Allocation
39   USE timing          ! Timing
40
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC   ssh_nxt    ! called by step.F90
45   PUBLIC   wzv        ! called by step.F90
46   PUBLIC   ssh_swp    ! called by step.F90
47
48   !! * Substitutions
49#  include "domzgr_substitute.h90"
50#  include "vectopt_loop_substitute.h90"
51   !!----------------------------------------------------------------------
52   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
53   !! $Id$
54   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
55   !!----------------------------------------------------------------------
56CONTAINS
57
58   SUBROUTINE ssh_nxt( kt )
59      !!----------------------------------------------------------------------
60      !!                ***  ROUTINE ssh_nxt  ***
61      !!                   
62      !! ** Purpose :   compute the after ssh (ssha)
63      !!
64      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
65      !!      is computed by integrating the horizontal divergence and multiply by
66      !!      by the time step.
67      !!
68      !! ** action  :   ssha    : after sea surface height
69      !!
70      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
71      !!----------------------------------------------------------------------
72      !
73      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv
74      INTEGER, INTENT(in) ::   kt                      ! time step
75      !
76      INTEGER             ::   ji, jj, jk   ! dummy loop indices
77      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
78      !!----------------------------------------------------------------------
79      !
80      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
81      !
82      CALL wrk_alloc( jpi, jpj, zhdiv ) 
83      !
84      IF( kt == nit000 ) THEN
85         !
86         IF(lwp) WRITE(numout,*)
87         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
88         IF(lwp) WRITE(numout,*) '~~~~~~~ '
89         !
90      ENDIF
91      !
92      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
93      !
94      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
95      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
96
97
98#if defined key_asminc
99      !                                                ! Include the IAU weighted SSH increment
100      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
101         CALL ssh_asm_inc( kt )
102#if defined key_vvl
103! Don't directly adjust ssh but change hdivn at all levels instead
104! In trasbc also add in the heat and salt content associated with these changes at each level 
105        DO jk = 1, jpkm1                                 
106                 hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) ) * ( e3t_0(:,:,jk) / fse3t_n(:,:,jk) ) * tmask(:,:,jk) 
107        END DO
108      ENDIF
109#endif
110#endif
111
112
113      !                                           !------------------------------!
114      !                                           !   After Sea Surface Height   !
115      !                                           !------------------------------!
116      zhdiv(:,:) = 0._wp
117      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
118        zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
119      END DO
120      !                                                ! Sea surface elevation time stepping
121      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
122      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
123      !
124      z1_rau0 = 0.5_wp * r1_rau0
125      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
126
127#if ! defined key_dynspg_ts
128      ! These lines are not necessary with time splitting since
129      ! boundary condition on sea level is set during ts loop
130#if defined key_agrif
131      CALL agrif_ssh( kt )
132#endif
133#if defined key_bdy
134      IF (lk_bdy) THEN
135         CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
136         CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
137      ENDIF
138#endif
139#endif
140
141
142      !                                           !------------------------------!
143      !                                           !           outputs            !
144      !                                           !------------------------------!
145      !
146      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
147      !
148      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
149      !
150      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
151      !
152   END SUBROUTINE ssh_nxt
153
154   
155   SUBROUTINE wzv( kt )
156      !!----------------------------------------------------------------------
157      !!                ***  ROUTINE wzv  ***
158      !!                   
159      !! ** Purpose :   compute the now vertical velocity
160      !!
161      !! ** Method  : - Using the incompressibility hypothesis, the vertical
162      !!      velocity is computed by integrating the horizontal divergence 
163      !!      from the bottom to the surface minus the scale factor evolution.
164      !!        The boundary conditions are w=0 at the bottom (no flux) and.
165      !!
166      !! ** action  :   wn      : now vertical velocity
167      !!
168      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
169      !!----------------------------------------------------------------------
170      !
171      INTEGER, INTENT(in) ::   kt           ! time step
172      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
173      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
174      !
175      INTEGER             ::   ji, jj, jk   ! dummy loop indices
176      REAL(wp)            ::   z1_2dt       ! local scalars
177      !!----------------------------------------------------------------------
178     
179      IF( nn_timing == 1 )  CALL timing_start('wzv')
180      !
181      IF( kt == nit000 ) THEN
182         !
183         IF(lwp) WRITE(numout,*)
184         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
185         IF(lwp) WRITE(numout,*) '~~~~~ '
186         !
187         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
188         !
189      ENDIF
190      !                                           !------------------------------!
191      !                                           !     Now Vertical Velocity    !
192      !                                           !------------------------------!
193      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
194      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
195      !
196      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
197         CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
198         !
199         DO jk = 1, jpkm1
200            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
201            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
202            DO jj = 2, jpjm1
203               DO ji = fs_2, fs_jpim1   ! vector opt.
204                  zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
205               END DO
206            END DO
207         END DO
208         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
209         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
210         !                             ! Same question holds for hdivn. Perhaps just for security
211         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
212            ! computation of w
213            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    &
214               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
215         END DO
216         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
217         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 
218      ELSE   ! z_star and linear free surface cases
219         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
220            ! computation of w
221            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   &
222               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
223         END DO
224      ENDIF
225
226#if defined key_bdy
227      IF (lk_bdy) THEN
228         DO jk = 1, jpkm1
229            wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
230         END DO
231      ENDIF
232#endif
233      !
234      IF( nn_timing == 1 )  CALL timing_stop('wzv')
235
236
237   END SUBROUTINE wzv
238
239   SUBROUTINE ssh_swp( kt )
240      !!----------------------------------------------------------------------
241      !!                    ***  ROUTINE ssh_nxt  ***
242      !!
243      !! ** Purpose :   achieve the sea surface  height time stepping by
244      !!              applying Asselin time filter and swapping the arrays
245      !!              ssha  already computed in ssh_nxt 
246      !!
247      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
248      !!              from the filter, see Leclair and Madec 2010) and swap :
249      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
250      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
251      !!                sshn = ssha
252      !!
253      !! ** action  : - sshb, sshn   : before & now sea surface height
254      !!                               ready for the next time step
255      !!
256      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
257      !!----------------------------------------------------------------------
258      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
259      !!----------------------------------------------------------------------
260      !
261      IF( nn_timing == 1 )  CALL timing_start('ssh_swp')
262      !
263      IF( kt == nit000 ) THEN
264         IF(lwp) WRITE(numout,*)
265         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
266         IF(lwp) WRITE(numout,*) '~~~~~~~ '
267      ENDIF
268
269# if defined key_dynspg_ts
270      IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN    !** Euler time-stepping: no filter
271# else
272      IF ( neuler == 0 .AND. kt == nit000 ) THEN   !** Euler time-stepping at first time-step : no filter
273#endif
274         sshb(:,:) = sshn(:,:)                           ! before <-- now
275         sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now)
276      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
277         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
278         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:)    - emp(:,:)    &
279                                &                                 - rnf_b(:,:)    + rnf(:,:)    &
280                                &                                 + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:)
281         sshn(:,:) = ssha(:,:)                           ! now <-- after
282      ENDIF
283      !
284      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
285      !
286      IF( nn_timing == 1 )  CALL timing_stop('ssh_swp')
287      !
288   END SUBROUTINE ssh_swp
289
290   !!======================================================================
291END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.