Changeset 11074
- Timestamp:
- 2019-06-04T17:46:46+02:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn2d.F90
r11071 r11074 14 14 !! bdy_ssh : Duplicate sea level across open boundaries 15 15 !!---------------------------------------------------------------------- 16 USE oce ! ocean dynamics and tracers16 USE oce, only : sshdta => spgu ! ocean dynamics and tracers 17 17 USE dom_oce ! ocean space and time domain 18 18 USE bdy_oce ! ocean open boundary conditions … … 169 169 170 170 INTEGER :: jb, igrd ! dummy loop indices 171 INTEGER :: ii, ij, iim1, iip1, ijm1, ijp1 ! 2D addresses 171 INTEGER :: ii, ij ! 2D addresses 172 INTEGER :: iiTrim, ijTrim ! T pts i/j-indice on the rim 173 INTEGER :: iiToce, ijToce, iiUoce, ijVoce ! T, U and V pts i/j-indice of the ocean next to the rim 172 174 REAL(wp), POINTER :: flagu, flagv ! short cuts 173 REAL(wp) :: zcorr ! Flather correction 174 REAL(wp) :: zforc ! temporary scalar 175 REAL(wp) :: zflag, z1_2 ! " " 175 REAL(wp) :: zfla ! Flather correction 176 REAL(wp) :: z1_2 ! 176 177 !!---------------------------------------------------------------------- 177 178 … … 180 181 ! ---------------------------------! 181 182 ! Flather boundary conditions :! 182 ! ---------------------------------! 183 184 !!! REPLACE spgu with nemo_wrk work space 185 186 ! Fill temporary array with ssh data (here spgu): 183 ! ---------------------------------! 184 185 ! Fill temporary array with ssh data (here we use spgu with the alias sshdta): 187 186 igrd = 1 188 s pgu(:,:) = 0.0187 sshdta(:,:) = 0.0 189 188 DO jb = 1, idx%nblenrim(igrd) 190 189 ii = idx%nbi(jb,igrd) 191 190 ij = idx%nbj(jb,igrd) 192 191 IF( ll_wd ) THEN 193 s pgu(ii, ij) = dta%ssh(jb) - ssh_ref192 sshdta(ii, ij) = dta%ssh(jb) - ssh_ref 194 193 ELSE 195 s pgu(ii, ij) = dta%ssh(jb)194 sshdta(ii, ij) = dta%ssh(jb) 196 195 ENDIF 197 196 END DO … … 201 200 ! ! I think we should rather use after ssh ? 202 201 DO jb = 1, idx%nblenrim(igrd) 203 ii 204 ij 202 ii = idx%nbi(jb,igrd) 203 ij = idx%nbj(jb,igrd) 205 204 flagu => idx%flagu(jb,igrd) 206 iim1 = ii + MAX( 0, INT( flagu ) ) ! T pts i-indice inside the boundary 207 iip1 = ii - MIN( 0, INT( flagu ) ) ! T pts i-indice outside the boundary 208 IF( iim1 > jpi .OR. iip1 > jpi ) CYCLE 209 ! 210 zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 211 212 ! jchanut tschanges: Set zflag to 0 below to revert to Flather scheme 213 ! Use characteristics method instead 214 zflag = ABS(flagu) 215 zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij) 216 pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 205 IF( flagu == 0. ) THEN 206 pua2d(ii,ij) = dta%u2d(jb) 207 ELSE ! T pts j-indice on the rim on the ocean next to the rim on T and U points 208 IF( flagu == 1. ) THEN ; iiTrim = ii ; iiToce = ii+1 ; iiUoce = ii+1 ; ENDIF 209 IF( flagu == -1. ) THEN ; iiTrim = ii+1 ; iiToce = ii ; iiUoce = ii-1 ; ENDIF 210 ! 211 ! Rare case : rim is parallel to the mpi subdomain border and on the halo : point will be received 212 IF( iiTrim > jpi .OR. iiToce > jpi .OR. iiUoce > jpi .OR. iiUoce < 1 ) CYCLE 213 ! 214 zfla = pua2d(iiUoce,ij) - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iiToce,ij) - sshdta(iiTrim,ij) ) 215 ! 216 ! jchanut tschanges, use characteristics method (Blayo et Debreu, 2005) : 217 ! mix Flather scheme with velocity of the ocean next to the rim 218 pua2d(ii,ij) = z1_2 * ( dta%u2d(jb) + zfla ) 219 END IF 217 220 END DO 218 221 ! … … 220 223 ! ! remember that flagv=-1 if normal velocity direction is outward 221 224 DO jb = 1, idx%nblenrim(igrd) 222 ii 223 ij 225 ii = idx%nbi(jb,igrd) 226 ij = idx%nbj(jb,igrd) 224 227 flagv => idx%flagv(jb,igrd) 225 ijm1 = ij + MAX( 0, INT( flagv ) ) ! T pts j-indice inside the boundary 226 ijp1 = ij - MIN( 0, INT( flagv ) ) ! T pts j-indice outside the boundary 227 IF( ijm1 > jpj .OR. ijp1 > jpj ) CYCLE 228 ! 229 zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 230 231 ! jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme 232 ! Use characteristics method instead 233 zflag = ABS(flagv) 234 zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1) 235 pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1) 228 IF( flagv == 0. ) THEN 229 pva2d(ii,ij) = dta%v2d(jb) 230 ELSE ! T pts j-indice on the rim on the ocean next to the rim on T and V points 231 IF( flagv == 1. ) THEN ; ijTrim = ij ; ijToce = ij+1 ; ijVoce = ij+1 ; ENDIF 232 IF( flagv == -1. ) THEN ; ijTrim = ij+1 ; ijToce = ij ; ijVoce = ij-1 ; ENDIF 233 ! 234 ! Rare case : rim is parallel to the mpi subdomain border and on the halo : point will be received 235 IF( ijTrim > jpj .OR. ijToce > jpj .OR. ijVoce > jpj .OR. ijVoce < 1 ) CYCLE 236 ! 237 zfla = pva2d(ii,ijVoce) - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii,ijToce) - sshdta(ii,ijTrim) ) 238 ! 239 ! jchanut tschanges, use characteristics method (Blayo et Debreu, 2005) : 240 ! mix Flather scheme with velocity of the ocean next to the rim 241 pva2d(ii,ij) = z1_2 * ( dta%v2d(jb) + zfla ) 242 END IF 236 243 END DO 237 244 !
Note: See TracChangeset
for help on using the changeset viewer.