source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 20 months ago

GMED 450 add flush after prints

File size: 13.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   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC   ssh_nxt    ! called by step.F90
45   PUBLIC   wzv        ! called by step.F90
46   PUBLIC   ssh_swp    ! 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   !!----------------------------------------------------------------------
56CONTAINS
57
58   SUBROUTINE ssh_nxt( kt )
59      !!----------------------------------------------------------------------
60      !!                ***  ROUTINE ssh_nxt  ***
61      !!                   
62      !! ** Purpose :   compute the after ssh (ssha)
63      !!
64      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
65      !!      is computed by integrating the horizontal divergence and multiply by
66      !!      by the time step.
67      !!
68      !! ** action  :   ssha    : after sea surface height
69      !!
70      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
71      !!----------------------------------------------------------------------
72      !
73      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv
74      INTEGER, INTENT(in) ::   kt                      ! time step
75      !
76      INTEGER             ::   jk                      ! dummy loop indices
77      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
78      !!----------------------------------------------------------------------
79      !
80      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
81      !
82      CALL wrk_alloc( jpi, jpj, zhdiv ) 
83      !
84      IF( kt == nit000 ) THEN
85         !
86         IF(lwp) WRITE(numout,*)
87         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
88         IF(lwp) WRITE(numout,*) '~~~~~~~ '
89         IF(lwp .AND. lflush) CALL flush(numout)
90         !
91      ENDIF
92      !
93      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
94      !
95      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
96      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
97
98
99#if defined key_asminc
100      !                                                ! Include the IAU weighted SSH increment
101      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
102         CALL ssh_asm_inc( kt )
103#if defined key_vvl
104! Don't directly adjust ssh but change hdivn at all levels instead
105! In trasbc also add in the heat and salt content associated with these changes at each level 
106        DO jk = 1, jpkm1                                 
107                 hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) ) * ( e3t_0(:,:,jk) / fse3t_n(:,:,jk) ) * tmask(:,:,jk) 
108        END DO
109      ENDIF
110#endif
111#endif
112
113
114      !                                           !------------------------------!
115      !                                           !   After Sea Surface Height   !
116      !                                           !------------------------------!
117      zhdiv(:,:) = 0._wp
118      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
119        zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
120      END DO
121      !                                                ! Sea surface elevation time stepping
122      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
123      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
124      !
125      z1_rau0 = 0.5_wp * r1_rau0
126      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
127
128#if ! defined key_dynspg_ts
129      ! These lines are not necessary with time splitting since
130      ! boundary condition on sea level is set during ts loop
131#if defined key_agrif
132      CALL agrif_ssh( kt )
133#endif
134#if defined key_bdy
135      IF (lk_bdy) THEN
136         CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
137         CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
138      ENDIF
139#endif
140#endif
141
142
143      !                                           !------------------------------!
144      !                                           !           outputs            !
145      !                                           !------------------------------!
146      !
147      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
148      !
149      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
150      !
151      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
152      !
153   END SUBROUTINE ssh_nxt
154
155   
156   SUBROUTINE wzv( kt )
157      !!----------------------------------------------------------------------
158      !!                ***  ROUTINE wzv  ***
159      !!                   
160      !! ** Purpose :   compute the now vertical velocity
161      !!
162      !! ** Method  : - Using the incompressibility hypothesis, the vertical
163      !!      velocity is computed by integrating the horizontal divergence 
164      !!      from the bottom to the surface minus the scale factor evolution.
165      !!        The boundary conditions are w=0 at the bottom (no flux) and.
166      !!
167      !! ** action  :   wn      : now vertical velocity
168      !!
169      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
170      !!----------------------------------------------------------------------
171      !
172      INTEGER, INTENT(in) ::   kt           ! time step
173      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
174      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
175      !
176      INTEGER             ::   ji, jj, jk   ! dummy loop indices
177      REAL(wp)            ::   z1_2dt       ! local scalars
178      !!----------------------------------------------------------------------
179     
180      IF( nn_timing == 1 )  CALL timing_start('wzv')
181      !
182      IF( kt == nit000 ) THEN
183         !
184         IF(lwp) WRITE(numout,*)
185         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
186         IF(lwp) WRITE(numout,*) '~~~~~ '
187         IF(lwp .AND. lflush) CALL flush(numout)
188         !
189         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
190         !
191      ENDIF
192      !                                           !------------------------------!
193      !                                           !     Now Vertical Velocity    !
194      !                                           !------------------------------!
195      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
196      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
197      !
198      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
199         CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
200         !
201         DO jk = 1, jpkm1
202            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
203            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
204            DO jj = 2, jpjm1
205               DO ji = fs_2, fs_jpim1   ! vector opt.
206                  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) )
207               END DO
208            END DO
209         END DO
210         CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
211         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
212         !                             ! Same question holds for hdivn. Perhaps just for security
213         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
214            ! computation of w
215            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    &
216               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
217         END DO
218         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
219         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 
220      ELSE   ! z_star and linear free surface cases
221         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
222            ! computation of w
223            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   &
224               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
225         END DO
226      ENDIF
227
228#if defined key_bdy
229      IF (lk_bdy) THEN
230         DO jk = 1, jpkm1
231            wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
232         END DO
233      ENDIF
234#endif
235      !
236      IF( nn_timing == 1 )  CALL timing_stop('wzv')
237
238
239   END SUBROUTINE wzv
240
241   SUBROUTINE ssh_swp( kt )
242      !!----------------------------------------------------------------------
243      !!                    ***  ROUTINE ssh_nxt  ***
244      !!
245      !! ** Purpose :   achieve the sea surface  height time stepping by
246      !!              applying Asselin time filter and swapping the arrays
247      !!              ssha  already computed in ssh_nxt 
248      !!
249      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
250      !!              from the filter, see Leclair and Madec 2010) and swap :
251      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
252      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
253      !!                sshn = ssha
254      !!
255      !! ** action  : - sshb, sshn   : before & now sea surface height
256      !!                               ready for the next time step
257      !!
258      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
259      !!----------------------------------------------------------------------
260      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
261      !!----------------------------------------------------------------------
262      !
263      IF( nn_timing == 1 )  CALL timing_start('ssh_swp')
264      !
265      IF( kt == nit000 ) THEN
266         IF(lwp) WRITE(numout,*)
267         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
268         IF(lwp) WRITE(numout,*) '~~~~~~~ '
269         IF(lwp .AND. lflush) CALL flush(numout)
270      ENDIF
271
272# if defined key_dynspg_ts
273      IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN    !** Euler time-stepping: no filter
274# else
275      IF ( neuler == 0 .AND. kt == nit000 ) THEN   !** Euler time-stepping at first time-step : no filter
276#endif
277         sshb(:,:) = sshn(:,:)                           ! before <-- now
278         sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now)
279      ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
280         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
281         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:)    - emp(:,:)    &
282                                &                                 - rnf_b(:,:)    + rnf(:,:)    &
283                                &                                 + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:)
284         sshn(:,:) = ssha(:,:)                           ! now <-- after
285      ENDIF
286      !
287      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
288      !
289      IF( nn_timing == 1 )  CALL timing_stop('ssh_swp')
290      !
291   END SUBROUTINE ssh_swp
292
293   !!======================================================================
294END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.