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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBC – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBC/obcvol.F90 @ 4400

Last change on this file since 4400 was 3432, checked in by trackstand2, 12 years ago

Merge branch 'ksection_partition'

  • Property svn:keywords set to Id
File size: 12.8 KB
RevLine 
[3]1MODULE obcvol
2   !!=================================================================================
3   !!                       ***  MODULE  obcvol  ***
[32]4   !! Ocean dynamic :  Volume constraint when OBC and Free surface are used
[3]5   !!=================================================================================
[1528]6#if   defined key_obc   &&   ! defined key_vvl
[3]7   !!---------------------------------------------------------------------------------
[1528]8   !!   'key_obc'               and   NOT                 open boundary conditions
9   !!   'key_vvl'                                         constant volume free surface
[3]10   !!---------------------------------------------------------------------------------
11   !! * Modules used
12   USE oce             ! ocean dynamics and tracers
13   USE dom_oce         ! ocean space and time domain
[1151]14   USE sbc_oce         ! ocean surface boundary conditions
[3]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
[367]24   PUBLIC obc_vol        ! routine called by dynspg_flt
[3]25
[3211]26   !! * Control permutation of array indices
27#  include "oce_ftrans.h90"
28#  include "dom_oce_ftrans.h90"
29#  include "sbc_oce_ftrans.h90"
30#  include "obc_oce_ftrans.h90"
31
[3]32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "obc_vectopt_loop_substitute.h90"
35   !!---------------------------------------------------------------------------------
[2528]36   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]37   !! $Id$
[2528]38   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]39   !!---------------------------------------------------------------------------------
40
41CONTAINS
42
43   SUBROUTINE obc_vol ( kt )
44      !!------------------------------------------------------------------------------
45      !!                      ***  ROUTINE obcvol  ***
46      !!
47      !! ** Purpose :
[367]48      !!      This routine is called in dynspg_flt to control
[3]49      !!      the volume of the system. A correction velocity is calculated
50      !!      to correct the total transport through the OBC.
51      !!      The total depth used is constant (H0) to be consistent with the
52      !!      linear free surface coded in OPA 8.2
53      !!
54      !! ** Method : 
55      !!      The correction velocity (zubtpecor here) is defined calculating
56      !!      the total transport through all open boundaries (trans_obc) minus
57      !!      the cumulate E-P flux (zCflxemp) divided by the total lateral
58      !!      surface (obcsurftot) of these OBC.
59      !!
60      !!      zubtpecor = [trans_obc - zCflxemp ]*(1./obcsurftot)
61      !!
62      !!      with zCflxemp => sum of (Evaporation minus Precipitation)
63      !!                       over all the domain in m3/s at each time step.
64      !!
65      !!      zCflxemp < 0 when precipitation dominate
66      !!      zCflxemp > 0 when evaporation dominate
67      !!
68      !!      There are 2 options (user's desiderata):
69      !!
70      !!         1/ The volume changes according to E-P, this is the default
71      !!            option. In this case the cumulate E-P flux are setting to
72      !!            zero (zCflxemp=0) to calculate the correction velocity. So
73      !!            it will only balance the flux through open boundaries.
74      !!            (set volemp to 0 in tne namelist for this option)
75      !!
76      !!         2/ The volume is constant even with E-P flux. In this case
77      !!            the correction velocity must balance both the flux
78      !!            through open boundaries and the ones through the free
79      !!            surface.
80      !!            (set volemp to 1 in tne namelist for this option)
81      !!
82      !! History :
[32]83      !!   8.5  !  02-10 (C. Talandier, A-M. Treguier) Original code
[3]84      !!----------------------------------------------------------------------------
85      !! * Arguments
86      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
87
88      !! * Local declarations
89      INTEGER ::   ji, jj, jk
90      REAL(wp) ::   zubtpecor
91      REAL(wp) ::   zCflxemp
92      REAL(wp) ::   ztransw, ztranse, ztransn, ztranss, ztranst
93      !!-----------------------------------------------------------------------------
94
95      IF( kt == nit000 ) THEN
96         IF(lwp) WRITE(numout,*)'        '
97         IF(lwp) WRITE(numout,*)'obc_vol : Correction of velocities along OBC'
98         IF(lwp) WRITE(numout,*)'~~~~~~~'
99         IF(lwp) WRITE(numout,*)'        '
100      END IF 
101
102      ! 1. Calculate the cumulate surface Flux zCflxemp (m3/s) over all the domain.
103      ! ---------------------------------------------------------------------------
104
[2528]105      zCflxemp = SUM ( ( emp(:,:)-rnf(:,:) )*obctmsk(:,:)* e1t(:,:) * e2t(:,:)  / rau0 ) 
[1151]106
[32]107      IF( lk_mpp )   CALL mpp_sum( zCflxemp )   ! sum over the global domain
[3]108
109      ! 2. Barotropic velocity for each open boundary
110      ! ---------------------------------------------
111
112      zubtpecor = 0.e0
113
[1151]114      ! ... East open boundary
115      IF( lp_obc_east ) THEN                      ! ... Total transport through the East OBC
[3211]116#if defined key_z_first
117         DO jj = 1, jpj
118            DO ji = fs_nie0, fs_nie1 ! Vector opt.
119               DO jk = 1, jpkm1
120#else
[1151]121         DO ji = fs_nie0, fs_nie1 ! Vector opt.
[3]122            DO jk = 1, jpkm1
123               DO jj = 1, jpj
[3211]124#endif
[1151]125                  zubtpecor = zubtpecor - ua(ji,jj,jk)*e2u(ji,jj)*fse3u(ji,jj,jk) * &
126             &     uemsk(jj,jk)*MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) )
[3]127               END DO
128            END DO
129         END DO
130      END IF 
131
[1151]132      ! ... West open boundary
133      IF( lp_obc_west ) THEN                      ! ... Total transport through the West OBC
[3211]134#if defined key_z_first
135         DO jj = 1, jpj
136            DO ji = fs_niw0, fs_niw1 ! Vector opt.
137               DO jk = 1, jpkm1
138#else
[1151]139         DO ji = fs_niw0, fs_niw1 ! Vector opt.
[3]140            DO jk = 1, jpkm1
141               DO jj = 1, jpj
[3211]142#endif
[1151]143                  zubtpecor = zubtpecor + ua(ji,jj,jk)*e2u(ji,jj)*fse3u(ji,jj,jk) * &
144             &    uwmsk(jj,jk) *MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) )
[3]145               END DO
146            END DO
147         END DO
[3211]148      ENDIF
[3]149
150      ! ... North open boundary
[78]151      IF( lp_obc_north ) THEN                     ! ... Total transport through the North OBC
[3211]152#if defined key_z_first
153         DO ji = 1, jpi
154            DO jj = fs_njn0, fs_njn1 ! Vector opt.
155               DO jk = 1, jpkm1
156#else
[3]157         DO jj = fs_njn0, fs_njn1 ! Vector opt.
158            DO jk = 1, jpkm1
159               DO ji = 1, jpi
[3211]160#endif
[1151]161                  zubtpecor = zubtpecor - va(ji,jj,jk)*e1v(ji,jj)*fse3v(ji,jj,jk) * &
162             &    vnmsk(ji,jk) * MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) )
[3]163               END DO
164            END DO
165         END DO
[3211]166      ENDIF
[3]167
168      ! ... South open boundary
[78]169      IF( lp_obc_south ) THEN                     ! ... Total transport through the South OBC
[3]170         DO jj = fs_njs0, fs_njs1 ! Vector opt.
[3211]171#if defined key_z_first
172            DO ji = 1, jpi
173               DO jk = 1, jpkm1
174#else
[3]175            DO jk = 1, jpkm1
176               DO ji = 1, jpi
[3211]177#endif
[1151]178                  zubtpecor = zubtpecor + va(ji,jj,jk)*e1v(ji,jj)*fse3v(ji,jj,jk) * &
179             &    vsmsk(ji,jk) * MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) )
[3]180               END DO
181            END DO
182         END DO
[3211]183      ENDIF
[3]184
[32]185      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain
[3]186
[32]187
[3]188      ! 3. The normal velocity correction
189      ! ---------------------------------
[32]190      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
[3]191         IF(lwp) WRITE(numout,*)'        '
192         IF(lwp) WRITE(numout,*)'obc_vol : time step :', kt
193         IF(lwp) WRITE(numout,*)'~~~~~~~ '
194         IF(lwp) WRITE(numout,*)'          cumulate flux EMP :', zCflxemp,' (m3/s)'
[1151]195         IF(lwp) WRITE(numout,*)'          lateral transport :',zubtpecor,'(m3/s)'
196         IF(lwp) WRITE(numout,*)'          net inflow        :',zubtpecor-zCflxemp,'(m3/s)'
197      ENDIF
198
199      zubtpecor = (zubtpecor - zCflxemp*volemp)*(1./obcsurftot)
200
201      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
[3]202         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC :',obcsurftot,'(m2)'
203         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor :',zubtpecor,'(m/s)'
204         IF(lwp) WRITE(numout,*)'        '
205      END IF 
206
207      ! 4. Correction of the total velocity on each open
[1151]208      !    boundary to respect the mass flux conservation
[3]209      ! -------------------------------------------------
210
[1151]211      ztranse = 0.e0   ; ztransw = 0.e0 ; ztransn = 0.e0 ; ztranss = 0.e0
212      ztranst = 0.e0  ! total
[3]213
[78]214      IF( lp_obc_west ) THEN
[3]215         ! ... correction of the west velocity
[3211]216#if defined key_z_first
217         DO jj = 1, jpj
218            DO ji = fs_niw0, fs_niw1 ! Vector opt.
219               DO jk = 1, jpkm1
220#else
[3]221         DO ji = fs_niw0, fs_niw1 ! Vector opt.
222            DO jk = 1, jpkm1
223               DO jj = 1, jpj
[3211]224#endif
[3]225                  ua(ji,jj,jk) = ua(ji,jj,jk) - zubtpecor*uwmsk(jj,jk)
[1151]226                  ztransw= ztransw + ua(ji,jj,jk)*fse3u(ji,jj,jk)*e2u(ji,jj)*uwmsk(jj,jk) * &
227             &    MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) )
[3]228               END DO
229            END DO
230         END DO
231
[32]232         IF( lk_mpp )   CALL mpp_sum( ztransw )   ! sum over the global domain
[3]233
[1151]234         IF( lwp .AND. MOD( kt, nwrite ) == 0)  WRITE(numout,*)'          West OB transport ztransw :', ztransw,'(m3/s)'
[3]235      END IF
236
[78]237      IF( lp_obc_east ) THEN
[3]238
239         ! ... correction of the east velocity
[3211]240#if defined key_z_first
241         DO jj = 1, jpj
242            DO ji = fs_nie0, fs_nie1 ! Vector opt.
243               DO jk = 1, jpkm1
244#else
[3]245         DO ji = fs_nie0, fs_nie1 ! Vector opt.
246            DO jk = 1, jpkm1
247               DO jj = 1, jpj
[3211]248#endif
[3]249                  ua(ji,jj,jk) = ua(ji,jj,jk) + zubtpecor*uemsk(jj,jk)
[1151]250                  ztranse= ztranse + ua(ji,jj,jk)*fse3u(ji,jj,jk)*e2u(ji,jj)*uemsk(jj,jk) * &
251            &     MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) )
[3]252               END DO
253            END DO
254         END DO
255
[32]256         IF( lk_mpp )   CALL mpp_sum( ztranse )   ! sum over the global domain
[3]257
[32]258         IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
[3]259            IF(lwp) WRITE(numout,*)'          East OB transport ztranse :', ztranse,'(m3/s)'
260         END IF
261
262      END IF
263
[78]264      IF( lp_obc_north ) THEN
[3]265
266         ! ... correction of the north velocity
267         DO jj = fs_njn0, fs_njn1 ! Vector opt.
[3211]268#if defined key_z_first
269            DO ji =  1, jpi
270               DO jk = 1, jpkm1
271#else
[3]272            DO jk = 1, jpkm1
273               DO ji =  1, jpi
[3211]274#endif
[3]275                  va(ji,jj,jk) = va(ji,jj,jk) + zubtpecor*vnmsk(ji,jk)
[1151]276                  ztransn= ztransn + va(ji,jj,jk)*fse3v(ji,jj,jk)*e1v(ji,jj)*vnmsk(ji,jk) * &
277           &      MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) )
[3]278               END DO
279            END DO
280         END DO
[32]281         IF( lk_mpp )   CALL mpp_sum( ztransn )   ! sum over the global domain
[3]282
[32]283         IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
[3]284            IF(lwp) WRITE(numout,*)'          North OB transport ztransn :', ztransn,'(m3/s)'
285         END IF
286
287      END IF
288
[78]289      IF( lp_obc_south ) THEN
[3]290
291         ! ... correction of the south velocity
292         DO jj = fs_njs0, fs_njs1 ! Vector opt.
[3211]293#if defined key_z_first
294            DO ji =  1, jpi
295               DO jk = 1, jpkm1
296#else
[3]297            DO jk = 1, jpkm1
298               DO ji =  1, jpi
[3211]299#endif
[3]300                  va(ji,jj,jk) = va(ji,jj,jk) - zubtpecor*vsmsk(ji,jk)
[1151]301                  ztranss= ztranss + va(ji,jj,jk)*fse3v(ji,jj,jk)*e1v(ji,jj)*vsmsk(ji,jk) * &
302            &     MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) )
[3]303               END DO
304            END DO
305         END DO
[32]306         IF( lk_mpp )   CALL mpp_sum( ztranss )   ! sum over the global domain
[3]307
[32]308         IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
[3]309            IF(lwp) WRITE(numout,*)'          South OB transport ztranss :', ztranss,'(m3/s)'
310         END IF
311
312      END IF 
313
314      ! 5. Check the cumulate transport through OBC
315      !    once barotropic velocities corrected
316      ! -------------------------------------------
317
318
[32]319      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
[1151]320         ztranst = ztransw - ztranse + ztranss - ztransn
[3]321         IF(lwp) WRITE(numout,*)'        '
322         IF(lwp) WRITE(numout,*)'          Cumulate transport ztranst =', ztranst,'(m3/s)'
[1151]323         IF(lwp) WRITE(numout,*)'          Balance  =', ztranst - zCflxemp ,'(m3/s)'
[3]324         IF(lwp) WRITE(numout,*)'        '
325      END IF
326
327   END SUBROUTINE obc_vol
328
329#else
330   !!---------------------------------------------------------------------------------
331   !!  Default option :                                                   Empty module
332   !!---------------------------------------------------------------------------------
333CONTAINS
[3432]334   SUBROUTINE obc_vol(kt)        ! Empty routine
335      !! * Arguments
336      INTEGER, INTENT( in ) ::   kt
[3]337   END SUBROUTINE obc_vol
338#endif
339
340   !!=================================================================================
341END MODULE obcvol
Note: See TracBrowser for help on using the repository browser.