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

source: branches/2012/dev_LOCEAN_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 3653

Last change on this file since 3653 was 3653, checked in by cetlod, 11 years ago

commit the changes from LOCEAN & UKMO merge, see ticket #1021
  • 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
186      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
187      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
188         DO jj = 1, jpjm1
189            DO ji = 1, jpim1      ! NO Vector Opt.
190               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
191                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
192                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
193               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
194                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
195                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
196            END DO
197         END DO
198         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
199      ENDIF
200     
201#if defined key_asminc
202      !                                                ! Include the IAU weighted SSH increment
203      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
204         CALL ssh_asm_inc( kt )
205         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
206      ENDIF
207#endif
208
209      !                                           !------------------------------!
210      !                                           !     Now Vertical Velocity    !
211      !                                           !------------------------------!
212      z1_2dt = 1.e0 / z2dt
213      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
214         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
215         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
216            &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
217            &                         * tmask(:,:,jk) * z1_2dt
218#if defined key_bdy
219         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
220#endif
221      END DO
222
223      !                                           !------------------------------!
224      !                                           !           outputs            !
225      !                                           !------------------------------!
226      CALL iom_put( "woce", wn                    )   ! vertical velocity
227      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
228      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
229      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
230         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
231         CALL wrk_alloc( jpi,jpj,jpk, z3d )
232         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
233         DO jk = 1, jpk
234            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
235         END DO
236         CALL iom_put( "w_masstr" , z3d                     ) 
237         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
238         CALL wrk_dealloc( jpi,jpj,jpk, z3d )
239      ENDIF
240      !
241      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
242      !
243      CALL wrk_dealloc( jpi, jpj, z2d, zhdiv ) 
244      !
245      IF( nn_timing == 1 )  CALL timing_stop('ssh_wzv')
246      !
247   END SUBROUTINE ssh_wzv
248
249
250   SUBROUTINE ssh_nxt( kt )
251      !!----------------------------------------------------------------------
252      !!                    ***  ROUTINE ssh_nxt  ***
253      !!
254      !! ** Purpose :   achieve the sea surface  height time stepping by
255      !!              applying Asselin time filter and swapping the arrays
256      !!              ssha  already computed in ssh_wzv 
257      !!
258      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
259      !!              from the filter, see Leclair and Madec 2010) and swap :
260      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
261      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
262      !!                sshn = ssha
263      !!
264      !! ** action  : - sshb, sshn   : before & now sea surface height
265      !!                               ready for the next time step
266      !!
267      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
268      !!----------------------------------------------------------------------
269      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
270      !!
271      INTEGER  ::   ji, jj   ! dummy loop indices
272      REAL(wp) ::   zec      ! temporary scalar
273      !!----------------------------------------------------------------------
274      !
275      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
276      !
277      IF( kt == nit000 ) THEN
278         IF(lwp) WRITE(numout,*)
279         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
280         IF(lwp) WRITE(numout,*) '~~~~~~~ '
281      ENDIF
282
283      !                       !--------------------------!
284      IF( lk_vvl ) THEN       !  Variable volume levels  !     (ssh at t-, u-, v, f-points)
285         !                    !--------------------------!
286         !
287         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
288            sshn  (:,:) = ssha  (:,:)                       ! now <-- after  (before already = now)
289            sshu_n(:,:) = sshu_a(:,:)
290            sshv_n(:,:) = sshv_a(:,:)
291            DO jj = 1, jpjm1                                ! ssh now at f-point
292               DO ji = 1, jpim1      ! NO Vector Opt.
293                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
294                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
295                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
296                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
297               END DO
298            END DO
299            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
300            !
301         ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
302            zec = atfp * rdt / rau0
303            DO jj = 1, jpj
304               DO ji = 1, jpi                               ! before <-- now filtered
305                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
306                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
307                  sshn  (ji,jj) = ssha  (ji,jj)             ! now <-- after
308                  sshu_n(ji,jj) = sshu_a(ji,jj)
309                  sshv_n(ji,jj) = sshv_a(ji,jj)
310               END DO
311            END DO
312            DO jj = 1, jpjm1                                ! ssh now at f-point
313               DO ji = 1, jpim1      ! NO Vector Opt.
314                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
315                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
316                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
317                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
318               END DO
319            END DO
320            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
321            !
322            DO jj = 1, jpjm1                                ! ssh before at u- & v-points
323               DO ji = 1, jpim1      ! NO Vector Opt.
324                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
325                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
326                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
327                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
328                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
329                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
330               END DO
331            END DO
332            CALL lbc_lnk( sshu_b, 'U', 1. )
333            CALL lbc_lnk( sshv_b, 'V', 1. )            !  Boundaries conditions
334            !
335         ENDIF
336         !                    !--------------------------!
337      ELSE                    !        fixed levels      !     (ssh at t-point only)
338         !                    !--------------------------!
339         !
340         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
341            sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now)
342            !
343         ELSE                                               ! Leap-Frog time-stepping: Asselin filter + swap
344            DO jj = 1, jpj
345               DO ji = 1, jpi                               ! before <-- now filtered
346                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
347                  sshn(ji,jj) = ssha(ji,jj)                 ! now <-- after
348               END DO
349            END DO
350         ENDIF
351         !
352      ENDIF
353      !
354      ! Update velocity at AGRIF zoom boundaries
355#if defined key_agrif
356      IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
357#endif
358      !
359      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
360      !
361      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
362      !
363   END SUBROUTINE ssh_nxt
364
365   !!======================================================================
366END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.