source: branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90 @ 5605

Last change on this file since 5605 was 5605, checked in by mathiot, 6 years ago

Marine glacier termini: initial commit

  • Property svn:keywords set to Id
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   USE traicw
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC bdy_vol        ! routine called by dynspg_flt.h90
32
33   !! * Substitutions
34#  include "domzgr_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
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 is called in dynspg_flt to control
47      !!      the volume of the system. 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      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:)+rdivisf*fwfisf(:,:) -qtot(:,:)) * bdytmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0
95      IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain
96
97      ! Transport through the unstructured open boundary
98      ! ------------------------------------------------
99      zubtpecor = 0.e0
100      DO ib_bdy = 1, nb_bdy
101         idx => idx_bdy(ib_bdy)
102
103         jgrd = 2                               ! cumulate u component contribution first
104         DO jb = 1, idx%nblenrim(jgrd)
105            DO jk = 1, jpkm1
106               ii = idx%nbi(jb,jgrd)
107               ij = idx%nbj(jb,jgrd)
108               zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
109            END DO
110         END DO
111         jgrd = 3                               ! then add v component contribution
112         DO jb = 1, idx%nblenrim(jgrd)
113            DO jk = 1, jpkm1
114               ii = idx%nbi(jb,jgrd)
115               ij = idx%nbj(jb,jgrd)
116               zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
117            END DO
118         END DO
119
120      END DO
121      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain
122
123      ! The normal velocity correction
124      ! ------------------------------
125      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp) / bdysurftot 
126      ELSE                   ;   zubtpecor =   zubtpecor             / bdysurftot
127      END IF
128
129      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation
130      ! -------------------------------------------------------------
131      ztranst = 0.e0
132      DO ib_bdy = 1, nb_bdy
133         idx => idx_bdy(ib_bdy)
134
135         jgrd = 2                               ! correct u component
136         DO jb = 1, idx%nblenrim(jgrd)
137            DO jk = 1, jpkm1
138               ii = idx%nbi(jb,jgrd)
139               ij = idx%nbj(jb,jgrd)
140               ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb,jgrd) * zubtpecor * umask(ii,ij,jk)
141               ztranst = ztranst + idx%flagu(jb,jgrd) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
142            END DO
143         END DO
144         jgrd = 3                              ! correct v component
145         DO jb = 1, idx%nblenrim(jgrd)
146            DO jk = 1, jpkm1
147               ii = idx%nbi(jb,jgrd)
148               ij = idx%nbj(jb,jgrd)
149               va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb,jgrd) * zubtpecor * vmask(ii,ij,jk)
150               ztranst = ztranst + idx%flagv(jb,jgrd) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk)
151            END DO
152         END DO
153
154      END DO
155      IF( lk_mpp )   CALL mpp_sum( ztranst )   ! sum over the global domain
156 
157      ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected
158      ! ------------------------------------------------------
159      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
160         IF(lwp) WRITE(numout,*)
161         IF(lwp) WRITE(numout,*)'bdy_vol : time step :', kt
162         IF(lwp) WRITE(numout,*)'~~~~~~~ '
163         IF(lwp) WRITE(numout,*)'          cumulate flux EMP             =', z_cflxemp  , ' (m3/s)'
164         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC  =', bdysurftot, '(m2)'
165         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor =', zubtpecor , '(m/s)'
166         IF(lwp) WRITE(numout,*)'          cumulated transport ztranst   =', ztranst   , '(m3/s)'
167      END IF 
168      !
169      IF( nn_timing == 1 ) CALL timing_stop('bdy_vol')
170      !
171      END IF ! ln_vol
172
173   END SUBROUTINE bdy_vol
174
175#else
176   !!----------------------------------------------------------------------
177   !!   Dummy module                   NO Unstruct Open Boundary Conditions
178   !!----------------------------------------------------------------------
179CONTAINS
180   SUBROUTINE bdy_vol( kt )        ! Empty routine
181      WRITE(*,*) 'bdy_vol: You should not have seen this print! error?', kt
182   END SUBROUTINE bdy_vol
183#endif
184
185   !!======================================================================
186END MODULE bdyvol
Note: See TracBrowser for help on using the repository browser.