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

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 2789

Last change on this file since 2789 was 2789, checked in by cetlod, 13 years ago

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