source: trunk/NEMO/OPA_SRC/DYN/dynspg_flt_jki.F90 @ 508

Last change on this file since 508 was 508, checked in by opalod, 15 years ago

nemo_v1_update_071:RB: add iom for restart and reorganization of restart

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 19.6 KB
Line 
1MODULE dynspg_flt_jki
2   !!======================================================================
3   !!                  ***  MODULE  dynspg_flt_jki  ***
4   !! Ocean dynamics:  surface pressure gradient trend
5   !!======================================================================
6#if ( defined key_dynspg_flt && defined key_mpp_omp )   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_dynspg_flt'                              filtered free surface
9   !!   'key_mpp_omp'                              j-k-i loop (vector opt.)
10   !!----------------------------------------------------------------------
11   !!----------------------------------------------------------------------
12   !!   dyn_spg_flt_jki : Update the momentum trend with the surface pressure
13   !!                     gradient for the free surf. constant volume case
14   !!                     with auto-tasking (j-slab) (no vectior opt.)
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE zdf_oce         ! ocean vertical physics
19   USE phycst          ! physical constant
20   USE ocesbc          ! Ocean Surface Boundary condition
21   USE flxrnf          ! ocean runoffs
22   USE sol_oce         ! ocean elliptic solver
23   USE solpcg          ! preconditionned conjugate gradient solver
24   USE solsor          ! Successive Over-relaxation solver
25   USE solfet          ! FETI solver
26   USE solsor_e        ! Successive Over-relaxation solver with MPP optimization
27   USE solver
28   USE obc_oce         ! Lateral open boundary condition
29   USE obcdyn          ! ocean open boundary condition (obc_dyn routines)
30   USE obcvol          ! ocean open boundary condition (obc_vol routines)
31   USE lib_mpp         ! distributed memory computing library
32   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
33   USE cla_dynspg      ! cross land advection
34   USE prtctl          ! Print control
35   USE in_out_manager  ! I/O manager
36   USE solmat          ! matrix construction for elliptic solvers
37   USE agrif_opa_interp
38   USE restart         ! only for lrst_oce
39   USE iom
40
41   IMPLICIT NONE
42   PRIVATE
43
44   !! * Accessibility
45   PUBLIC dyn_spg_flt_jki        ! routine called by step.F90
46
47   !! * Substitutions
48#  include "domzgr_substitute.h90"
49   !!----------------------------------------------------------------------
50   !!   OPA 9.0 , LOCEAN-IPSL (2005)
51   !! $Header$
52   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54
55CONTAINS
56
57   SUBROUTINE dyn_spg_flt_jki( kt, kindic )
58      !!----------------------------------------------------------------------
59      !!                  ***  routine dyn_spg_flt_jki  ***
60      !!
61      !! ** Purpose :   Compute the now trend due to the surface pressure
62      !!      gradient for free surface formulation with a constant ocean
63      !!      volume case, add it to the general trend of momentum equation.
64      !!
65      !! ** Method  :   Free surface formulation. The surface pressure gradient
66      !!      is given by:
67      !!         spgu = 1/rau0 d/dx(ps) =  1/e1u di( etn + rnu btda )
68      !!         spgv = 1/rau0 d/dy(ps) =  1/e2v dj( etn + rnu btda )
69      !!      where etn is the free surface elevation and btda is the after
70      !!      of the free surface elevation
71      !!       -1- compute the after sea surface elevation from the cinematic
72      !!      surface boundary condition:
73      !!              zssha = sshb + 2 rdt ( wn - emp )
74      !!           Time filter applied on now sea surface elevation to avoid
75      !!      the divergence of two consecutive time-steps and swap of free
76      !!      surface arrays to start the next time step:
77      !!              sshb = sshn + atfp * [ sshb + zssha - 2 sshn ]
78      !!              sshn = zssha
79      !!       -2- 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]
86      !!            spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
87      !!                 - grav 2 rdt hv /e2v dj[sshn + emp]
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      !!       -3- compute and add the free surface pressure gradient inclu-
96      !!      ding the additional force used to stabilize the equation.
97      !!      several slabs used ('key-autotasking')
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 the solver doesn t converge
105      !!
106      INTEGER  ::   ji, jj, jk               ! dummy loop indices
107      REAL(wp) ::   &              ! temporary scalars
108         z2dt, z2dtg, zraur, znugdt, znurau,   &
109         zssha, zgcb, zbtd, ztdgu, ztdgv, zgwgt
110      !!----------------------------------------------------------------------
111
112      IF( kt == nit000 ) THEN
113         IF(lwp) WRITE(numout,*)
114         IF(lwp) WRITE(numout,*) 'dyn_spg_flt_jki : surface pressure gradient trend'
115         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~  (free surface constant volume, autotasking case)'
116
117         ! set to zero free surface specific arrays
118         spgu(:,:) = 0.e0                     ! surface pressure gradient (i-direction)
119         
120         spgv(:,:) = 0.e0                     ! surface pressure gradient (j-direction)
121         
122         CALL solver_init( nit000 )           ! Elliptic solver initialisation
123
124         ! read filtered free surface arrays in restart file
125         CALL flt_rst( nit000, 'READ' )       ! read or initialize the following fields:
126         !                                    ! gcx, gcxb, sshb, sshn
127      ENDIF
128
129      ! 0. Local constant initialization
130      ! --------------------------------
131      ! time step: leap-frog
132      z2dt = 2. * rdt
133      ! time step: Euler if restart from rest
134      IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
135      IF( neuler == 0 .AND. kt == nit000+1 ) CALL sol_mat(kt)
136      ! coefficients
137      z2dtg  = grav * z2dt
138      zraur  = 1. / rauw
139      znugdt =  rnu * grav * z2dt
140      znurau =  znugdt * zraur
141
142      !                                                ! ===============
143      DO jj = 2, jpjm1                                 !  Vertical slab
144         !                                             ! ===============
145         ! Surface pressure gradient (now)
146         DO ji = 2, jpim1
147            spgu(ji,jj) = - grav * ( sshn(ji+1,jj  ) - sshn(ji,jj) ) / e1u(ji,jj)
148            spgv(ji,jj) = - grav * ( sshn(ji  ,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
149         END DO 
150         !  Add the surface pressure trend to the general trend
151         DO jk = 1, jpkm1
152            DO ji = 2, jpim1
153               ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
154               va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
155            END DO
156         END DO
157         ! Evaluate the masked next velocity (effect of the additional force not included)
158         DO jk = 1, jpkm1
159            DO ji = 2, jpim1
160               ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk)
161               va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk)
162            END DO
163         END DO
164         !                                             ! ===============
165      END DO                                           !   End of slab
166      !                                                ! ===============
167
168#if defined key_obc
169         ! Update velocities on each open boundary with the radiation algorithm
170         CALL obc_dyn( kt )
171         ! Correction of the barotropic componant velocity to control the volume of the system
172         CALL obc_vol( kt )
173#endif
174#if defined key_agrif
175      ! Update velocities on each coarse/fine interfaces
176
177      CALL Agrif_dyn( kt )
178#endif
179#if defined key_orca_r2
180      IF( n_cla == 1 )   CALL dyn_spg_cla( kt )      ! Cross Land Advection (Update (ua,va))
181#endif
182
183      !                                                ! ===============
184      DO jj = 2, jpjm1                                 !  Vertical slab
185         !                                             ! ===============
186
187         ! 2. compute the next vertically averaged velocity
188         ! ------------------------------------------------
189         !     (effect of the additional force not included)
190         ! initialize to zero
191         DO ji = 2, jpim1
192            spgu(ji,jj) = 0.e0
193            spgv(ji,jj) = 0.e0
194         END DO
195         ! vertical sum
196         DO jk = 1, jpk
197            DO ji = 2, jpim1
198               spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk)
199               spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk)
200            END DO
201         END DO
202         ! transport: multiplied by the horizontal scale factor
203         DO ji = 2, jpim1
204            spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
205            spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
206         END DO
207         !                                             ! ===============
208      END DO                                           !   End of slab
209      !                                                ! ===============
210
211      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
212     
213      ! Boundary conditions on (spgu,spgv)
214      CALL lbc_lnk( spgu, 'U', -1. )
215      CALL lbc_lnk( spgv, 'V', -1. )
216
217      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
218
219      ! 3. Right hand side of the elliptic equation and first guess
220      ! -----------------------------------------------------------
221      DO jj = 2, jpjm1
222         DO ji = 2, jpim1
223            ! Divergence of the after vertically averaged velocity
224            zgcb =  spgu(ji,jj) - spgu(ji-1,jj)   &
225                  + spgv(ji,jj) - spgv(ji,jj-1)
226            gcb(ji,jj) = gcdprc(ji,jj) * zgcb
227            ! First guess of the after barotropic transport divergence
228            zbtd = gcx(ji,jj)
229            gcx (ji,jj) = 2. * zbtd   - gcxb(ji,jj)
230            gcxb(ji,jj) =      zbtd
231         END DO
232      END DO
233      ! applied the lateral boundary conditions
234      IF( nsolv == 4)   CALL lbc_lnk_e( gcb, c_solver_pt, 1. ) 
235
236#if defined key_agrif
237      IF( .NOT. AGRIF_ROOT() ) THEN
238         ! add contribution of gradient of after barotropic transport divergence
239        IF( (nbondi == -1) .OR. (nbondi == 2) ) gcb(3,:) = gcb(3,:) &
240           &            -znugdt * z2dt*laplacu(2,:)*gcdprc(3,:)*hu(2,:)*e2u(2,:)
241        IF( (nbondi ==  1) .OR. (nbondi == 2) )  gcb(nlci-2,:) = gcb(nlci-2,:) &
242           &           +znugdt * z2dt*laplacu(nlci-2,:)*gcdprc(nlci-2,:)*hu(nlci-2,:)*e2u(nlci-2,:)
243        IF( (nbondj == -1) .OR. (nbondj == 2) ) gcb(:,3) = gcb(:,3) &
244           &           -znugdt * z2dt*laplacv(:,2)*gcdprc(:,3)*hv(:,2)*e1v(:,2)
245        IF( (nbondj ==  1) .OR. (nbondj == 2) )  gcb(:,nlcj-2) = gcb(:,nlcj-2) &
246           &           +znugdt * z2dt*laplacv(:,nlcj-2)*gcdprc(:,nlcj-2)*hv(:,nlcj-2)*e1v(:,nlcj-2)
247      ENDIF
248#endif
249
250      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
251     
252      ! 4. Relative precision (computation on one processor)
253      ! ---------------------
254      rnorme =0.
255      DO jj = 1, jpj
256         DO ji = 1, jpi
257            zgwgt  = gcdmat(ji,jj) * gcb(ji,jj)
258            rnorme = rnorme + gcb(ji,jj) * zgwgt * bmask(ji,jj)
259         END DO
260      END DO
261      IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
262
263      epsr = eps * eps * rnorme
264      ncut = 0
265      ! if rnorme is 0, the solution is 0, the solver isn't called
266      IF( rnorme == 0.e0 ) THEN
267         gcx(:,:) = 0.e0
268         res   = 0.e0
269         niter = 0
270         ncut  = 999
271      ENDIF
272     
273      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
274     
275      ! 5. Evaluate the next transport divergence
276      ! -----------------------------------------
277      !    Iterarive solver for the elliptic equation (except IF sol.=0)
278      !    (output in gcx with boundary conditions applied)
279      kindic = 0
280      IF( ncut == 0 ) THEN
281         IF( nsolv == 1 ) THEN         ! diagonal preconditioned conjuguate gradient
282            CALL sol_pcg( kindic )
283         ELSEIF( nsolv == 2 ) THEN     ! successive-over-relaxation
284            CALL sol_sor( kindic )
285         ELSEIF( nsolv == 3 ) THEN     ! FETI solver
286            CALL sol_fet( kindic )
287         ELSEIF( nsolv == 4 ) THEN     ! successive-over-relaxation with extra outer halo
288            CALL sol_sor_e( kindic )
289         ELSE                          ! e r r o r in nsolv namelist parameter
290            WRITE(ctmp1,*) ' ~~~~~~~~~~~~~~~~                not = ', nsolv
291            CALL ctl_stop( ' dyn_spg_flt_jki : e r r o r, nsolv = 1, 2, 3 or 4', ctmp1 )
292         ENDIF
293      ENDIF
294
295      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
296
297      !                                                ! ===============
298      DO jj = 2, jpjm1                                 !  Vertical slab
299         !                                             ! ===============
300         
301         ! 6. Transport divergence gradient multiplied by z2dt
302         ! -----------------------------------------------====
303         DO ji = 2, jpim1
304            ! trend of Transport divergence gradient
305            ztdgu = znugdt * (gcx(ji+1,jj  ) - gcx(ji,jj) ) / e1u(ji,jj)
306            ztdgv = znugdt * (gcx(ji  ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj)
307            ! multiplied by z2dt
308#if defined key_obc
309            ! caution : grad D = 0 along open boundaries
310            spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
311            spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
312#else
313            spgu(ji,jj) = z2dt * ztdgu
314            spgv(ji,jj) = z2dt * ztdgv
315#endif
316         END DO
317         
318#if defined key_agrif     
319      IF( .NOT. Agrif_Root() ) THEN
320         ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
321         IF( (nbondi == -1) .OR. (nbondi == 2) ) spgu(2,:) = znugdt * z2dt * laplacu(2,:) * umask(2,:,1)
322         IF( (nbondi ==  1) .OR. (nbondi == 2) ) spgu(nlci-2,:) = znugdt * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
323         IF( (nbondj == -1) .OR. (nbondj == 2) ) spgv(:,2) = znugdt * z2dt * laplacv(:,2) * vmask(:,2,1)
324         IF( (nbondj ==  1) .OR. (nbondj == 2) ) spgv(:,nlcj-2) = znugdt * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
325      ENDIF
326#endif
327     
328         ! 7.  Add the trends multiplied by z2dt to the after velocity
329         ! -----------------------------------------------------------
330         !     ( c a u t i o n : (ua,va) here are the after velocity not the
331         !                       trend, the leap-frog time stepping will not
332         !                       be done in dynnxt.F routine)
333         DO jk = 1, jpkm1
334            DO ji = 2, jpim1
335               ua(ji,jj,jk) = (ua(ji,jj,jk) + spgu(ji,jj)) * umask(ji,jj,jk)
336               va(ji,jj,jk) = (va(ji,jj,jk) + spgv(ji,jj)) * vmask(ji,jj,jk)
337            END DO
338         END DO
339
340         ! 8. Sea surface elevation time stepping
341         ! --------------------------------------
342         IF( neuler == 0 .AND. kt == nit000 ) THEN         ! Euler (forward) time stepping, no time filter
343            DO ji = 1, jpi
344               ! after free surface elevation
345               zssha = sshb(ji,jj) + rdt * ( wn(ji,jj,1) - emp(ji,jj) * zraur ) * tmask(ji,jj,1)
346               ! swap of arrays
347               sshb(ji,jj) = sshn(ji,jj)
348               sshn(ji,jj) = zssha
349            END DO
350         ELSE                                              ! Leap-frog time stepping and time filter
351            DO ji = 1, jpi
352               ! after free surface elevation
353               zssha = sshb(ji,jj) + z2dt * ( wn(ji,jj,1) - emp(ji,jj) * zraur ) * tmask(ji,jj,1)
354               ! time filter and array swap
355               sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj)
356               sshn(ji,jj) = zssha
357            END DO
358         ENDIF
359         !                                             ! ===============
360      END DO                                           !   End of slab
361      !                                                ! ===============
362
363      !Boundary conditions on sshn
364      CALL lbc_lnk( sshn, 'T', 1. )
365
366      ! write filtered free surface arrays in restart file
367      ! --------------------------------------------------
368      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
369
370      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
371         CALL prt_ctl( tab3d_1=ua  , clinfo1=' spg  - Ua : ', mask1=umask,   &
372            &          tab3d_2=va  , clinfo2='        Va : ', mask2=vmask )
373         CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg  - ssh: ', mask1=tmask)
374      ENDIF
375      !
376   END SUBROUTINE dyn_spg_flt_jki
377
378   SUBROUTINE flt_rst( kt, cdrw )
379     !!---------------------------------------------------------------------
380     !!                   ***  ROUTINE ts_rst  ***
381     !!
382     !! ** Purpose : Read or write filtered free surface arrays in restart file
383     !!----------------------------------------------------------------------
384     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
385     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
386     !!----------------------------------------------------------------------
387
388     IF( TRIM(cdrw) == 'READ' ) THEN
389        IF( iom_varid( numror, 'gcx' ) > 0 ) THEN
390! Caution : extra-hallow
391! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
392           CALL iom_get( numror, jpdom_local, 'gcx' , gcx (1:jpi,1:jpj) )
393           CALL iom_get( numror, jpdom_local, 'gcxb', gcxb(1:jpi,1:jpj) )
394           CALL iom_get( numror, jpdom_local, 'sshb', sshb(:,:)         )
395           CALL iom_get( numror, jpdom_local, 'sshn', sshn(:,:)         )
396           IF( neuler == 0 ) THEN
397              sshb(:,:) = sshn(:,:)
398              gcxb(:,:) = gcx (:,:)
399           ENDIF
400        ELSE
401           gcx (:,:) = 0.e0
402           gcxb(:,:) = 0.e0
403           sshb(:,:) = 0.e0
404           sshn(:,:) = 0.e0
405        ENDIF
406     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
407! Caution : extra-hallow
408! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
409        CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx( 1:jpi,1:jpj) )
410        CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
411        CALL iom_rstput( kt, nitrst, numrow, 'sshb', sshb(:,:)         )
412        CALL iom_rstput( kt, nitrst, numrow, 'sshn', sshn(:,:)         )
413     ENDIF
414     !
415   END SUBROUTINE flt_rst
416
417#else
418   !!----------------------------------------------------------------------
419   !!   Default case :                                         Empty module
420   !!----------------------------------------------------------------------
421CONTAINS
422   SUBROUTINE dyn_spg_flt_jki( kt, kindic )      ! Empty module
423      WRITE(*,*) 'dyn_spg_flt_jki: You should not have seen this print! error?', kt, kindic
424   END SUBROUTINE dyn_spg_flt_jki
425#endif
426   
427   !!======================================================================
428END MODULE dynspg_flt_jki
Note: See TracBrowser for help on using the repository browser.