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