source: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90 @ 5930

Last change on this file since 5930 was 5930, checked in by jchanut, 5 years ago

#1620 Merge free surface simplification into trunk

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