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

source: trunk/NEMO/OPA_SRC/OBC/obcvol.F90 @ 78

Last change on this file since 78 was 78, checked in by opalod, 20 years ago

CT : UPDATE052 : change logical lpXXXobc to lp_obc_XXX for Open Boundaries case
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.9 KB
Line 
1MODULE obcvol
2   !!=================================================================================
3   !!                       ***  MODULE  obcvol  ***
4   !! Ocean dynamic :  Volume constraint when OBC and Free surface are used
5   !!=================================================================================
6#if   defined key_obc   &&   defined key_dynspg_fsc
7   !!---------------------------------------------------------------------------------
8   !!   'key_obc'               and                           open boundary conditions
9   !!   'key_dynspg_fsc'                                  constant volume free surface
10   !!---------------------------------------------------------------------------------
11   !! * Modules used
12   USE oce             ! ocean dynamics and tracers
13   USE dom_oce         ! ocean space and time domain
14   USE phycst          ! physical constants
15   USE obc_oce         ! ocean open boundary conditions
16   USE lib_mpp         ! for mppsum
17   USE in_out_manager  ! I/O manager
18   USE ocesbc          ! ocean surface boundary conditions
19
20   IMPLICIT NONE
21   PRIVATE
22
23   !! * Accessibility
24   PUBLIC obc_vol        ! routine called by dynspg_fsc.h90
25
26   !! * Substitutions
27#  include "domzgr_substitute.h90"
28#  include "obc_vectopt_loop_substitute.h90"
29   !!---------------------------------------------------------------------------------
30   !!   OPA 9.0 , LODYC-IPSL  (2003)
31   !!---------------------------------------------------------------------------------
32
33CONTAINS
34
35   SUBROUTINE obc_vol ( kt )
36      !!------------------------------------------------------------------------------
37      !!                      ***  ROUTINE obcvol  ***
38      !!
39      !! ** Purpose :
40      !!      This routine is called in dynspg_fsc to control
41      !!      the volume of the system. A correction velocity is calculated
42      !!      to correct the total transport through the OBC.
43      !!      The total depth used is constant (H0) to be consistent with the
44      !!      linear free surface coded in OPA 8.2
45      !!
46      !! ** Method : 
47      !!      The correction velocity (zubtpecor here) is defined calculating
48      !!      the total transport through all open boundaries (trans_obc) minus
49      !!      the cumulate E-P flux (zCflxemp) divided by the total lateral
50      !!      surface (obcsurftot) of these OBC.
51      !!
52      !!      zubtpecor = [trans_obc - zCflxemp ]*(1./obcsurftot)
53      !!
54      !!      with zCflxemp => sum of (Evaporation minus Precipitation)
55      !!                       over all the domain in m3/s at each time step.
56      !!
57      !!      zCflxemp < 0 when precipitation dominate
58      !!      zCflxemp > 0 when evaporation dominate
59      !!
60      !!      There are 2 options (user's desiderata):
61      !!
62      !!         1/ The volume changes according to E-P, this is the default
63      !!            option. In this case the cumulate E-P flux are setting to
64      !!            zero (zCflxemp=0) to calculate the correction velocity. So
65      !!            it will only balance the flux through open boundaries.
66      !!            (set volemp to 0 in tne namelist for this option)
67      !!
68      !!         2/ The volume is constant even with E-P flux. In this case
69      !!            the correction velocity must balance both the flux
70      !!            through open boundaries and the ones through the free
71      !!            surface.
72      !!            (set volemp to 1 in tne namelist for this option)
73      !!
74      !! History :
75      !!   8.5  !  02-10 (C. Talandier, A-M. Treguier) Original code
76      !!----------------------------------------------------------------------------
77      !! * Arguments
78      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
79
80      !! * Local declarations
81      INTEGER ::   ji, jj, jk
82      REAL(wp) ::   zubtpecor
83      REAL(wp) ::   zCflxemp
84      REAL(wp) ::   ztransw, ztranse, ztransn, ztranss, ztranst
85      !!-----------------------------------------------------------------------------
86
87      IF( kt == nit000 ) THEN
88         IF(lwp) WRITE(numout,*)'        '
89         IF(lwp) WRITE(numout,*)'obc_vol : Correction of velocities along OBC'
90         IF(lwp) WRITE(numout,*)'~~~~~~~'
91         IF(lwp) WRITE(numout,*)'        '
92      END IF 
93
94      ! 1. Calculate the cumulate surface Flux zCflxemp (m3/s) over all the domain.
95      ! ---------------------------------------------------------------------------
96 
97      zCflxemp = 0.e0
98
99      DO jj = 1, jpj
100         DO ji = 1, jpi
101            zCflxemp = zCflxemp + ( (emp(ji,jj)*obctmsk(ji,jj) )/rauw)*e1v(ji,jj)*e2u(ji,jj)
102         END DO
103      END DO
104      IF( lk_mpp )   CALL mpp_sum( zCflxemp )   ! sum over the global domain
105
106      ! 2. Barotropic velocity for each open boundary
107      ! ---------------------------------------------
108
109      zubtpecor = 0.e0
110
111      ! ... West open boundary
112      IF( lp_obc_west ) THEN                      ! ... Total transport through the West OBC
113         DO ji = fs_niw0, fs_niw1 ! Vector opt.
114            DO jk = 1, jpkm1
115               DO jj = 1, jpj
116                  zubtpecor = zubtpecor + ua(ji,jj,jk)*e2u(ji,jj)*fse3u(ji,jj,jk) * uwmsk(jj,jk)
117               END DO
118            END DO
119         END DO
120      END IF 
121
122      ! ... East open boundary
123      IF( lp_obc_east ) THEN                      ! ... Total transport through the East OBC
124         DO ji = fs_nie0, fs_nie1 ! Vector opt.
125            DO jk = 1, jpkm1
126               DO jj = 1, jpj
127                  zubtpecor = zubtpecor - ua(ji,jj,jk)*e2u(ji,jj)*fse3u(ji,jj,jk) * uemsk(jj,jk)
128               END DO
129            END DO
130         END DO
131      END IF 
132
133      ! ... North open boundary
134      IF( lp_obc_north ) THEN                     ! ... Total transport through the North OBC
135         DO jj = fs_njn0, fs_njn1 ! Vector opt.
136            DO jk = 1, jpkm1
137               DO ji = 1, jpi
138                  zubtpecor = zubtpecor - va(ji,jj,jk)*e1v(ji,jj)*fse3v(ji,jj,jk) * vnmsk(ji,jk)
139               END DO
140            END DO
141         END DO
142      END IF 
143
144      ! ... South open boundary
145      IF( lp_obc_south ) THEN                     ! ... Total transport through the South OBC
146         DO jj = fs_njs0, fs_njs1 ! Vector opt.
147            DO jk = 1, jpkm1
148               DO ji = 1, jpi
149                  zubtpecor = zubtpecor + va(ji,jj,jk)*e1v(ji,jj)*fse3v(ji,jj,jk) * vsmsk(ji,jk)
150               END DO
151            END DO
152         END DO
153      END IF
154
155      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain
156
157
158      ! 3. The normal velocity correction
159      ! ---------------------------------
160
161      zubtpecor = (zubtpecor - zCflxemp*volemp)*(1./obcsurftot)
162
163      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
164         IF(lwp) WRITE(numout,*)'        '
165         IF(lwp) WRITE(numout,*)'obc_vol : time step :', kt
166         IF(lwp) WRITE(numout,*)'~~~~~~~ '
167         IF(lwp) WRITE(numout,*)'        '
168         IF(lwp) WRITE(numout,*)'          cumulate flux EMP :', zCflxemp,' (m3/s)'
169         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC :',obcsurftot,'(m2)'
170         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor :',zubtpecor,'(m/s)'
171         IF(lwp) WRITE(numout,*)'        '
172      END IF 
173
174      ! 4. Correction of the total velocity on each open
175      !    boundary torespect the mass flux conservation
176      ! -------------------------------------------------
177
178      ztransw = 0.e0
179      ztranse = 0.e0
180      ztransn = 0.e0
181      ztranss = 0.e0
182      ztranst = 0.e0
183
184      IF( lp_obc_west ) THEN
185
186         ! ... correction of the west velocity
187         DO ji = fs_niw0, fs_niw1 ! Vector opt.
188            DO jk = 1, jpkm1
189               DO jj = 1, jpj
190                  ua(ji,jj,jk) = ua(ji,jj,jk) - zubtpecor*uwmsk(jj,jk)
191                  ztransw= ztransw + ua(ji,jj,jk)*fse3u(ji,jj,jk)*e2u(ji,jj)*uwmsk(jj,jk)
192               END DO
193            END DO
194         END DO
195
196         IF( lk_mpp )   CALL mpp_sum( ztransw )   ! sum over the global domain
197
198         IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
199            IF(lwp) WRITE(numout,*)'          West OB transport ztransw :', ztransw,'(m3/s)'
200         END IF
201
202      END IF
203
204      IF( lp_obc_east ) THEN
205
206         ! ... correction of the east velocity
207         DO ji = fs_nie0, fs_nie1 ! Vector opt.
208            DO jk = 1, jpkm1
209               DO jj = 1, jpj
210                  ua(ji,jj,jk) = ua(ji,jj,jk) + zubtpecor*uemsk(jj,jk)
211                  ztranse= ztranse + ua(ji,jj,jk)*fse3u(ji,jj,jk)*e2u(ji,jj)*uemsk(jj,jk)
212               END DO
213            END DO
214         END DO
215
216         IF( lk_mpp )   CALL mpp_sum( ztranse )   ! sum over the global domain
217
218         IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
219            IF(lwp) WRITE(numout,*)'          East OB transport ztranse :', ztranse,'(m3/s)'
220         END IF
221
222      END IF
223
224      IF( lp_obc_north ) THEN
225
226         ! ... correction of the north velocity
227         DO jj = fs_njn0, fs_njn1 ! Vector opt.
228            DO jk = 1, jpkm1
229               DO ji =  1, jpi
230                  va(ji,jj,jk) = va(ji,jj,jk) + zubtpecor*vnmsk(ji,jk)
231                  ztransn= ztransn + va(ji,jj,jk)*fse3v(ji,jj,jk)*e1v(ji,jj)*vnmsk(ji,jk)
232               END DO
233            END DO
234         END DO
235         IF( lk_mpp )   CALL mpp_sum( ztransn )   ! sum over the global domain
236
237         IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
238            IF(lwp) WRITE(numout,*)'          North OB transport ztransn :', ztransn,'(m3/s)'
239         END IF
240
241      END IF
242
243      IF( lp_obc_south ) THEN
244
245         ! ... correction of the south velocity
246         DO jj = fs_njs0, fs_njs1 ! Vector opt.
247            DO jk = 1, jpkm1
248               DO ji =  1, jpi
249                  va(ji,jj,jk) = va(ji,jj,jk) - zubtpecor*vsmsk(ji,jk)
250                  ztranss= ztranss + va(ji,jj,jk)*fse3v(ji,jj,jk)*e1v(ji,jj)*vsmsk(ji,jk)
251               END DO
252            END DO
253         END DO
254         IF( lk_mpp )   CALL mpp_sum( ztranss )   ! sum over the global domain
255
256         IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
257            IF(lwp) WRITE(numout,*)'          South OB transport ztranss :', ztranss,'(m3/s)'
258         END IF
259
260      END IF 
261
262      ! 5. Check the cumulate transport through OBC
263      !    once barotropic velocities corrected
264      ! -------------------------------------------
265
266      ztranst = ztransw - ztranse + ztranss - ztransn
267
268      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
269         IF(lwp) WRITE(numout,*)'        '
270         IF(lwp) WRITE(numout,*)'          Cumulate transport ztranst =', ztranst,'(m3/s)'
271         IF(lwp) WRITE(numout,*)'        '
272      END IF
273
274   END SUBROUTINE obc_vol
275
276#else
277   !!---------------------------------------------------------------------------------
278   !!  Default option :                                                   Empty module
279   !!---------------------------------------------------------------------------------
280CONTAINS
281   SUBROUTINE obc_vol        ! Empty routine
282   END SUBROUTINE obc_vol
283#endif
284
285   !!=================================================================================
286END MODULE obcvol
Note: See TracBrowser for help on using the repository browser.