source: branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90 @ 5591

Last change on this file since 5591 was 5477, checked in by cguiavarch, 5 years ago

Clear svn keywords from UKMO/dev_r5107_hadgem3_cplseq

File size: 8.6 KB
Line 
1MODULE bdyvol
2   !!======================================================================
3   !!                       ***  MODULE  bdyvol  ***
4   !! Ocean dynamic :  Volume constraint when unstructured boundary
5   !!                  and filtered free surface are used
6   !!======================================================================
7   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code
8   !!             -   !  2006-01  (J. Chanut) Bug correction
9   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
10   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
11   !!----------------------------------------------------------------------
12#if   defined key_bdy   &&   defined key_dynspg_flt
13   !!----------------------------------------------------------------------
14   !!   'key_bdy'            AND      unstructured open boundary conditions
15   !!   'key_dynspg_flt'                              filtered free surface
16   !!----------------------------------------------------------------------
17   USE timing          ! Timing
18   USE oce             ! ocean dynamics and tracers
19   USE sbcisf          ! ice shelf
20   USE dom_oce         ! ocean space and time domain
21   USE phycst          ! physical constants
22   USE bdy_oce         ! ocean open boundary conditions
23   USE lib_mpp         ! for mppsum
24   USE in_out_manager  ! I/O manager
25   USE sbc_oce         ! ocean surface boundary conditions
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC bdy_vol        ! routine called by dynspg_flt.h90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
36   !! $Id$
37   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE bdy_vol( kt )
42      !!----------------------------------------------------------------------
43      !!                      ***  ROUTINE bdyvol  ***
44      !!
45      !! ** Purpose :   This routine is called in dynspg_flt to control
46      !!      the volume of the system. A correction velocity is calculated
47      !!      to correct the total transport through the unstructured OBC.
48      !!      The total depth used is constant (H0) to be consistent with the
49      !!      linear free surface coded in OPA 8.2
50      !!
51      !! ** Method  :   The correction velocity (zubtpecor here) is defined calculating
52      !!      the total transport through all open boundaries (trans_bdy) minus
53      !!      the cumulate E-P flux (z_cflxemp) divided by the total lateral
54      !!      surface (bdysurftot) of the unstructured boundary.
55      !!         zubtpecor = [trans_bdy - z_cflxemp ]*(1./bdysurftot)
56      !!      with z_cflxemp => sum of (Evaporation minus Precipitation)
57      !!                       over all the domain in m3/s at each time step.
58      !!      z_cflxemp < 0 when precipitation dominate
59      !!      z_cflxemp > 0 when evaporation dominate
60      !!
61      !!      There are 2 options (user's desiderata):
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 (z_cflxemp=0) to calculate the correction velocity. So
65      !!            it will only balance the flux through open boundaries.
66      !!            (set nn_volctl to 0 in tne namelist for this option)
67      !!         2/ The volume is constant even with E-P flux. In this case
68      !!            the correction velocity must balance both the flux
69      !!            through open boundaries and the ones through the free
70      !!            surface.
71      !!            (set nn_volctl to 1 in tne namelist for this option)
72      !!----------------------------------------------------------------------
73      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
74      !!
75      INTEGER  ::   ji, jj, jk, jb, jgrd
76      INTEGER  ::   ib_bdy, ii, ij
77      REAL(wp) ::   zubtpecor, z_cflxemp, ztranst
78      TYPE(OBC_INDEX), POINTER :: idx
79      !!-----------------------------------------------------------------------------
80
81      IF( nn_timing == 1 ) CALL timing_start('bdy_vol')
82
83      IF( ln_vol ) THEN
84
85      IF( kt == nit000 ) THEN
86         IF(lwp) WRITE(numout,*)
87         IF(lwp) WRITE(numout,*)'bdy_vol : Correction of velocities along unstructured OBC'
88         IF(lwp) WRITE(numout,*)'~~~~~~~'
89      END IF 
90
91      ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain
92      ! -----------------------------------------------------------------------
93      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:)+rdivisf*fwfisf(:,:) ) * bdytmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0
94      IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain
95
96      ! Transport through the unstructured open boundary
97      ! ------------------------------------------------
98      zubtpecor = 0.e0
99      DO ib_bdy = 1, nb_bdy
100         idx => idx_bdy(ib_bdy)
101
102         jgrd = 2                               ! cumulate u component contribution first
103         DO jb = 1, idx%nblenrim(jgrd)
104            DO jk = 1, jpkm1
105               ii = idx%nbi(jb,jgrd)
106               ij = idx%nbj(jb,jgrd)
107               zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
108            END DO
109         END DO
110         jgrd = 3                               ! then add v component contribution
111         DO jb = 1, idx%nblenrim(jgrd)
112            DO jk = 1, jpkm1
113               ii = idx%nbi(jb,jgrd)
114               ij = idx%nbj(jb,jgrd)
115               zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
116            END DO
117         END DO
118
119      END DO
120      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain
121
122      ! The normal velocity correction
123      ! ------------------------------
124      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp) / bdysurftot 
125      ELSE                   ;   zubtpecor =   zubtpecor             / bdysurftot
126      END IF
127
128      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation
129      ! -------------------------------------------------------------
130      ztranst = 0.e0
131      DO ib_bdy = 1, nb_bdy
132         idx => idx_bdy(ib_bdy)
133
134         jgrd = 2                               ! correct u component
135         DO jb = 1, idx%nblenrim(jgrd)
136            DO jk = 1, jpkm1
137               ii = idx%nbi(jb,jgrd)
138               ij = idx%nbj(jb,jgrd)
139               ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb,jgrd) * zubtpecor * umask(ii,ij,jk)
140               ztranst = ztranst + idx%flagu(jb,jgrd) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
141            END DO
142         END DO
143         jgrd = 3                              ! correct v component
144         DO jb = 1, idx%nblenrim(jgrd)
145            DO jk = 1, jpkm1
146               ii = idx%nbi(jb,jgrd)
147               ij = idx%nbj(jb,jgrd)
148               va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb,jgrd) * zubtpecor * vmask(ii,ij,jk)
149               ztranst = ztranst + idx%flagv(jb,jgrd) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk)
150            END DO
151         END DO
152
153      END DO
154      IF( lk_mpp )   CALL mpp_sum( ztranst )   ! sum over the global domain
155 
156      ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected
157      ! ------------------------------------------------------
158      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
159         IF(lwp) WRITE(numout,*)
160         IF(lwp) WRITE(numout,*)'bdy_vol : time step :', kt
161         IF(lwp) WRITE(numout,*)'~~~~~~~ '
162         IF(lwp) WRITE(numout,*)'          cumulate flux EMP             =', z_cflxemp  , ' (m3/s)'
163         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC  =', bdysurftot, '(m2)'
164         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor =', zubtpecor , '(m/s)'
165         IF(lwp) WRITE(numout,*)'          cumulated transport ztranst   =', ztranst   , '(m3/s)'
166      END IF 
167      !
168      IF( nn_timing == 1 ) CALL timing_stop('bdy_vol')
169      !
170      END IF ! ln_vol
171
172   END SUBROUTINE bdy_vol
173
174#else
175   !!----------------------------------------------------------------------
176   !!   Dummy module                   NO Unstruct Open Boundary Conditions
177   !!----------------------------------------------------------------------
178CONTAINS
179   SUBROUTINE bdy_vol( kt )        ! Empty routine
180      WRITE(*,*) 'bdy_vol: You should not have seen this print! error?', kt
181   END SUBROUTINE bdy_vol
182#endif
183
184   !!======================================================================
185END MODULE bdyvol
Note: See TracBrowser for help on using the repository browser.