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

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 2797

Last change on this file since 2797 was 2797, checked in by davestorkey, 13 years ago

Delete BDY module and first implementation of new OBC module.

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