- Timestamp:
- 2011-07-11T12:53:56+02:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcvol.F90
r2528 r2797 1 1 MODULE obcvol 2 !!====================================================================== ===========2 !!====================================================================== 3 3 !! *** MODULE obcvol *** 4 !! Ocean dynamic : Volume constraint when OBC and Free surface are used 5 !!================================================================================= 6 #if defined key_obc && ! defined key_vvl 7 !!--------------------------------------------------------------------------------- 8 !! 'key_obc' and NOT open boundary conditions 9 !! 'key_vvl' constant volume free surface 10 !!--------------------------------------------------------------------------------- 11 !! * Modules used 12 USE oce ! ocean dynamics and tracers 13 USE dom_oce ! ocean space and time domain 14 USE sbc_oce ! ocean surface boundary conditions 15 USE phycst ! physical constants 16 USE obc_oce ! ocean open boundary conditions 17 USE lib_mpp ! for mppsum 18 USE in_out_manager ! I/O manager 19 20 IMPLICIT NONE 21 PRIVATE 22 23 !! * Accessibility 24 PUBLIC obc_vol ! routine called by dynspg_flt 25 26 !! * Substitutions 27 # include "domzgr_substitute.h90" 28 # include "obc_vectopt_loop_substitute.h90" 29 !!--------------------------------------------------------------------------------- 30 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 31 !! $Id$ 32 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 33 !!--------------------------------------------------------------------------------- 34 4 !! Ocean dynamic : Volume constraint when unstructured boundary 5 !! and 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 !!---------------------------------------------------------------------- 11 #if defined key_obc && defined key_dynspg_flt 12 !!$ !!---------------------------------------------------------------------- 13 !!$ !! 'key_obc' AND unstructured open boundary conditions 14 !!$ !! 'key_dynspg_flt' filtered free surface 15 !!$ !!---------------------------------------------------------------------- 16 !!$ USE oce ! ocean dynamics and tracers 17 !!$ USE dom_oce ! ocean space and time domain 18 !!$ USE phycst ! physical constants 19 !!$ USE obc_oce ! ocean open boundary conditions 20 !!$ USE lib_mpp ! for mppsum 21 !!$ USE in_out_manager ! I/O manager 22 !!$ USE sbc_oce ! ocean surface boundary conditions 23 !!$ 24 !!$ IMPLICIT NONE 25 !!$ PRIVATE 26 !!$ 27 !!$ PUBLIC obc_vol ! routine called by dynspg_flt.h90 28 !!$ 29 !!$ !! * Substitutions 30 !!$# include "domzgr_substitute.h90" 31 !!$ !!---------------------------------------------------------------------- 32 !!$ !! NEMO/OPA 3.3 , NEMO Consortium (2010) 33 !!$ !! $Id$ 34 !!$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 35 !!$ !!---------------------------------------------------------------------- 36 !!$CONTAINS 37 !!$ 38 !!$ SUBROUTINE obc_vol( kt ) 39 !!$ !!---------------------------------------------------------------------- 40 !!$ !! *** ROUTINE obcvol *** 41 !!$ !! 42 !!$ !! ** Purpose : This routine is called in dynspg_flt to control 43 !!$ !! the volume of the system. A correction velocity is calculated 44 !!$ !! to correct the total transport through the unstructured OBC. 45 !!$ !! The total depth used is constant (H0) to be consistent with the 46 !!$ !! linear free surface coded in OPA 8.2 47 !!$ !! 48 !!$ !! ** Method : The correction velocity (zubtpecor here) is defined calculating 49 !!$ !! the total transport through all open boundaries (trans_obc) minus 50 !!$ !! the cumulate E-P flux (z_cflxemp) divided by the total lateral 51 !!$ !! surface (obcsurftot) of the unstructured boundary. 52 !!$ !! zubtpecor = [trans_obc - z_cflxemp ]*(1./obcsurftot) 53 !!$ !! with z_cflxemp => sum of (Evaporation minus Precipitation) 54 !!$ !! over all the domain in m3/s at each time step. 55 !!$ !! z_cflxemp < 0 when precipitation dominate 56 !!$ !! z_cflxemp > 0 when evaporation dominate 57 !!$ !! 58 !!$ !! There are 2 options (user's desiderata): 59 !!$ !! 1/ The volume changes according to E-P, this is the default 60 !!$ !! option. In this case the cumulate E-P flux are setting to 61 !!$ !! zero (z_cflxemp=0) to calculate the correction velocity. So 62 !!$ !! it will only balance the flux through open boundaries. 63 !!$ !! (set nn_volctl to 0 in tne namelist for this option) 64 !!$ !! 2/ The volume is constant even with E-P flux. In this case 65 !!$ !! the correction velocity must balance both the flux 66 !!$ !! through open boundaries and the ones through the free 67 !!$ !! surface. 68 !!$ !! (set nn_volctl to 1 in tne namelist for this option) 69 !!$ !!---------------------------------------------------------------------- 70 !!$ INTEGER, INTENT( in ) :: kt ! ocean time-step index 71 !!$ !! 72 !!$ INTEGER :: ji, jj, jk, jb, jgrd 73 !!$ INTEGER :: ii, ij 74 !!$ REAL(wp) :: zubtpecor, z_cflxemp, ztranst 75 !!$ !!----------------------------------------------------------------------------- 76 !!$ 77 !!$ IF( ln_vol ) THEN 78 !!$ 79 !!$ IF( kt == nit000 ) THEN 80 !!$ IF(lwp) WRITE(numout,*) 81 !!$ IF(lwp) WRITE(numout,*)'obc_vol : Correction of velocities along unstructured OBC' 82 !!$ IF(lwp) WRITE(numout,*)'~~~~~~~' 83 !!$ END IF 84 !!$ 85 !!$ ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain 86 !!$ ! ----------------------------------------------------------------------- 87 !!$ z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:) ) * obctmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0 88 !!$ IF( lk_mpp ) CALL mpp_sum( z_cflxemp ) ! sum over the global domain 89 !!$ 90 !!$ ! Transport through the unstructured open boundary 91 !!$ ! ------------------------------------------------ 92 !!$ zubtpecor = 0.e0 93 !!$ jgrd = 2 ! cumulate u component contribution first 94 !!$ DO jb = 1, nblenrim(jgrd) 95 !!$ DO jk = 1, jpkm1 96 !!$ ii = nbi(jb,jgrd) 97 !!$ ij = nbj(jb,jgrd) 98 !!$ zubtpecor = zubtpecor + flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 99 !!$ END DO 100 !!$ END DO 101 !!$ jgrd = 3 ! then add v component contribution 102 !!$ DO jb = 1, nblenrim(jgrd) 103 !!$ DO jk = 1, jpkm1 104 !!$ ii = nbi(jb,jgrd) 105 !!$ ij = nbj(jb,jgrd) 106 !!$ zubtpecor = zubtpecor + flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 107 !!$ END DO 108 !!$ END DO 109 !!$ IF( lk_mpp ) CALL mpp_sum( zubtpecor ) ! sum over the global domain 110 !!$ 111 !!$ ! The normal velocity correction 112 !!$ ! ------------------------------ 113 !!$ IF( nn_volctl==1 ) THEN ; zubtpecor = ( zubtpecor - z_cflxemp) / obcsurftot 114 !!$ ELSE ; zubtpecor = zubtpecor / obcsurftot 115 !!$ END IF 116 !!$ 117 !!$ ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation 118 !!$ ! ------------------------------------------------------------- 119 !!$ ztranst = 0.e0 120 !!$ jgrd = 2 ! correct u component 121 !!$ DO jb = 1, nblenrim(jgrd) 122 !!$ DO jk = 1, jpkm1 123 !!$ ii = nbi(jb,jgrd) 124 !!$ ij = nbj(jb,jgrd) 125 !!$ ua(ii,ij,jk) = ua(ii,ij,jk) - flagu(jb) * zubtpecor * umask(ii,ij,jk) 126 !!$ ztranst = ztranst + flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 127 !!$ END DO 128 !!$ END DO 129 !!$ jgrd = 3 ! correct v component 130 !!$ DO jb = 1, nblenrim(jgrd) 131 !!$ DO jk = 1, jpkm1 132 !!$ ii = nbi(jb,jgrd) 133 !!$ ij = nbj(jb,jgrd) 134 !!$ va(ii,ij,jk) = va(ii,ij,jk) -flagv(jb) * zubtpecor * vmask(ii,ij,jk) 135 !!$ ztranst = ztranst + flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 136 !!$ END DO 137 !!$ END DO 138 !!$ IF( lk_mpp ) CALL mpp_sum( ztranst ) ! sum over the global domain 139 !!$ 140 !!$ ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected 141 !!$ ! ------------------------------------------------------ 142 !!$ IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 143 !!$ IF(lwp) WRITE(numout,*) 144 !!$ IF(lwp) WRITE(numout,*)'obc_vol : time step :', kt 145 !!$ IF(lwp) WRITE(numout,*)'~~~~~~~ ' 146 !!$ IF(lwp) WRITE(numout,*)' cumulate flux EMP =', z_cflxemp , ' (m3/s)' 147 !!$ IF(lwp) WRITE(numout,*)' total lateral surface of OBC =', obcsurftot, '(m2)' 148 !!$ IF(lwp) WRITE(numout,*)' correction velocity zubtpecor =', zubtpecor , '(m/s)' 149 !!$ IF(lwp) WRITE(numout,*)' cumulated transport ztranst =', ztranst , '(m3/s)' 150 !!$ END IF 151 !!$ ! 152 !!$ END IF ! ln_vol 153 !!$ 154 !!$ END SUBROUTINE obc_vol 155 !!$ 156 !!$#else 157 !!---------------------------------------------------------------------- 158 !! Dummy module NO Unstruct Open Boundary Conditions 159 !!---------------------------------------------------------------------- 35 160 CONTAINS 36 37 SUBROUTINE obc_vol ( kt ) 38 !!------------------------------------------------------------------------------ 39 !! *** ROUTINE obcvol *** 40 !! 41 !! ** Purpose : 42 !! This routine is called in dynspg_flt to control 43 !! the volume of the system. A correction velocity is calculated 44 !! to correct the total transport through the OBC. 45 !! The total depth used is constant (H0) to be consistent with the 46 !! linear free surface coded in OPA 8.2 47 !! 48 !! ** Method : 49 !! The correction velocity (zubtpecor here) is defined calculating 50 !! the total transport through all open boundaries (trans_obc) minus 51 !! the cumulate E-P flux (zCflxemp) divided by the total lateral 52 !! surface (obcsurftot) of these OBC. 53 !! 54 !! zubtpecor = [trans_obc - zCflxemp ]*(1./obcsurftot) 55 !! 56 !! with zCflxemp => sum of (Evaporation minus Precipitation) 57 !! over all the domain in m3/s at each time step. 58 !! 59 !! zCflxemp < 0 when precipitation dominate 60 !! zCflxemp > 0 when evaporation dominate 61 !! 62 !! There are 2 options (user's desiderata): 63 !! 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 (zCflxemp=0) to calculate the correction velocity. So 67 !! it will only balance the flux through open boundaries. 68 !! (set volemp to 0 in tne namelist for this option) 69 !! 70 !! 2/ The volume is constant even with E-P flux. In this case 71 !! the correction velocity must balance both the flux 72 !! through open boundaries and the ones through the free 73 !! surface. 74 !! (set volemp to 1 in tne namelist for this option) 75 !! 76 !! History : 77 !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Original code 78 !!---------------------------------------------------------------------------- 79 !! * Arguments 80 INTEGER, INTENT( in ) :: kt ! ocean time-step index 81 82 !! * Local declarations 83 INTEGER :: ji, jj, jk 84 REAL(wp) :: zubtpecor 85 REAL(wp) :: zCflxemp 86 REAL(wp) :: ztransw, ztranse, ztransn, ztranss, ztranst 87 !!----------------------------------------------------------------------------- 88 89 IF( kt == nit000 ) THEN 90 IF(lwp) WRITE(numout,*)' ' 91 IF(lwp) WRITE(numout,*)'obc_vol : Correction of velocities along OBC' 92 IF(lwp) WRITE(numout,*)'~~~~~~~' 93 IF(lwp) WRITE(numout,*)' ' 94 END IF 95 96 ! 1. Calculate the cumulate surface Flux zCflxemp (m3/s) over all the domain. 97 ! --------------------------------------------------------------------------- 98 99 zCflxemp = SUM ( ( emp(:,:)-rnf(:,:) )*obctmsk(:,:)* e1t(:,:) * e2t(:,:) / rau0 ) 100 101 IF( lk_mpp ) CALL mpp_sum( zCflxemp ) ! sum over the global domain 102 103 ! 2. Barotropic velocity for each open boundary 104 ! --------------------------------------------- 105 106 zubtpecor = 0.e0 107 108 ! ... East open boundary 109 IF( lp_obc_east ) THEN ! ... Total transport through the East OBC 110 DO ji = fs_nie0, fs_nie1 ! Vector opt. 111 DO jk = 1, jpkm1 112 DO jj = 1, jpj 113 zubtpecor = zubtpecor - ua(ji,jj,jk)*e2u(ji,jj)*fse3u(ji,jj,jk) * & 114 & uemsk(jj,jk)*MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 115 END DO 116 END DO 117 END DO 118 END IF 119 120 ! ... West open boundary 121 IF( lp_obc_west ) THEN ! ... Total transport through the West OBC 122 DO ji = fs_niw0, fs_niw1 ! Vector opt. 123 DO jk = 1, jpkm1 124 DO jj = 1, jpj 125 zubtpecor = zubtpecor + ua(ji,jj,jk)*e2u(ji,jj)*fse3u(ji,jj,jk) * & 126 & uwmsk(jj,jk) *MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 127 END DO 128 END DO 129 END DO 130 ENDIF 131 132 ! ... North open boundary 133 IF( lp_obc_north ) THEN ! ... Total transport through the North OBC 134 DO jj = fs_njn0, fs_njn1 ! Vector opt. 135 DO jk = 1, jpkm1 136 DO ji = 1, jpi 137 zubtpecor = zubtpecor - va(ji,jj,jk)*e1v(ji,jj)*fse3v(ji,jj,jk) * & 138 & vnmsk(ji,jk) * MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) ) 139 END DO 140 END DO 141 END DO 142 ENDIF 143 144 ! ... South open boundary 145 IF( lp_obc_south ) THEN ! ... Total transport through the South OBC 146 DO jj = fs_njs0, fs_njs1 ! Vector opt. 147 DO jk = 1, jpkm1 148 DO ji = 1, jpi 149 zubtpecor = zubtpecor + va(ji,jj,jk)*e1v(ji,jj)*fse3v(ji,jj,jk) * & 150 & vsmsk(ji,jk) * MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) ) 151 END DO 152 END DO 153 END DO 154 ENDIF 155 156 IF( lk_mpp ) CALL mpp_sum( zubtpecor ) ! sum over the global domain 157 158 159 ! 3. The normal velocity correction 160 ! --------------------------------- 161 IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 162 IF(lwp) WRITE(numout,*)' ' 163 IF(lwp) WRITE(numout,*)'obc_vol : time step :', kt 164 IF(lwp) WRITE(numout,*)'~~~~~~~ ' 165 IF(lwp) WRITE(numout,*)' cumulate flux EMP :', zCflxemp,' (m3/s)' 166 IF(lwp) WRITE(numout,*)' lateral transport :',zubtpecor,'(m3/s)' 167 IF(lwp) WRITE(numout,*)' net inflow :',zubtpecor-zCflxemp,'(m3/s)' 168 ENDIF 169 170 zubtpecor = (zubtpecor - zCflxemp*volemp)*(1./obcsurftot) 171 172 IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 173 IF(lwp) WRITE(numout,*)' total lateral surface of OBC :',obcsurftot,'(m2)' 174 IF(lwp) WRITE(numout,*)' correction velocity zubtpecor :',zubtpecor,'(m/s)' 175 IF(lwp) WRITE(numout,*)' ' 176 END IF 177 178 ! 4. Correction of the total velocity on each open 179 ! boundary to respect the mass flux conservation 180 ! ------------------------------------------------- 181 182 ztranse = 0.e0 ; ztransw = 0.e0 ; ztransn = 0.e0 ; ztranss = 0.e0 183 ztranst = 0.e0 ! total 184 185 IF( lp_obc_west ) THEN 186 ! ... correction of the west velocity 187 DO ji = fs_niw0, fs_niw1 ! Vector opt. 188 DO jk = 1, jpkm1 189 DO jj = 1, jpj 190 ua(ji,jj,jk) = ua(ji,jj,jk) - zubtpecor*uwmsk(jj,jk) 191 ztransw= ztransw + ua(ji,jj,jk)*fse3u(ji,jj,jk)*e2u(ji,jj)*uwmsk(jj,jk) * & 192 & MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 193 END DO 194 END DO 195 END DO 196 197 IF( lk_mpp ) CALL mpp_sum( ztransw ) ! sum over the global domain 198 199 IF( lwp .AND. MOD( kt, nwrite ) == 0) WRITE(numout,*)' West OB transport ztransw :', ztransw,'(m3/s)' 200 END IF 201 202 IF( lp_obc_east ) THEN 203 204 ! ... correction of the east velocity 205 DO ji = fs_nie0, fs_nie1 ! Vector opt. 206 DO jk = 1, jpkm1 207 DO jj = 1, jpj 208 ua(ji,jj,jk) = ua(ji,jj,jk) + zubtpecor*uemsk(jj,jk) 209 ztranse= ztranse + ua(ji,jj,jk)*fse3u(ji,jj,jk)*e2u(ji,jj)*uemsk(jj,jk) * & 210 & MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 211 END DO 212 END DO 213 END DO 214 215 IF( lk_mpp ) CALL mpp_sum( ztranse ) ! sum over the global domain 216 217 IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 218 IF(lwp) WRITE(numout,*)' East OB transport ztranse :', ztranse,'(m3/s)' 219 END IF 220 221 END IF 222 223 IF( lp_obc_north ) THEN 224 225 ! ... correction of the north velocity 226 DO jj = fs_njn0, fs_njn1 ! Vector opt. 227 DO jk = 1, jpkm1 228 DO ji = 1, jpi 229 va(ji,jj,jk) = va(ji,jj,jk) + zubtpecor*vnmsk(ji,jk) 230 ztransn= ztransn + va(ji,jj,jk)*fse3v(ji,jj,jk)*e1v(ji,jj)*vnmsk(ji,jk) * & 231 & MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) ) 232 END DO 233 END DO 234 END DO 235 IF( lk_mpp ) CALL mpp_sum( ztransn ) ! sum over the global domain 236 237 IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 238 IF(lwp) WRITE(numout,*)' North OB transport ztransn :', ztransn,'(m3/s)' 239 END IF 240 241 END IF 242 243 IF( lp_obc_south ) THEN 244 245 ! ... correction of the south velocity 246 DO jj = fs_njs0, fs_njs1 ! Vector opt. 247 DO jk = 1, jpkm1 248 DO ji = 1, jpi 249 va(ji,jj,jk) = va(ji,jj,jk) - zubtpecor*vsmsk(ji,jk) 250 ztranss= ztranss + va(ji,jj,jk)*fse3v(ji,jj,jk)*e1v(ji,jj)*vsmsk(ji,jk) * & 251 & MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) ) 252 END DO 253 END DO 254 END DO 255 IF( lk_mpp ) CALL mpp_sum( ztranss ) ! sum over the global domain 256 257 IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 258 IF(lwp) WRITE(numout,*)' South OB transport ztranss :', ztranss,'(m3/s)' 259 END IF 260 261 END IF 262 263 ! 5. Check the cumulate transport through OBC 264 ! once barotropic velocities corrected 265 ! ------------------------------------------- 266 267 268 IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 269 ztranst = ztransw - ztranse + ztranss - ztransn 270 IF(lwp) WRITE(numout,*)' ' 271 IF(lwp) WRITE(numout,*)' Cumulate transport ztranst =', ztranst,'(m3/s)' 272 IF(lwp) WRITE(numout,*)' Balance =', ztranst - zCflxemp ,'(m3/s)' 273 IF(lwp) WRITE(numout,*)' ' 274 END IF 275 276 END SUBROUTINE obc_vol 277 278 #else 279 !!--------------------------------------------------------------------------------- 280 !! Default option : Empty module 281 !!--------------------------------------------------------------------------------- 282 CONTAINS 283 SUBROUTINE obc_vol ! Empty routine 161 SUBROUTINE obc_vol( kt ) ! Empty routine 162 WRITE(*,*) 'obc_vol: You should not have seen this print! error?', kt 284 163 END SUBROUTINE obc_vol 285 164 #endif 286 165 287 !!====================================================================== ===========166 !!====================================================================== 288 167 END MODULE obcvol
Note: See TracChangeset
for help on using the changeset viewer.
