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.
dynspg_flt.F90 in branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 2797

Last change on this file since 2797 was 2797, checked in by davestorkey, 13 years ago

Delete BDY module and first implementation of new OBC module.

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