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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBC/obcvol.F90 @ 2281

Last change on this file since 2281 was 2281, checked in by smasson, 14 years ago

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