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/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 3186

Last change on this file since 3186 was 3186, checked in by smasson, 12 years ago

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