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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 2636

Last change on this file since 2636 was 2636, checked in by gm, 13 years ago

dynamic mem: #785 ; move ctl_stop & warn in lib_mpp to avoid a circular dependency + ctl_stop improvment
  • Property svn:keywords set to Id
File size: 18.4 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   !!----------------------------------------------------------------------
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      USE oce, ONLY :   z3d => ta   ! use ta as 3D workspace
79      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
80      USE wrk_nemo, ONLY: zhdiv => wrk_2d_1, z2d => wrk_2d_2
81      !!
82      INTEGER, INTENT(in) ::   kt   ! time step
83      !!
84      INTEGER  ::   ji, jj, jk      ! dummy loop indices
85      REAL(wp) ::   zcoefu, zcoefv, zcoeff      ! temporary scalars
86      REAL(wp) ::   z2dt, z1_2dt, z1_rau0       ! temporary scalars
87      !!----------------------------------------------------------------------
88
89      IF( wrk_in_use(2, 1,2) ) THEN
90         CALL ctl_stop('ssh_wzv: requested workspace arrays unavailable')   ;   RETURN
91      ENDIF
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.e0                   ! 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.e0 - umask(:,:,1) )
153         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) )
154         !
155      ENDIF
156      !
157      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
158      !
159      z2dt = 2. * 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.e0
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(:,:) )  ) &
174      &                      * tmask(:,:,1)
175
176#if defined key_agrif
177      CALL agrif_ssh(kt)
178#endif
179#if defined key_obc
180      IF( Agrif_Root() ) THEN
181         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
182         CALL lbc_lnk( ssha, 'T', 1. )                ! absolutly compulsory !! (jmm)
183      ENDIF
184#endif
185#if defined key_bdy
186      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
187      CALL lbc_lnk( ssha, 'T', 1. ) 
188#endif
189
190      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
191      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
192         DO jj = 1, jpjm1
193            DO ji = 1, jpim1      ! NO Vector Opt.
194               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
195                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
196                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
197               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
198                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
199                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
200            END DO
201         END DO
202         ! Boundaries conditions
203         CALL lbc_lnk( sshu_a, 'U', 1. )
204         CALL lbc_lnk( sshv_a, 'V', 1. )
205      ENDIF
206! Include the IAU weighted SSH increment
207#if defined key_asminc
208      IF( ( lk_asminc ).AND.( ln_sshinc ).AND.( ln_asmiau ) ) THEN
209         CALL ssh_asm_inc( kt )
210         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
211      ENDIF
212#endif
213
214      !                                           !------------------------------!
215      !                                           !     Now Vertical Velocity    !
216      !                                           !------------------------------!
217      z1_2dt = 1.e0 / z2dt
218      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
219         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
220         wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
221            &                      - (  fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
222            &                         * tmask(:,:,jk) * z1_2dt
223#if defined key_bdy
224         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
225#endif
226      END DO
227
228      !                                           !------------------------------!
229      !                                           !           outputs            !
230      !                                           !------------------------------!
231      CALL iom_put( "woce", wn                    )   ! vertical velocity
232      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
233      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
234      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
235         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
236         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:)
237         DO jk = 1, jpk
238            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
239         END DO
240         CALL iom_put( "w_masstr" , z3d                     ) 
241         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
242      ENDIF
243      !
244      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
245      !
246      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('ssh_wzv: failed to release workspace arrays')
247      !
248   END SUBROUTINE ssh_wzv
249
250
251   SUBROUTINE ssh_nxt( kt )
252      !!----------------------------------------------------------------------
253      !!                    ***  ROUTINE ssh_nxt  ***
254      !!
255      !! ** Purpose :   achieve the sea surface  height time stepping by
256      !!              applying Asselin time filter and swapping the arrays
257      !!              ssha  already computed in ssh_wzv 
258      !!
259      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
260      !!              from the filter, see Leclair and Madec 2010) and swap :
261      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
262      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
263      !!                sshn = ssha
264      !!
265      !! ** action  : - sshb, sshn   : before & now sea surface height
266      !!                               ready for the next time step
267      !!
268      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
269      !!----------------------------------------------------------------------
270      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
271      !!
272      INTEGER  ::   ji, jj   ! dummy loop indices
273      REAL(wp) ::   zec      ! temporary scalar
274      !!----------------------------------------------------------------------
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  !
284         !                    !--------------------------!
285         !
286         ! ssh at t-, u-, v, f-points
287         !===========================
288         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
289            sshn  (:,:) = ssha  (:,:)                        ! now <-- after  (before already = now)
290            sshu_n(:,:) = sshu_a(:,:)
291            sshv_n(:,:) = sshv_a(:,:)
292            DO jj = 1, jpjm1
293               DO ji = 1, jpim1      ! NO Vector Opt.
294                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
295                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
296                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
297                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
298               END DO
299            END DO
300            ! Boundaries conditions
301            CALL lbc_lnk( sshf_n, 'F', 1. )
302         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
303            zec = atfp * rdt / rau0
304            DO jj = 1, jpj
305               DO ji = 1, jpi                                ! before <-- now filtered
306                  sshb  (ji,jj) = sshn  (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )   &
307                     &                          - zec  * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1)
308                  sshn  (ji,jj) = ssha  (ji,jj)              ! now <-- after
309                  sshu_n(ji,jj) = sshu_a(ji,jj)
310                  sshv_n(ji,jj) = sshv_a(ji,jj)
311               END DO
312            END DO
313            DO jj = 1, jpjm1
314               DO ji = 1, jpim1      ! NO Vector Opt.
315                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                 &
316                     &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
317                     &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
318                     &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
319               END DO
320            END DO
321            ! Boundaries conditions
322            CALL lbc_lnk( sshf_n, 'F', 1. )
323            DO jj = 1, jpjm1
324               DO ji = 1, jpim1      ! NO Vector Opt.
325                  sshu_b(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
326                     &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
327                     &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
328                  sshv_b(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
329                     &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
330                     &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
331               END DO
332            END DO
333            ! Boundaries conditions
334            CALL lbc_lnk( sshu_b, 'U', 1. )
335            CALL lbc_lnk( sshv_b, 'V', 1. )
336         ENDIF
337         !                    !--------------------------!
338      ELSE                    !        fixed levels      !
339         !                    !--------------------------!
340         !
341         ! ssh at t-point only
342         !====================
343         IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler time-stepping at first time-step : no filter
344            sshn(:,:) = ssha(:,:)                            ! now <-- after  (before already = now)
345            !
346         ELSE                                           ! Leap-Frog time-stepping: Asselin filter + swap
347            DO jj = 1, jpj
348               DO ji = 1, jpi                                ! before <-- now filtered
349                  sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) )
350                  sshn(ji,jj) = ssha(ji,jj)                  ! now <-- after
351               END DO
352            END DO
353         ENDIF
354         !
355      ENDIF
356      !
357      ! Update velocity at AGRIF zoom boundaries
358#if defined key_agrif
359      IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
360#endif
361      !
362      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
363      !
364   END SUBROUTINE ssh_nxt
365
366   !!======================================================================
367END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.