Changeset 11831 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/BDY/bdydyn2d.F90
- Timestamp:
- 2019-10-29T18:14:49+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/BDY/bdydyn2d.F90
r11072 r11831 14 14 !! bdy_ssh : Duplicate sea level across open boundaries 15 15 !!---------------------------------------------------------------------- 16 USE oce ! ocean dynamics and tracers17 16 USE dom_oce ! ocean space and time domain 18 17 USE bdy_oce ! ocean open boundary conditions … … 50 49 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh 51 50 !! 52 INTEGER :: ib_bdy ! Loop counter 53 54 DO ib_bdy=1, nb_bdy 55 56 SELECT CASE( cn_dyn2d(ib_bdy) ) 57 CASE('none') 58 CYCLE 59 CASE('frs') 60 CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, pua2d, pva2d ) 61 CASE('flather') 62 CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, pua2d, pva2d, pssh, phur, phvr ) 63 CASE('orlanski') 64 CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, & 65 & pua2d, pva2d, pub2d, pvb2d, ll_npo=.false.) 66 CASE('orlanski_npo') 67 CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, & 68 & pua2d, pva2d, pub2d, pvb2d, ll_npo=.true. ) 69 CASE DEFAULT 70 CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' ) 71 END SELECT 72 ENDDO 73 51 INTEGER :: ib_bdy, ir ! BDY set index, rim index 52 LOGICAL :: llrim0 ! indicate if rim 0 is treated 53 LOGICAL, DIMENSION(4) :: llsend2, llrecv2, llsend3, llrecv3 ! indicate how communications are to be carried out 54 55 llsend2(:) = .false. ; llrecv2(:) = .false. 56 llsend3(:) = .false. ; llrecv3(:) = .false. 57 DO ir = 1, 0, -1 ! treat rim 1 before rim 0 58 IF( ir == 0 ) THEN ; llrim0 = .TRUE. 59 ELSE ; llrim0 = .FALSE. 60 END IF 61 DO ib_bdy=1, nb_bdy 62 SELECT CASE( cn_dyn2d(ib_bdy) ) 63 CASE('none') 64 CYCLE 65 CASE('frs') ! treat the whole boundary at once 66 IF( llrim0 ) CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, pua2d, pva2d ) 67 CASE('flather') 68 CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, pua2d, pva2d, pssh, phur, phvr, llrim0 ) 69 CASE('orlanski') 70 CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, & 71 & pua2d, pva2d, pub2d, pvb2d, llrim0, ll_npo=.false. ) 72 CASE('orlanski_npo') 73 CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, & 74 & pua2d, pva2d, pub2d, pvb2d, llrim0, ll_npo=.true. ) 75 CASE DEFAULT 76 CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' ) 77 END SELECT 78 ENDDO 79 ! 80 IF( nn_hls > 1 .AND. ir == 1 ) CYCLE ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0 81 IF( nn_hls == 1 ) THEN 82 llsend2(:) = .false. ; llrecv2(:) = .false. 83 llsend3(:) = .false. ; llrecv3(:) = .false. 84 END IF 85 DO ib_bdy=1, nb_bdy 86 SELECT CASE( cn_dyn2d(ib_bdy) ) 87 CASE('flather') 88 llsend2(1:2) = llsend2(1:2) .OR. lsend_bdyint(ib_bdy,2,1:2,ir) ! west/east, U points 89 llsend2(1) = llsend2(1) .OR. lsend_bdyext(ib_bdy,2,1,ir) ! neighbour might search point towards its east bdy 90 llrecv2(1:2) = llrecv2(1:2) .OR. lrecv_bdyint(ib_bdy,2,1:2,ir) ! west/east, U points 91 llrecv2(2) = llrecv2(2) .OR. lrecv_bdyext(ib_bdy,2,2,ir) ! might search point towards bdy on the east 92 llsend3(3:4) = llsend3(3:4) .OR. lsend_bdyint(ib_bdy,3,3:4,ir) ! north/south, V points 93 llsend3(3) = llsend3(3) .OR. lsend_bdyext(ib_bdy,3,3,ir) ! neighbour might search point towards its north bdy 94 llrecv3(3:4) = llrecv3(3:4) .OR. lrecv_bdyint(ib_bdy,3,3:4,ir) ! north/south, V points 95 llrecv3(4) = llrecv3(4) .OR. lrecv_bdyext(ib_bdy,3,4,ir) ! might search point towards bdy on the north 96 CASE('orlanski', 'orlanski_npo') 97 llsend2(:) = llsend2(:) .OR. lsend_bdy(ib_bdy,2,:,ir) ! possibly every direction, U points 98 llrecv2(:) = llrecv2(:) .OR. lrecv_bdy(ib_bdy,2,:,ir) ! possibly every direction, U points 99 llsend3(:) = llsend3(:) .OR. lsend_bdy(ib_bdy,3,:,ir) ! possibly every direction, V points 100 llrecv3(:) = llrecv3(:) .OR. lrecv_bdy(ib_bdy,3,:,ir) ! possibly every direction, V points 101 END SELECT 102 END DO 103 IF( ANY(llsend2) .OR. ANY(llrecv2) ) THEN ! if need to send/recv in at least one direction 104 CALL lbc_lnk( 'bdydyn2d', pua2d, 'U', -1., kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 ) 105 END IF 106 IF( ANY(llsend3) .OR. ANY(llrecv3) ) THEN ! if need to send/recv in at least one direction 107 CALL lbc_lnk( 'bdydyn2d', pva2d, 'V', -1., kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 ) 108 END IF 109 ! 110 END DO ! ir 111 ! 74 112 END SUBROUTINE bdy_dyn2d 75 113 … … 90 128 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d 91 129 !! 92 INTEGER :: jb , jk! dummy loop indices130 INTEGER :: jb ! dummy loop indices 93 131 INTEGER :: ii, ij, igrd ! local integers 94 132 REAL(wp) :: zwgt ! boundary weight … … 110 148 pva2d(ii,ij) = ( pva2d(ii,ij) + zwgt * ( dta%v2d(jb) - pva2d(ii,ij) ) ) * vmask(ii,ij,1) 111 149 END DO 112 CALL lbc_bdy_lnk( 'bdydyn2d', pua2d, 'U', -1., ib_bdy )113 CALL lbc_bdy_lnk( 'bdydyn2d', pva2d, 'V', -1., ib_bdy) ! Boundary points should be updated114 150 ! 115 151 END SUBROUTINE bdy_dyn2d_frs 116 152 117 153 118 SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy, pua2d, pva2d, pssh, phur, phvr )154 SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy, pua2d, pva2d, pssh, phur, phvr, llrim0 ) 119 155 !!---------------------------------------------------------------------- 120 156 !! *** SUBROUTINE bdy_dyn2d_fla *** … … 139 175 INTEGER, INTENT(in) :: ib_bdy ! BDY set index 140 176 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d 141 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh, phur, phvr 142 177 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh, phur, phvr 178 LOGICAL , INTENT(in) :: llrim0 ! indicate if rim 0 is treated 179 INTEGER :: ibeg, iend ! length of rim to be treated (rim 0 or rim 1) 143 180 INTEGER :: jb, igrd ! dummy loop indices 144 INTEGER :: ii, ij, iim1, iip1, ijm1, ijp1 ! 2D addresses 145 REAL(wp), POINTER :: flagu, flagv ! short cuts 146 REAL(wp) :: zcorr ! Flather correction 147 REAL(wp) :: zforc ! temporary scalar 148 REAL(wp) :: zflag, z1_2 ! " " 181 INTEGER :: ii, ij ! 2D addresses 182 INTEGER :: iiTrim, ijTrim ! T pts i/j-indice on the rim 183 INTEGER :: iiToce, ijToce, iiUoce, ijVoce ! T, U and V pts i/j-indice of the ocean next to the rim 184 REAL(wp) :: flagu, flagv ! short cuts 185 REAL(wp) :: zfla ! Flather correction 186 REAL(wp) :: z1_2 ! 187 REAL(wp), DIMENSION(jpi,jpj) :: sshdta ! 2D version of dta%ssh 149 188 !!---------------------------------------------------------------------- 150 189 … … 153 192 ! ---------------------------------! 154 193 ! Flather boundary conditions :! 155 ! ---------------------------------! 156 157 !!! REPLACE spgu with nemo_wrk work space 158 159 ! Fill temporary array with ssh data (here spgu): 194 ! ---------------------------------! 195 196 ! Fill temporary array with ssh data (here we use spgu with the alias sshdta): 160 197 igrd = 1 161 spgu(:,:) = 0.0 162 DO jb = 1, idx%nblenrim(igrd) 198 IF( llrim0 ) THEN ; ibeg = 1 ; iend = idx%nblenrim0(igrd) 199 ELSE ; ibeg = idx%nblenrim0(igrd)+1 ; iend = idx%nblenrim(igrd) 200 END IF 201 ! 202 DO jb = ibeg, iend 163 203 ii = idx%nbi(jb,igrd) 164 204 ij = idx%nbj(jb,igrd) 165 IF( ll_wd ) THEN 166 spgu(ii, ij) = dta%ssh(jb) - ssh_ref 167 ELSE 168 spgu(ii, ij) = dta%ssh(jb) 205 IF( ll_wd ) THEN ; sshdta(ii, ij) = dta%ssh(jb) - ssh_ref 206 ELSE ; sshdta(ii, ij) = dta%ssh(jb) 169 207 ENDIF 170 208 END DO 171 172 CALL lbc_bdy_lnk( 'bdydyn2d', spgu(:,:), 'T', 1., ib_bdy ) 173 ! 174 igrd = 2 ! Flather bc on u-velocity; 209 ! 210 igrd = 2 ! Flather bc on u-velocity 175 211 ! ! remember that flagu=-1 if normal velocity direction is outward 176 212 ! ! I think we should rather use after ssh ? 177 DO jb = 1, idx%nblenrim(igrd) 178 ii = idx%nbi(jb,igrd) 179 ij = idx%nbj(jb,igrd) 180 flagu => idx%flagu(jb,igrd) 181 iim1 = ii + MAX( 0, INT( flagu ) ) ! T pts i-indice inside the boundary 182 iip1 = ii - MIN( 0, INT( flagu ) ) ! T pts i-indice outside the boundary 183 ! 184 zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 185 186 ! jchanut tschanges: Set zflag to 0 below to revert to Flather scheme 187 ! Use characteristics method instead 188 zflag = ABS(flagu) 189 zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(ii+NINT(flagu),ij) 190 pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 213 IF( llrim0 ) THEN ; ibeg = 1 ; iend = idx%nblenrim0(igrd) 214 ELSE ; ibeg = idx%nblenrim0(igrd)+1 ; iend = idx%nblenrim(igrd) 215 END IF 216 DO jb = ibeg, iend 217 ii = idx%nbi(jb,igrd) 218 ij = idx%nbj(jb,igrd) 219 flagu = idx%flagu(jb,igrd) 220 IF( flagu == 0. ) THEN 221 pua2d(ii,ij) = dta%u2d(jb) 222 ELSE ! T pts j-indice on the rim on the ocean next to the rim on T and U points 223 IF( flagu == 1. ) THEN ; iiTrim = ii ; iiToce = ii+1 ; iiUoce = ii+1 ; ENDIF 224 IF( flagu == -1. ) THEN ; iiTrim = ii+1 ; iiToce = ii ; iiUoce = ii-1 ; ENDIF 225 ! 226 ! Rare case : rim is parallel to the mpi subdomain border and on the halo : point will be received 227 IF( iiTrim > jpi .OR. iiToce > jpi .OR. iiUoce > jpi .OR. iiUoce < 1 ) CYCLE 228 ! 229 zfla = dta%u2d(jb) - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iiToce,ij) - sshdta(iiTrim,ij) ) 230 ! 231 ! jchanut tschanges, use characteristics method (Blayo et Debreu, 2005) : 232 ! mix Flather scheme with velocity of the ocean next to the rim 233 pua2d(ii,ij) = z1_2 * ( pua2d(iiUoce,ij) + zfla ) 234 END IF 191 235 END DO 192 236 ! 193 237 igrd = 3 ! Flather bc on v-velocity 194 238 ! ! remember that flagv=-1 if normal velocity direction is outward 195 DO jb = 1, idx%nblenrim(igrd) 196 ii = idx%nbi(jb,igrd) 197 ij = idx%nbj(jb,igrd) 198 flagv => idx%flagv(jb,igrd) 199 ijm1 = ij + MAX( 0, INT( flagv ) ) ! T pts j-indice inside the boundary 200 ijp1 = ij - MIN( 0, INT( flagv ) ) ! T pts j-indice outside the boundary 201 ! 202 zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 203 204 ! jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme 205 ! Use characteristics method instead 206 zflag = ABS(flagv) 207 zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ij+NINT(flagv)) 208 pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1) 209 END DO 210 CALL lbc_bdy_lnk( 'bdydyn2d', pua2d, 'U', -1., ib_bdy ) ! Boundary points should be updated 211 CALL lbc_bdy_lnk( 'bdydyn2d', pva2d, 'V', -1., ib_bdy ) ! 239 IF( llrim0 ) THEN ; ibeg = 1 ; iend = idx%nblenrim0(igrd) 240 ELSE ; ibeg = idx%nblenrim0(igrd)+1 ; iend = idx%nblenrim(igrd) 241 END IF 242 DO jb = ibeg, iend 243 ii = idx%nbi(jb,igrd) 244 ij = idx%nbj(jb,igrd) 245 flagv = idx%flagv(jb,igrd) 246 IF( flagv == 0. ) THEN 247 pva2d(ii,ij) = dta%v2d(jb) 248 ELSE ! T pts j-indice on the rim on the ocean next to the rim on T and V points 249 IF( flagv == 1. ) THEN ; ijTrim = ij ; ijToce = ij+1 ; ijVoce = ij+1 ; ENDIF 250 IF( flagv == -1. ) THEN ; ijTrim = ij+1 ; ijToce = ij ; ijVoce = ij-1 ; ENDIF 251 ! 252 ! Rare case : rim is parallel to the mpi subdomain border and on the halo : point will be received 253 IF( ijTrim > jpj .OR. ijToce > jpj .OR. ijVoce > jpj .OR. ijVoce < 1 ) CYCLE 254 ! 255 zfla = dta%v2d(jb) - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii,ijToce) - sshdta(ii,ijTrim) ) 256 ! 257 ! jchanut tschanges, use characteristics method (Blayo et Debreu, 2005) : 258 ! mix Flather scheme with velocity of the ocean next to the rim 259 pva2d(ii,ij) = z1_2 * ( pva2d(ii,ijVoce) + zfla ) 260 END IF 261 END DO 212 262 ! 213 263 END SUBROUTINE bdy_dyn2d_fla 214 264 215 265 216 SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, pua2d, pva2d, pub2d, pvb2d, ll _npo )266 SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, pua2d, pva2d, pub2d, pvb2d, llrim0, ll_npo ) 217 267 !!---------------------------------------------------------------------- 218 268 !! *** SUBROUTINE bdy_dyn2d_orlanski *** … … 231 281 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub2d, pvb2d 232 282 LOGICAL, INTENT(in) :: ll_npo ! flag for NPO version 233 283 LOGICAL, INTENT(in) :: llrim0 ! indicate if rim 0 is treated 234 284 INTEGER :: ib, igrd ! dummy loop indices 235 285 INTEGER :: ii, ij, iibm1, ijbm1 ! indices … … 238 288 igrd = 2 ! Orlanski bc on u-velocity; 239 289 ! 240 CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, ll _npo )290 CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, llrim0, ll_npo ) 241 291 242 292 igrd = 3 ! Orlanski bc on v-velocity 243 293 ! 244 CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo ) 245 ! 246 CALL lbc_bdy_lnk( 'bdydyn2d', pua2d, 'U', -1., ib_bdy ) ! Boundary points should be updated 247 CALL lbc_bdy_lnk( 'bdydyn2d', pva2d, 'V', -1., ib_bdy ) ! 294 CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, llrim0, ll_npo ) 248 295 ! 249 296 END SUBROUTINE bdy_dyn2d_orlanski … … 257 304 !! 258 305 !!---------------------------------------------------------------------- 259 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zssh ! Sea level 260 !! 261 INTEGER :: ib_bdy, ib, igrd ! local integers 262 INTEGER :: ii, ij, zcoef, zcoef1, zcoef2, ip, jp ! " " 263 264 igrd = 1 ! Everything is at T-points here 265 266 DO ib_bdy = 1, nb_bdy 267 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 268 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 269 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 270 ! Set gradient direction: 271 zcoef1 = bdytmask(ii-1,ij ) + bdytmask(ii+1,ij ) 272 zcoef2 = bdytmask(ii ,ij-1) + bdytmask(ii ,ij+1) 273 IF ( zcoef1+zcoef2 == 0 ) THEN ! corner 274 zcoef = bdytmask(ii-1,ij-1) + bdytmask(ii+1,ij+1) + bdytmask(ii+1,ij-1) + bdytmask(ii-1,ij+1) 275 zssh(ii,ij) = zssh( ii-1, ij-1 ) * bdytmask( ii-1, ij-1) + & 276 & zssh( ii+1, ij+1 ) * bdytmask( ii+1, ij+1) + & 277 & zssh( ii+1, ij-1 ) * bdytmask( ii+1, ij-1) + & 278 & zssh( ii-1, ij+1 ) * bdytmask( ii-1, ij+1) 279 zssh(ii,ij) = ( zssh(ii,ij) / MAX( 1, zcoef) ) * tmask(ii,ij,1) 280 ELSE 281 ip = bdytmask(ii+1,ij ) - bdytmask(ii-1,ij ) 282 jp = bdytmask(ii ,ij+1) - bdytmask(ii ,ij-1) 283 zssh(ii,ij) = zssh(ii+ip,ij+jp) * tmask(ii+ip,ij+jp,1) 284 ENDIF 306 REAL(wp), DIMENSION(jpi,jpj,1), INTENT(inout) :: zssh ! Sea level, need 3 dimensions to be used by bdy_nmn 307 !! 308 INTEGER :: ib_bdy, ir ! bdy index, rim index 309 INTEGER :: ibeg, iend ! length of rim to be treated (rim 0 or rim 1) 310 LOGICAL :: llrim0 ! indicate if rim 0 is treated 311 LOGICAL, DIMENSION(4) :: llsend1, llrecv1 ! indicate how communications are to be carried out 312 !!---------------------------------------------------------------------- 313 llsend1(:) = .false. ; llrecv1(:) = .false. 314 DO ir = 1, 0, -1 ! treat rim 1 before rim 0 315 IF( nn_hls == 1 ) THEN ; llsend1(:) = .false. ; llrecv1(:) = .false. ; END IF 316 IF( ir == 0 ) THEN ; llrim0 = .TRUE. 317 ELSE ; llrim0 = .FALSE. 318 END IF 319 DO ib_bdy = 1, nb_bdy 320 CALL bdy_nmn( idx_bdy(ib_bdy), 1, zssh, llrim0 ) ! zssh is masked 321 llsend1(:) = llsend1(:) .OR. lsend_bdyint(ib_bdy,1,:,ir) ! possibly every direction, T points 322 llrecv1(:) = llrecv1(:) .OR. lrecv_bdyint(ib_bdy,1,:,ir) ! possibly every direction, T points 285 323 END DO 286 287 ! Boundary points should be updated 288 CALL lbc_bdy_lnk( 'bdydyn2d', zssh(:,:), 'T', 1., ib_bdy ) 289 END DO 290 324 IF( nn_hls > 1 .AND. ir == 1 ) CYCLE ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0 325 IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN ! if need to send/recv in at least one direction 326 CALL lbc_lnk( 'bdydyn2d', zssh(:,:,1), 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 327 END IF 328 END DO 329 ! 291 330 END SUBROUTINE bdy_ssh 292 331
Note: See TracChangeset
for help on using the changeset viewer.