Changeset 1125 for trunk/NEMO/OPA_SRC/BDY/bdyvol.F90
- Timestamp:
- 2008-06-23T11:05:02+02:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/BDY/bdyvol.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/BDY/bdyvol.F90
r911 r1125 1 1 MODULE bdyvol 2 !!====================================================================== ===========2 !!====================================================================== 3 3 !! *** MODULE bdyvol *** 4 4 !! Ocean dynamic : Volume constraint when unstructured boundary 5 5 !! and Free surface are used 6 !!================================================================================= 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 !!---------------------------------------------------------------------- 7 11 #if defined key_bdy && defined key_dynspg_flt 8 !!--------------------------------------------------------------------------------- 9 !! 'key_bdy' and unstructured open boundary conditions 10 !! 'key_dynspg_flt' constant volume free surface 11 !!--------------------------------------------------------------------------------- 12 !! * Modules used 12 !!---------------------------------------------------------------------- 13 !! 'key_bdy' and unstructured open boundary conditions 14 !! 'key_dynspg_flt' filtered free surface 15 !!---------------------------------------------------------------------- 13 16 USE oce ! ocean dynamics and tracers 14 17 USE dom_oce ! ocean space and time domain … … 22 25 PRIVATE 23 26 24 !! * Accessibility25 27 PUBLIC bdy_vol ! routine called by dynspg_flt.h90 26 28 27 29 !! * Substitutions 28 30 # include "domzgr_substitute.h90" 29 !!--------------------------------------------------------------------------------- 30 !! OPA 9.0 , LODYC-IPSL (2003) 31 !!--------------------------------------------------------------------------------- 31 !!---------------------------------------------------------------------- 32 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 33 !! $Id: $ 34 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 35 !!---------------------------------------------------------------------- 32 36 33 37 CONTAINS 34 38 35 SUBROUTINE bdy_vol 36 !!---------------------------------------------------------------------- --------39 SUBROUTINE bdy_vol( kt ) 40 !!---------------------------------------------------------------------- 37 41 !! *** ROUTINE bdyvol *** 38 42 !! 39 !! ** Purpose : 40 !! This routine is called in dynspg_flt to control 43 !! ** Purpose : This routine is called in dynspg_flt to control 41 44 !! the volume of the system. A correction velocity is calculated 42 45 !! to correct the total transport through the unstructured OBC. … … 44 47 !! linear free surface coded in OPA 8.2 45 48 !! 46 !! ** Method : 47 !! The correction velocity (zubtpecor here) is defined calculating 49 !! ** Method : The correction velocity (zubtpecor here) is defined calculating 48 50 !! the total transport through all open boundaries (trans_bdy) minus 49 !! the cumulate E-P flux (z Cflxemp) divided by the total lateral51 !! the cumulate E-P flux (z_cflxemp) divided by the total lateral 50 52 !! surface (bdysurftot) of the unstructured boundary. 51 !! 52 !! zubtpecor = [trans_bdy - zCflxemp ]*(1./bdysurftot) 53 !! 54 !! with zCflxemp => sum of (Evaporation minus Precipitation) 53 !! zubtpecor = [trans_bdy - z_cflxemp ]*(1./bdysurftot) 54 !! with z_cflxemp => sum of (Evaporation minus Precipitation) 55 55 !! over all the domain in m3/s at each time step. 56 !! 57 !! zCflxemp < 0 when precipitation dominate 58 !! zCflxemp > 0 when evaporation dominate 56 !! z_cflxemp < 0 when precipitation dominate 57 !! z_cflxemp > 0 when evaporation dominate 59 58 !! 60 59 !! There are 2 options (user's desiderata): 61 !!62 60 !! 1/ The volume changes according to E-P, this is the default 63 61 !! option. In this case the cumulate E-P flux are setting to 64 !! zero (z Cflxemp=0) to calculate the correction velocity. So62 !! zero (z_cflxemp=0) to calculate the correction velocity. So 65 63 !! it will only balance the flux through open boundaries. 66 64 !! (set volbdy to 0 in tne namelist for this option) 67 !!68 65 !! 2/ The volume is constant even with E-P flux. In this case 69 66 !! the correction velocity must balance both the flux … … 71 68 !! surface. 72 69 !! (set volbdy to 1 in tne namelist for this option) 70 !!---------------------------------------------------------------------- 71 INTEGER, INTENT( in ) :: kt ! ocean time-step index 73 72 !! 74 !! History : 75 !! 8.5 ! 05-01 (J. Chanut, A. Sellar) Original code 76 !! ! 06-01 (J. Chanut) Bug correction 77 !!---------------------------------------------------------------------------- 78 !! * Arguments 79 INTEGER, INTENT( in ) :: kt ! ocean time-step index 80 81 !! * Local declarations 82 INTEGER :: ji,jj,jb, jgrd, jk 83 REAL(wp) :: zubtpecor 84 REAL(wp) :: zCflxemp 85 REAL(wp) :: ztranst 73 INTEGER :: ji, jj, jk, jb, jgrd 74 INTEGER :: ii, ij 75 REAL(wp) :: zubtpecor, z_cflxemp, ztranst 86 76 !!----------------------------------------------------------------------------- 87 77 88 78 IF( kt == nit000 ) THEN 89 IF(lwp) WRITE(numout,*) ' '79 IF(lwp) WRITE(numout,*) 90 80 IF(lwp) WRITE(numout,*)'bdy_vol : Correction of velocities along unstructured OBC' 91 81 IF(lwp) WRITE(numout,*)'~~~~~~~' 92 IF(lwp) WRITE(numout,*)' '93 82 END IF 94 83 95 ! 1. Calculate the cumulate surface Flux zCflxemp (m3/s) over all the domain. 96 ! --------------------------------------------------------------------------- 97 98 zCflxemp = 0.e0 99 84 ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain 85 ! ----------------------------------------------------------------------- 86 z_cflxemp = 0.e0 100 87 DO jj = 1, jpj 101 88 DO ji = 1, jpi 102 zCflxemp = zCflxemp + ( (emp(ji,jj)*bdytmask(ji,jj)*tmask_i(ji,jj) )/rauw) & 103 *e1v(ji,jj)*e2u(ji,jj) 89 z_cflxemp = z_cflxemp + emp(ji,jj) * bdytmask(ji,jj) * tmask_i(ji,jj) / rauw * e1v(ji,jj) * e2u(ji,jj) 104 90 END DO 105 91 END DO 106 IF( lk_mpp ) CALL mpp_sum( z Cflxemp ) ! sum over the global domain92 IF( lk_mpp ) CALL mpp_sum( z_cflxemp ) ! sum over the global domain 107 93 108 ! 2. Barotropic velocity through the unstructured open boundary 109 ! ------------------------------------------------------------- 110 94 ! Barotropic velocity through the unstructured open boundary 95 ! ---------------------------------------------------------- 111 96 zubtpecor = 0.e0 112 113 jgrd = 2 ! cumulate u component contribution first 97 jgrd = 2 ! cumulate u component contribution first 114 98 DO jb = 1, nblenrim(jgrd) 115 DO jk = 1, jpkm1116 zubtpecor = zubtpecor + flagu(jb) * ua (nbi(jb,jgrd), nbj(jb,jgrd), jk) &117 * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) &118 * fse3u(nbi(jb,jgrd), nbj(jb,jgrd),jk)119 END DO99 DO jk = 1, jpkm1 100 ii = nbi(jb,jgrd) 101 ij = nbj(jb,jgrd) 102 zubtpecor = zubtpecor + flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 103 END DO 120 104 END DO 121 122 jgrd = 3 ! then add v component contribution 123 DO jb = 1, nblenrim(jgrd) 124 DO jk = 1, jpkm1 125 zubtpecor = zubtpecor + flagv(jb) * va (nbi(jb,jgrd), nbj(jb,jgrd), jk) & 126 * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) & 127 * fse3v(nbi(jb,jgrd), nbj(jb,jgrd), jk) 128 END DO 105 jgrd = 3 ! then add v component contribution 106 DO jb = 1, nblenrim(jgrd) 107 DO jk = 1, jpkm1 108 ii = nbi(jb,jgrd) 109 ij = nbj(jb,jgrd) 110 zubtpecor = zubtpecor + flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 111 END DO 129 112 END DO 130 131 113 IF( lk_mpp ) CALL mpp_sum( zubtpecor ) ! sum over the global domain 132 114 133 134 ! 3. The normal velocity correction 135 ! --------------------------------- 136 137 IF (volbdy==1) THEN 138 zubtpecor = (zubtpecor - zCflxemp)*(1./bdysurftot) 139 ELSE 140 zubtpecor = zubtpecor*(1./bdysurftot) 115 ! The normal velocity correction 116 ! ------------------------------ 117 IF (volbdy==1) THEN ; zubtpecor = ( zubtpecor - z_cflxemp) / bdysurftot 118 ELSE ; zubtpecor = zubtpecor / bdysurftot 141 119 END IF 142 120 121 ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation 122 ! ------------------------------------------------------------- 123 ztranst = 0.e0 124 jgrd = 2 ! correct u component 125 DO jb = 1, nblenrim(jgrd) 126 DO jk = 1, jpkm1 127 ii = nbi(jb,jgrd) 128 ij = nbj(jb,jgrd) 129 ua(ii,ij,jk) = ua(ii,ij,jk) - flagu(jb) * zubtpecor * umask(ii,ij,jk) 130 ztranst = ztranst + flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 131 END DO 132 END DO 133 jgrd = 3 ! correct v component 134 DO jb = 1, nblenrim(jgrd) 135 DO jk = 1, jpkm1 136 ii = nbi(jb,jgrd) 137 ij = nbj(jb,jgrd) 138 va(ii,ij,jk) = va(ii,ij,jk) -flagv(jb) * zubtpecor * vmask(ii,ij,jk) 139 ztranst = ztranst + flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 140 END DO 141 END DO 142 IF( lk_mpp ) CALL mpp_sum( ztranst ) ! sum over the global domain 143 144 ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected 145 ! ------------------------------------------------------ 146 143 147 IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 144 IF(lwp) WRITE(numout,*) ' '148 IF(lwp) WRITE(numout,*) 145 149 IF(lwp) WRITE(numout,*)'bdy_vol : time step :', kt 146 150 IF(lwp) WRITE(numout,*)'~~~~~~~ ' 147 IF(lwp) WRITE(numout,*)' ' 148 IF(lwp) WRITE(numout,*)' cumulate flux EMP :', zCflxemp,' (m3/s)' 149 IF(lwp) WRITE(numout,*)' total lateral surface of OBC :',bdysurftot,'(m2)' 150 IF(lwp) WRITE(numout,*)' correction velocity zubtpecor :',zubtpecor,'(m/s)' 151 IF(lwp) WRITE(numout,*)' ' 151 IF(lwp) WRITE(numout,*)' cumulate flux EMP =', z_cflxemp , ' (m3/s)' 152 IF(lwp) WRITE(numout,*)' total lateral surface of OBC =', bdysurftot, '(m2)' 153 IF(lwp) WRITE(numout,*)' correction velocity zubtpecor =', zubtpecor , '(m/s)' 154 IF(lwp) WRITE(numout,*)' cumulated transport ztranst =', ztranst , '(m3/s)' 152 155 END IF 153 154 ! 4. Correction of the total velocity on the unstructured 155 ! boundary to respect the mass flux conservation 156 ! ------------------------------------------------------- 157 158 ztranst = 0.e0 159 160 jgrd = 2 ! correct u component 161 DO jb = 1, nblenrim(jgrd) 162 DO jk = 1, jpkm1 163 ua(nbi(jb, jgrd), nbj(jb, jgrd), jk) = ua(nbi(jb, jgrd), nbj(jb, jgrd), jk) & 164 -flagu(jb) * zubtpecor * umask(nbi(jb, jgrd), nbj(jb, jgrd), jk) 165 ztranst = ztranst + flagu(jb) * ua (nbi(jb,jgrd), nbj(jb,jgrd), jk) & 166 * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) & 167 * fse3u(nbi(jb,jgrd), nbj(jb,jgrd), jk) 168 END DO 169 END DO 170 171 jgrd = 3 ! correct v component 172 DO jb = 1, nblenrim(jgrd) 173 DO jk = 1, jpkm1 174 va(nbi(jb, jgrd), nbj(jb, jgrd), jk) = va(nbi(jb, jgrd), nbj(jb, jgrd), jk) & 175 -flagv(jb) * zubtpecor * vmask(nbi(jb, jgrd), nbj(jb, jgrd), jk) 176 ztranst = ztranst + flagv(jb) * va (nbi(jb,jgrd), nbj(jb,jgrd), jk) & 177 * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) & 178 * fse3v(nbi(jb,jgrd), nbj(jb,jgrd), jk) 179 END DO 180 END DO 181 182 IF( lk_mpp ) CALL mpp_sum( ztranst ) ! sum over the global domain 183 184 ! 5. Check the cumulate transport through unstructured OBC 185 ! once barotropic velocities corrected 186 ! -------------------------------------------------------- 187 188 IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN 189 IF(lwp) WRITE(numout,*)' ' 190 IF(lwp) WRITE(numout,*)' Cumulate transport ztranst =', ztranst,'(m3/s)' 191 IF(lwp) WRITE(numout,*)' ' 192 END IF 193 156 ! 194 157 END SUBROUTINE bdy_vol 195 158 196 159 #else 197 !!---------------------------------------------------------------------- -----------198 !! Default option : Empty module199 !!---------------------------------------------------------------------- -----------160 !!---------------------------------------------------------------------- 161 !! Dummy module NO Unstruct Open Boundary Conditions 162 !!---------------------------------------------------------------------- 200 163 CONTAINS 201 SUBROUTINE bdy_vol ! Empty routine 164 SUBROUTINE bdy_vol( kt ) ! Empty routine 165 WRITE(*,*) 'bdy_vol: You should not have seen this print! error?', kt 202 166 END SUBROUTINE bdy_vol 203 167 #endif 204 168 205 !!====================================================================== ===========169 !!====================================================================== 206 170 END MODULE bdyvol
Note: See TracChangeset
for help on using the changeset viewer.
