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/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 4147

Last change on this file since 4147 was 3764, checked in by smasson, 11 years ago

dev_MERGE_2012: report bugfixes done in the trunk from r3555 to r3763 into dev_MERGE_2012
  • Property svn:keywords set to Id
File size: 18.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   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   ssh_wzv        : after ssh & now vertical velocity
14   !!   ssh_nxt        : filter ans swap the ssh arrays
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers variables
17   USE dom_oce         ! ocean space and time domain variables
18   USE sbc_oce         ! surface boundary condition: ocean
19   USE domvvl          ! Variable volume
20   USE divcur          ! hor. divergence and curl      (div & cur routines)
21   USE iom             ! I/O library
22   USE in_out_manager  ! I/O manager
23   USE prtctl          ! Print control
24   USE phycst
25   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
26   USE lib_mpp         ! MPP library
27   USE obc_par         ! open boundary cond. parameter
28   USE obc_oce
29   USE bdy_oce
30   USE diaar5, ONLY:   lk_diaar5
31   USE iom
32   USE sbcrnf, ONLY: h_rnf, nk_rnf   ! River runoff
33#if defined key_agrif
34   USE agrif_opa_update
35   USE agrif_opa_interp
36#endif
37#if defined key_asminc   
38   USE asminc          ! Assimilation increment
39#endif
40   USE wrk_nemo        ! Memory Allocation
41   USE timing          ! Timing
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC   ssh_wzv    ! called by step.F90
47   PUBLIC   ssh_nxt    ! called by step.F90
48
49   !! * Substitutions
50#  include "domzgr_substitute.h90"
51#  include "vectopt_loop_substitute.h90"
52   !!----------------------------------------------------------------------
53   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
54   !! $Id$
55   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   SUBROUTINE ssh_wzv( kt ) 
60      !!----------------------------------------------------------------------
61      !!                ***  ROUTINE ssh_wzv  ***
62      !!                   
63      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
64      !!              and update the now vertical coordinate (lk_vvl=T).
65      !!
66      !! ** Method  : - Using the incompressibility hypothesis, the vertical
67      !!      velocity is computed by integrating the horizontal divergence 
68      !!      from the bottom to the surface minus the scale factor evolution.
69      !!        The boundary conditions are w=0 at the bottom (no flux) and.
70      !!
71      !! ** action  :   ssha    : after sea surface height
72      !!                wn      : now vertical velocity
73      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T)
74      !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points
75      !!
76      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
77      !!----------------------------------------------------------------------
78      INTEGER, INTENT(in) ::   kt   ! time step
79      !
80      INTEGER  ::   ji, jj, jk   ! dummy loop indices
81      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
82      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d, zhdiv
83      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d
84      !!----------------------------------------------------------------------
85      !
86      IF( nn_timing == 1 )  CALL timing_start('ssh_wzv')
87      !
88      CALL wrk_alloc( jpi, jpj, z2d, zhdiv ) 
89      !
90      IF( kt == nit000 ) THEN
91         !
92         IF(lwp) WRITE(numout,*)
93         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
94         IF(lwp) WRITE(numout,*) '~~~~~~~ '
95         !
96         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
97         !
98         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
99            DO jj = 1, jpjm1
100               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
101                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
102                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
103                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
104                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
105                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
106                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
107                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
108                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
109                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
110                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
111                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
112               END DO
113            END DO
114            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
115            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
116            DO jj = 1, jpjm1
117               DO ji = 1, jpim1      ! NO Vector Opt.
118                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
119                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
120                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
121                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
122               END DO
123            END DO
124            CALL lbc_lnk( sshf_n, 'F', 1. )
125         ENDIF
126         !
127      ENDIF
128
129      !                                           !------------------------------------------!
130      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
131         !                                        !------------------------------------------!
132         DO jk = 1, jpkm1
133            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
134            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
135            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
136            !
137            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
138            fse3u (:,:,jk) = fse3u_n (:,:,jk)
139            fse3v (:,:,jk) = fse3v_n (:,:,jk)
140            fse3f (:,:,jk) = fse3f_n (:,:,jk)
141            fse3w (:,:,jk) = fse3w_n (:,:,jk)
142            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
143            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
144         END DO
145         !
146         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
147         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
148         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
149         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )
150         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) )
151         !
152      ENDIF
153      !
154      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
155      !
156      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
157      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
158
159      !                                           !------------------------------!
160      !                                           !   After Sea Surface Height   !
161      !                                           !------------------------------!
162      zhdiv(:,:) = 0._wp
163      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
164        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
165      END DO
166      !                                                ! Sea surface elevation time stepping
167      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
168      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
169      z1_rau0 = 0.5 / rau0
170      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
171
172#if defined key_agrif
173      CALL agrif_ssh( kt )
174#endif
175#if defined key_obc
176      IF( Agrif_Root() ) THEN
177         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
178         CALL lbc_lnk( ssha, 'T', 1. )                 ! absolutly compulsory !! (jmm)
179      ENDIF
180#endif
181#if defined key_bdy
182      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
183      CALL lbc_lnk( ssha, 'T', 1. )                    ! absolutly compulsory !! (jmm)
184#endif
185#if defined key_asminc
186      !                                                ! Include the IAU weighted SSH increment
187      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
188         CALL ssh_asm_inc( kt )
189         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
190      ENDIF
191#endif
192      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
193      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
194         DO jj = 1, jpjm1
195            DO ji = 1, jpim1      ! NO Vector Opt.
196               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
197                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
198                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
199               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
200                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
201                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
202            END DO
203         END DO
204         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
205      ENDIF
206
207      !                                           !------------------------------!
208      !                                           !     Now Vertical Velocity    !
209      !                                           !------------------------------!
210      z1_2dt = 1.e0 / z2dt
211      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
212         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
213         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
214            &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
215            &                         * tmask(:,:,jk) * z1_2dt
216#if defined key_bdy
217         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
218#endif
219      END DO
220
221      !                                           !------------------------------!
222      !                                           !           outputs            !
223      !                                           !------------------------------!
224      CALL iom_put( "woce", wn                    )   ! vertical velocity
225      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
226      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
227      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
228         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
229         CALL wrk_alloc( jpi,jpj,jpk, z3d )
230         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
231         DO jk = 1, jpk
232            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
233         END DO
234         CALL iom_put( "w_masstr" , z3d                     ) 
235         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
236         CALL wrk_dealloc( jpi,jpj,jpk, z3d )
237      ENDIF
238      !
239      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
240      !
241      CALL wrk_dealloc( jpi, jpj, z2d, zhdiv ) 
242      !
243      IF( nn_timing == 1 )  CALL timing_stop('ssh_wzv')
244      !
245   END SUBROUTINE ssh_wzv
246
247
248   SUBROUTINE ssh_nxt( kt )
249      !!----------------------------------------------------------------------
250      !!                    ***  ROUTINE ssh_nxt  ***
251      !!
252      !! ** Purpose :   achieve the sea surface  height time stepping by
253      !!              applying Asselin time filter and swapping the arrays
254      !!              ssha  already computed in ssh_wzv 
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      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
259      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
260      !!                sshn = ssha
261      !!
262      !! ** action  : - sshb, sshn   : before & now sea surface height
263      !!                               ready for the next time step
264      !!
265      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
266      !!----------------------------------------------------------------------
267      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
268      !!
269      INTEGER  ::   ji, jj   ! dummy loop indices
270      REAL(wp) ::   zec      ! temporary scalar
271      !!----------------------------------------------------------------------
272      !
273      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
274      !
275      IF( kt == nit000 ) THEN
276         IF(lwp) WRITE(numout,*)
277         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
278         IF(lwp) WRITE(numout,*) '~~~~~~~ '
279      ENDIF
280
281      !                       !--------------------------!
282      IF( lk_vvl ) THEN       !  Variable volume levels  !     (ssh at t-, u-, v, f-points)
283         !                    !--------------------------!
284         !
285         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
286            sshn  (:,:) = ssha  (:,:)                       ! now <-- after  (before already = now)
287            sshu_n(:,:) = sshu_a(:,:)
288            sshv_n(:,:) = sshv_a(:,:)
289            DO jj = 1, jpjm1                                ! ssh now at f-point
290               DO ji = 1, jpim1      ! NO Vector Opt.
291                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
292                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
293                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
294                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
295               END DO
296            END DO
297            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
298            !
299         ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
300            zec = atfp * rdt / rau0
301            DO jj = 1, jpj
302               DO ji = 1, jpi                               ! before <-- now filtered
303                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
304                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
305                  sshn  (ji,jj) = ssha  (ji,jj)             ! now <-- after
306                  sshu_n(ji,jj) = sshu_a(ji,jj)
307                  sshv_n(ji,jj) = sshv_a(ji,jj)
308               END DO
309            END DO
310            DO jj = 1, jpjm1                                ! ssh now at f-point
311               DO ji = 1, jpim1      ! NO Vector Opt.
312                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
313                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
314                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
315                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
316               END DO
317            END DO
318            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
319            !
320            DO jj = 1, jpjm1                                ! ssh before at u- & v-points
321               DO ji = 1, jpim1      ! NO Vector Opt.
322                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
323                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
324                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
325                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
326                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
327                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
328               END DO
329            END DO
330            CALL lbc_lnk( sshu_b, 'U', 1. )
331            CALL lbc_lnk( sshv_b, 'V', 1. )            !  Boundaries conditions
332            !
333         ENDIF
334         !                    !--------------------------!
335      ELSE                    !        fixed levels      !     (ssh at t-point only)
336         !                    !--------------------------!
337         !
338         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
339            sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now)
340            !
341         ELSE                                               ! Leap-Frog time-stepping: Asselin filter + swap
342            DO jj = 1, jpj
343               DO ji = 1, jpi                               ! before <-- now filtered
344                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
345                  sshn(ji,jj) = ssha(ji,jj)                 ! now <-- after
346               END DO
347            END DO
348         ENDIF
349         !
350      ENDIF
351      !
352      ! Update velocity at AGRIF zoom boundaries
353#if defined key_agrif
354      IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
355#endif
356      !
357      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
358      !
359      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
360      !
361   END SUBROUTINE ssh_nxt
362
363   !!======================================================================
364END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.