source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 12 months ago

The Dr Hook changes from my perl code.

File size: 14.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   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   ssh_nxt        : after ssh
15   !!   ssh_swp        : filter ans swap the ssh arrays
16   !!   wzv            : compute now vertical velocity
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers variables
19   USE dom_oce         ! ocean space and time domain variables
20   USE sbc_oce         ! surface boundary condition: ocean
21   USE domvvl          ! Variable volume
22   USE divcur          ! hor. divergence and curl      (div & cur routines)
23   USE restart         ! only for lrst_oce
24   USE in_out_manager  ! I/O manager
25   USE prtctl          ! Print control
26   USE phycst
27   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
28   USE lib_mpp         ! MPP library
29   USE bdy_oce
30   USE bdy_par         
31   USE bdydyn2d        ! bdy_ssh routine
32#if defined key_agrif
33   USE agrif_opa_interp
34#endif
35#if defined key_asminc   
36   USE asminc          ! Assimilation increment
37#endif
38   USE wrk_nemo        ! Memory Allocation
39   USE timing          ! Timing
40
41   USE yomhook, ONLY: lhook, dr_hook
42   USE parkind1, ONLY: jprb, jpim
43
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC   ssh_nxt    ! called by step.F90
48   PUBLIC   wzv        ! called by step.F90
49   PUBLIC   ssh_swp    ! called by step.F90
50
51   !! * Substitutions
52#  include "domzgr_substitute.h90"
53#  include "vectopt_loop_substitute.h90"
54   !!----------------------------------------------------------------------
55   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
56   !! $Id$
57   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE ssh_nxt( kt )
62      !!----------------------------------------------------------------------
63      !!                ***  ROUTINE ssh_nxt  ***
64      !!                   
65      !! ** Purpose :   compute the after ssh (ssha)
66      !!
67      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
68      !!      is computed by integrating the horizontal divergence and multiply by
69      !!      by the time step.
70      !!
71      !! ** action  :   ssha    : after sea surface height
72      !!
73      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
74      !!----------------------------------------------------------------------
75      !
76      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv
77      INTEGER, INTENT(in) ::   kt                      ! time step
78      !
79      INTEGER             ::   jk                      ! dummy loop indices
80      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
81      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
82      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
83      REAL(KIND=jprb)               :: zhook_handle
84
85      CHARACTER(LEN=*), PARAMETER :: RoutineName='SSH_NXT'
86
87      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
88
89      !!----------------------------------------------------------------------
90      !
91      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
92      !
93      CALL wrk_alloc( jpi, jpj, zhdiv ) 
94      !
95      IF( kt == nit000 ) THEN
96         !
97         IF(lwp) WRITE(numout,*)
98         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
99         IF(lwp) WRITE(numout,*) '~~~~~~~ '
100         !
101      ENDIF
102      !
103      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
104      !
105      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
106      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
107
108
109#if defined key_asminc
110      !                                                ! Include the IAU weighted SSH increment
111      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
112         CALL ssh_asm_inc( kt )
113#if defined key_vvl
114! Don't directly adjust ssh but change hdivn at all levels instead
115! In trasbc also add in the heat and salt content associated with these changes at each level 
116        DO jk = 1, jpkm1                                 
117                 hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) ) * ( e3t_0(:,:,jk) / fse3t_n(:,:,jk) ) * tmask(:,:,jk) 
118        END DO
119      ENDIF
120#endif
121#endif
122
123
124      !                                           !------------------------------!
125      !                                           !   After Sea Surface Height   !
126      !                                           !------------------------------!
127      zhdiv(:,:) = 0._wp
128      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
129        zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
130      END DO
131      !                                                ! Sea surface elevation time stepping
132      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
133      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
134      !
135      z1_rau0 = 0.5_wp * r1_rau0
136      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
137
138#if ! defined key_dynspg_ts
139      ! These lines are not necessary with time splitting since
140      ! boundary condition on sea level is set during ts loop
141#if defined key_agrif
142      CALL agrif_ssh( kt )
143#endif
144#if defined key_bdy
145      IF (lk_bdy) THEN
146         CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
147         CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
148      ENDIF
149#endif
150#endif
151
152
153      !                                           !------------------------------!
154      !                                           !           outputs            !
155      !                                           !------------------------------!
156      !
157      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
158      !
159      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
160      !
161      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
162      !
163      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
164   END SUBROUTINE ssh_nxt
165
166   
167   SUBROUTINE wzv( kt )
168      !!----------------------------------------------------------------------
169      !!                ***  ROUTINE wzv  ***
170      !!                   
171      !! ** Purpose :   compute the now vertical velocity
172      !!
173      !! ** Method  : - Using the incompressibility hypothesis, the vertical
174      !!      velocity is computed by integrating the horizontal divergence 
175      !!      from the bottom to the surface minus the scale factor evolution.
176      !!        The boundary conditions are w=0 at the bottom (no flux) and.
177      !!
178      !! ** action  :   wn      : now vertical velocity
179      !!
180      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
181      !!----------------------------------------------------------------------
182      !
183      INTEGER, INTENT(in) ::   kt           ! time step
184      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
185      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
186      !
187      INTEGER             ::   ji, jj, jk   ! dummy loop indices
188      REAL(wp)            ::   z1_2dt       ! local scalars
189      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
190      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
191      REAL(KIND=jprb)               :: zhook_handle
192
193      CHARACTER(LEN=*), PARAMETER :: RoutineName='WZV'
194
195      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
196
197      !!----------------------------------------------------------------------
198     
199      IF( nn_timing == 1 )  CALL timing_start('wzv')
200      !
201      IF( kt == nit000 ) THEN
202         !
203         IF(lwp) WRITE(numout,*)
204         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
205         IF(lwp) WRITE(numout,*) '~~~~~ '
206         !
207         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
208         !
209      ENDIF
210      !                                           !------------------------------!
211      !                                           !     Now Vertical Velocity    !
212      !                                           !------------------------------!
213      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
214      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
215      !
216      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
217         CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
218         !
219         DO jk = 1, jpkm1
220            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
221            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
222            DO jj = 2, jpjm1
223               DO ji = fs_2, fs_jpim1   ! vector opt.
224                  zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
225               END DO
226            END DO
227         END DO
228         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
229         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
230         !                             ! Same question holds for hdivn. Perhaps just for security
231         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
232            ! computation of w
233            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    &
234               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
235         END DO
236         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
237         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 
238      ELSE   ! z_star and linear free surface cases
239         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
240            ! computation of w
241            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   &
242               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
243         END DO
244      ENDIF
245
246#if defined key_bdy
247      IF (lk_bdy) THEN
248         DO jk = 1, jpkm1
249            wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
250         END DO
251      ENDIF
252#endif
253      !
254      IF( nn_timing == 1 )  CALL timing_stop('wzv')
255
256
257      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
258   END SUBROUTINE wzv
259
260   SUBROUTINE ssh_swp( kt )
261      !!----------------------------------------------------------------------
262      !!                    ***  ROUTINE ssh_nxt  ***
263      !!
264      !! ** Purpose :   achieve the sea surface  height time stepping by
265      !!              applying Asselin time filter and swapping the arrays
266      !!              ssha  already computed in ssh_nxt 
267      !!
268      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
269      !!              from the filter, see Leclair and Madec 2010) and swap :
270      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
271      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
272      !!                sshn = ssha
273      !!
274      !! ** action  : - sshb, sshn   : before & now sea surface height
275      !!                               ready for the next time step
276      !!
277      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
278      !!----------------------------------------------------------------------
279      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
280      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
281      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
282      REAL(KIND=jprb)               :: zhook_handle
283
284      CHARACTER(LEN=*), PARAMETER :: RoutineName='SSH_SWP'
285
286      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
287
288      !!----------------------------------------------------------------------
289      !
290      IF( nn_timing == 1 )  CALL timing_start('ssh_swp')
291      !
292      IF( kt == nit000 ) THEN
293         IF(lwp) WRITE(numout,*)
294         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
295         IF(lwp) WRITE(numout,*) '~~~~~~~ '
296      ENDIF
297
298# if defined key_dynspg_ts
299      IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN    !** Euler time-stepping: no filter
300# else
301      IF ( neuler == 0 .AND. kt == nit000 ) THEN   !** Euler time-stepping at first time-step : no filter
302#endif
303         sshb(:,:) = sshn(:,:)                           ! before <-- now
304         sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now)
305      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
306         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
307         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:)    - emp(:,:)    &
308                                &                                 - rnf_b(:,:)    + rnf(:,:)    &
309                                &                                 + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:)
310         sshn(:,:) = ssha(:,:)                           ! now <-- after
311      ENDIF
312      !
313      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
314      !
315      IF( nn_timing == 1 )  CALL timing_stop('ssh_swp')
316      !
317      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
318   END SUBROUTINE ssh_swp
319
320   !!======================================================================
321END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.