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

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 4616

Last change on this file since 4616 was 4616, checked in by gm, 10 years ago

#1260 : see the associated wiki page for explanation

  • Property svn:keywords set to Id
File size: 18.0 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) ) * r1_e1u(ji,jj)
169               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) * r1_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      ! vertical sum
195      DO jj = 2, jpjm1
196         DO ji = fs_2, fs_jpim1   ! vector opt.
197            spgu(ji,jj) = fse3u_a(ji,jj,1) * ua(ji,jj,1)
198            spgv(ji,jj) = fse3v_a(ji,jj,1) * va(ji,jj,1)
199         END DO
200      END DO
201      DO jk = 2, jpkm1
202         DO jj = 2, jpjm1
203            DO ji = 2, jpim1
204               spgu(ji,jj) = spgu(ji,jj) + fse3u_a(ji,jj,jk) * ua(ji,jj,jk)
205               spgv(ji,jj) = spgv(ji,jj) + fse3v_a(ji,jj,jk) * va(ji,jj,jk)
206            END DO
207         END DO
208      END DO
209
210      ! transport: multiplied by the horizontal scale factor
211      DO jj = 2, jpjm1
212         DO ji = fs_2, fs_jpim1   ! vector opt.
213            spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
214            spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
215         END DO
216      END DO
217      CALL lbc_lnk( spgu, 'U', -1. )       ! lateral boundary conditions
218      CALL lbc_lnk( spgv, 'V', -1. )
219
220      IF( lk_vvl ) CALL sol_mat( kt )      ! build the matrix at kt (vvl case only)
221
222      ! Right hand side of the elliptic equation and first guess
223      ! --------------------------------------------------------
224      DO jj = 2, jpjm1
225         DO ji = fs_2, fs_jpim1   ! vector opt.
226            ! Divergence of the after vertically averaged velocity
227            zgcb =  spgu(ji,jj) - spgu(ji-1,jj)   &
228                  + spgv(ji,jj) - spgv(ji,jj-1)
229            gcb(ji,jj) = gcdprc(ji,jj) * zgcb
230            ! First guess of the after barotropic transport divergence
231            zbtd = gcx(ji,jj)
232            gcx (ji,jj) = 2. * zbtd   - gcxb(ji,jj)
233            gcxb(ji,jj) =      zbtd
234         END DO
235      END DO
236      ! applied the lateral boundary conditions
237      IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 )   CALL lbc_lnk_e( gcb, c_solver_pt, 1., jpr2di, jpr2dj )   
238
239#if defined key_agrif
240      IF( .NOT. AGRIF_ROOT() ) THEN
241         ! add contribution of gradient of after barotropic transport divergence
242         IF( nbondi == -1 .OR. nbondi == 2 )   gcb(3     ,:) =   &
243            &    gcb(3     ,:) - z2dtg * z2dt * laplacu(2     ,:) * gcdprc(3     ,:) * hu(2     ,:) * e2u(2     ,:)
244         IF( nbondi ==  1 .OR. nbondi == 2 )   gcb(nlci-2,:) =   &
245            &    gcb(nlci-2,:) + z2dtg * z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:)
246         IF( nbondj == -1 .OR. nbondj == 2 )   gcb(:     ,3) =   &
247            &    gcb(:,3     ) - z2dtg * z2dt * laplacv(:,2     ) * gcdprc(:,3     ) * hv(:,2     ) * e1v(:,2     )
248         IF( nbondj ==  1 .OR. nbondj == 2 )   gcb(:,nlcj-2) =   &
249            &    gcb(:,nlcj-2) + z2dtg * z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2)
250      ENDIF
251#endif
252
253
254      ! Relative precision (computation on one processor)
255      ! ------------------
256      rnorme =0.e0
257      rnorme = GLOB_SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) )
258
259      epsr = eps * eps * rnorme
260      ncut = 0
261      ! if rnorme is 0, the solution is 0, the solver is not called
262      IF( rnorme == 0._wp ) THEN
263         gcx(:,:) = 0._wp
264         res   = 0._wp
265         niter = 0
266         ncut  = 999
267      ENDIF
268
269      ! Evaluate the next transport divergence
270      ! --------------------------------------
271      !    Iterarive solver for the elliptic equation (except IF sol.=0)
272      !    (output in gcx with boundary conditions applied)
273      kindic = 0
274      IF( ncut == 0 ) THEN
275         IF    ( nn_solv == 1 ) THEN   ;   CALL sol_pcg( kindic )      ! diagonal preconditioned conjuguate gradient
276         ELSEIF( nn_solv == 2 ) THEN   ;   CALL sol_sor( kindic )      ! successive-over-relaxation
277         ENDIF
278      ENDIF
279
280      ! Transport divergence gradient multiplied by z2dt
281      ! --------------------------------------------====
282      DO jj = 2, jpjm1
283         DO ji = fs_2, fs_jpim1   ! vector opt.
284            ! trend of Transport divergence gradient
285            ztdgu = z2dtg * (gcx(ji+1,jj  ) - gcx(ji,jj) ) * r1_e1u(ji,jj)
286            ztdgv = z2dtg * (gcx(ji  ,jj+1) - gcx(ji,jj) ) * r1_e2v(ji,jj)
287            ! multiplied by z2dt
288#if defined key_bdy
289            IF(lk_bdy) THEN
290            ! caution : grad D = 0 along open boundaries
291               spgu(ji,jj) = z2dt * ztdgu * bdyumask(ji,jj)
292               spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj)
293            ELSE
294               spgu(ji,jj) = z2dt * ztdgu
295               spgv(ji,jj) = z2dt * ztdgv
296            ENDIF
297#else
298            spgu(ji,jj) = z2dt * ztdgu
299            spgv(ji,jj) = z2dt * ztdgv
300#endif
301         END DO
302      END DO
303
304#if defined key_agrif     
305      IF( .NOT. Agrif_Root() ) THEN
306         ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
307         IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2     ,:) = z2dtg * z2dt * laplacu(2     ,:) * umask(2     ,:,1)
308         IF( nbondi ==  1 .OR. nbondi == 2 ) spgu(nlci-2,:) = z2dtg * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
309         IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2     ) = z2dtg * z2dt * laplacv(:,2     ) * vmask(:     ,2,1)
310         IF( nbondj ==  1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = z2dtg * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
311      ENDIF
312#endif     
313      ! Add the trends multiplied by z2dt to the after velocity
314      ! -------------------------------------------------------
315      !     ( c a u t i o n : (ua,va) here are the after velocity not the
316      !                       trend, the leap-frog time stepping will not
317      !                       be done in dynnxt.F90 routine)
318      DO jk = 1, jpkm1
319         DO jj = 2, jpjm1
320            DO ji = fs_2, fs_jpim1   ! vector opt.
321               ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj) ) * umask(ji,jj,jk)
322               va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj) ) * vmask(ji,jj,jk)
323            END DO
324         END DO
325      END DO
326
327      ! write filtered free surface arrays in restart file
328      ! --------------------------------------------------
329      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
330      !
331      !
332      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_flt')
333      !
334   END SUBROUTINE dyn_spg_flt
335
336
337   SUBROUTINE flt_rst( kt, cdrw )
338      !!---------------------------------------------------------------------
339      !!                   ***  ROUTINE ts_rst  ***
340      !!
341      !! ** Purpose : Read or write filtered free surface arrays in restart file
342      !!----------------------------------------------------------------------
343      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
344      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
345      !!----------------------------------------------------------------------
346      !
347      IF( TRIM(cdrw) == 'READ' ) THEN
348         IF( iom_varid( numror, 'gcx', ldstop = .FALSE. ) > 0 ) THEN
349! Caution : extra-hallow
350! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
351            CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) )
352            CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) )
353            IF( neuler == 0 )   gcxb(:,:) = gcx (:,:)
354         ELSE
355            gcx (:,:) = 0.e0
356            gcxb(:,:) = 0.e0
357         ENDIF
358      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
359! Caution : extra-hallow
360! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
361         CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )
362         CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
363      ENDIF
364      !
365   END SUBROUTINE flt_rst
366
367#else
368   !!----------------------------------------------------------------------
369   !!   Default case :   Empty module   No standart free surface cst volume
370   !!----------------------------------------------------------------------
371CONTAINS
372   SUBROUTINE dyn_spg_flt( kt, kindic )       ! Empty routine
373      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt, kindic
374   END SUBROUTINE dyn_spg_flt
375   SUBROUTINE flt_rst    ( kt, cdrw )         ! Empty routine
376      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
377      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
378      WRITE(*,*) 'flt_rst: You should not have seen this print! error?', kt, cdrw
379   END SUBROUTINE flt_rst
380#endif
381   
382   !!======================================================================
383END MODULE dynspg_flt
Note: See TracBrowser for help on using the repository browser.