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_exp.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynspg_exp.F90 @ 363

Last change on this file since 363 was 359, checked in by opalod, 18 years ago

nemo_v1_update_033 : RB + CT : Add new surface pressure gradient algorithms
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 8.8 KB
Line 
1MODULE dynspg_exp
2   !!======================================================================
3   !!                   ***  MODULE  dynspg_exp  ***
4   !! Ocean dynamics:  surface pressure gradient trend
5   !!======================================================================
6#if ( defined key_dynspg_exp && ! defined key_autotasking ) ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_dynspg_exp'       free sfce cst vol. without filter nor ts
9   !!   NOT 'key_autotasking'                      k-j-i loop (vector opt.)
10   !!----------------------------------------------------------------------
11   !!   dyn_spg_exp  : update the momentum trend with the surface
12   !!                      pressure gradient in the free surface constant 
13   !!                      volume case with vector optimization
14   !!----------------------------------------------------------------------
15   !! * Modules used
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE in_out_manager  ! I/O manager
19   USE phycst          ! physical constants
20   USE ocesbc          ! ocean surface boundary condition
21   USE obc_oce         ! Lateral open boundary condition
22   USE lib_mpp         ! distributed memory computing library
23   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
24   USE prtctl          ! Print control
25
26   IMPLICIT NONE
27   PRIVATE
28
29   !! * Accessibility
30   PUBLIC dyn_spg_exp  ! routine called by step.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !!  OPA 9.0 , LOCEAN-IPSL (2005)
37   !! $Header$
38   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
39   !!----------------------------------------------------------------------
40
41
42CONTAINS
43
44   SUBROUTINE dyn_spg_exp( kt )
45      !!----------------------------------------------------------------------
46      !!                  ***  routine dyn_spg_exp  ***
47      !!
48      !! ** Purpose :   Compute the now trend due to the surface pressure
49      !!      gradient in case of explicit free surface formulation and
50      !!      add it to the general trend of momentum equation. Compute
51      !!      the free surface.
52      !!
53      !! ** Method  :   Explicit free surface formulation. The surface pressure
54      !!      gradient is given by:
55      !!         spgu = 1/rau0 d/dx(ps) =  g/e1u di( sshn )
56      !!         spgv = 1/rau0 d/dy(ps) =  g/e2v dj( sshn )
57      !!      -1- Compute the now surface pressure gradient
58      !!      -2- Add it to the general trend
59      !!      -3- Compute the horizontal divergence of velocities
60      !!      - the now divergence is given by :
61      !!         zhdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
62      !!      - integrate the horizontal divergence from the bottom
63      !!         to the surface
64      !!      - apply lateral boundary conditions on zhdivn
65      !!      -4- Estimate the after sea surface elevation from the kinematic
66      !!         surface boundary condition:
67      !!         zssha = sshb - 2 rdt ( zhdiv + emp )
68      !!      - Time filter applied on now sea surface elevation to avoid
69      !!         the divergence of two consecutive time-steps and swap of free
70      !!         surface arrays to start the next time step:
71      !!         sshb = sshn + atfp * [ sshb + zssha - 2 sshn ]
72      !!         sshn = zssha
73      !!      - apply lateral boundary conditions on sshn
74      !!
75      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
76      !!
77      !! References :
78      !!
79      !! History :
80      !!   9.0  !  05-11  (V. Garnier, G. Madec, L. Bessieres) Original code
81      !!
82      !!---------------------------------------------------------------------
83      !! * Arguments
84      INTEGER, INTENT( in )  ::   kt         ! ocean time-step index
85
86      !! * Local declarations
87      INTEGER  ::   ji, jj, jk               ! dummy loop indices
88      REAL(wp) ::   z2dt, zraur, zssha       ! temporary scalars
89      REAL(wp), DIMENSION(jpi,jpj)    ::  &  ! temporary arrays
90         &         zhdiv
91      !!----------------------------------------------------------------------
92
93      IF( kt == nit000 ) THEN
94         IF(lwp) WRITE(numout,*)
95         IF(lwp) WRITE(numout,*) 'dyn_spg_exp : surface pressure gradient trend'
96         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (explicit free surface)'
97
98         ! set to zero free surface specific arrays
99         spgu(:,:) = 0.e0                     ! surface pressure gradient (i-direction)
100         spgv(:,:) = 0.e0                     ! surface pressure gradient (j-direction)
101      ENDIF
102
103      ! 0. Local constant initialization
104      ! --------------------------------
105      z2dt = 2. * rdt                                       ! time step: leap-frog
106      IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest
107      zraur = 1.e0 / rauw
108
109
110      ! 1. Surface pressure gradient (now)
111      ! ----------------------------
112      DO jj = 2, jpjm1
113         DO ji = fs_2, fs_jpim1   ! vector opt.
114            spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
115            spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
116         END DO
117      END DO
118
119      ! 2. Add the surface pressure trend to the general trend
120      ! ------------------------------------------------------
121      DO jk = 1, jpkm1
122         DO jj = 2, jpjm1
123            DO ji = fs_2, fs_jpim1   ! vector opt.
124               ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
125               va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
126            END DO
127         END DO
128      END DO
129     
130      ! 3. Vertical integration of horizontal divergence of velocities
131      ! --------------------------------
132      zhdiv(:,:) = 0.e0
133      DO jk = jpkm1, 1, -1
134         DO jj = 2, jpjm1
135            DO ji = fs_2, fs_jpim1   ! vector opt.
136               zhdiv(ji,jj) = zhdiv(ji,jj) + (  e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)      &
137                  &                           - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)      &
138                  &                           + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)      &
139                  &                           - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  )   &
140                  &                        / ( e1t(ji,jj) * e2t(ji,jj) )
141            END DO
142         END DO
143      END DO
144
145#if defined key_obc
146      ! open boundaries (div must be zero behind the open boundary)
147      !  mpp remark: The zeroing of zhdiv can probably be extended to 1->jpi/jpj for the correct row/column
148      IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0.e0      ! east
149      IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0.e0      ! west
150      IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
151      IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0.e0      ! south
152#endif
153
154
155      ! 4. Sea surface elevation time stepping
156      ! --------------------------------------
157      ! Euler (forward) time stepping, no time filter
158      IF( neuler == 0 .AND. kt == nit000 ) THEN
159         DO jj = 1, jpj
160            DO ji = 1, jpi
161               ! after free surface elevation
162               zssha = sshb(ji,jj) - rdt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1)
163               ! swap of arrays
164               sshb(ji,jj) = sshn(ji,jj)
165               sshn(ji,jj) = zssha
166            END DO
167         END DO
168      ELSE
169         ! Leap-frog time stepping and time filter
170         DO jj = 1, jpj
171            DO ji = 1, jpi
172               ! after free surface elevation
173               zssha = sshb(ji,jj) - z2dt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1)
174               ! time filter and array swap
175               sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj)
176               sshn(ji,jj) = zssha
177            END DO
178         END DO
179      ENDIF
180
181      !Boundary conditions on sshn
182      CALL lbc_lnk( sshn, 'T', 1. )
183 
184      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
185         CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask)
186      ENDIF
187     
188   END SUBROUTINE dyn_spg_exp
189
190#else
191   !!----------------------------------------------------------------------
192   !!   Default case :   Empty module   No standart explicit free surface
193   !!----------------------------------------------------------------------
194CONTAINS
195   SUBROUTINE dyn_spg_exp( kt )       ! Empty routine
196      WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt
197   END SUBROUTINE dyn_spg_exp
198#endif
199   
200   !!======================================================================
201END MODULE dynspg_exp
Note: See TracBrowser for help on using the repository browser.