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

Last change on this file since 5836 was 5836, checked in by cetlod, 5 years ago

merge the simplification branch onto the trunk, see ticket #1612

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