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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 2305

Last change on this file since 2305 was 2305, checked in by davestorkey, 14 years ago

Updates to BDY - DOCTOR standards

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