source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 4328

Last change on this file since 4328 was 4328, checked in by davestorkey, 7 years ago

Remove OBC module at NEMO 3.6. See ticket #1189.

  • Property svn:keywords set to Id
File size: 18.3 KB
Line 
1MODULE dynspg_flt
2   !!======================================================================
3   !!                   ***  MODULE  dynspg_flt  ***
4   !! Ocean dynamics:  surface pressure gradient trend
5   !!======================================================================
6   !! History    OPA  !  1998-05  (G. Roullet)  free surface
7   !!                 !  1998-10  (G. Madec, M. Imbard)  release 8.2
8   !!   NEMO     O.1  !  2002-08  (G. Madec)  F90: Free form and module
9   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
10   !!            1.0  !  2004-08  (C. Talandier) New trends organization
11   !!             -   !  2005-11  (V. Garnier) Surface pressure gradient organization
12   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
13   !!             -   !  2006-08  (J.Chanut, A.Sellar) Calls to BDY routines.
14   !!            3.2  !  2009-03  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module
15   !!----------------------------------------------------------------------
16#if defined key_dynspg_flt   ||   defined key_esopa 
17   !!----------------------------------------------------------------------
18   !!   'key_dynspg_flt'                              filtered free surface
19   !!----------------------------------------------------------------------
20   !!   dyn_spg_flt  : update the momentum trend with the surface pressure gradient in the filtered free surface case
21   !!   flt_rst      : read/write the time-splitting restart fields in the ocean restart file
22   !!----------------------------------------------------------------------
23   USE oce             ! ocean dynamics and tracers
24   USE dom_oce         ! ocean space and time domain
25   USE zdf_oce         ! ocean vertical physics
26   USE sbc_oce         ! surface boundary condition: ocean
27   USE bdy_oce         ! Lateral open boundary condition
28   USE sol_oce         ! ocean elliptic solver
29   USE phycst          ! physical constants
30   USE domvvl          ! variable volume
31   USE dynadv          ! advection
32   USE solmat          ! matrix construction for elliptic solvers
33   USE solpcg          ! preconditionned conjugate gradient solver
34   USE solsor          ! Successive Over-relaxation solver
35   USE bdydyn          ! ocean open boundary condition on dynamics
36   USE bdyvol          ! ocean open boundary condition (bdy_vol routine)
37   USE cla             ! cross land advection
38   USE in_out_manager  ! I/O manager
39   USE lib_mpp         ! distributed memory computing library
40   USE wrk_nemo        ! Memory Allocation
41   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
42   USE prtctl          ! Print control
43   USE iom
44   USE lib_fortran
45#if defined key_agrif
46   USE agrif_opa_interp
47#endif
48   USE timing          ! Timing
49
50   IMPLICIT NONE
51   PRIVATE
52
53   PUBLIC   dyn_spg_flt  ! routine called by step.F90
54   PUBLIC   flt_rst      ! routine called by istate.F90
55
56   !! * Substitutions
57#  include "domzgr_substitute.h90"
58#  include "vectopt_loop_substitute.h90"
59   !!----------------------------------------------------------------------
60   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
61   !! $Id$
62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE dyn_spg_flt( kt, kindic )
67      !!----------------------------------------------------------------------
68      !!                  ***  routine dyn_spg_flt  ***
69      !!
70      !! ** Purpose :   Compute the now trend due to the surface pressure
71      !!      gradient in case of filtered free surface formulation  and add
72      !!      it to the general trend of momentum equation.
73      !!
74      !! ** Method  :   Filtered free surface formulation. The surface
75      !!      pressure gradient is given by:
76      !!         spgu = 1/rau0 d/dx(ps) =  1/e1u di( sshn + btda )
77      !!         spgv = 1/rau0 d/dy(ps) =  1/e2v dj( sshn + btda )
78      !!      where sshn is the free surface elevation and btda is the after
79      !!      time derivative of the free surface elevation
80      !!       -1- evaluate the surface presure trend (including the addi-
81      !!      tional force) in three steps:
82      !!        a- compute the right hand side of the elliptic equation:
83      !!            gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ]
84      !!         where (spgu,spgv) are given by:
85      !!            spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ]
86      !!                 - grav 2 rdt hu /e1u di[sshn + (emp-rnf)]
87      !!            spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
88      !!                 - grav 2 rdt hv /e2v dj[sshn + (emp-rnf)]
89      !!         and define the first guess from previous computation :
90      !!            zbtd = btda
91      !!            btda = 2 zbtd - btdb
92      !!            btdb = zbtd
93      !!        b- compute the relative accuracy to be reached by the
94      !!         iterative solver
95      !!        c- apply the solver by a call to sol... routine
96      !!       -2- compute and add the free surface pressure gradient inclu-
97      !!      ding the additional force used to stabilize the equation.
98      !!
99      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
100      !!
101      !! References : Roullet and Madec 1999, JGR.
102      !!---------------------------------------------------------------------
103      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
104      INTEGER, INTENT(  out) ::   kindic   ! solver convergence flag (<0 if not converge)
105      !!                                   
106      INTEGER  ::   ji, jj, jk   ! dummy loop indices
107      REAL(wp) ::   z2dt, z2dtg, zgcb, zbtd, ztdgu, ztdgv   ! local scalars
108      !!----------------------------------------------------------------------
109      !
110      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_flt')
111      !
112      !
113      IF( kt == nit000 ) THEN
114         IF(lwp) WRITE(numout,*)
115         IF(lwp) WRITE(numout,*) 'dyn_spg_flt : surface pressure gradient trend'
116         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (free surface constant volume case)'
117       
118         ! set to zero free surface specific arrays
119         spgu(:,:) = 0._wp                     ! surface pressure gradient (i-direction)
120         spgv(:,:) = 0._wp                     ! surface pressure gradient (j-direction)
121
122         ! read filtered free surface arrays in restart file
123         ! when using agrif, sshn, gcx have to be read in istate
124         IF(.NOT. lk_agrif)   CALL flt_rst( nit000, 'READ' )      ! read or initialize the following fields:
125         !                                                        ! gcx, gcxb
126      ENDIF
127
128      ! Local constant initialization
129      z2dt = 2. * rdt                                             ! time step: leap-frog
130      IF( neuler == 0 .AND. kt == nit000   )   z2dt = rdt         ! time step: Euler if restart from rest
131      IF( neuler == 0 .AND. kt == nit000+1 )   CALL sol_mat( kt )
132      z2dtg  = grav * z2dt
133
134      ! Evaluate the masked next velocity (effect of the additional force not included)
135      ! --------------------------------- 
136      IF( lk_vvl ) THEN          ! variable volume  (surface pressure gradient already included in dyn_hpg)
137         !
138         IF( ln_dynadv_vec ) THEN      ! vector form : applied on velocity
139            DO jk = 1, jpkm1
140               DO jj = 2, jpjm1
141                  DO ji = fs_2, fs_jpim1   ! vector opt.
142                     ua(ji,jj,jk) = (  ub(ji,jj,jk) + z2dt * ua(ji,jj,jk)  ) * umask(ji,jj,jk)
143                     va(ji,jj,jk) = (  vb(ji,jj,jk) + z2dt * va(ji,jj,jk)  ) * vmask(ji,jj,jk)
144                  END DO
145               END DO
146            END DO
147            !
148         ELSE                          ! flux form : applied on thickness weighted velocity
149            DO jk = 1, jpkm1
150               DO jj = 2, jpjm1
151                  DO ji = fs_2, fs_jpim1   ! vector opt.
152                     ua(ji,jj,jk) = (        ub(ji,jj,jk) * fse3u_b(ji,jj,jk)      &
153                        &           + z2dt * ua(ji,jj,jk) * fse3u_n(ji,jj,jk)  )   &
154                        &         / fse3u_a(ji,jj,jk) * umask(ji,jj,jk)
155                     va(ji,jj,jk) = (        vb(ji,jj,jk) * fse3v_b(ji,jj,jk)      &
156                        &           + z2dt * va(ji,jj,jk) * fse3v_n(ji,jj,jk)  )   &
157                        &         / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk)
158                 END DO
159               END DO
160            END DO
161            !
162         ENDIF
163         !
164      ELSE                       ! fixed volume  (add the surface pressure gradient + unweighted time stepping)
165         !
166         DO jj = 2, jpjm1              ! Surface pressure gradient (now)
167            DO ji = fs_2, fs_jpim1   ! vector opt.
168               spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
169               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
170            END DO
171         END DO
172         DO jk = 1, jpkm1              ! unweighted time stepping
173            DO jj = 2, jpjm1
174               DO ji = fs_2, fs_jpim1   ! vector opt.
175                  ua(ji,jj,jk) = (  ub(ji,jj,jk) + z2dt * ( ua(ji,jj,jk) + spgu(ji,jj) )  ) * umask(ji,jj,jk)
176                  va(ji,jj,jk) = (  vb(ji,jj,jk) + z2dt * ( va(ji,jj,jk) + spgv(ji,jj) )  ) * vmask(ji,jj,jk)
177               END DO
178            END DO
179         END DO
180         !
181      ENDIF
182
183#if defined key_bdy
184      IF( lk_bdy ) CALL bdy_dyn( kt )   ! Update velocities on each open boundary
185      IF( lk_bdy ) CALL bdy_vol( kt )   ! Correction of the barotropic component velocity to control the volume of the system
186#endif
187#if defined key_agrif
188      CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces
189#endif
190      IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_dynspg( kt )      ! Cross Land Advection (update (ua,va))
191
192      ! compute the next vertically averaged velocity (effect of the additional force not included)
193      ! ---------------------------------------------
194      DO jj = 2, jpjm1
195         DO ji = fs_2, fs_jpim1   ! vector opt.
196            spgu(ji,jj) = 0._wp
197            spgv(ji,jj) = 0._wp
198         END DO
199      END DO
200
201      ! vertical sum
202!CDIR NOLOOPCHG
203      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
204         DO jk = 1, jpkm1
205            DO ji = 1, jpij
206               spgu(ji,1) = spgu(ji,1) + fse3u_a(ji,1,jk) * ua(ji,1,jk)
207               spgv(ji,1) = spgv(ji,1) + fse3v_a(ji,1,jk) * va(ji,1,jk)
208            END DO
209         END DO
210      ELSE                        ! No  vector opt.
211         DO jk = 1, jpkm1
212            DO jj = 2, jpjm1
213               DO ji = 2, jpim1
214                  spgu(ji,jj) = spgu(ji,jj) + fse3u_a(ji,jj,jk) * ua(ji,jj,jk)
215                  spgv(ji,jj) = spgv(ji,jj) + fse3v_a(ji,jj,jk) * va(ji,jj,jk)
216               END DO
217            END DO
218         END DO
219      ENDIF
220
221      ! transport: multiplied by the horizontal scale factor
222      DO jj = 2, jpjm1
223         DO ji = fs_2, fs_jpim1   ! vector opt.
224            spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
225            spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
226         END DO
227      END DO
228      CALL lbc_lnk( spgu, 'U', -1. )       ! lateral boundary conditions
229      CALL lbc_lnk( spgv, 'V', -1. )
230
231      IF( lk_vvl ) CALL sol_mat( kt )      ! build the matrix at kt (vvl case only)
232
233      ! Right hand side of the elliptic equation and first guess
234      ! --------------------------------------------------------
235      DO jj = 2, jpjm1
236         DO ji = fs_2, fs_jpim1   ! vector opt.
237            ! Divergence of the after vertically averaged velocity
238            zgcb =  spgu(ji,jj) - spgu(ji-1,jj)   &
239                  + spgv(ji,jj) - spgv(ji,jj-1)
240            gcb(ji,jj) = gcdprc(ji,jj) * zgcb
241            ! First guess of the after barotropic transport divergence
242            zbtd = gcx(ji,jj)
243            gcx (ji,jj) = 2. * zbtd   - gcxb(ji,jj)
244            gcxb(ji,jj) =      zbtd
245         END DO
246      END DO
247      ! applied the lateral boundary conditions
248      IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 )   CALL lbc_lnk_e( gcb, c_solver_pt, 1., jpr2di, jpr2dj )   
249
250#if defined key_agrif
251      IF( .NOT. AGRIF_ROOT() ) THEN
252         ! add contribution of gradient of after barotropic transport divergence
253         IF( nbondi == -1 .OR. nbondi == 2 )   gcb(3     ,:) =   &
254            &    gcb(3     ,:) - z2dtg * z2dt * laplacu(2     ,:) * gcdprc(3     ,:) * hu(2     ,:) * e2u(2     ,:)
255         IF( nbondi ==  1 .OR. nbondi == 2 )   gcb(nlci-2,:) =   &
256            &    gcb(nlci-2,:) + z2dtg * z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:)
257         IF( nbondj == -1 .OR. nbondj == 2 )   gcb(:     ,3) =   &
258            &    gcb(:,3     ) - z2dtg * z2dt * laplacv(:,2     ) * gcdprc(:,3     ) * hv(:,2     ) * e1v(:,2     )
259         IF( nbondj ==  1 .OR. nbondj == 2 )   gcb(:,nlcj-2) =   &
260            &    gcb(:,nlcj-2) + z2dtg * z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2)
261      ENDIF
262#endif
263
264
265      ! Relative precision (computation on one processor)
266      ! ------------------
267      rnorme =0.e0
268      rnorme = GLOB_SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) )
269
270      epsr = eps * eps * rnorme
271      ncut = 0
272      ! if rnorme is 0, the solution is 0, the solver is not called
273      IF( rnorme == 0._wp ) THEN
274         gcx(:,:) = 0._wp
275         res   = 0._wp
276         niter = 0
277         ncut  = 999
278      ENDIF
279
280      ! Evaluate the next transport divergence
281      ! --------------------------------------
282      !    Iterarive solver for the elliptic equation (except IF sol.=0)
283      !    (output in gcx with boundary conditions applied)
284      kindic = 0
285      IF( ncut == 0 ) THEN
286         IF    ( nn_solv == 1 ) THEN   ;   CALL sol_pcg( kindic )      ! diagonal preconditioned conjuguate gradient
287         ELSEIF( nn_solv == 2 ) THEN   ;   CALL sol_sor( kindic )      ! successive-over-relaxation
288         ENDIF
289      ENDIF
290
291      ! Transport divergence gradient multiplied by z2dt
292      ! --------------------------------------------====
293      DO jj = 2, jpjm1
294         DO ji = fs_2, fs_jpim1   ! vector opt.
295            ! trend of Transport divergence gradient
296            ztdgu = z2dtg * (gcx(ji+1,jj  ) - gcx(ji,jj) ) / e1u(ji,jj)
297            ztdgv = z2dtg * (gcx(ji  ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj)
298            ! multiplied by z2dt
299#if defined key_bdy
300            IF(lk_bdy) THEN
301            ! caution : grad D = 0 along open boundaries
302               spgu(ji,jj) = z2dt * ztdgu * bdyumask(ji,jj)
303               spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj)
304            ELSE
305               spgu(ji,jj) = z2dt * ztdgu
306               spgv(ji,jj) = z2dt * ztdgv
307            ENDIF
308#else
309            spgu(ji,jj) = z2dt * ztdgu
310            spgv(ji,jj) = z2dt * ztdgv
311#endif
312         END DO
313      END DO
314
315#if defined key_agrif     
316      IF( .NOT. Agrif_Root() ) THEN
317         ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
318         IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2     ,:) = z2dtg * z2dt * laplacu(2     ,:) * umask(2     ,:,1)
319         IF( nbondi ==  1 .OR. nbondi == 2 ) spgu(nlci-2,:) = z2dtg * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
320         IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2     ) = z2dtg * z2dt * laplacv(:,2     ) * vmask(:     ,2,1)
321         IF( nbondj ==  1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = z2dtg * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
322      ENDIF
323#endif     
324      ! Add the trends multiplied by z2dt to the after velocity
325      ! -------------------------------------------------------
326      !     ( c a u t i o n : (ua,va) here are the after velocity not the
327      !                       trend, the leap-frog time stepping will not
328      !                       be done in dynnxt.F90 routine)
329      DO jk = 1, jpkm1
330         DO jj = 2, jpjm1
331            DO ji = fs_2, fs_jpim1   ! vector opt.
332               ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj) ) * umask(ji,jj,jk)
333               va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj) ) * vmask(ji,jj,jk)
334            END DO
335         END DO
336      END DO
337
338      ! write filtered free surface arrays in restart file
339      ! --------------------------------------------------
340      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
341      !
342      !
343      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_flt')
344      !
345   END SUBROUTINE dyn_spg_flt
346
347
348   SUBROUTINE flt_rst( kt, cdrw )
349      !!---------------------------------------------------------------------
350      !!                   ***  ROUTINE ts_rst  ***
351      !!
352      !! ** Purpose : Read or write filtered free surface arrays in restart file
353      !!----------------------------------------------------------------------
354      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
355      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
356      !!----------------------------------------------------------------------
357      !
358      IF( TRIM(cdrw) == 'READ' ) THEN
359         IF( iom_varid( numror, 'gcx', ldstop = .FALSE. ) > 0 ) THEN
360! Caution : extra-hallow
361! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
362            CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) )
363            CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) )
364            IF( neuler == 0 )   gcxb(:,:) = gcx (:,:)
365         ELSE
366            gcx (:,:) = 0.e0
367            gcxb(:,:) = 0.e0
368         ENDIF
369      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
370! Caution : extra-hallow
371! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
372         CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )
373         CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
374      ENDIF
375      !
376   END SUBROUTINE flt_rst
377
378#else
379   !!----------------------------------------------------------------------
380   !!   Default case :   Empty module   No standart free surface cst volume
381   !!----------------------------------------------------------------------
382CONTAINS
383   SUBROUTINE dyn_spg_flt( kt, kindic )       ! Empty routine
384      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt, kindic
385   END SUBROUTINE dyn_spg_flt
386   SUBROUTINE flt_rst    ( kt, cdrw )         ! Empty routine
387      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
388      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
389      WRITE(*,*) 'flt_rst: You should not have seen this print! error?', kt, cdrw
390   END SUBROUTINE flt_rst
391#endif
392   
393   !!======================================================================
394END MODULE dynspg_flt
Note: See TracBrowser for help on using the repository browser.