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

source: branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 6736

Last change on this file since 6736 was 6736, checked in by jamesharle, 8 years ago

FASTNEt code modifications

  • 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      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      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d, zhdiv
84      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d
85      !!----------------------------------------------------------------------
86      !
87      IF( nn_timing == 1 )  CALL timing_start('ssh_wzv')
88      !
89      CALL wrk_alloc( jpi, jpj, z2d, zhdiv ) 
90      !
91      IF( kt == nit000 ) THEN
92         !
93         IF(lwp) WRITE(numout,*)
94         IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity '
95         IF(lwp) WRITE(numout,*) '~~~~~~~ '
96         !
97         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
98         !
99         IF( lk_vvl ) THEN                    ! before and now Sea SSH at u-, v-, f-points (vvl case only)
100            DO jj = 1, jpjm1
101               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
102                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
103                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
104                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
105                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
106                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
107                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
108                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
109                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
110                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
111                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
112                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
113               END DO
114            END DO
115            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
116            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
117            DO jj = 1, jpjm1
118               DO ji = 1, jpim1      ! NO Vector Opt.
119                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
120                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
121                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
122                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
123               END DO
124            END DO
125            CALL lbc_lnk( sshf_n, 'F', 1. )
126         ENDIF
127         !
128      ENDIF
129
130      !                                           !------------------------------------------!
131      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
132         !                                        !------------------------------------------!
133         DO jk = 1, jpkm1
134            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
135            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
136            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
137            !
138            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
139            fse3u (:,:,jk) = fse3u_n (:,:,jk)
140            fse3v (:,:,jk) = fse3v_n (:,:,jk)
141            fse3f (:,:,jk) = fse3f_n (:,:,jk)
142            fse3w (:,:,jk) = fse3w_n (:,:,jk)
143            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
144            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
145         END DO
146         !
147         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
148         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
149         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
150         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )
151         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) )
152         !
153      ENDIF
154      !
155      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
156      !
157      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
158      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
159
160      !                                           !------------------------------!
161      !                                           !   After Sea Surface Height   !
162      !                                           !------------------------------!
163      zhdiv(:,:) = 0._wp
164      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
165        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
166      END DO
167      !                                                ! Sea surface elevation time stepping
168      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
169      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
170      z1_rau0 = 0.5 / rau0
171      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
172
173#if defined key_agrif
174      CALL agrif_ssh( kt )
175#endif
176#if defined key_obc
177      IF( Agrif_Root() ) THEN
178         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
179         CALL lbc_lnk( ssha, 'T', 1. )                 ! absolutly compulsory !! (jmm)
180      ENDIF
181#endif
182#if defined key_bdy
183      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
184      CALL lbc_lnk( ssha, 'T', 1. )                    ! absolutly compulsory !! (jmm)
185#endif
186#if defined key_asminc
187      !                                                ! Include the IAU weighted SSH increment
188      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
189         CALL ssh_asm_inc( kt )
190         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
191      ENDIF
192#endif
193      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
194      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
195         DO jj = 1, jpjm1
196            DO ji = 1, jpim1      ! NO Vector Opt.
197               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
198                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
199                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
200               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
201                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
202                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
203            END DO
204         END DO
205         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
206      ENDIF
207
208      !                                           !------------------------------!
209      !                                           !     Now Vertical Velocity    !
210      !                                           !------------------------------!
211      z1_2dt = 1.e0 / z2dt
212      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
213         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
214         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
215            &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
216            &                         * tmask(:,:,jk) * z1_2dt
217#if defined key_bdy
218         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
219#endif
220      END DO
221
222      !                                           !------------------------------!
223      !                                           !           outputs            !
224      !                                           !------------------------------!
225      CALL iom_put( "woce", wn                    )   ! vertical velocity
226      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
227      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
228      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
229         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
230         CALL wrk_alloc( jpi,jpj,jpk, z3d )
231         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
232         DO jk = 1, jpk
233            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
234         END DO
235         CALL iom_put( "w_masstr" , z3d                     ) 
236         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
237         CALL wrk_dealloc( jpi,jpj,jpk, z3d )
238      ENDIF
239      !
240      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
241      !
242      CALL wrk_dealloc( jpi, jpj, z2d, zhdiv ) 
243      !
244      IF( nn_timing == 1 )  CALL timing_stop('ssh_wzv')
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( nn_timing == 1 )  CALL timing_start('ssh_nxt')
275      !
276      IF( kt == nit000 ) THEN
277         IF(lwp) WRITE(numout,*)
278         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)'
279         IF(lwp) WRITE(numout,*) '~~~~~~~ '
280      ENDIF
281
282      !                       !--------------------------!
283      IF( lk_vvl ) THEN       !  Variable volume levels  !     (ssh at t-, u-, v, f-points)
284         !                    !--------------------------!
285         !
286         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
287            sshn  (:,:) = ssha  (:,:)                       ! now <-- after  (before already = now)
288            sshu_n(:,:) = sshu_a(:,:)
289            sshv_n(:,:) = sshv_a(:,:)
290            DO jj = 1, jpjm1                                ! ssh now at f-point
291               DO ji = 1, jpim1      ! NO Vector Opt.
292                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
293                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
294                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
295                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
296               END DO
297            END DO
298            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
299            !
300         ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
301            zec = atfp * rdt / rau0
302            DO jj = 1, jpj
303               DO ji = 1, jpi                               ! before <-- now filtered
304                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
305                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
306                  sshn  (ji,jj) = ssha  (ji,jj)             ! now <-- after
307                  sshu_n(ji,jj) = sshu_a(ji,jj)
308                  sshv_n(ji,jj) = sshv_a(ji,jj)
309               END DO
310            END DO
311            DO jj = 1, jpjm1                                ! ssh now at f-point
312               DO ji = 1, jpim1      ! NO Vector Opt.
313                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
314                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
315                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
316                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
317               END DO
318            END DO
319            CALL lbc_lnk( sshf_n, 'F', 1. )                 ! Boundaries conditions
320            !
321            DO jj = 1, jpjm1                                ! ssh before at u- & v-points
322               DO ji = 1, jpim1      ! NO Vector Opt.
323                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
324                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
325                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
326                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
327                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
328                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
329               END DO
330            END DO
331            CALL lbc_lnk( sshu_b, 'U', 1. )
332            CALL lbc_lnk( sshv_b, 'V', 1. )            !  Boundaries conditions
333            !
334         ENDIF
335         !                    !--------------------------!
336      ELSE                    !        fixed levels      !     (ssh at t-point only)
337         !                    !--------------------------!
338         !
339         IF( neuler == 0 .AND. kt == nit000 ) THEN    !** Euler time-stepping at first time-step : no filter
340            sshn(:,:) = ssha(:,:)                           ! now <-- after  (before already = now)
341            !
342         ELSE                                               ! Leap-Frog time-stepping: Asselin filter + swap
343            DO jj = 1, jpj
344               DO ji = 1, jpi                               ! before <-- now filtered
345                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
346                  sshn(ji,jj) = ssha(ji,jj)                 ! now <-- after
347               END DO
348            END DO
349         ENDIF
350         !
351      ENDIF
352      !
353      ! Update velocity at AGRIF zoom boundaries
354#if defined key_agrif
355      IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
356#endif
357      !
358      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
359      !
360      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
361      !
362   END SUBROUTINE ssh_nxt
363
364   !!======================================================================
365END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.