Changeset 3294 for trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
- Timestamp:
- 2012-01-28T17:44:18+01:00 (12 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
r2715 r3294 10 10 !! 3.3 ! 2010-09 (E.O'Dea) updates for Shelf configurations 11 11 !! 3.3 ! 2010-09 (D.Storkey) add ice boundary conditions 12 !! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge 12 13 !!---------------------------------------------------------------------- 13 14 #if defined key_bdy … … 17 18 !! bdy_init : Initialization of unstructured open boundaries 18 19 !!---------------------------------------------------------------------- 20 USE timing ! Timing 19 21 USE oce ! ocean dynamics and tracers variables 20 22 USE dom_oce ! ocean space and time domain 21 USE obc_par ! ocean open boundary conditions22 23 USE bdy_oce ! unstructured open boundary conditions 23 USE bdydta, ONLY: bdy_dta_alloc ! open boundary data24 USE bdytides ! tides at open boundaries initialization (tide_init routine)25 24 USE in_out_manager ! I/O units 26 25 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 31 30 PRIVATE 32 31 33 PUBLIC bdy_init ! routine called by opa.F9032 PUBLIC bdy_init ! routine called in nemo_init 34 33 35 34 !!---------------------------------------------------------------------- … … 52 51 !! ** Input : bdy_init.nc, input file for unstructured open boundaries 53 52 !!---------------------------------------------------------------------- 54 INTEGER :: ii, ij, ik, igrd, ib, ir ! dummy loop indices 55 INTEGER :: icount, icountr, ib_len, ibr_max ! local integers 56 INTEGER :: iw, ie, is, in, inum, id_dummy ! - - 57 INTEGER :: igrd_start, igrd_end ! - - 58 REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars 59 INTEGER, DIMENSION (2) :: kdimsz 60 INTEGER, DIMENSION(jpbdta, jpbgrd) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta 61 INTEGER, DIMENSION(jpbdta, jpbgrd) :: nbrdta ! Discrete distance from rim points 62 REAL(wp), DIMENSION(jpidta,jpjdta) :: zmask ! global domain mask 63 REAL(wp), DIMENSION(jpbdta,1) :: zdta ! temporary array 64 CHARACTER(LEN=80),DIMENSION(6) :: clfile 53 ! namelist variables 54 !------------------- 55 INTEGER, PARAMETER :: jp_nseg = 100 56 INTEGER :: nbdysege, nbdysegw, nbdysegn, nbdysegs 57 INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft 58 INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft 59 INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft 60 INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft 61 62 ! local variables 63 !------------------- 64 INTEGER :: ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices 65 INTEGER :: icount, icountr, ibr_max, ilen1, ibm1 ! local integers 66 INTEGER :: iw, ie, is, in, inum, id_dummy ! - - 67 INTEGER :: igrd_start, igrd_end, jpbdta ! - - 68 INTEGER, POINTER :: nbi, nbj, nbr ! short cuts 69 REAL , POINTER :: flagu, flagv ! - - 70 REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars 71 INTEGER, DIMENSION (2) :: kdimsz 72 INTEGER, DIMENSION(jpbgrd,jp_bdy) :: nblendta ! Length of index arrays 73 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta 74 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points 75 REAL(wp), DIMENSION(jpidta,jpjdta) :: zmask ! global domain mask 76 CHARACTER(LEN=80),DIMENSION(jpbgrd) :: clfile 77 CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid 65 78 !! 66 NAMELIST/nambdy/cn_mask, cn_dta_frs_T, cn_dta_frs_U, cn_dta_frs_V, & 67 & cn_dta_fla_T, cn_dta_fla_U, cn_dta_fla_V, & 68 & ln_tides, ln_clim, ln_vol, ln_mask, & 69 & ln_dyn_fla, ln_dyn_frs, ln_tra_frs,ln_ice_frs, & 70 & nn_dtactl, nn_rimwidth, nn_volctl 79 NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file, & 80 & ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, & 81 & nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta, & 82 #if defined key_lim2 83 & nn_ice_lim2, nn_ice_lim2_dta, & 84 #endif 85 & ln_vol, nn_volctl, nn_rimwidth 86 !! 87 NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft, & 88 nbdysegw, jpiwob, jpjwdt, jpjwft, & 89 nbdysegn, jpjnob, jpindt, jpinft, & 90 nbdysegs, jpjsob, jpisdt, jpisft 91 71 92 !!---------------------------------------------------------------------- 72 93 94 IF( nn_timing == 1 ) CALL timing_start('bdy_init') 95 96 IF( bdy_oce_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate oce arrays' ) 97 73 98 IF(lwp) WRITE(numout,*) 74 IF(lwp) WRITE(numout,*) 'bdy_init : initialization of unstructuredopen boundaries'99 IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries' 75 100 IF(lwp) WRITE(numout,*) '~~~~~~~~' 76 101 ! 77 ! ! allocate bdy_oce arrays78 IF( bdy_oce_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate oce arrays' )79 IF( bdy_dta_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate dta arrays' )80 102 81 103 IF( jperio /= 0 ) CALL ctl_stop( 'Cyclic or symmetric,', & 82 & ' and unstructured open boundary condition are not compatible' ) 83 84 IF( lk_obc ) CALL ctl_stop( 'Straight open boundaries,', & 85 & ' and unstructured open boundaries are not compatible' ) 86 87 ! --------------------------- 88 REWIND( numnam ) ! Read namelist parameters 104 & ' and general open boundary condition are not compatible' ) 105 106 cgrid= (/'t','u','v'/) 107 108 ! ----------------------------------------- 109 ! Initialise and read namelist parameters 110 ! ----------------------------------------- 111 112 nb_bdy = 0 113 ln_coords_file(:) = .false. 114 cn_coords_file(:) = '' 115 ln_mask_file = .false. 116 cn_mask_file(:) = '' 117 nn_dyn2d(:) = 0 118 nn_dyn2d_dta(:) = -1 ! uninitialised flag 119 nn_dyn3d(:) = 0 120 nn_dyn3d_dta(:) = -1 ! uninitialised flag 121 nn_tra(:) = 0 122 nn_tra_dta(:) = -1 ! uninitialised flag 123 #if defined key_lim2 124 nn_ice_lim2(:) = 0 125 nn_ice_lim2_dta(:)= -1 ! uninitialised flag 126 #endif 127 ln_vol = .false. 128 nn_volctl = -1 ! uninitialised flag 129 nn_rimwidth(:) = -1 ! uninitialised flag 130 131 REWIND( numnam ) 89 132 READ ( numnam, nambdy ) 133 134 ! ----------------------------------------- 135 ! Check and write out namelist parameters 136 ! ----------------------------------------- 90 137 91 138 ! ! control prints 92 139 IF(lwp) WRITE(numout,*) ' nambdy' 93 140 94 ! ! check type of data used (nn_dtactl value) 95 IF(lwp) WRITE(numout,*) 'nn_dtactl =', nn_dtactl 96 IF(lwp) WRITE(numout,*) 97 SELECT CASE( nn_dtactl ) ! 98 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data' 99 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file' 100 CASE DEFAULT ; CALL ctl_stop( 'nn_dtactl must be 0 or 1' ) 101 END SELECT 102 103 IF(lwp) WRITE(numout,*) 104 IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nn_rimwidth = ', nn_rimwidth 105 106 IF(lwp) WRITE(numout,*) 107 IF(lwp) WRITE(numout,*) ' nn_volctl = ', nn_volctl 108 109 IF( ln_vol ) THEN ! check volume conservation (nn_volctl value) 110 SELECT CASE ( nn_volctl ) 141 IF( nb_bdy .eq. 0 ) THEN 142 IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.' 143 ELSE 144 IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_bdy 145 ENDIF 146 147 DO ib_bdy = 1,nb_bdy 148 IF(lwp) WRITE(numout,*) ' ' 149 IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------' 150 151 IF( ln_coords_file(ib_bdy) ) THEN 152 IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy)) 153 ELSE 154 IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.' 155 ENDIF 156 IF(lwp) WRITE(numout,*) 157 158 IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution: ' 159 SELECT CASE( nn_dyn2d(ib_bdy) ) 160 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 161 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 162 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Flather radiation condition' 163 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) 164 END SELECT 165 IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN 166 SELECT CASE( nn_dyn2d_dta(ib_bdy) ) ! 167 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data' 168 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file' 169 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' tidal harmonic forcing taken from file' 170 CASE( 3 ) ; IF(lwp) WRITE(numout,*) ' boundary data AND tidal harmonic forcing taken from files' 171 CASE DEFAULT ; CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 172 END SELECT 173 ENDIF 174 IF(lwp) WRITE(numout,*) 175 176 IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities: ' 177 SELECT CASE( nn_dyn3d(ib_bdy) ) 178 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 179 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 180 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn3d' ) 181 END SELECT 182 IF( nn_dyn3d(ib_bdy) .gt. 0 ) THEN 183 SELECT CASE( nn_dyn3d_dta(ib_bdy) ) ! 184 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data' 185 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file' 186 CASE DEFAULT ; CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' ) 187 END SELECT 188 ENDIF 189 IF(lwp) WRITE(numout,*) 190 191 IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity: ' 192 SELECT CASE( nn_tra(ib_bdy) ) 193 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 194 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 195 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_tra' ) 196 END SELECT 197 IF( nn_tra(ib_bdy) .gt. 0 ) THEN 198 SELECT CASE( nn_tra_dta(ib_bdy) ) ! 199 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data' 200 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file' 201 CASE DEFAULT ; CALL ctl_stop( 'nn_tra_dta must be 0 or 1' ) 202 END SELECT 203 ENDIF 204 IF(lwp) WRITE(numout,*) 205 206 #if defined key_lim2 207 IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice: ' 208 SELECT CASE( nn_ice_lim2(ib_bdy) ) 209 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 210 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 211 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_tra' ) 212 END SELECT 213 IF( nn_ice_lim2(ib_bdy) .gt. 0 ) THEN 214 SELECT CASE( nn_ice_lim2_dta(ib_bdy) ) ! 215 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data' 216 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file' 217 CASE DEFAULT ; CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' ) 218 END SELECT 219 ENDIF 220 IF(lwp) WRITE(numout,*) 221 #endif 222 223 IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_bdy) 224 IF(lwp) WRITE(numout,*) 225 226 ENDDO 227 228 IF( ln_vol ) THEN ! check volume conservation (nn_volctl value) 229 IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 230 IF(lwp) WRITE(numout,*) 231 SELECT CASE ( nn_volctl ) 111 232 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' The total volume will be constant' 112 233 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' The total volume will vary according to the surface E-P flux' 113 234 CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 114 END SELECT 115 IF(lwp) WRITE(numout,*) 116 ELSE 117 IF(lwp) WRITE(numout,*) 'No volume correction with unstructured open boundaries' 118 IF(lwp) WRITE(numout,*) 119 ENDIF 120 121 IF( ln_tides ) THEN 122 IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing at unstructured open boundaries' 123 IF(lwp) WRITE(numout,*) 124 ENDIF 125 126 IF( ln_dyn_fla ) THEN 127 IF(lwp) WRITE(numout,*) 'Flather condition on U, V at unstructured open boundaries' 128 IF(lwp) WRITE(numout,*) 129 ENDIF 130 131 IF( ln_dyn_frs ) THEN 132 IF(lwp) WRITE(numout,*) 'FRS condition on U and V at unstructured open boundaries' 133 IF(lwp) WRITE(numout,*) 134 ENDIF 135 136 IF( ln_tra_frs ) THEN 137 IF(lwp) WRITE(numout,*) 'FRS condition on T & S fields at unstructured open boundaries' 138 IF(lwp) WRITE(numout,*) 139 ENDIF 140 141 IF( ln_ice_frs ) THEN 142 IF(lwp) WRITE(numout,*) 'FRS condition on ice fields at unstructured open boundaries' 143 IF(lwp) WRITE(numout,*) 144 ENDIF 145 146 IF( ln_tides ) CALL tide_init ! Read tides namelist 147 148 149 ! Read arrays defining unstructured open boundaries 235 END SELECT 236 IF(lwp) WRITE(numout,*) 237 ELSE 238 IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 239 IF(lwp) WRITE(numout,*) 240 ENDIF 241 150 242 ! ------------------------------------------------- 243 ! Initialise indices arrays for open boundaries 244 ! ------------------------------------------------- 245 246 ! Work out global dimensions of boundary data 247 ! --------------------------------------------- 248 REWIND( numnam ) 249 DO ib_bdy = 1, nb_bdy 250 251 jpbdta = 1 252 IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 253 254 ! No REWIND here because may need to read more than one nambdy_index namelist. 255 READ ( numnam, nambdy_index ) 256 257 ! Automatic boundary definition: if nbdysegX = -1 258 ! set boundary to whole side of model domain. 259 IF( nbdysege == -1 ) THEN 260 nbdysege = 1 261 jpieob(1) = jpiglo - 1 262 jpjedt(1) = 2 263 jpjeft(1) = jpjglo - 1 264 ENDIF 265 IF( nbdysegw == -1 ) THEN 266 nbdysegw = 1 267 jpiwob(1) = 2 268 jpjwdt(1) = 2 269 jpjwft(1) = jpjglo - 1 270 ENDIF 271 IF( nbdysegn == -1 ) THEN 272 nbdysegn = 1 273 jpjnob(1) = jpjglo - 1 274 jpindt(1) = 2 275 jpinft(1) = jpiglo - 1 276 ENDIF 277 IF( nbdysegs == -1 ) THEN 278 nbdysegs = 1 279 jpjsob(1) = 2 280 jpisdt(1) = 2 281 jpisft(1) = jpiglo - 1 282 ENDIF 283 284 nblendta(:,ib_bdy) = 0 285 DO iseg = 1, nbdysege 286 igrd = 1 287 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1 288 igrd = 2 289 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1 290 igrd = 3 291 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) 292 ENDDO 293 DO iseg = 1, nbdysegw 294 igrd = 1 295 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1 296 igrd = 2 297 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1 298 igrd = 3 299 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) 300 ENDDO 301 DO iseg = 1, nbdysegn 302 igrd = 1 303 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1 304 igrd = 2 305 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) 306 igrd = 3 307 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1 308 ENDDO 309 DO iseg = 1, nbdysegs 310 igrd = 1 311 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1 312 igrd = 2 313 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) 314 igrd = 3 315 nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1 316 ENDDO 317 318 nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy) 319 jpbdta = MAXVAL(nblendta(:,ib_bdy)) 320 321 322 ELSE ! Read size of arrays in boundary coordinates file. 323 324 325 CALL iom_open( cn_coords_file(ib_bdy), inum ) 326 jpbdta = 1 327 DO igrd = 1, jpbgrd 328 id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 329 nblendta(igrd,ib_bdy) = kdimsz(1) 330 jpbdta = MAX(jpbdta, kdimsz(1)) 331 ENDDO 332 333 ENDIF 334 335 ENDDO ! ib_bdy 336 337 ! Allocate arrays 338 !--------------- 339 ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), & 340 & nbrdta(jpbdta, jpbgrd, nb_bdy) ) 341 342 ALLOCATE( dta_global(jpbdta, 1, jpk) ) 343 344 ! Calculate global boundary index arrays or read in from file 345 !------------------------------------------------------------ 346 REWIND( numnam ) 347 DO ib_bdy = 1, nb_bdy 348 349 IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters 350 351 ! No REWIND here because may need to read more than one nambdy_index namelist. 352 READ ( numnam, nambdy_index ) 353 354 ! Automatic boundary definition: if nbdysegX = -1 355 ! set boundary to whole side of model domain. 356 IF( nbdysege == -1 ) THEN 357 nbdysege = 1 358 jpieob(1) = jpiglo - 1 359 jpjedt(1) = 2 360 jpjeft(1) = jpjglo - 1 361 ENDIF 362 IF( nbdysegw == -1 ) THEN 363 nbdysegw = 1 364 jpiwob(1) = 2 365 jpjwdt(1) = 2 366 jpjwft(1) = jpjglo - 1 367 ENDIF 368 IF( nbdysegn == -1 ) THEN 369 nbdysegn = 1 370 jpjnob(1) = jpjglo - 1 371 jpindt(1) = 2 372 jpinft(1) = jpiglo - 1 373 ENDIF 374 IF( nbdysegs == -1 ) THEN 375 nbdysegs = 1 376 jpjsob(1) = 2 377 jpisdt(1) = 2 378 jpisft(1) = jpiglo - 1 379 ENDIF 380 381 ! ------------ T points ------------- 382 igrd = 1 383 icount = 0 384 DO ir = 1, nn_rimwidth(ib_bdy) 385 ! east 386 DO iseg = 1, nbdysege 387 DO ij = jpjedt(iseg), jpjeft(iseg) 388 icount = icount + 1 389 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 390 nbjdta(icount, igrd, ib_bdy) = ij 391 nbrdta(icount, igrd, ib_bdy) = ir 392 ENDDO 393 ENDDO 394 ! west 395 DO iseg = 1, nbdysegw 396 DO ij = jpjwdt(iseg), jpjwft(iseg) 397 icount = icount + 1 398 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 399 nbjdta(icount, igrd, ib_bdy) = ij 400 nbrdta(icount, igrd, ib_bdy) = ir 401 ENDDO 402 ENDDO 403 ! north 404 DO iseg = 1, nbdysegn 405 DO ii = jpindt(iseg), jpinft(iseg) 406 icount = icount + 1 407 nbidta(icount, igrd, ib_bdy) = ii 408 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 409 nbrdta(icount, igrd, ib_bdy) = ir 410 ENDDO 411 ENDDO 412 ! south 413 DO iseg = 1, nbdysegs 414 DO ii = jpisdt(iseg), jpisft(iseg) 415 icount = icount + 1 416 nbidta(icount, igrd, ib_bdy) = ii 417 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1 418 nbrdta(icount, igrd, ib_bdy) = ir 419 ENDDO 420 ENDDO 421 ENDDO 422 423 ! ------------ U points ------------- 424 igrd = 2 425 icount = 0 426 DO ir = 1, nn_rimwidth(ib_bdy) 427 ! east 428 DO iseg = 1, nbdysege 429 DO ij = jpjedt(iseg), jpjeft(iseg) 430 icount = icount + 1 431 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir 432 nbjdta(icount, igrd, ib_bdy) = ij 433 nbrdta(icount, igrd, ib_bdy) = ir 434 ENDDO 435 ENDDO 436 ! west 437 DO iseg = 1, nbdysegw 438 DO ij = jpjwdt(iseg), jpjwft(iseg) 439 icount = icount + 1 440 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 441 nbjdta(icount, igrd, ib_bdy) = ij 442 nbrdta(icount, igrd, ib_bdy) = ir 443 ENDDO 444 ENDDO 445 ! north 446 DO iseg = 1, nbdysegn 447 DO ii = jpindt(iseg), jpinft(iseg) - 1 448 icount = icount + 1 449 nbidta(icount, igrd, ib_bdy) = ii 450 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 451 nbrdta(icount, igrd, ib_bdy) = ir 452 ENDDO 453 ENDDO 454 ! south 455 DO iseg = 1, nbdysegs 456 DO ii = jpisdt(iseg), jpisft(iseg) - 1 457 icount = icount + 1 458 nbidta(icount, igrd, ib_bdy) = ii 459 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1 460 nbrdta(icount, igrd, ib_bdy) = ir 461 ENDDO 462 ENDDO 463 ENDDO 464 465 ! ------------ V points ------------- 466 igrd = 3 467 icount = 0 468 DO ir = 1, nn_rimwidth(ib_bdy) 469 ! east 470 DO iseg = 1, nbdysege 471 DO ij = jpjedt(iseg), jpjeft(iseg) - 1 472 icount = icount + 1 473 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 474 nbjdta(icount, igrd, ib_bdy) = ij 475 nbrdta(icount, igrd, ib_bdy) = ir 476 ENDDO 477 ENDDO 478 ! west 479 DO iseg = 1, nbdysegw 480 DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 481 icount = icount + 1 482 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 483 nbjdta(icount, igrd, ib_bdy) = ij 484 nbrdta(icount, igrd, ib_bdy) = ir 485 ENDDO 486 ENDDO 487 ! north 488 DO iseg = 1, nbdysegn 489 DO ii = jpindt(iseg), jpinft(iseg) 490 icount = icount + 1 491 nbidta(icount, igrd, ib_bdy) = ii 492 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir 493 nbrdta(icount, igrd, ib_bdy) = ir 494 ENDDO 495 ENDDO 496 ! south 497 DO iseg = 1, nbdysegs 498 DO ii = jpisdt(iseg), jpisft(iseg) 499 icount = icount + 1 500 nbidta(icount, igrd, ib_bdy) = ii 501 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1 502 nbrdta(icount, igrd, ib_bdy) = ir 503 ENDDO 504 ENDDO 505 ENDDO 506 507 ELSE ! Read global index arrays from boundary coordinates file. 508 509 DO igrd = 1, jpbgrd 510 CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 511 DO ii = 1,nblendta(igrd,ib_bdy) 512 nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 513 END DO 514 CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 515 DO ii = 1,nblendta(igrd,ib_bdy) 516 nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 517 END DO 518 CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 519 DO ii = 1,nblendta(igrd,ib_bdy) 520 nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 521 END DO 522 523 ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) ) 524 IF(lwp) WRITE(numout,*) 525 IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max 526 IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy) 527 IF (ibr_max < nn_rimwidth(ib_bdy)) & 528 CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 529 530 END DO 531 CALL iom_close( inum ) 532 533 ENDIF 534 535 ENDDO 536 537 ! Work out dimensions of boundary data on each processor 538 ! ------------------------------------------------------ 539 540 iw = mig(1) + 1 ! if monotasking and no zoom, iw=2 541 ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1 542 is = mjg(1) + 1 ! if monotasking and no zoom, is=2 543 in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1 544 545 DO ib_bdy = 1, nb_bdy 546 DO igrd = 1, jpbgrd 547 icount = 0 548 icountr = 0 549 idx_bdy(ib_bdy)%nblen(igrd) = 0 550 idx_bdy(ib_bdy)%nblenrim(igrd) = 0 551 DO ib = 1, nblendta(igrd,ib_bdy) 552 ! check that data is in correct order in file 553 ibm1 = MAX(1,ib-1) 554 IF(lwp) THEN ! Since all procs read global data only need to do this check on one proc... 555 IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN 556 CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', & 557 'A utility for re-ordering boundary coordinates and data files exists in the TOOLS/OBC directory') 558 ENDIF 559 ENDIF 560 ! check if point is in local domain 561 IF( nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND. & 562 & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in ) THEN 563 ! 564 icount = icount + 1 565 ! 566 IF( nbrdta(ib,igrd,ib_bdy) == 1 ) icountr = icountr+1 567 ENDIF 568 ENDDO 569 idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc 570 idx_bdy(ib_bdy)%nblen (igrd) = icount !: length of boundary data on each proc 571 ENDDO ! igrd 572 573 ! Allocate index arrays for this boundary set 574 !-------------------------------------------- 575 ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:)) 576 ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) ) 577 ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) ) 578 ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) ) 579 ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) ) 580 ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) ) 581 ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) ) 582 ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) ) 583 584 ! Dispatch mapping indices and discrete distances on each processor 585 ! ----------------------------------------------------------------- 586 587 DO igrd = 1, jpbgrd 588 icount = 0 589 ! Loop on rimwidth to ensure outermost points come first in the local arrays. 590 DO ir=1, nn_rimwidth(ib_bdy) 591 DO ib = 1, nblendta(igrd,ib_bdy) 592 ! check if point is in local domain and equals ir 593 IF( nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND. & 594 & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND. & 595 & nbrdta(ib,igrd,ib_bdy) == ir ) THEN 596 ! 597 icount = icount + 1 598 idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1 599 idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 600 idx_bdy(ib_bdy)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_bdy) 601 idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib 602 ENDIF 603 ENDDO 604 ENDDO 605 ENDDO 606 607 ! Compute rim weights for FRS scheme 608 ! ---------------------------------- 609 DO igrd = 1, jpbgrd 610 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 611 nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 612 idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ! tanh formulation 613 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2 ! quadratic 614 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth) ! linear 615 END DO 616 END DO 617 618 ENDDO 619 620 ! ------------------------------------------------------ 621 ! Initialise masks and find normal/tangential directions 622 ! ------------------------------------------------------ 151 623 152 624 ! Read global 2D mask at T-points: bdytmask 153 ! *****************************************625 ! ----------------------------------------- 154 626 ! bdytmask = 1 on the computational domain AND on open boundaries 155 627 ! = 0 elsewhere … … 158 630 zmask( : ,:) = 0.e0 159 631 zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0 160 ELSE IF( ln_mask ) THEN161 CALL iom_open( cn_mask , inum )632 ELSE IF( ln_mask_file ) THEN 633 CALL iom_open( cn_mask_file, inum ) 162 634 CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) ) 163 635 CALL iom_close( inum ) … … 184 656 185 657 186 ! Read discrete distance and mapping indices187 ! ******************************************188 nbidta(:,:) = 0.e0189 nbjdta(:,:) = 0.e0190 nbrdta(:,:) = 0.e0191 192 IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN193 icount = 0194 DO ir = 1, nn_rimwidth ! Define west boundary (from ii=2 to ii=1+nn_rimwidth):195 DO ij = 3, jpjglo-2196 icount = icount + 1197 nbidta(icount,:) = ir + 1 + (jpizoom-1)198 nbjdta(icount,:) = ij + (jpjzoom-1)199 nbrdta(icount,:) = ir200 END DO201 END DO202 !203 DO ir = 1, nn_rimwidth ! Define east boundary (from ii=jpiglo-1 to ii=jpiglo-nn_rimwidth):204 DO ij=3,jpjglo-2205 icount = icount + 1206 nbidta(icount,:) = jpiglo-ir + (jpizoom-1)207 nbidta(icount,2) = jpiglo-ir-1 + (jpizoom-1) ! special case for u points208 nbjdta(icount,:) = ij + (jpjzoom-1)209 nbrdta(icount,:) = ir210 END DO211 END DO212 !213 ELSE ! Read indices and distances in unstructured boundary data files214 !215 IF( ln_tides ) THEN ! Read tides input files for preference in case there are no bdydata files216 clfile(4) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_T.nc'217 clfile(5) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_U.nc'218 clfile(6) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_V.nc'219 ENDIF220 IF( ln_dyn_fla .AND. .NOT. ln_tides ) THEN221 clfile(4) = cn_dta_fla_T222 clfile(5) = cn_dta_fla_U223 clfile(6) = cn_dta_fla_V224 ENDIF225 226 IF( ln_tra_frs ) THEN227 clfile(1) = cn_dta_frs_T228 IF( .NOT. ln_dyn_frs ) THEN229 clfile(2) = cn_dta_frs_T ! Dummy read re read T file for sake of 6 files230 clfile(3) = cn_dta_frs_T !231 ENDIF232 ENDIF233 IF( ln_dyn_frs ) THEN234 IF( .NOT. ln_tra_frs ) clfile(1) = cn_dta_frs_U ! Dummy Read235 clfile(2) = cn_dta_frs_U236 clfile(3) = cn_dta_frs_V237 ENDIF238 239 ! ! how many files are we to read in?240 IF(ln_tides .OR. ln_dyn_fla) igrd_start = 4241 !242 IF(ln_tra_frs ) THEN ; igrd_start = 1243 ELSEIF(ln_dyn_frs) THEN ; igrd_start = 2244 ENDIF245 !246 IF( ln_tra_frs ) igrd_end = 1247 !248 IF(ln_dyn_fla .OR. ln_tides) THEN ; igrd_end = 6249 ELSEIF( ln_dyn_frs ) THEN ; igrd_end = 3250 ENDIF251 252 DO igrd = igrd_start, igrd_end253 CALL iom_open( clfile(igrd), inum )254 id_dummy = iom_varid( inum, 'nbidta', kdimsz=kdimsz )255 IF(lwp) WRITE(numout,*) 'kdimsz : ',kdimsz256 ib_len = kdimsz(1)257 IF( ib_len > jpbdta) CALL ctl_stop( 'Boundary data array in file too long.', &258 & 'File :', TRIM(clfile(igrd)),'increase parameter jpbdta.' )259 260 CALL iom_get( inum, jpdom_unknown, 'nbidta', zdta(1:ib_len,:) )261 DO ii = 1,ib_len262 nbidta(ii,igrd) = INT( zdta(ii,1) )263 END DO264 CALL iom_get( inum, jpdom_unknown, 'nbjdta', zdta(1:ib_len,:) )265 DO ii = 1,ib_len266 nbjdta(ii,igrd) = INT( zdta(ii,1) )267 END DO268 CALL iom_get( inum, jpdom_unknown, 'nbrdta', zdta(1:ib_len,:) )269 DO ii = 1,ib_len270 nbrdta(ii,igrd) = INT( zdta(ii,1) )271 END DO272 CALL iom_close( inum )273 274 IF( igrd < 4) THEN ! Check that rimwidth in file is big enough for Frs case(barotropic is one):275 ibr_max = MAXVAL( nbrdta(:,igrd) )276 IF(lwp) WRITE(numout,*)277 IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max278 IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth279 IF (ibr_max < nn_rimwidth) CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file' )280 ENDIF !Check igrd < 4281 !282 END DO283 !284 ENDIF285 286 ! Dispatch mapping indices and discrete distances on each processor287 ! *****************************************************************288 289 iw = mig(1) + 1 ! if monotasking and no zoom, iw=2290 ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1291 is = mjg(1) + 1 ! if monotasking and no zoom, is=2292 in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1293 294 DO igrd = igrd_start, igrd_end295 icount = 0296 icountr = 0297 nblen (igrd) = 0298 nblenrim(igrd) = 0299 nblendta(igrd) = 0300 DO ir=1, nn_rimwidth301 DO ib = 1, jpbdta302 ! check if point is in local domain and equals ir303 IF( nbidta(ib,igrd) >= iw .AND. nbidta(ib,igrd) <= ie .AND. &304 & nbjdta(ib,igrd) >= is .AND. nbjdta(ib,igrd) <= in .AND. &305 & nbrdta(ib,igrd) == ir ) THEN306 !307 icount = icount + 1308 !309 IF( ir == 1 ) icountr = icountr+1310 IF (icount > jpbdim) THEN311 IF(lwp) WRITE(numout,*) 'bdy_ini: jpbdim too small'312 nstop = nstop + 1313 ELSE314 nbi(icount, igrd) = nbidta(ib,igrd)- mig(1)+1315 nbj(icount, igrd) = nbjdta(ib,igrd)- mjg(1)+1316 nbr(icount, igrd) = nbrdta(ib,igrd)317 nbmap(icount,igrd) = ib318 ENDIF319 ENDIF320 END DO321 END DO322 nblenrim(igrd) = icountr !: length of rim boundary data on each proc323 nblen (igrd) = icount !: length of boundary data on each proc324 END DO325 326 ! Compute rim weights327 ! -------------------328 DO igrd = igrd_start, igrd_end329 DO ib = 1, nblen(igrd)330 nbw(ib,igrd) = 1.- TANH( FLOAT( nbr(ib,igrd) - 1 ) *0.5 ) ! tanh formulation331 ! nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr(ib,igrd))/FLOAT(nn_rimwidth))**2 ! quadratic332 ! nbw(ib,igrd) = FLOAT(nn_rimwidth+1-nbr(ib,igrd))/FLOAT(nn_rimwidth) ! linear333 END DO334 END DO335 336 658 ! Mask corrections 337 659 ! ---------------- … … 361 683 ! bdy masks and bmask are now set to zero on boundary points: 362 684 igrd = 1 ! In the free surface case, bmask is at T-points 363 DO ib = 1, nblenrim(igrd) 364 bmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 365 END DO 685 DO ib_bdy = 1, nb_bdy 686 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 687 bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0 688 ENDDO 689 ENDDO 366 690 ! 367 691 igrd = 1 368 DO ib = 1, nblenrim(igrd) 369 bdytmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 370 END DO 692 DO ib_bdy = 1, nb_bdy 693 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 694 bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0 695 ENDDO 696 ENDDO 371 697 ! 372 698 igrd = 2 373 DO ib = 1, nblenrim(igrd) 374 bdyumask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 375 END DO 699 DO ib_bdy = 1, nb_bdy 700 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 701 bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0 702 ENDDO 703 ENDDO 376 704 ! 377 705 igrd = 3 378 DO ib = 1, nblenrim(igrd) 379 bdyvmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0 380 END DO 706 DO ib_bdy = 1, nb_bdy 707 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 708 bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0 709 ENDDO 710 ENDDO 381 711 382 712 ! Lateral boundary conditions … … 384 714 CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) 385 715 386 IF( ln_vol .OR. ln_dyn_fla ) THEN ! Indices and directions of rim velocity components 387 ! 716 DO ib_bdy = 1, nb_bdy ! Indices and directions of rim velocity components 717 718 idx_bdy(ib_bdy)%flagu(:) = 0.e0 719 idx_bdy(ib_bdy)%flagv(:) = 0.e0 720 icount = 0 721 388 722 !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward 389 723 !flagu = 0 : u is tangential 390 724 !flagu = 1 : u is normal to the boundary and is direction is inward 391 icount = 0 392 flagu(:) = 0.e0 393 725 394 726 igrd = 2 ! u-component 395 DO ib = 1, nblenrim(igrd) 396 zefl=bdytmask(nbi(ib,igrd) , nbj(ib,igrd)) 397 zwfl=bdytmask(nbi(ib,igrd)+1, nbj(ib,igrd)) 398 IF( zefl + zwfl ==2 ) THEN 399 icount = icount +1 727 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 728 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 729 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 730 zefl = bdytmask(nbi ,nbj) 731 zwfl = bdytmask(nbi+1,nbj) 732 IF( zefl + zwfl == 2 ) THEN 733 icount = icount + 1 400 734 ELSE 401 flagu(ib)=-zefl+zwfl735 idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl 402 736 ENDIF 403 737 END DO … … 406 740 !flagv = 0 : u is tangential 407 741 !flagv = 1 : u is normal to the boundary and is direction is inward 408 flagv(:) = 0.e0409 742 410 743 igrd = 3 ! v-component 411 DO ib = 1, nblenrim(igrd) 412 znfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd)) 413 zsfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd)+1) 414 IF( znfl + zsfl ==2 ) THEN 744 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 745 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 746 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 747 znfl = bdytmask(nbi,nbj ) 748 zsfl = bdytmask(nbi,nbj+1) 749 IF( znfl + zsfl == 2 ) THEN 415 750 icount = icount + 1 416 751 ELSE 417 flagv(ib) = -znfl + zsfl752 idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl 418 753 END IF 419 754 END DO … … 422 757 IF(lwp) WRITE(numout,*) 423 758 IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,', & 424 ' are not boundary points. Check nbi, nbj, indices .'759 ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy 425 760 IF(lwp) WRITE(numout,*) ' ========== ' 426 761 IF(lwp) WRITE(numout,*) … … 428 763 ENDIF 429 764 430 END IF765 ENDDO 431 766 432 767 ! Compute total lateral surface for volume correction: … … 435 770 IF( ln_vol ) THEN 436 771 igrd = 2 ! Lateral surface at U-points 437 DO ib = 1, nblenrim(igrd) 438 bdysurftot = bdysurftot + hu (nbi(ib,igrd) ,nbj(ib,igrd)) & 439 & * e2u (nbi(ib,igrd) ,nbj(ib,igrd)) * ABS( flagu(ib) ) & 440 & * tmask_i(nbi(ib,igrd) ,nbj(ib,igrd)) & 441 & * tmask_i(nbi(ib,igrd)+1,nbj(ib,igrd)) 442 END DO 772 DO ib_bdy = 1, nb_bdy 773 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 774 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 775 nbj => idx_bdy(ib_bdy)%nbi(ib,igrd) 776 flagu => idx_bdy(ib_bdy)%flagu(ib) 777 bdysurftot = bdysurftot + hu (nbi , nbj) & 778 & * e2u (nbi , nbj) * ABS( flagu ) & 779 & * tmask_i(nbi , nbj) & 780 & * tmask_i(nbi+1, nbj) 781 ENDDO 782 ENDDO 443 783 444 784 igrd=3 ! Add lateral surface at V-points 445 DO ib = 1, nblenrim(igrd) 446 bdysurftot = bdysurftot + hv (nbi(ib,igrd),nbj(ib,igrd) ) & 447 & * e1v (nbi(ib,igrd),nbj(ib,igrd) ) * ABS( flagv(ib) ) & 448 & * tmask_i(nbi(ib,igrd),nbj(ib,igrd) ) & 449 & * tmask_i(nbi(ib,igrd),nbj(ib,igrd)+1) 450 END DO 785 DO ib_bdy = 1, nb_bdy 786 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 787 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 788 nbj => idx_bdy(ib_bdy)%nbi(ib,igrd) 789 flagv => idx_bdy(ib_bdy)%flagv(ib) 790 bdysurftot = bdysurftot + hv (nbi, nbj ) & 791 & * e1v (nbi, nbj ) * ABS( flagv ) & 792 & * tmask_i(nbi, nbj ) & 793 & * tmask_i(nbi, nbj+1) 794 ENDDO 795 ENDDO 451 796 ! 452 797 IF( lk_mpp ) CALL mpp_sum( bdysurftot ) ! sum over the global domain 453 798 END IF 454 455 ! Initialise bdy data arrays456 ! --------------------------457 tbdy(:,:) = 0.e0458 sbdy(:,:) = 0.e0459 ubdy(:,:) = 0.e0460 vbdy(:,:) = 0.e0461 sshbdy(:) = 0.e0462 ubtbdy(:) = 0.e0463 vbtbdy(:) = 0.e0464 #if defined key_lim2465 frld_bdy(:) = 0.e0466 hicif_bdy(:) = 0.e0467 hsnif_bdy(:) = 0.e0468 #endif469 470 ! Read in tidal constituents and adjust for model start time471 ! ----------------------------------------------------------472 IF( ln_tides ) CALL tide_data473 799 ! 800 ! Tidy up 801 !-------- 802 DEALLOCATE(nbidta, nbjdta, nbrdta) 803 804 IF( nn_timing == 1 ) CALL timing_stop('bdy_init') 805 474 806 END SUBROUTINE bdy_init 475 807 476 808 #else 477 809 !!--------------------------------------------------------------------------------- 478 !! Dummy module NO unstructuredopen boundaries810 !! Dummy module NO open boundaries 479 811 !!--------------------------------------------------------------------------------- 480 812 CONTAINS
Note: See TracChangeset
for help on using the changeset viewer.
