source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

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