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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk

  • Property svn:keywords set to Id
File size: 18.9 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 obcdyn          ! ocean open boundary condition (obc_dyn routines)
36   USE obcvol          ! ocean open boundary condition (obc_vol routines)
37   USE bdy_oce         ! Unstructured open boundaries condition
38   USE bdydyn          ! Unstructured open boundaries condition (bdy_dyn routine)
39   USE bdyvol          ! Unstructured open boundaries condition (bdy_vol routine)
40   USE cla             ! cross land advection
41   USE in_out_manager  ! I/O manager
42   USE lib_mpp         ! distributed memory computing library
43   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
44   USE prtctl          ! Print control
45   USE iom
46   USE restart         ! only for lrst_oce
47   USE lib_fortran
48#if defined key_agrif
49   USE agrif_opa_interp
50#endif
51
52   IMPLICIT NONE
53   PRIVATE
54
55   PUBLIC   dyn_spg_flt  ! routine called by step.F90
56   PUBLIC   flt_rst      ! routine called by istate.F90
57
58   !! * Substitutions
59#  include "domzgr_substitute.h90"
60#  include "vectopt_loop_substitute.h90"
61   !!----------------------------------------------------------------------
62   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
63   !! $Id$
64   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
65   !!----------------------------------------------------------------------
66CONTAINS
67
68   SUBROUTINE dyn_spg_flt( kt, kindic )
69      !!----------------------------------------------------------------------
70      !!                  ***  routine dyn_spg_flt  ***
71      !!
72      !! ** Purpose :   Compute the now trend due to the surface pressure
73      !!      gradient in case of filtered free surface formulation  and add
74      !!      it to the general trend of momentum equation.
75      !!
76      !! ** Method  :   Filtered free surface formulation. The surface
77      !!      pressure gradient is given by:
78      !!         spgu = 1/rau0 d/dx(ps) =  1/e1u di( sshn + btda )
79      !!         spgv = 1/rau0 d/dy(ps) =  1/e2v dj( sshn + btda )
80      !!      where sshn is the free surface elevation and btda is the after
81      !!      time derivative of the free surface elevation
82      !!       -1- evaluate the surface presure trend (including the addi-
83      !!      tional force) in three steps:
84      !!        a- compute the right hand side of the elliptic equation:
85      !!            gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ]
86      !!         where (spgu,spgv) are given by:
87      !!            spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ]
88      !!                 - grav 2 rdt hu /e1u di[sshn + (emp-rnf)]
89      !!            spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
90      !!                 - grav 2 rdt hv /e2v dj[sshn + (emp-rnf)]
91      !!         and define the first guess from previous computation :
92      !!            zbtd = btda
93      !!            btda = 2 zbtd - btdb
94      !!            btdb = zbtd
95      !!        b- compute the relative accuracy to be reached by the
96      !!         iterative solver
97      !!        c- apply the solver by a call to sol... routine
98      !!       -2- compute and add the free surface pressure gradient inclu-
99      !!      ding the additional force used to stabilize the equation.
100      !!
101      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
102      !!
103      !! References : Roullet and Madec 1999, JGR.
104      !!---------------------------------------------------------------------
105      USE oce, ONLY:   zub   => ta , zvb   => sa   ! (ta,sa) used as workspace
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      !!----------------------------------------------------------------------
113      !
114      IF( kt == nit000 ) THEN
115         IF(lwp) WRITE(numout,*)
116         IF(lwp) WRITE(numout,*) 'dyn_spg_flt : surface pressure gradient trend'
117         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (free surface constant volume case)'
118       
119         ! set to zero free surface specific arrays
120         spgu(:,:) = 0._wp                     ! surface pressure gradient (i-direction)
121         spgv(:,:) = 0._wp                     ! surface pressure gradient (j-direction)
122
123         ! read filtered free surface arrays in restart file
124         ! when using agrif, sshn, gcx have to be read in istate
125         IF(.NOT. lk_agrif)   CALL flt_rst( nit000, 'READ' )      ! read or initialize the following fields:
126         !                                                        ! gcx, gcxb
127      ENDIF
128
129      ! Local constant initialization
130      z2dt = 2. * rdt                                             ! time step: leap-frog
131      IF( neuler == 0 .AND. kt == nit000   )   z2dt = rdt         ! time step: Euler if restart from rest
132      IF( neuler == 0 .AND. kt == nit000+1 )   CALL sol_mat( kt )
133      z2dtg  = grav * z2dt
134
135      ! Evaluate the masked next velocity (effect of the additional force not included)
136      ! --------------------------------- 
137      IF( lk_vvl ) THEN          ! variable volume  (surface pressure gradient already included in dyn_hpg)
138         !
139         IF( ln_dynadv_vec ) THEN      ! vector form : applied on velocity
140            DO jk = 1, jpkm1
141               DO jj = 2, jpjm1
142                  DO ji = fs_2, fs_jpim1   ! vector opt.
143                     ua(ji,jj,jk) = (  ub(ji,jj,jk) + z2dt * ua(ji,jj,jk)  ) * umask(ji,jj,jk)
144                     va(ji,jj,jk) = (  vb(ji,jj,jk) + z2dt * va(ji,jj,jk)  ) * vmask(ji,jj,jk)
145                  END DO
146               END DO
147            END DO
148            !
149         ELSE                          ! flux form : applied on thickness weighted velocity
150            DO jk = 1, jpkm1
151               DO jj = 2, jpjm1
152                  DO ji = fs_2, fs_jpim1   ! vector opt.
153                     ua(ji,jj,jk) = (        ub(ji,jj,jk) * fse3u_b(ji,jj,jk)      &
154                        &           + z2dt * ua(ji,jj,jk) * fse3u_n(ji,jj,jk)  )   &
155                        &         / fse3u_a(ji,jj,jk) * umask(ji,jj,jk)
156                     va(ji,jj,jk) = (        vb(ji,jj,jk) * fse3v_b(ji,jj,jk)      &
157                        &           + z2dt * va(ji,jj,jk) * fse3v_n(ji,jj,jk)  )   &
158                        &         / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk)
159                 END DO
160               END DO
161            END DO
162            !
163         ENDIF
164         !
165      ELSE                       ! fixed volume  (add the surface pressure gradient + unweighted time stepping)
166         !
167         DO jj = 2, jpjm1              ! Surface pressure gradient (now)
168            DO ji = fs_2, fs_jpim1   ! vector opt.
169               spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
170               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
171            END DO
172         END DO
173         DO jk = 1, jpkm1              ! unweighted time stepping
174            DO jj = 2, jpjm1
175               DO ji = fs_2, fs_jpim1   ! vector opt.
176                  ua(ji,jj,jk) = (  ub(ji,jj,jk) + z2dt * ( ua(ji,jj,jk) + spgu(ji,jj) )  ) * umask(ji,jj,jk)
177                  va(ji,jj,jk) = (  vb(ji,jj,jk) + z2dt * ( va(ji,jj,jk) + spgv(ji,jj) )  ) * vmask(ji,jj,jk)
178               END DO
179            END DO
180         END DO
181         !
182      ENDIF
183
184#if defined key_obc
185      CALL obc_dyn( kt )      ! Update velocities on each open boundary with the radiation algorithm
186      CALL obc_vol( kt )      ! Correction of the barotropic componant velocity to control the volume of the system
187#endif
188#if defined key_bdy
189      CALL bdy_dyn_frs( kt )       ! Update velocities on unstructured boundary using the Flow Relaxation Scheme
190      CALL bdy_vol( kt )           ! Correction of the barotropic component velocity to control the volume of the system
191#endif
192#if defined key_agrif
193      CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces
194#endif
195      IF( nn_cla == 1 )   CALL cla_dynspg( kt )      ! Cross Land Advection (update (ua,va))
196
197      ! compute the next vertically averaged velocity (effect of the additional force not included)
198      ! ---------------------------------------------
199      DO jj = 2, jpjm1
200         DO ji = fs_2, fs_jpim1   ! vector opt.
201            spgu(ji,jj) = 0._wp
202            spgv(ji,jj) = 0._wp
203         END DO
204      END DO
205
206      ! vertical sum
207!CDIR NOLOOPCHG
208      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
209         DO jk = 1, jpkm1
210            DO ji = 1, jpij
211               spgu(ji,1) = spgu(ji,1) + fse3u(ji,1,jk) * ua(ji,1,jk)
212               spgv(ji,1) = spgv(ji,1) + fse3v(ji,1,jk) * va(ji,1,jk)
213            END DO
214         END DO
215      ELSE                        ! No  vector opt.
216         DO jk = 1, jpkm1
217            DO jj = 2, jpjm1
218               DO ji = 2, jpim1
219                  spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk)
220                  spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk)
221               END DO
222            END DO
223         END DO
224      ENDIF
225
226      ! transport: multiplied by the horizontal scale factor
227      DO jj = 2, jpjm1
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. )   
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_obc
305            ! caution : grad D = 0 along open boundaries
306            spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
307            spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
308#elif defined key_bdy
309            ! caution : grad D = 0 along open boundaries
310            ! Remark: The filtering force could be reduced here in the FRS zone
311            !         by multiplying spgu/spgv by (1-alpha) ?? 
312            spgu(ji,jj) = z2dt * ztdgu * bdyumask(ji,jj)
313            spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj)           
314#else
315            spgu(ji,jj) = z2dt * ztdgu
316            spgv(ji,jj) = z2dt * ztdgv
317#endif
318         END DO
319      END DO
320
321#if defined key_agrif     
322      IF( .NOT. Agrif_Root() ) THEN
323         ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
324         IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2     ,:) = z2dtg * z2dt * laplacu(2     ,:) * umask(2     ,:,1)
325         IF( nbondi ==  1 .OR. nbondi == 2 ) spgu(nlci-2,:) = z2dtg * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
326         IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2     ) = z2dtg * z2dt * laplacv(:,2     ) * vmask(:     ,2,1)
327         IF( nbondj ==  1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = z2dtg * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
328      ENDIF
329#endif     
330      ! Add the trends multiplied by z2dt to the after velocity
331      ! -------------------------------------------------------
332      !     ( c a u t i o n : (ua,va) here are the after velocity not the
333      !                       trend, the leap-frog time stepping will not
334      !                       be done in dynnxt.F90 routine)
335      DO jk = 1, jpkm1
336         DO jj = 2, jpjm1
337            DO ji = fs_2, fs_jpim1   ! vector opt.
338               ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj) ) * umask(ji,jj,jk)
339               va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj) ) * vmask(ji,jj,jk)
340            END DO
341         END DO
342      END DO
343
344      ! write filtered free surface arrays in restart file
345      ! --------------------------------------------------
346      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
347      !
348   END SUBROUTINE dyn_spg_flt
349
350
351   SUBROUTINE flt_rst( kt, cdrw )
352      !!---------------------------------------------------------------------
353      !!                   ***  ROUTINE ts_rst  ***
354      !!
355      !! ** Purpose : Read or write filtered free surface arrays in restart file
356      !!----------------------------------------------------------------------
357      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
358      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
359      !!----------------------------------------------------------------------
360      !
361      IF( TRIM(cdrw) == 'READ' ) THEN
362         IF( iom_varid( numror, 'gcx', ldstop = .FALSE. ) > 0 ) THEN
363! Caution : extra-hallow
364! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
365            CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) )
366            CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) )
367            IF( neuler == 0 )   gcxb(:,:) = gcx (:,:)
368         ELSE
369            gcx (:,:) = 0.e0
370            gcxb(:,:) = 0.e0
371         ENDIF
372      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
373! Caution : extra-hallow
374! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
375         CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )
376         CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
377      ENDIF
378      !
379   END SUBROUTINE flt_rst
380
381#else
382   !!----------------------------------------------------------------------
383   !!   Default case :   Empty module   No standart free surface cst volume
384   !!----------------------------------------------------------------------
385CONTAINS
386   SUBROUTINE dyn_spg_flt( kt, kindic )       ! Empty routine
387      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt, kindic
388   END SUBROUTINE dyn_spg_flt
389   SUBROUTINE flt_rst    ( kt, cdrw )         ! Empty routine
390      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
391      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
392      WRITE(*,*) 'flt_rst: You should not have seen this print! error?', kt, cdrw
393   END SUBROUTINE flt_rst
394#endif
395   
396   !!======================================================================
397END MODULE dynspg_flt
Note: See TracBrowser for help on using the repository browser.