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.
bdyvol.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

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

Last change on this file since 4244 was 3294, checked in by rblod, 12 years ago

Merge of 3.4beta into the trunk

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