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 branches/2017/dev_r8789_sbc/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2017/dev_r8789_sbc/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 8795

Last change on this file since 8795 was 8795, checked in by charris, 7 years ago

Fix typo.

  • Property svn:keywords set to Id
File size: 13.7 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
25   USE bdydyn2d       ! bdy_ssh routine
26#if defined key_agrif
27   USE agrif_opa_interp
28#endif
29#if defined key_asminc   
30   USE   asminc       ! Assimilation increment
31#endif
32   !
33   USE in_out_manager ! I/O manager
34   USE restart        ! only for lrst_oce
35   USE prtctl         ! Print control
36   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
37   USE lib_mpp        ! MPP library
38   USE wrk_nemo       ! Memory Allocation
39   USE timing         ! Timing
40   USE wet_dry         ! Wetting/Drying flux limting
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   ssh_nxt    ! called by step.F90
46   PUBLIC   wzv        ! called by step.F90
47   PUBLIC   ssh_swp    ! called by step.F90
48
49   !! * Substitutions
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      INTEGER, INTENT(in) ::   kt   ! time step
73      !
74      INTEGER  ::   jk            ! dummy loop indice
75      REAL(wp) ::   z2dt, zcoef   ! local scalars
76      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhdiv   ! 2D workspace
77      !!----------------------------------------------------------------------
78      !
79      IF( nn_timing == 1 )   CALL timing_start('ssh_nxt')
80      !
81      CALL wrk_alloc( jpi,jpj,   zhdiv ) 
82      !
83      IF( kt == nit000 ) THEN
84         IF(lwp) WRITE(numout,*)
85         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
86         IF(lwp) WRITE(numout,*) '~~~~~~~ '
87      ENDIF
88      !
89      z2dt = 2._wp * rdt                          ! set time step size (Euler/Leapfrog)
90      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
91      zcoef = 0.5_wp * r1_rau0
92
93      !                                           !------------------------------!
94      !                                           !   After Sea Surface Height   !
95      !                                           !------------------------------!
96      IF(ln_wd) THEN
97         CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)
98      ENDIF
99
100      CALL div_hor( kt )                               ! Horizontal divergence
101      !
102#if defined key_asminc
103      !                                                ! Include the IAU weighted SSH increment
104      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
105         CALL ssh_asm_inc( kt )
106         IF ( ln_linssh ) THEN
107            ! Adjust sea surface height directly
108            ! (note that qns/sfx are not modified to take into account heat/salt content added as effect would
109            ! just need to be removed by concentration / dilution effect in trasbc)
110            ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
111         ELSE
112            ! Don't directly adjust sea surface height but change horizontal divergence at all levels instead
113            ! (in trasbc the heat and salt content associated with these changes are also taken into account)
114            DO jk = 1, jpkm1                                 
115               hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) )   &
116                                               * ( e3t_0(:,:,jk) / e3t_n(:,:,jk) ) * tmask(:,:,jk) 
117            END DO
118            CALL lbc_lnk( hdivn, 'T', 1. )
119         ENDIF
120      ENDIF
121      !
122#endif
123      !
124      zhdiv(:,:) = 0._wp
125      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
126        zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
127      END DO
128      !                                                ! Sea surface elevation time stepping
129      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
130      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
131      !
132      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
133
134      IF ( .NOT.ln_dynspg_ts ) THEN
135         ! These lines are not necessary with time splitting since
136         ! boundary condition on sea level is set during ts loop
137# if defined key_agrif
138         CALL agrif_ssh( kt )
139# endif
140         IF( ln_bdy ) THEN
141            CALL lbc_lnk( ssha, 'T', 1. )    ! Not sure that's necessary
142            CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries
143         ENDIF
144      ENDIF
145      !                                           !------------------------------!
146      !                                           !           outputs            !
147      !                                           !------------------------------!
148      !
149      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
150      !
151      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
152      !
153      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
154      !
155   END SUBROUTINE ssh_nxt
156
157   
158   SUBROUTINE wzv( kt )
159      !!----------------------------------------------------------------------
160      !!                ***  ROUTINE wzv  ***
161      !!                   
162      !! ** Purpose :   compute the now vertical velocity
163      !!
164      !! ** Method  : - Using the incompressibility hypothesis, the vertical
165      !!      velocity is computed by integrating the horizontal divergence 
166      !!      from the bottom to the surface minus the scale factor evolution.
167      !!        The boundary conditions are w=0 at the bottom (no flux) and.
168      !!
169      !! ** action  :   wn      : now vertical velocity
170      !!
171      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
172      !!----------------------------------------------------------------------
173      INTEGER, INTENT(in) ::   kt   ! time step
174      !
175      INTEGER  ::   ji, jj, jk   ! dummy loop indices
176      REAL(wp) ::   z1_2dt       ! local scalars
177      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
178      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
179      !!----------------------------------------------------------------------
180      !
181      IF( nn_timing == 1 )   CALL timing_start('wzv')
182      !
183      IF( kt == nit000 ) THEN
184         IF(lwp) WRITE(numout,*)
185         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
186         IF(lwp) WRITE(numout,*) '~~~~~ '
187         !
188         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
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_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) )
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) - (  e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)    &
214               &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_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) - (  e3t_n(:,:,jk) * hdivn(:,:,jk)                 &
222               &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )  ) * tmask(:,:,jk)
223         END DO
224      ENDIF
225
226      IF( ln_bdy ) THEN
227         DO jk = 1, jpkm1
228            wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
229         END DO
230      ENDIF
231      !
232      IF( nn_timing == 1 )  CALL timing_stop('wzv')
233      !
234   END SUBROUTINE wzv
235
236
237   SUBROUTINE ssh_swp( kt )
238      !!----------------------------------------------------------------------
239      !!                    ***  ROUTINE ssh_nxt  ***
240      !!
241      !! ** Purpose :   achieve the sea surface  height time stepping by
242      !!              applying Asselin time filter and swapping the arrays
243      !!              ssha  already computed in ssh_nxt 
244      !!
245      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
246      !!              from the filter, see Leclair and Madec 2010) and swap :
247      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
248      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
249      !!                sshn = ssha
250      !!
251      !! ** action  : - sshb, sshn   : before & now sea surface height
252      !!                               ready for the next time step
253      !!
254      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
255      !!----------------------------------------------------------------------
256      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
257      !
258      REAL(wp) ::   zcoef   ! local scalar
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      !              !==  Euler time-stepping: no filter, just swap  ==!
269      IF(  ( neuler == 0 .AND. kt == nit000 ) .OR.    &
270         & ( ln_bt_fw    .AND. ln_dynspg_ts )      ) THEN
271         sshb(:,:) = sshn(:,:)                              ! before <-- now
272         sshn(:,:) = ssha(:,:)                              ! now    <-- after  (before already = now)
273         !
274      ELSE           !==  Leap-Frog time-stepping: Asselin filter + swap  ==!
275         !                                                  ! before <-- now filtered
276         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )
277         IF( .NOT.ln_linssh ) THEN                          ! before <-- with forcing removed
278            zcoef = atfp * rdt * r1_rau0
279            sshb(:,:) = sshb(:,:) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
280               &                             -    rnf_b(:,:) + rnf   (:,:)   &
281               &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:)
282         ENDIF
283         sshn(:,:) = ssha(:,:)                              ! now <-- after
284      ENDIF
285      !
286      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
287      !
288      IF( nn_timing == 1 )   CALL timing_stop('ssh_swp')
289      !
290   END SUBROUTINE ssh_swp
291
292   !!======================================================================
293END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.