- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/BDY/bdyini.F90 (modified) (42 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/BDY/bdyini.F90
r10983 r13463 19 19 USE oce ! ocean dynamics and tracers variables 20 20 USE dom_oce ! ocean space and time domain 21 USE sbc_oce , ONLY: nn_ice 21 22 USE bdy_oce ! unstructured open boundary conditions 22 23 USE bdydta ! open boundary cond. setting (bdy_dta_init routine) 23 24 USE bdytides ! open boundary cond. setting (bdytide_init routine) 24 USE sbctide ! Tidal forcing or not25 USE phycst 25 USE tide_mod, ONLY: ln_tide ! tidal forcing 26 USE phycst , ONLY: rday 26 27 ! 27 28 USE in_out_manager ! I/O units … … 33 34 PRIVATE 34 35 35 PUBLIC bdy_init ! routine called in nemo_init 36 PUBLIC bdy_init ! routine called in nemo_init 37 PUBLIC find_neib ! routine called in bdy_nmn 36 38 37 39 INTEGER, PARAMETER :: jp_nseg = 100 ! 38 INTEGER, PARAMETER :: nrimmax = 20 ! maximum rimwidth in structured39 ! open boundary data files40 40 ! Straight open boundary segment parameters: 41 41 INTEGER :: nbdysege, nbdysegw, nbdysegn, nbdysegs … … 68 68 & ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 69 69 & cn_ice, nn_ice_dta, & 70 & rn_ice_tem, rn_ice_sal, rn_ice_age, & 71 & ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 70 & ln_vol, nn_volctl, nn_rimwidth 72 71 ! 73 72 INTEGER :: ios ! Local integer output status for namelist read … … 77 76 ! Read namelist parameters 78 77 ! ------------------------ 79 REWIND( numnam_ref ) ! Namelist nambdy in reference namelist :Unstructured open boundaries80 78 READ ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901) 81 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp ) 82 REWIND( numnam_cfg ) ! Namelist nambdy in configuration namelist :Unstructured open boundaries 79 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist' ) 80 ! make sur that all elements of the namelist variables have a default definition from namelist_ref 81 ln_coords_file (2:jp_bdy) = ln_coords_file (1) 82 cn_coords_file (2:jp_bdy) = cn_coords_file (1) 83 cn_dyn2d (2:jp_bdy) = cn_dyn2d (1) 84 nn_dyn2d_dta (2:jp_bdy) = nn_dyn2d_dta (1) 85 cn_dyn3d (2:jp_bdy) = cn_dyn3d (1) 86 nn_dyn3d_dta (2:jp_bdy) = nn_dyn3d_dta (1) 87 cn_tra (2:jp_bdy) = cn_tra (1) 88 nn_tra_dta (2:jp_bdy) = nn_tra_dta (1) 89 ln_tra_dmp (2:jp_bdy) = ln_tra_dmp (1) 90 ln_dyn3d_dmp (2:jp_bdy) = ln_dyn3d_dmp (1) 91 rn_time_dmp (2:jp_bdy) = rn_time_dmp (1) 92 rn_time_dmp_out(2:jp_bdy) = rn_time_dmp_out(1) 93 cn_ice (2:jp_bdy) = cn_ice (1) 94 nn_ice_dta (2:jp_bdy) = nn_ice_dta (1) 83 95 READ ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 ) 84 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist' , lwp)96 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist' ) 85 97 IF(lwm) WRITE ( numond, nambdy ) 86 98 87 99 IF( .NOT. Agrif_Root() ) ln_bdy = .FALSE. ! forced for Agrif children 100 101 IF( nb_bdy == 0 ) ln_bdy = .FALSE. 88 102 89 103 ! ----------------------------------------- … … 96 110 ! 97 111 ! Open boundaries definition (arrays and masks) 98 CALL bdy_segs 112 CALL bdy_def 113 IF( ln_meshmask ) CALL bdy_meshwri() 99 114 ! 100 115 ! Open boundaries initialisation of external data arrays … … 114 129 115 130 116 SUBROUTINE bdy_ segs131 SUBROUTINE bdy_def 117 132 !!---------------------------------------------------------------------- 118 133 !! *** ROUTINE bdy_init *** … … 125 140 !! ** Input : bdy_init.nc, input file for unstructured open boundaries 126 141 !!---------------------------------------------------------------------- 127 INTEGER :: ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices 128 INTEGER :: icount, icountr, ibr_max, ilen1, ibm1 ! local integers 142 INTEGER :: ib_bdy, ii, ij, igrd, ib, ir, iseg ! dummy loop indices 143 INTEGER :: icount, icountr, icountr0, ibr_max ! local integers 144 INTEGER :: ilen1 ! - - 129 145 INTEGER :: iwe, ies, iso, ino, inum, id_dummy ! - - 130 INTEGER :: igrd_start, igrd_end, jpbdta ! - - 131 INTEGER :: jpbdtau, jpbdtas ! - - 146 INTEGER :: jpbdta ! - - 132 147 INTEGER :: ib_bdy1, ib_bdy2, ib1, ib2 ! - - 133 INTEGER :: i_offset, j_offset ! - - 134 INTEGER , POINTER :: nbi, nbj, nbr ! short cuts 135 REAL(wp), POINTER, DIMENSION(:,:) :: pmask ! pointer to 2D mask fields 136 REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars 137 INTEGER, DIMENSION (2) :: kdimsz 138 INTEGER, DIMENSION(jpbgrd,jp_bdy) :: nblendta ! Length of index arrays 139 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta 140 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points 141 CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid 142 INTEGER :: com_east, com_west, com_south, com_north, jpk_max ! Flags for boundaries sending 143 INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b ! Flags for boundaries receiving 144 INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4) ! Arrays for neighbours coordinates 145 REAL(wp), TARGET, DIMENSION(jpi,jpj) :: zfmask ! temporary fmask array excluding coastal boundary condition (shlat) 146 !! 147 CHARACTER(LEN=1) :: ctypebdy ! - - 148 INTEGER :: nbdyind, nbdybeg, nbdyend 149 !! 150 NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 151 INTEGER :: ios ! Local integer output status for namelist read 148 INTEGER :: ii1, ii2, ii3, ij1, ij2, ij3 ! - - 149 INTEGER :: iibe, ijbe, iibi, ijbi ! - - 150 INTEGER :: flagu, flagv ! short cuts 151 INTEGER :: nbdyind, nbdybeg, nbdyend 152 INTEGER , DIMENSION(4) :: kdimsz 153 INTEGER , DIMENSION(jpbgrd,jp_bdy) :: nblendta ! Length of index arrays 154 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta 155 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points 156 CHARACTER(LEN=1) , DIMENSION(jpbgrd) :: cgrid 157 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zz_read ! work space for 2D global boundary data 158 REAL(wp), POINTER , DIMENSION(:,:) :: zmask ! pointer to 2D mask fields 159 REAL(wp) , DIMENSION(jpi,jpj) :: zfmask ! temporary fmask array excluding coastal boundary condition (shlat) 160 REAL(wp) , DIMENSION(jpi,jpj) :: ztmask, zumask, zvmask ! temporary u/v mask array 152 161 !!---------------------------------------------------------------------- 153 162 ! … … 160 169 & ' and general open boundary condition are not compatible' ) 161 170 162 IF( nb_bdy == 0 ) THEN 163 IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.' 164 ELSE 165 IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy 171 IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy 172 173 DO ib_bdy = 1,nb_bdy 174 175 IF(lwp) THEN 176 WRITE(numout,*) ' ' 177 WRITE(numout,*) '------ Open boundary data set ',ib_bdy,' ------' 178 IF( ln_coords_file(ib_bdy) ) THEN 179 WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy)) 180 ELSE 181 WRITE(numout,*) 'Boundary defined in namelist.' 182 ENDIF 183 WRITE(numout,*) 184 ENDIF 185 186 ! barotropic bdy 187 !---------------- 188 IF(lwp) THEN 189 WRITE(numout,*) 'Boundary conditions for barotropic solution: ' 190 SELECT CASE( cn_dyn2d(ib_bdy) ) 191 CASE( 'none' ) ; WRITE(numout,*) ' no open boundary condition' 192 CASE( 'frs' ) ; WRITE(numout,*) ' Flow Relaxation Scheme' 193 CASE( 'flather' ) ; WRITE(numout,*) ' Flather radiation condition' 194 CASE( 'orlanski' ) ; WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging' 195 CASE( 'orlanski_npo' ) ; WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging' 196 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn2d' ) 197 END SELECT 198 ENDIF 199 200 dta_bdy(ib_bdy)%lneed_ssh = cn_dyn2d(ib_bdy) == 'flather' 201 dta_bdy(ib_bdy)%lneed_dyn2d = cn_dyn2d(ib_bdy) /= 'none' 202 203 IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn2d ) THEN 204 SELECT CASE( nn_dyn2d_dta(ib_bdy) ) ! 205 CASE( 0 ) ; WRITE(numout,*) ' initial state used for bdy data' 206 CASE( 1 ) ; WRITE(numout,*) ' boundary data taken from file' 207 CASE( 2 ) ; WRITE(numout,*) ' tidal harmonic forcing taken from file' 208 CASE( 3 ) ; WRITE(numout,*) ' boundary data AND tidal harmonic forcing taken from files' 209 CASE DEFAULT ; CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 210 END SELECT 211 ENDIF 212 IF ( dta_bdy(ib_bdy)%lneed_dyn2d .AND. nn_dyn2d_dta(ib_bdy) .GE. 2 .AND. .NOT.ln_tide ) THEN 213 CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' ) 214 ENDIF 215 IF(lwp) WRITE(numout,*) 216 217 ! baroclinic bdy 218 !---------------- 219 IF(lwp) THEN 220 WRITE(numout,*) 'Boundary conditions for baroclinic velocities: ' 221 SELECT CASE( cn_dyn3d(ib_bdy) ) 222 CASE('none') ; WRITE(numout,*) ' no open boundary condition' 223 CASE('frs') ; WRITE(numout,*) ' Flow Relaxation Scheme' 224 CASE('specified') ; WRITE(numout,*) ' Specified value' 225 CASE('neumann') ; WRITE(numout,*) ' Neumann conditions' 226 CASE('zerograd') ; WRITE(numout,*) ' Zero gradient for baroclinic velocities' 227 CASE('zero') ; WRITE(numout,*) ' Zero baroclinic velocities (runoff case)' 228 CASE('orlanski') ; WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging' 229 CASE('orlanski_npo') ; WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging' 230 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn3d' ) 231 END SELECT 232 ENDIF 233 234 dta_bdy(ib_bdy)%lneed_dyn3d = cn_dyn3d(ib_bdy) == 'frs' .OR. cn_dyn3d(ib_bdy) == 'specified' & 235 & .OR. cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' 236 237 IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn3d ) THEN 238 SELECT CASE( nn_dyn3d_dta(ib_bdy) ) ! 239 CASE( 0 ) ; WRITE(numout,*) ' initial state used for bdy data' 240 CASE( 1 ) ; WRITE(numout,*) ' boundary data taken from file' 241 CASE DEFAULT ; CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' ) 242 END SELECT 243 END IF 244 245 IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 246 IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN 247 IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.' 248 ln_dyn3d_dmp(ib_bdy) = .false. 249 ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN 250 CALL ctl_stop( 'Use FRS OR relaxation' ) 251 ELSE 252 IF(lwp) WRITE(numout,*) ' + baroclinic velocities relaxation zone' 253 IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days' 254 IF(rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 255 dta_bdy(ib_bdy)%lneed_dyn3d = .TRUE. 256 ENDIF 257 ELSE 258 IF(lwp) WRITE(numout,*) ' NO relaxation on baroclinic velocities' 259 ENDIF 260 IF(lwp) WRITE(numout,*) 261 262 ! tra bdy 263 !---------------- 264 IF(lwp) THEN 265 WRITE(numout,*) 'Boundary conditions for temperature and salinity: ' 266 SELECT CASE( cn_tra(ib_bdy) ) 267 CASE('none') ; WRITE(numout,*) ' no open boundary condition' 268 CASE('frs') ; WRITE(numout,*) ' Flow Relaxation Scheme' 269 CASE('specified') ; WRITE(numout,*) ' Specified value' 270 CASE('neumann') ; WRITE(numout,*) ' Neumann conditions' 271 CASE('runoff') ; WRITE(numout,*) ' Runoff conditions : Neumann for T and specified to 0.1 for salinity' 272 CASE('orlanski') ; WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging' 273 CASE('orlanski_npo') ; WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging' 274 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_tra' ) 275 END SELECT 276 ENDIF 277 278 dta_bdy(ib_bdy)%lneed_tra = cn_tra(ib_bdy) == 'frs' .OR. cn_tra(ib_bdy) == 'specified' & 279 & .OR. cn_tra(ib_bdy) == 'orlanski' .OR. cn_tra(ib_bdy) == 'orlanski_npo' 280 281 IF( lwp .AND. dta_bdy(ib_bdy)%lneed_tra ) THEN 282 SELECT CASE( nn_tra_dta(ib_bdy) ) ! 283 CASE( 0 ) ; WRITE(numout,*) ' initial state used for bdy data' 284 CASE( 1 ) ; WRITE(numout,*) ' boundary data taken from file' 285 CASE DEFAULT ; CALL ctl_stop( 'nn_tra_dta must be 0 or 1' ) 286 END SELECT 287 ENDIF 288 289 IF ( ln_tra_dmp(ib_bdy) ) THEN 290 IF ( cn_tra(ib_bdy) == 'none' ) THEN 291 IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 292 ln_tra_dmp(ib_bdy) = .false. 293 ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN 294 CALL ctl_stop( 'Use FRS OR relaxation' ) 295 ELSE 296 IF(lwp) WRITE(numout,*) ' + T/S relaxation zone' 297 IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days' 298 IF(lwp) WRITE(numout,*) ' Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days' 299 IF(lwp.AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 300 dta_bdy(ib_bdy)%lneed_tra = .TRUE. 301 ENDIF 302 ELSE 303 IF(lwp) WRITE(numout,*) ' NO T/S relaxation' 304 ENDIF 305 IF(lwp) WRITE(numout,*) 306 307 #if defined key_si3 308 IF(lwp) THEN 309 WRITE(numout,*) 'Boundary conditions for sea ice: ' 310 SELECT CASE( cn_ice(ib_bdy) ) 311 CASE('none') ; WRITE(numout,*) ' no open boundary condition' 312 CASE('frs') ; WRITE(numout,*) ' Flow Relaxation Scheme' 313 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_ice' ) 314 END SELECT 315 ENDIF 316 317 dta_bdy(ib_bdy)%lneed_ice = cn_ice(ib_bdy) /= 'none' 318 319 IF( dta_bdy(ib_bdy)%lneed_ice .AND. nn_ice /= 2 ) THEN 320 WRITE(ctmp1,*) 'bdy number ', ib_bdy,', needs ice model but nn_ice = ', nn_ice 321 CALL ctl_stop( ctmp1 ) 322 ENDIF 323 324 IF( lwp .AND. dta_bdy(ib_bdy)%lneed_ice ) THEN 325 SELECT CASE( nn_ice_dta(ib_bdy) ) ! 326 CASE( 0 ) ; WRITE(numout,*) ' initial state used for bdy data' 327 CASE( 1 ) ; WRITE(numout,*) ' boundary data taken from file' 328 CASE DEFAULT ; CALL ctl_stop( 'nn_ice_dta must be 0 or 1' ) 329 END SELECT 330 ENDIF 331 #else 332 dta_bdy(ib_bdy)%lneed_ice = .FALSE. 333 #endif 334 ! 335 IF(lwp) WRITE(numout,*) 336 IF(lwp) WRITE(numout,*) ' Width of relaxation zone = ', nn_rimwidth(ib_bdy) 337 IF(lwp) WRITE(numout,*) 338 ! 339 END DO ! nb_bdy 340 341 IF( lwp ) THEN 342 IF( ln_vol ) THEN ! check volume conservation (nn_volctl value) 343 WRITE(numout,*) 'Volume correction applied at open boundaries' 344 WRITE(numout,*) 345 SELECT CASE ( nn_volctl ) 346 CASE( 1 ) ; WRITE(numout,*) ' The total volume will be constant' 347 CASE( 0 ) ; WRITE(numout,*) ' The total volume will vary according to the surface E-P flux' 348 CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 349 END SELECT 350 WRITE(numout,*) 351 ! 352 ! sanity check if used with tides 353 IF( ln_tide ) THEN 354 WRITE(numout,*) ' The total volume correction is not working with tides. ' 355 WRITE(numout,*) ' Set ln_vol to .FALSE. ' 356 WRITE(numout,*) ' or ' 357 WRITE(numout,*) ' equilibriate your bdy input files ' 358 CALL ctl_stop( 'The total volume correction is not working with tides.' ) 359 END IF 360 ELSE 361 WRITE(numout,*) 'No volume correction applied at open boundaries' 362 WRITE(numout,*) 363 ENDIF 166 364 ENDIF 167 168 DO ib_bdy = 1,nb_bdy169 IF(lwp) WRITE(numout,*) ' '170 IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------'171 172 IF( ln_coords_file(ib_bdy) ) THEN173 IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))174 ELSE175 IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.'176 ENDIF177 IF(lwp) WRITE(numout,*)178 179 IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution: '180 SELECT CASE( cn_dyn2d(ib_bdy) )181 CASE( 'none' )182 IF(lwp) WRITE(numout,*) ' no open boundary condition'183 dta_bdy(ib_bdy)%ll_ssh = .false.184 dta_bdy(ib_bdy)%ll_u2d = .false.185 dta_bdy(ib_bdy)%ll_v2d = .false.186 CASE( 'frs' )187 IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'188 dta_bdy(ib_bdy)%ll_ssh = .false.189 dta_bdy(ib_bdy)%ll_u2d = .true.190 dta_bdy(ib_bdy)%ll_v2d = .true.191 CASE( 'flather' )192 IF(lwp) WRITE(numout,*) ' Flather radiation condition'193 dta_bdy(ib_bdy)%ll_ssh = .true.194 dta_bdy(ib_bdy)%ll_u2d = .true.195 dta_bdy(ib_bdy)%ll_v2d = .true.196 CASE( 'orlanski' )197 IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'198 dta_bdy(ib_bdy)%ll_ssh = .false.199 dta_bdy(ib_bdy)%ll_u2d = .true.200 dta_bdy(ib_bdy)%ll_v2d = .true.201 CASE( 'orlanski_npo' )202 IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'203 dta_bdy(ib_bdy)%ll_ssh = .false.204 dta_bdy(ib_bdy)%ll_u2d = .true.205 dta_bdy(ib_bdy)%ll_v2d = .true.206 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn2d' )207 END SELECT208 IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN209 SELECT CASE( nn_dyn2d_dta(ib_bdy) ) !210 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'211 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'212 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' tidal harmonic forcing taken from file'213 CASE( 3 ) ; IF(lwp) WRITE(numout,*) ' boundary data AND tidal harmonic forcing taken from files'214 CASE DEFAULT ; CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )215 END SELECT216 IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.ln_tide)) THEN217 CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' )218 ENDIF219 ENDIF220 IF(lwp) WRITE(numout,*)221 222 IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities: '223 SELECT CASE( cn_dyn3d(ib_bdy) )224 CASE('none')225 IF(lwp) WRITE(numout,*) ' no open boundary condition'226 dta_bdy(ib_bdy)%ll_u3d = .false.227 dta_bdy(ib_bdy)%ll_v3d = .false.228 CASE('frs')229 IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'230 dta_bdy(ib_bdy)%ll_u3d = .true.231 dta_bdy(ib_bdy)%ll_v3d = .true.232 CASE('specified')233 IF(lwp) WRITE(numout,*) ' Specified value'234 dta_bdy(ib_bdy)%ll_u3d = .true.235 dta_bdy(ib_bdy)%ll_v3d = .true.236 CASE('neumann')237 IF(lwp) WRITE(numout,*) ' Neumann conditions'238 dta_bdy(ib_bdy)%ll_u3d = .false.239 dta_bdy(ib_bdy)%ll_v3d = .false.240 CASE('zerograd')241 IF(lwp) WRITE(numout,*) ' Zero gradient for baroclinic velocities'242 dta_bdy(ib_bdy)%ll_u3d = .false.243 dta_bdy(ib_bdy)%ll_v3d = .false.244 CASE('zero')245 IF(lwp) WRITE(numout,*) ' Zero baroclinic velocities (runoff case)'246 dta_bdy(ib_bdy)%ll_u3d = .false.247 dta_bdy(ib_bdy)%ll_v3d = .false.248 CASE('orlanski')249 IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'250 dta_bdy(ib_bdy)%ll_u3d = .true.251 dta_bdy(ib_bdy)%ll_v3d = .true.252 CASE('orlanski_npo')253 IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'254 dta_bdy(ib_bdy)%ll_u3d = .true.255 dta_bdy(ib_bdy)%ll_v3d = .true.256 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn3d' )257 END SELECT258 IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN259 SELECT CASE( nn_dyn3d_dta(ib_bdy) ) !260 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'261 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'262 CASE DEFAULT ; CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )263 END SELECT264 ENDIF265 266 IF ( ln_dyn3d_dmp(ib_bdy) ) THEN267 IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN268 IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.'269 ln_dyn3d_dmp(ib_bdy)=.false.270 ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN271 CALL ctl_stop( 'Use FRS OR relaxation' )272 ELSE273 IF(lwp) WRITE(numout,*) ' + baroclinic velocities relaxation zone'274 IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days'275 IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )276 dta_bdy(ib_bdy)%ll_u3d = .true.277 dta_bdy(ib_bdy)%ll_v3d = .true.278 ENDIF279 ELSE280 IF(lwp) WRITE(numout,*) ' NO relaxation on baroclinic velocities'281 ENDIF282 IF(lwp) WRITE(numout,*)283 284 IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity: '285 SELECT CASE( cn_tra(ib_bdy) )286 CASE('none')287 IF(lwp) WRITE(numout,*) ' no open boundary condition'288 dta_bdy(ib_bdy)%ll_tem = .false.289 dta_bdy(ib_bdy)%ll_sal = .false.290 CASE('frs')291 IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'292 dta_bdy(ib_bdy)%ll_tem = .true.293 dta_bdy(ib_bdy)%ll_sal = .true.294 CASE('specified')295 IF(lwp) WRITE(numout,*) ' Specified value'296 dta_bdy(ib_bdy)%ll_tem = .true.297 dta_bdy(ib_bdy)%ll_sal = .true.298 CASE('neumann')299 IF(lwp) WRITE(numout,*) ' Neumann conditions'300 dta_bdy(ib_bdy)%ll_tem = .false.301 dta_bdy(ib_bdy)%ll_sal = .false.302 CASE('runoff')303 IF(lwp) WRITE(numout,*) ' Runoff conditions : Neumann for T and specified to 0.1 for salinity'304 dta_bdy(ib_bdy)%ll_tem = .false.305 dta_bdy(ib_bdy)%ll_sal = .false.306 CASE('orlanski')307 IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'308 dta_bdy(ib_bdy)%ll_tem = .true.309 dta_bdy(ib_bdy)%ll_sal = .true.310 CASE('orlanski_npo')311 IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'312 dta_bdy(ib_bdy)%ll_tem = .true.313 dta_bdy(ib_bdy)%ll_sal = .true.314 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_tra' )315 END SELECT316 IF( cn_tra(ib_bdy) /= 'none' ) THEN317 SELECT CASE( nn_tra_dta(ib_bdy) ) !318 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'319 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'320 CASE DEFAULT ; CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )321 END SELECT322 ENDIF323 324 IF ( ln_tra_dmp(ib_bdy) ) THEN325 IF ( cn_tra(ib_bdy) == 'none' ) THEN326 IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'327 ln_tra_dmp(ib_bdy)=.false.328 ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN329 CALL ctl_stop( 'Use FRS OR relaxation' )330 ELSE331 IF(lwp) WRITE(numout,*) ' + T/S relaxation zone'332 IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days'333 IF(lwp) WRITE(numout,*) ' Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days'334 IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )335 dta_bdy(ib_bdy)%ll_tem = .true.336 dta_bdy(ib_bdy)%ll_sal = .true.337 ENDIF338 ELSE339 IF(lwp) WRITE(numout,*) ' NO T/S relaxation'340 ENDIF341 IF(lwp) WRITE(numout,*)342 343 #if defined key_si3344 IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice: '345 SELECT CASE( cn_ice(ib_bdy) )346 CASE('none')347 IF(lwp) WRITE(numout,*) ' no open boundary condition'348 dta_bdy(ib_bdy)%ll_a_i = .false.349 dta_bdy(ib_bdy)%ll_h_i = .false.350 dta_bdy(ib_bdy)%ll_h_s = .false.351 CASE('frs')352 IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'353 dta_bdy(ib_bdy)%ll_a_i = .true.354 dta_bdy(ib_bdy)%ll_h_i = .true.355 dta_bdy(ib_bdy)%ll_h_s = .true.356 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_ice' )357 END SELECT358 IF( cn_ice(ib_bdy) /= 'none' ) THEN359 SELECT CASE( nn_ice_dta(ib_bdy) ) !360 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data'361 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file'362 CASE DEFAULT ; CALL ctl_stop( 'nn_ice_dta must be 0 or 1' )363 END SELECT364 ENDIF365 IF(lwp) WRITE(numout,*)366 IF(lwp) WRITE(numout,*) ' tem of bdy sea-ice = ', rn_ice_tem(ib_bdy)367 IF(lwp) WRITE(numout,*) ' sal of bdy sea-ice = ', rn_ice_sal(ib_bdy)368 IF(lwp) WRITE(numout,*) ' age of bdy sea-ice = ', rn_ice_age(ib_bdy)369 #endif370 371 IF(lwp) WRITE(numout,*) ' Width of relaxation zone = ', nn_rimwidth(ib_bdy)372 IF(lwp) WRITE(numout,*)373 !374 END DO375 376 IF( nb_bdy > 0 ) THEN377 IF( ln_vol ) THEN ! check volume conservation (nn_volctl value)378 IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'379 IF(lwp) WRITE(numout,*)380 SELECT CASE ( nn_volctl )381 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' The total volume will be constant'382 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' The total volume will vary according to the surface E-P flux'383 CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' )384 END SELECT385 IF(lwp) WRITE(numout,*)386 !387 ! sanity check if used with tides388 IF( ln_tide ) THEN389 IF(lwp) WRITE(numout,*) ' The total volume correction is not working with tides. '390 IF(lwp) WRITE(numout,*) ' Set ln_vol to .FALSE. '391 IF(lwp) WRITE(numout,*) ' or '392 IF(lwp) WRITE(numout,*) ' equilibriate your bdy input files '393 CALL ctl_stop( 'The total volume correction is not working with tides.' )394 END IF395 ELSE396 IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'397 IF(lwp) WRITE(numout,*)398 ENDIF399 IF( nb_jpk_bdy(ib_bdy) > 0 ) THEN400 IF(lwp) WRITE(numout,*) '*** open boundary will be interpolate in the vertical onto the native grid ***'401 ELSE402 IF(lwp) WRITE(numout,*) '*** open boundary will be read straight onto the native grid without vertical interpolation ***'403 ENDIF404 ENDIF405 365 406 366 ! ------------------------------------------------- 407 367 ! Initialise indices arrays for open boundaries 408 368 ! ------------------------------------------------- 409 410 ! Work out global dimensions of boundary data411 ! ---------------------------------------------412 REWIND( numnam_cfg )413 369 414 370 nblendta(:,:) = 0 … … 417 373 nbdysegn = 0 418 374 nbdysegs = 0 419 icount = 0 ! count user defined segments 420 ! Dimensions below are used to allocate arrays to read external data 421 jpbdtas = 1 ! Maximum size of boundary data (structured case) 422 jpbdtau = 1 ! Maximum size of boundary data (unstructured case) 423 375 376 ! Define all boundaries 377 ! --------------------- 424 378 DO ib_bdy = 1, nb_bdy 425 426 IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 427 428 icount = icount + 1 429 ! No REWIND here because may need to read more than one nambdy_index namelist. 430 ! Read only namelist_cfg to avoid unseccessfull overwrite 431 ! keep full control of the configuration namelist 432 READ ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 ) 433 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_index in configuration namelist', lwp ) 434 IF(lwm) WRITE ( numond, nambdy_index ) 435 436 SELECT CASE ( TRIM(ctypebdy) ) 437 CASE( 'N' ) 438 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 439 nbdyind = jpjglo - 2 ! set boundary to whole side of model domain. 440 nbdybeg = 2 441 nbdyend = jpiglo - 1 442 ENDIF 443 nbdysegn = nbdysegn + 1 444 npckgn(nbdysegn) = ib_bdy ! Save bdy package number 445 jpjnob(nbdysegn) = nbdyind 446 jpindt(nbdysegn) = nbdybeg 447 jpinft(nbdysegn) = nbdyend 448 ! 449 CASE( 'S' ) 450 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 451 nbdyind = 2 ! set boundary to whole side of model domain. 452 nbdybeg = 2 453 nbdyend = jpiglo - 1 454 ENDIF 455 nbdysegs = nbdysegs + 1 456 npckgs(nbdysegs) = ib_bdy ! Save bdy package number 457 jpjsob(nbdysegs) = nbdyind 458 jpisdt(nbdysegs) = nbdybeg 459 jpisft(nbdysegs) = nbdyend 460 ! 461 CASE( 'E' ) 462 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 463 nbdyind = jpiglo - 2 ! set boundary to whole side of model domain. 464 nbdybeg = 2 465 nbdyend = jpjglo - 1 466 ENDIF 467 nbdysege = nbdysege + 1 468 npckge(nbdysege) = ib_bdy ! Save bdy package number 469 jpieob(nbdysege) = nbdyind 470 jpjedt(nbdysege) = nbdybeg 471 jpjeft(nbdysege) = nbdyend 472 ! 473 CASE( 'W' ) 474 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 475 nbdyind = 2 ! set boundary to whole side of model domain. 476 nbdybeg = 2 477 nbdyend = jpjglo - 1 478 ENDIF 479 nbdysegw = nbdysegw + 1 480 npckgw(nbdysegw) = ib_bdy ! Save bdy package number 481 jpiwob(nbdysegw) = nbdyind 482 jpjwdt(nbdysegw) = nbdybeg 483 jpjwft(nbdysegw) = nbdyend 484 ! 485 CASE DEFAULT ; CALL ctl_stop( 'ctypebdy must be N, S, E or W' ) 486 END SELECT 487 488 ! For simplicity we assume that in case of straight bdy, arrays have the same length 489 ! (even if it is true that last tangential velocity points 490 ! are useless). This simplifies a little bit boundary data format (and agrees with format 491 ! used so far in obc package) 492 493 nblendta(1:jpbgrd,ib_bdy) = (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy) 494 jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1)) 495 IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) & 496 & CALL ctl_stop( 'rimwidth must be lower than nrimmax' ) 497 498 ELSE ! Read size of arrays in boundary coordinates file. 379 ! 380 IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! build bdy coordinates with segments defined in namelist 381 382 CALL bdy_read_seg( ib_bdy, nblendta(:,ib_bdy) ) 383 384 ELSE ! Read size of arrays in boundary coordinates file. 385 499 386 CALL iom_open( cn_coords_file(ib_bdy), inum ) 500 387 DO igrd = 1, jpbgrd 501 388 id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 502 389 nblendta(igrd,ib_bdy) = MAXVAL(kdimsz) 503 jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz))504 390 END DO 505 391 CALL iom_close( inum ) 506 ! 507 ENDIF 392 ENDIF 508 393 ! 509 394 END DO ! ib_bdy 510 395 511 IF (nb_bdy>0) THEN512 jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))513 514 ! Allocate arrays515 !---------------516 ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), &517 & nbrdta(jpbdta, jpbgrd, nb_bdy) )518 519 jpk_max = MAXVAL(nb_jpk_bdy)520 jpk_max = MAX(jpk_max, jpk)521 522 ALLOCATE( dta_global(jpbdtau, 1, jpk_max) )523 ALLOCATE( dta_global_z(jpbdtau, 1, jpk_max) ) ! needed ?? TODO524 ALLOCATE( dta_global_dz(jpbdtau, 1, jpk_max) )! needed ?? TODO525 526 IF ( icount>0 ) THEN527 ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk_max) )528 ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk_max) ) ! needed ?? TODO529 ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk_max) )! needed ?? TODO530 ENDIF531 !532 ENDIF533 534 396 ! Now look for crossings in user (namelist) defined open boundary segments: 535 !-------------------------------------------------------------------------- 536 IF( icount>0 ) CALL bdy_ctl_seg 537 397 IF( nbdysege > 0 .OR. nbdysegw > 0 .OR. nbdysegn > 0 .OR. nbdysegs > 0) CALL bdy_ctl_seg 398 399 ! Allocate arrays 400 !--------------- 401 jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 402 ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), nbrdta(jpbdta, jpbgrd, nb_bdy) ) 403 nbrdta(:,:,:) = 0 ! initialize nbrdta as it may not be completely defined for each bdy 404 538 405 ! Calculate global boundary index arrays or read in from file 539 406 !------------------------------------------------------------ … … 543 410 IF( ln_coords_file(ib_bdy) ) THEN 544 411 ! 412 ALLOCATE( zz_read( MAXVAL(nblendta), 1 ) ) 545 413 CALL iom_open( cn_coords_file(ib_bdy), inum ) 414 ! 546 415 DO igrd = 1, jpbgrd 547 CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )416 CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) ) 548 417 DO ii = 1,nblendta(igrd,ib_bdy) 549 nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )418 nbidta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) + nn_hls 550 419 END DO 551 CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )420 CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) ) 552 421 DO ii = 1,nblendta(igrd,ib_bdy) 553 nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )422 nbjdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) + nn_hls 554 423 END DO 555 CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )424 CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) ) 556 425 DO ii = 1,nblendta(igrd,ib_bdy) 557 nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )426 nbrdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) 558 427 END DO 559 428 ! … … 563 432 IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy) 564 433 IF (ibr_max < nn_rimwidth(ib_bdy)) & 565 CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 566 END DO 434 CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 435 END DO 436 ! 567 437 CALL iom_close( inum ) 438 DEALLOCATE( zz_read ) 568 439 ! 569 ENDIF 570 ! 571 END DO 572 440 ENDIF 441 ! 442 END DO 443 573 444 ! 2. Now fill indices corresponding to straight open boundary arrays: 574 ! East 575 !----- 576 DO iseg = 1, nbdysege 577 ib_bdy = npckge(iseg) 578 ! 579 ! ------------ T points ------------- 580 igrd=1 581 icount=0 582 DO ir = 1, nn_rimwidth(ib_bdy) 583 DO ij = jpjedt(iseg), jpjeft(iseg) 584 icount = icount + 1 585 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 586 nbjdta(icount, igrd, ib_bdy) = ij 587 nbrdta(icount, igrd, ib_bdy) = ir 588 ENDDO 589 ENDDO 590 ! 591 ! ------------ U points ------------- 592 igrd=2 593 icount=0 594 DO ir = 1, nn_rimwidth(ib_bdy) 595 DO ij = jpjedt(iseg), jpjeft(iseg) 596 icount = icount + 1 597 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir 598 nbjdta(icount, igrd, ib_bdy) = ij 599 nbrdta(icount, igrd, ib_bdy) = ir 600 ENDDO 601 ENDDO 602 ! 603 ! ------------ V points ------------- 604 igrd=3 605 icount=0 606 DO ir = 1, nn_rimwidth(ib_bdy) 607 ! DO ij = jpjedt(iseg), jpjeft(iseg) - 1 608 DO ij = jpjedt(iseg), jpjeft(iseg) 609 icount = icount + 1 610 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 611 nbjdta(icount, igrd, ib_bdy) = ij 612 nbrdta(icount, igrd, ib_bdy) = ir 613 ENDDO 614 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 615 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 616 ENDDO 617 ENDDO 618 ! 619 ! West 620 !----- 621 DO iseg = 1, nbdysegw 622 ib_bdy = npckgw(iseg) 623 ! 624 ! ------------ T points ------------- 625 igrd=1 626 icount=0 627 DO ir = 1, nn_rimwidth(ib_bdy) 628 DO ij = jpjwdt(iseg), jpjwft(iseg) 629 icount = icount + 1 630 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 631 nbjdta(icount, igrd, ib_bdy) = ij 632 nbrdta(icount, igrd, ib_bdy) = ir 633 ENDDO 634 ENDDO 635 ! 636 ! ------------ U points ------------- 637 igrd=2 638 icount=0 639 DO ir = 1, nn_rimwidth(ib_bdy) 640 DO ij = jpjwdt(iseg), jpjwft(iseg) 641 icount = icount + 1 642 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 643 nbjdta(icount, igrd, ib_bdy) = ij 644 nbrdta(icount, igrd, ib_bdy) = ir 645 ENDDO 646 ENDDO 647 ! 648 ! ------------ V points ------------- 649 igrd=3 650 icount=0 651 DO ir = 1, nn_rimwidth(ib_bdy) 652 ! DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 653 DO ij = jpjwdt(iseg), jpjwft(iseg) 654 icount = icount + 1 655 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 656 nbjdta(icount, igrd, ib_bdy) = ij 657 nbrdta(icount, igrd, ib_bdy) = ir 658 ENDDO 659 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 660 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 661 ENDDO 662 ENDDO 663 ! 664 ! North 665 !----- 666 DO iseg = 1, nbdysegn 667 ib_bdy = npckgn(iseg) 668 ! 669 ! ------------ T points ------------- 670 igrd=1 671 icount=0 672 DO ir = 1, nn_rimwidth(ib_bdy) 673 DO ii = jpindt(iseg), jpinft(iseg) 674 icount = icount + 1 675 nbidta(icount, igrd, ib_bdy) = ii 676 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 677 nbrdta(icount, igrd, ib_bdy) = ir 678 ENDDO 679 ENDDO 680 ! 681 ! ------------ U points ------------- 682 igrd=2 683 icount=0 684 DO ir = 1, nn_rimwidth(ib_bdy) 685 ! DO ii = jpindt(iseg), jpinft(iseg) - 1 686 DO ii = jpindt(iseg), jpinft(iseg) 687 icount = icount + 1 688 nbidta(icount, igrd, ib_bdy) = ii 689 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 690 nbrdta(icount, igrd, ib_bdy) = ir 691 ENDDO 692 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 693 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 694 ENDDO 695 ! 696 ! ------------ V points ------------- 697 igrd=3 698 icount=0 699 DO ir = 1, nn_rimwidth(ib_bdy) 700 DO ii = jpindt(iseg), jpinft(iseg) 701 icount = icount + 1 702 nbidta(icount, igrd, ib_bdy) = ii 703 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir 704 nbrdta(icount, igrd, ib_bdy) = ir 705 ENDDO 706 ENDDO 707 ENDDO 708 ! 709 ! South 710 !----- 711 DO iseg = 1, nbdysegs 712 ib_bdy = npckgs(iseg) 713 ! 714 ! ------------ T points ------------- 715 igrd=1 716 icount=0 717 DO ir = 1, nn_rimwidth(ib_bdy) 718 DO ii = jpisdt(iseg), jpisft(iseg) 719 icount = icount + 1 720 nbidta(icount, igrd, ib_bdy) = ii 721 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 722 nbrdta(icount, igrd, ib_bdy) = ir 723 ENDDO 724 ENDDO 725 ! 726 ! ------------ U points ------------- 727 igrd=2 728 icount=0 729 DO ir = 1, nn_rimwidth(ib_bdy) 730 ! DO ii = jpisdt(iseg), jpisft(iseg) - 1 731 DO ii = jpisdt(iseg), jpisft(iseg) 732 icount = icount + 1 733 nbidta(icount, igrd, ib_bdy) = ii 734 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 735 nbrdta(icount, igrd, ib_bdy) = ir 736 ENDDO 737 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 738 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 739 ENDDO 740 ! 741 ! ------------ V points ------------- 742 igrd=3 743 icount=0 744 DO ir = 1, nn_rimwidth(ib_bdy) 745 DO ii = jpisdt(iseg), jpisft(iseg) 746 icount = icount + 1 747 nbidta(icount, igrd, ib_bdy) = ii 748 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 749 nbrdta(icount, igrd, ib_bdy) = ir 750 ENDDO 751 ENDDO 752 ENDDO 445 CALL bdy_coords_seg( nbidta, nbjdta, nbrdta ) 753 446 754 447 ! Deal with duplicated points … … 764 457 DO ib2 = 1, nblendta(igrd,ib_bdy2) 765 458 IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. & 766 & (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN767 ! IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &768 ! & nbidta(ib1, igrd, ib_bdy1), &769 ! & nbjdta(ib2, igrd, ib_bdy2)459 & (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN 460 ! IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', & 461 ! & nbidta(ib1, igrd, ib_bdy1), & 462 ! & nbjdta(ib2, igrd, ib_bdy2) 770 463 ! keep only points with the lowest distance to boundary: 771 464 IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN 772 nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2773 nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2465 nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2 466 nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2 774 467 ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN 775 nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1776 nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1777 ! Arbitrary choice if distances are the same:468 nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 469 nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 470 ! Arbitrary choice if distances are the same: 778 471 ELSE 779 nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1780 nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1472 nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 473 nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 781 474 ENDIF 782 475 END IF … … 787 480 END DO 788 481 END DO 789 790 ! Work out dimensions of boundary data on each processor 791 ! ------------------------------------------------------ 792 793 ! Rather assume that boundary data indices are given on global domain 794 ! TO BE DISCUSSED ? 795 ! iw = mig(1) + 1 ! if monotasking and no zoom, iw=2 796 ! ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1 797 ! is = mjg(1) + 1 ! if monotasking and no zoom, is=2 798 ! in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1 799 iwe = mig(1) - 1 + 2 ! if monotasking and no zoom, iw=2 800 ies = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1 801 iso = mjg(1) - 1 + 2 ! if monotasking and no zoom, is=2 802 ino = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1 803 804 ALLOCATE( nbondi_bdy(nb_bdy)) 805 ALLOCATE( nbondj_bdy(nb_bdy)) 806 nbondi_bdy(:)=2 807 nbondj_bdy(:)=2 808 ALLOCATE( nbondi_bdy_b(nb_bdy)) 809 ALLOCATE( nbondj_bdy_b(nb_bdy)) 810 nbondi_bdy_b(:)=2 811 nbondj_bdy_b(:)=2 812 813 ! Work out dimensions of boundary data on each neighbour process 814 IF(nbondi == 0) THEN 815 iw_b(1) = 1 + nimppt(nowe+1) 816 ie_b(1) = 1 + nimppt(nowe+1)+nlcit(nowe+1)-3 817 is_b(1) = 1 + njmppt(nowe+1) 818 in_b(1) = 1 + njmppt(nowe+1)+nlcjt(nowe+1)-3 819 820 iw_b(2) = 1 + nimppt(noea+1) 821 ie_b(2) = 1 + nimppt(noea+1)+nlcit(noea+1)-3 822 is_b(2) = 1 + njmppt(noea+1) 823 in_b(2) = 1 + njmppt(noea+1)+nlcjt(noea+1)-3 824 ELSEIF(nbondi == 1) THEN 825 iw_b(1) = 1 + nimppt(nowe+1) 826 ie_b(1) = 1 + nimppt(nowe+1)+nlcit(nowe+1)-3 827 is_b(1) = 1 + njmppt(nowe+1) 828 in_b(1) = 1 + njmppt(nowe+1)+nlcjt(nowe+1)-3 829 ELSEIF(nbondi == -1) THEN 830 iw_b(2) = 1 + nimppt(noea+1) 831 ie_b(2) = 1 + nimppt(noea+1)+nlcit(noea+1)-3 832 is_b(2) = 1 + njmppt(noea+1) 833 in_b(2) = 1 + njmppt(noea+1)+nlcjt(noea+1)-3 834 ENDIF 835 836 IF(nbondj == 0) THEN 837 iw_b(3) = 1 + nimppt(noso+1) 838 ie_b(3) = 1 + nimppt(noso+1)+nlcit(noso+1)-3 839 is_b(3) = 1 + njmppt(noso+1) 840 in_b(3) = 1 + njmppt(noso+1)+nlcjt(noso+1)-3 841 842 iw_b(4) = 1 + nimppt(nono+1) 843 ie_b(4) = 1 + nimppt(nono+1)+nlcit(nono+1)-3 844 is_b(4) = 1 + njmppt(nono+1) 845 in_b(4) = 1 + njmppt(nono+1)+nlcjt(nono+1)-3 846 ELSEIF(nbondj == 1) THEN 847 iw_b(3) = 1 + nimppt(noso+1) 848 ie_b(3) = 1 + nimppt(noso+1)+nlcit(noso+1)-3 849 is_b(3) = 1 + njmppt(noso+1) 850 in_b(3) = 1 + njmppt(noso+1)+nlcjt(noso+1)-3 851 ELSEIF(nbondj == -1) THEN 852 iw_b(4) = 1 + nimppt(nono+1) 853 ie_b(4) = 1 + nimppt(nono+1)+nlcit(nono+1)-3 854 is_b(4) = 1 + njmppt(nono+1) 855 in_b(4) = 1 + njmppt(nono+1)+nlcjt(nono+1)-3 856 ENDIF 857 482 ! 483 ! Find lenght of boundaries and rim on local mpi domain 484 !------------------------------------------------------ 485 ! 486 iwe = mig(1) 487 ies = mig(jpi) 488 iso = mjg(1) 489 ino = mjg(jpj) 490 ! 858 491 DO ib_bdy = 1, nb_bdy 859 492 DO igrd = 1, jpbgrd 860 icount = 0 861 icountr = 0 862 idx_bdy(ib_bdy)%nblen(igrd) = 0 863 idx_bdy(ib_bdy)%nblenrim(igrd) = 0 493 icount = 0 ! initialization of local bdy length 494 icountr = 0 ! initialization of local rim 0 and rim 1 bdy length 495 icountr0 = 0 ! initialization of local rim 0 bdy length 496 idx_bdy(ib_bdy)%nblen(igrd) = 0 497 idx_bdy(ib_bdy)%nblenrim(igrd) = 0 498 idx_bdy(ib_bdy)%nblenrim0(igrd) = 0 864 499 DO ib = 1, nblendta(igrd,ib_bdy) 865 500 ! check that data is in correct order in file 866 ibm1 = MAX(1,ib-1) 867 IF(lwp) THEN ! Since all procs read global data only need to do this check on one proc... 868 IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN 501 IF( ib > 1 ) THEN 502 IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ib-1,igrd,ib_bdy) ) THEN 869 503 CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', & 870 871 872 ENDIF 504 & ' in order of distance from edge nbr A utility for re-ordering ', & 505 & ' boundary coordinates and data files exists in the TOOLS/OBC directory') 506 ENDIF 873 507 ENDIF 874 508 ! check if point is in local domain … … 876 510 & nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino ) THEN 877 511 ! 878 icount = icount 879 !880 IF( nbrdta(ib,igrd,ib_bdy) == 1 ) icountr = icountr+1512 icount = icount + 1 513 IF( nbrdta(ib,igrd,ib_bdy) == 1 .OR. nbrdta(ib,igrd,ib_bdy) == 0 ) icountr = icountr + 1 514 IF( nbrdta(ib,igrd,ib_bdy) == 0 ) icountr0 = icountr0 + 1 881 515 ENDIF 882 516 END DO 883 idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc 884 idx_bdy(ib_bdy)%nblen (igrd) = icount !: length of boundary data on each proc 885 END DO ! igrd 517 idx_bdy(ib_bdy)%nblen (igrd) = icount !: length of boundary data on each proc 518 idx_bdy(ib_bdy)%nblenrim (igrd) = icountr !: length of rim 0 and rim 1 boundary data on each proc 519 idx_bdy(ib_bdy)%nblenrim0(igrd) = icountr0 !: length of rim 0 boundary data on each proc 520 END DO ! igrd 886 521 887 522 ! Allocate index arrays for this boundary set … … 893 528 & idx_bdy(ib_bdy)%nbd (ilen1,jpbgrd) , & 894 529 & idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) , & 530 & idx_bdy(ib_bdy)%ntreat(ilen1,jpbgrd) , & 895 531 & idx_bdy(ib_bdy)%nbmap (ilen1,jpbgrd) , & 896 532 & idx_bdy(ib_bdy)%nbw (ilen1,jpbgrd) , & … … 900 536 ! Dispatch mapping indices and discrete distances on each processor 901 537 ! ----------------------------------------------------------------- 902 903 com_east = 0904 com_west = 0905 com_south = 0906 com_north = 0907 908 com_east_b = 0909 com_west_b = 0910 com_south_b = 0911 com_north_b = 0912 913 538 DO igrd = 1, jpbgrd 914 539 icount = 0 915 ! Loop on rimwidth to ensure outermost points come first in the local arrays.916 DO ir =1, nn_rimwidth(ib_bdy)540 ! Outer loop on rimwidth to ensure outermost points come first in the local arrays. 541 DO ir = 0, nn_rimwidth(ib_bdy) 917 542 DO ib = 1, nblendta(igrd,ib_bdy) 918 543 ! check if point is in local domain and equals ir … … 922 547 ! 923 548 icount = icount + 1 924 925 ! Rather assume that boundary data indices are given on global domain 926 ! TO BE DISCUSSED ? 927 ! idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1 928 ! idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 929 idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1 930 idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 931 ! check if point has to be sent 932 ii = idx_bdy(ib_bdy)%nbi(icount,igrd) 933 ij = idx_bdy(ib_bdy)%nbj(icount,igrd) 934 if((com_east .ne. 1) .and. (ii == (nlci-1)) .and. (nbondi .le. 0)) then 935 com_east = 1 936 elseif((com_west .ne. 1) .and. (ii == 2) .and. (nbondi .ge. 0) .and. (nbondi .ne. 2)) then 937 com_west = 1 938 endif 939 if((com_south .ne. 1) .and. (ij == 2) .and. (nbondj .ge. 0) .and. (nbondj .ne. 2)) then 940 com_south = 1 941 elseif((com_north .ne. 1) .and. (ij == (nlcj-1)) .and. (nbondj .le. 0)) then 942 com_north = 1 943 endif 549 idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1 ! global to local indexes 550 idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 ! global to local indexes 944 551 idx_bdy(ib_bdy)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_bdy) 945 552 idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib 946 553 ENDIF 947 ! check if point has to be received from a neighbour 948 IF(nbondi == 0) THEN 949 IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND. & 950 & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND. & 951 & nbrdta(ib,igrd,ib_bdy) == ir ) THEN 952 ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2 953 if( ii == (nlcit(nowe+1)-1) ) then 954 ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2 955 if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then 956 com_south = 1 957 elseif((ij == nlcjt(nowe+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then 958 com_north = 1 959 endif 960 com_west_b = 1 961 endif 962 ENDIF 963 IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND. & 964 & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND. & 965 & nbrdta(ib,igrd,ib_bdy) == ir ) THEN 966 ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2 967 if( ii == 2 ) then 968 ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2 969 if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then 970 com_south = 1 971 elseif((ij == nlcjt(noea+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then 972 com_north = 1 973 endif 974 com_east_b = 1 975 endif 976 ENDIF 977 ELSEIF(nbondi == 1) THEN 978 IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND. & 979 & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND. & 980 & nbrdta(ib,igrd,ib_bdy) == ir ) THEN 981 ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2 982 if( ii == (nlcit(nowe+1)-1) ) then 983 ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2 984 if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then 985 com_south = 1 986 elseif((ij == nlcjt(nowe+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then 987 com_north = 1 988 endif 989 com_west_b = 1 990 endif 991 ENDIF 992 ELSEIF(nbondi == -1) THEN 993 IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND. & 994 & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND. & 995 & nbrdta(ib,igrd,ib_bdy) == ir ) THEN 996 ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2 997 if( ii == 2 ) then 998 ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2 999 if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then 1000 com_south = 1 1001 elseif((ij == nlcjt(noea+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then 1002 com_north = 1 1003 endif 1004 com_east_b = 1 1005 endif 1006 ENDIF 1007 ENDIF 1008 IF(nbondj == 0) THEN 1009 IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1 & 1010 & .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. & 1011 & nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN 1012 com_north_b = 1 1013 ENDIF 1014 IF(com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 & 1015 &.OR. nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. & 1016 & nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN 1017 com_south_b = 1 1018 ENDIF 1019 IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND. & 1020 & nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND. & 1021 & nbrdta(ib,igrd,ib_bdy) == ir ) THEN 1022 ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2 1023 if((com_south_b .ne. 1) .and. (ij == (nlcjt(noso+1)-1))) then 1024 com_south_b = 1 1025 endif 1026 ENDIF 1027 IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND. & 1028 & nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND. & 1029 & nbrdta(ib,igrd,ib_bdy) == ir ) THEN 1030 ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2 1031 if((com_north_b .ne. 1) .and. (ij == 2)) then 1032 com_north_b = 1 1033 endif 1034 ENDIF 1035 ELSEIF(nbondj == 1) THEN 1036 IF( com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 .OR. & 1037 & nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. & 1038 & nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN 1039 com_south_b = 1 1040 ENDIF 1041 IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND. & 1042 & nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND. & 1043 & nbrdta(ib,igrd,ib_bdy) == ir ) THEN 1044 ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2 1045 if((com_south_b .ne. 1) .and. (ij == (nlcjt(noso+1)-1))) then 1046 com_south_b = 1 1047 endif 1048 ENDIF 1049 ELSEIF(nbondj == -1) THEN 1050 IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1 & 1051 & .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. & 1052 & nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN 1053 com_north_b = 1 1054 ENDIF 1055 IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND. & 1056 & nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND. & 1057 & nbrdta(ib,igrd,ib_bdy) == ir ) THEN 1058 ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2 1059 if((com_north_b .ne. 1) .and. (ij == 2)) then 1060 com_north_b = 1 1061 endif 1062 ENDIF 1063 ENDIF 1064 ENDDO 1065 ENDDO 1066 ENDDO 1067 1068 ! definition of the i- and j- direction local boundaries arrays used for sending the boundaries 1069 IF( (com_east == 1) .and. (com_west == 1) ) THEN ; nbondi_bdy(ib_bdy) = 0 1070 ELSEIF( (com_east == 1) .and. (com_west == 0) ) THEN ; nbondi_bdy(ib_bdy) = -1 1071 ELSEIF( (com_east == 0) .and. (com_west == 1) ) THEN ; nbondi_bdy(ib_bdy) = 1 1072 ENDIF 1073 IF( (com_north == 1) .and. (com_south == 1) ) THEN ; nbondj_bdy(ib_bdy) = 0 1074 ELSEIF( (com_north == 1) .and. (com_south == 0) ) THEN ; nbondj_bdy(ib_bdy) = -1 1075 ELSEIF( (com_north == 0) .and. (com_south == 1) ) THEN ; nbondj_bdy(ib_bdy) = 1 1076 ENDIF 1077 1078 ! definition of the i- and j- direction local boundaries arrays used for receiving the boundaries 1079 IF( (com_east_b == 1) .and. (com_west_b == 1) ) THEN ; nbondi_bdy_b(ib_bdy) = 0 1080 ELSEIF( (com_east_b == 1) .and. (com_west_b == 0) ) THEN ; nbondi_bdy_b(ib_bdy) = -1 1081 ELSEIF( (com_east_b == 0) .and. (com_west_b == 1) ) THEN ; nbondi_bdy_b(ib_bdy) = 1 1082 ENDIF 1083 IF( (com_north_b == 1) .and. (com_south_b == 1) ) THEN ; nbondj_bdy_b(ib_bdy) = 0 1084 ELSEIF( (com_north_b == 1) .and. (com_south_b == 0) ) THEN ; nbondj_bdy_b(ib_bdy) = -1 1085 ELSEIF( (com_north_b == 0) .and. (com_south_b == 1) ) THEN ; nbondj_bdy_b(ib_bdy) = 1 1086 ENDIF 554 END DO 555 END DO 556 END DO ! igrd 557 558 END DO ! ib_bdy 559 560 ! Initialize array indicating communications in bdy 561 ! ------------------------------------------------- 562 ALLOCATE( lsend_bdy(nb_bdy,jpbgrd,4,0:1), lrecv_bdy(nb_bdy,jpbgrd,4,0:1) ) 563 lsend_bdy(:,:,:,:) = .false. 564 lrecv_bdy(:,:,:,:) = .false. 565 566 DO ib_bdy = 1, nb_bdy 567 DO igrd = 1, jpbgrd 568 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! only the rim triggers communications, see bdy routines 569 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 570 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 571 IF( ib .LE. idx_bdy(ib_bdy)%nblenrim0(igrd) ) THEN ; ir = 0 572 ELSE ; ir = 1 573 END IF 574 ! 575 ! check if point has to be sent to a neighbour 576 ! W neighbour and on the inner left side 577 IF( ii == 2 .and. (nbondi == 0 .or. nbondi == 1) ) lsend_bdy(ib_bdy,igrd,1,ir) = .true. 578 ! E neighbour and on the inner right side 579 IF( ii == jpi-1 .and. (nbondi == 0 .or. nbondi == -1) ) lsend_bdy(ib_bdy,igrd,2,ir) = .true. 580 ! S neighbour and on the inner down side 581 IF( ij == 2 .and. (nbondj == 0 .or. nbondj == 1) ) lsend_bdy(ib_bdy,igrd,3,ir) = .true. 582 ! N neighbour and on the inner up side 583 IF( ij == jpj-1 .and. (nbondj == 0 .or. nbondj == -1) ) lsend_bdy(ib_bdy,igrd,4,ir) = .true. 584 ! 585 ! check if point has to be received from a neighbour 586 ! W neighbour and on the outter left side 587 IF( ii == 1 .and. (nbondi == 0 .or. nbondi == 1) ) lrecv_bdy(ib_bdy,igrd,1,ir) = .true. 588 ! E neighbour and on the outter right side 589 IF( ii == jpi .and. (nbondi == 0 .or. nbondi == -1) ) lrecv_bdy(ib_bdy,igrd,2,ir) = .true. 590 ! S neighbour and on the outter down side 591 IF( ij == 1 .and. (nbondj == 0 .or. nbondj == 1) ) lrecv_bdy(ib_bdy,igrd,3,ir) = .true. 592 ! N neighbour and on the outter up side 593 IF( ij == jpj .and. (nbondj == 0 .or. nbondj == -1) ) lrecv_bdy(ib_bdy,igrd,4,ir) = .true. 594 ! 595 END DO 596 END DO ! igrd 1087 597 1088 598 ! Compute rim weights for FRS scheme … … 1090 600 DO igrd = 1, jpbgrd 1091 601 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 1092 nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)1093 idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( nbr - 1 ) *0.5 ) ! tanh formulation1094 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic1095 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)) ! linear1096 END DO 1097 END DO 602 ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) ) ! both rim 0 and rim 1 have the same weights 603 idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( ir - 1 ) *0.5 ) ! tanh formulation 604 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic 605 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)) ! linear 606 END DO 607 END DO 1098 608 1099 609 ! Compute damping coefficients … … 1101 611 DO igrd = 1, jpbgrd 1102 612 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 1103 nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)613 ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) ) ! both rim 0 and rim 1 have the same damping coefficients 1104 614 idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 1105 & *(REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic615 & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic 1106 616 idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) & 1107 & *(REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic1108 END DO 1109 END DO 1110 1111 END DO 617 & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic 618 END DO 619 END DO 620 621 END DO ! ib_bdy 1112 622 1113 623 ! ------------------------------------------------------ 1114 624 ! Initialise masks and find normal/tangential directions 1115 625 ! ------------------------------------------------------ 626 627 ! ------------------------------------------ 628 ! handle rim0, do as if rim 1 was free ocean 629 ! ------------------------------------------ 630 631 ztmask(:,:) = tmask(:,:,1) ; zumask(:,:) = umask(:,:,1) ; zvmask(:,:) = vmask(:,:,1) 632 ! For the flagu/flagv calculation below we require a version of fmask without 633 ! the land boundary condition (shlat) included: 634 DO ij = 1, jpjm1 635 DO ii = 1, jpim1 636 zfmask(ii,ij) = ztmask(ii,ij ) * ztmask(ii+1,ij ) & 637 & * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 638 END DO 639 END DO 640 CALL lbc_lnk( 'bdyini', zfmask, 'F', 1.0_wp ) 1116 641 1117 642 ! Read global 2D mask at T-points: bdytmask … … 1119 644 ! bdytmask = 1 on the computational domain AND on open boundaries 1120 645 ! = 0 elsewhere 1121 646 1122 647 bdytmask(:,:) = ssmask(:,:) 1123 648 1124 649 ! Derive mask on U and V grid from mask on T grid 1125 1126 bdyumask(:,:) = 0._wp1127 bdyvmask(:,:) = 0._wp1128 650 DO ij = 1, jpjm1 1129 651 DO ii = 1, jpim1 1130 bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1, ij)652 bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1,ij ) 1131 653 bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii ,ij+1) 1132 654 END DO 1133 655 END DO 1134 CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1. , bdyvmask, 'V', 1. ) ! Lateral boundary cond. 1135 1136 ! bdy masks are now set to zero on boundary points: 1137 ! 1138 igrd = 1 656 CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1.0_wp , bdyvmask, 'V', 1.0_wp ) ! Lateral boundary cond. 657 658 ! bdy masks are now set to zero on rim 0 points: 1139 659 DO ib_bdy = 1, nb_bdy 1140 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1141 bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp 1142 END DO 1143 END DO 1144 ! 1145 igrd = 2 660 DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1) ! extent of rim 0 661 bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp 662 END DO 663 DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2) ! extent of rim 0 664 bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp 665 END DO 666 DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3) ! extent of rim 0 667 bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp 668 END DO 669 END DO 670 671 CALL bdy_rim_treat( zumask, zvmask, zfmask, .true. ) ! compute flagu, flagv, ntreat on rim 0 672 673 ! ------------------------------------ 674 ! handle rim1, do as if rim 0 was land 675 ! ------------------------------------ 676 677 ! z[tuv]mask are now set to zero on rim 0 points: 1146 678 DO ib_bdy = 1, nb_bdy 1147 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1148 bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp 1149 END DO 1150 END DO 1151 ! 1152 igrd = 3 679 DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1) ! extent of rim 0 680 ztmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp 681 END DO 682 DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2) ! extent of rim 0 683 zumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp 684 END DO 685 DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3) ! extent of rim 0 686 zvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp 687 END DO 688 END DO 689 690 ! Recompute zfmask 691 DO ij = 1, jpjm1 692 DO ii = 1, jpim1 693 zfmask(ii,ij) = ztmask(ii,ij ) * ztmask(ii+1,ij ) & 694 & * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 695 END DO 696 END DO 697 CALL lbc_lnk( 'bdyini', zfmask, 'F', 1.0_wp ) 698 699 ! bdy masks are now set to zero on rim1 points: 1153 700 DO ib_bdy = 1, nb_bdy 1154 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1155 bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp 1156 END DO 1157 END DO 1158 1159 ! For the flagu/flagv calculation below we require a version of fmask without 1160 ! the land boundary condition (shlat) included: 1161 zfmask(:,:) = 0 1162 DO ij = 2, jpjm1 1163 DO ii = 2, jpim1 1164 zfmask(ii,ij) = tmask(ii,ij ,1) * tmask(ii+1,ij ,1) & 1165 & * tmask(ii,ij+1,1) * tmask(ii+1,ij+1,1) 1166 END DO 1167 END DO 1168 1169 ! Lateral boundary conditions 1170 CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. ) 1171 CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1. , bdyvmask, 'V', 1., bdytmask, 'T', 1. ) 701 DO ib = idx_bdy(ib_bdy)%nblenrim0(1) + 1, idx_bdy(ib_bdy)%nblenrim(1) ! extent of rim 1 702 bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp 703 END DO 704 DO ib = idx_bdy(ib_bdy)%nblenrim0(2) + 1, idx_bdy(ib_bdy)%nblenrim(2) ! extent of rim 1 705 bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp 706 END DO 707 DO ib = idx_bdy(ib_bdy)%nblenrim0(3) + 1, idx_bdy(ib_bdy)%nblenrim(3) ! extent of rim 1 708 bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp 709 END DO 710 END DO 711 712 CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. ) ! compute flagu, flagv, ntreat on rim 1 713 ! 714 ! Check which boundaries might need communication 715 ALLOCATE( lsend_bdyint(nb_bdy,jpbgrd,4,0:1), lrecv_bdyint(nb_bdy,jpbgrd,4,0:1) ) 716 lsend_bdyint(:,:,:,:) = .false. 717 lrecv_bdyint(:,:,:,:) = .false. 718 ALLOCATE( lsend_bdyext(nb_bdy,jpbgrd,4,0:1), lrecv_bdyext(nb_bdy,jpbgrd,4,0:1) ) 719 lsend_bdyext(:,:,:,:) = .false. 720 lrecv_bdyext(:,:,:,:) = .false. 721 ! 722 DO igrd = 1, jpbgrd 723 DO ib_bdy = 1, nb_bdy 724 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 725 IF( idx_bdy(ib_bdy)%ntreat(ib,igrd) == -1 ) CYCLE 726 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 727 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 728 ir = idx_bdy(ib_bdy)%nbr(ib,igrd) 729 flagu = NINT(idx_bdy(ib_bdy)%flagu(ib,igrd)) 730 flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd)) 731 iibe = ii - flagu ! neighbouring point towards the exterior of the computational domain 732 ijbe = ij - flagv 733 iibi = ii + flagu ! neighbouring point towards the interior of the computational domain 734 ijbi = ij + flagv 735 CALL find_neib( ii, ij, idx_bdy(ib_bdy)%ntreat(ib,igrd), ii1, ij1, ii2, ij2, ii3, ij3 ) ! free ocean neighbours 736 ! 737 ! search neighbour in the west/east direction 738 ! Rim is on the halo and computed ocean is towards exterior of mpi domain 739 ! <-- (o exterior) --> 740 ! (1) o|x OR (2) x|o 741 ! |___ ___| 742 IF( iibi == 0 .OR. ii1 == 0 .OR. ii2 == 0 .OR. ii3 == 0 ) lrecv_bdyint(ib_bdy,igrd,1,ir) = .true. 743 IF( iibi == jpi+1 .OR. ii1 == jpi+1 .OR. ii2 == jpi+1 .OR. ii3 == jpi+1 ) lrecv_bdyint(ib_bdy,igrd,2,ir) = .true. 744 IF( iibe == 0 ) lrecv_bdyext(ib_bdy,igrd,1,ir) = .true. 745 IF( iibe == jpi+1 ) lrecv_bdyext(ib_bdy,igrd,2,ir) = .true. 746 ! Check if neighbour has its rim parallel to its mpi subdomain border and located next to its halo 747 ! :¨¨¨¨¨|¨¨--> | | <--¨¨|¨¨¨¨¨: 748 ! : | x:o | neighbour limited by ... would need o | o:x | : 749 ! :.....|_._:_____| (1) W neighbour E neighbour (2) |_____:_._|.....: 750 IF( ii == 2 .AND. ( nbondi == 1 .OR. nbondi == 0 ) .AND. & 751 & ( iibi == 3 .OR. ii1 == 3 .OR. ii2 == 3 .OR. ii3 == 3 ) ) lsend_bdyint(ib_bdy,igrd,1,ir)=.true. 752 IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. & 753 & ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) ) lsend_bdyint(ib_bdy,igrd,2,ir)=.true. 754 IF( ii == 2 .AND. ( nbondi == 1 .OR. nbondi == 0 ) .AND. iibe == 3 ) lsend_bdyext(ib_bdy,igrd,1,ir)=.true. 755 IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2 ) lsend_bdyext(ib_bdy,igrd,2,ir)=.true. 756 ! 757 ! search neighbour in the north/south direction 758 ! Rim is on the halo and computed ocean is towards exterior of mpi domain 759 !(3) | | ^ ___o___ 760 ! | |___x___| OR | | x | 761 ! v o (4) | | 762 IF( ijbi == 0 .OR. ij1 == 0 .OR. ij2 == 0 .OR. ij3 == 0 ) lrecv_bdyint(ib_bdy,igrd,3,ir) = .true. 763 IF( ijbi == jpj+1 .OR. ij1 == jpj+1 .OR. ij2 == jpj+1 .OR. ij3 == jpj+1 ) lrecv_bdyint(ib_bdy,igrd,4,ir) = .true. 764 IF( ijbe == 0 ) lrecv_bdyext(ib_bdy,igrd,3,ir) = .true. 765 IF( ijbe == jpj+1 ) lrecv_bdyext(ib_bdy,igrd,4,ir) = .true. 766 ! Check if neighbour has its rim parallel to its mpi subdomain _________ border and next to its halo 767 ! ^ | o | : : 768 ! | |¨¨¨¨x¨¨¨¨| neighbour limited by ... would need o | |....x....| 769 ! :_________: (3) S neighbour N neighbour (4) v | o | 770 IF( ij == 2 .AND. ( nbondj == 1 .OR. nbondj == 0 ) .AND. & 771 & ( ijbi == 3 .OR. ij1 == 3 .OR. ij2 == 3 .OR. ij3 == 3 ) ) lsend_bdyint(ib_bdy,igrd,3,ir)=.true. 772 IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. & 773 & ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) ) lsend_bdyint(ib_bdy,igrd,4,ir)=.true. 774 IF( ij == 2 .AND. ( nbondj == 1 .OR. nbondj == 0 ) .AND. ijbe == 3 ) lsend_bdyext(ib_bdy,igrd,3,ir)=.true. 775 IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2 ) lsend_bdyext(ib_bdy,igrd,4,ir)=.true. 776 END DO 777 END DO 778 END DO 779 780 DO ib_bdy = 1,nb_bdy 781 IF( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. & 782 & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. & 783 & cn_tra(ib_bdy) == 'orlanski' .OR. cn_tra(ib_bdy) == 'orlanski_npo' ) THEN 784 DO igrd = 1, jpbgrd 785 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 786 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 787 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 788 IF( mig(ii) > 2 .AND. mig(ii) < jpiglo-2 .AND. mjg(ij) > 2 .AND. mjg(ij) < jpjglo-2 ) THEN 789 WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain' 790 CALL ctl_stop( ctmp1 ) 791 END IF 792 END DO 793 END DO 794 END IF 795 END DO 796 ! 797 DEALLOCATE( nbidta, nbjdta, nbrdta ) 798 ! 799 END SUBROUTINE bdy_def 800 801 802 SUBROUTINE bdy_rim_treat( pumask, pvmask, pfmask, lrim0 ) 803 !!---------------------------------------------------------------------- 804 !! *** ROUTINE bdy_rim_treat *** 805 !! 806 !! ** Purpose : Initialize structures ( flagu, flagv, ntreat ) indicating how rim points 807 !! are to be handled in the boundary condition treatment 808 !! 809 !! ** Method : - to handle rim 0 zmasks must indicate ocean points (set at one on rim 0 and rim 1 and interior) 810 !! and bdymasks must be set at 0 on rim 0 (set at one on rim 1 and interior) 811 !! (as if rim 1 was free ocean) 812 !! - to handle rim 1 zmasks must be set at 0 on rim 0 (set at one on rim 1 and interior) 813 !! and bdymasks must indicate free ocean points (set at one on interior) 814 !! (as if rim 0 was land) 815 !! - we can then check in which direction the interior of the computational domain is with the difference 816 !! mask array values on both sides to compute flagu and flagv 817 !! - and look at the ocean neighbours to compute ntreat 818 !!---------------------------------------------------------------------- 819 REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in ) :: pfmask ! temporary fmask excluding coastal boundary condition (shlat) 820 REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in ) :: pumask, pvmask ! temporary t/u/v mask array 821 LOGICAL , INTENT (in ) :: lrim0 ! .true. -> rim 0 .false. -> rim 1 822 INTEGER :: ib_bdy, ii, ij, igrd, ib, icount ! dummy loop indices 823 INTEGER :: i_offset, j_offset, inn ! local integer 824 INTEGER :: ibeg, iend ! local integer 825 LOGICAL :: llnon, llson, llean, llwen ! local logicals indicating the presence of a ocean neighbour 826 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! pointer to 2D mask fields 827 REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars 828 CHARACTER(LEN=1), DIMENSION(jpbgrd) :: cgrid 829 REAL(wp) , DIMENSION(jpi,jpj) :: ztmp 830 !!---------------------------------------------------------------------- 831 832 cgrid = (/'t','u','v'/) 833 1172 834 DO ib_bdy = 1, nb_bdy ! Indices and directions of rim velocity components 1173 1174 idx_bdy(ib_bdy)%flagu(:,:) = 0._wp1175 idx_bdy(ib_bdy)%flagv(:,:) = 0._wp1176 icount = 01177 835 1178 836 ! Calculate relationship of U direction to the local orientation of the boundary … … 1180 838 ! flagu = 0 : u is tangential 1181 839 ! flagu = 1 : u is normal to the boundary and is direction is inward 1182 1183 840 DO igrd = 1, jpbgrd 1184 841 SELECT CASE( igrd ) 1185 CASE( 1 ) ; pmask => umask (:,:,1); i_offset = 01186 CASE( 2 ) ; pmask => bdytmask(:,:); i_offset = 11187 CASE( 3 ) ; pmask => zfmask (:,:); i_offset = 0842 CASE( 1 ) ; zmask => pumask ; i_offset = 0 843 CASE( 2 ) ; zmask => bdytmask ; i_offset = 1 844 CASE( 3 ) ; zmask => pfmask ; i_offset = 0 1188 845 END SELECT 1189 846 icount = 0 1190 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1191 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 1192 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 1193 zefl = pmask(nbi+i_offset-1,nbj) 1194 zwfl = pmask(nbi+i_offset,nbj) 847 ztmp(:,:) = -999._wp 848 IF( lrim0 ) THEN ! extent of rim 0 849 ibeg = 1 ; iend = idx_bdy(ib_bdy)%nblenrim0(igrd) 850 ELSE ! extent of rim 1 851 ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1 ; iend = idx_bdy(ib_bdy)%nblenrim(igrd) 852 END IF 853 DO ib = ibeg, iend 854 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 855 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 856 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE 857 zwfl = zmask(ii+i_offset-1,ij) 858 zefl = zmask(ii+i_offset ,ij) 1195 859 ! This error check only works if you are using the bdyXmask arrays 1196 IF( i_offset == 1 .and. zefl + zwfl == 2 ) THEN860 IF( i_offset == 1 .and. zefl + zwfl == 2. ) THEN 1197 861 icount = icount + 1 1198 IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig( nbi),mjg(nbj)862 IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij) 1199 863 ELSE 1200 idx_bdy(ib_bdy)%flagu(ib,igrd) = -zefl + zwfl864 ztmp(ii,ij) = -zwfl + zefl 1201 865 ENDIF 1202 866 END DO 1203 867 IF( icount /= 0 ) THEN 1204 WRITE(ctmp1,*) ' E R R O R :Some ',cgrid(igrd),' grid points,', &868 WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,', & 1205 869 ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy 1206 WRITE(ctmp2,*) ' ========== ' 1207 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 870 CALL ctl_stop( ctmp1 ) 1208 871 ENDIF 872 SELECT CASE( igrd ) 873 CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp ) 874 CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp ) 875 CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp ) 876 END SELECT 877 DO ib = ibeg, iend 878 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 879 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 880 idx_bdy(ib_bdy)%flagu(ib,igrd) = ztmp(ii,ij) 881 END DO 1209 882 END DO 1210 883 … … 1213 886 ! flagv = 0 : v is tangential 1214 887 ! flagv = 1 : v is normal to the boundary and is direction is inward 1215 1216 888 DO igrd = 1, jpbgrd 1217 889 SELECT CASE( igrd ) 1218 CASE( 1 ) ; pmask => vmask (:,:,1); j_offset = 01219 CASE( 2 ) ; pmask => zfmask(:,:); j_offset = 01220 CASE( 3 ) ; pmask => bdytmask; j_offset = 1890 CASE( 1 ) ; zmask => pvmask ; j_offset = 0 891 CASE( 2 ) ; zmask => pfmask ; j_offset = 0 892 CASE( 3 ) ; zmask => bdytmask ; j_offset = 1 1221 893 END SELECT 1222 894 icount = 0 1223 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1224 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 1225 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 1226 znfl = pmask(nbi,nbj+j_offset-1) 1227 zsfl = pmask(nbi,nbj+j_offset ) 895 ztmp(:,:) = -999._wp 896 IF( lrim0 ) THEN ! extent of rim 0 897 ibeg = 1 ; iend = idx_bdy(ib_bdy)%nblenrim0(igrd) 898 ELSE ! extent of rim 1 899 ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1 ; iend = idx_bdy(ib_bdy)%nblenrim(igrd) 900 END IF 901 DO ib = ibeg, iend 902 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 903 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 904 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE 905 zsfl = zmask(ii,ij+j_offset-1) 906 znfl = zmask(ii,ij+j_offset ) 1228 907 ! This error check only works if you are using the bdyXmask arrays 1229 IF( j_offset == 1 .and. znfl + zsfl == 2 ) THEN1230 IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig( nbi),mjg(nbj)908 IF( j_offset == 1 .and. znfl + zsfl == 2. ) THEN 909 IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij) 1231 910 icount = icount + 1 1232 911 ELSE 1233 idx_bdy(ib_bdy)%flagv(ib,igrd) = -znfl + zsfl912 ztmp(ii,ij) = -zsfl + znfl 1234 913 END IF 1235 914 END DO 1236 915 IF( icount /= 0 ) THEN 1237 WRITE(ctmp1,*) ' E R R O R :Some ',cgrid(igrd),' grid points,', &916 WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,', & 1238 917 ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy 1239 WRITE(ctmp2,*) ' ========== ' 1240 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1241 ENDIF 1242 END DO 1243 ! 1244 END DO 1245 ! 1246 ! Tidy up 1247 !-------- 1248 IF( nb_bdy>0 ) DEALLOCATE( nbidta, nbjdta, nbrdta ) 1249 ! 1250 END SUBROUTINE bdy_segs 1251 918 CALL ctl_stop( ctmp1 ) 919 ENDIF 920 SELECT CASE( igrd ) 921 CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp ) 922 CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp ) 923 CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp ) 924 END SELECT 925 DO ib = ibeg, iend 926 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 927 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 928 idx_bdy(ib_bdy)%flagv(ib,igrd) = ztmp(ii,ij) 929 END DO 930 END DO 931 ! 932 END DO ! ib_bdy 933 934 DO ib_bdy = 1, nb_bdy 935 DO igrd = 1, jpbgrd 936 SELECT CASE( igrd ) 937 CASE( 1 ) ; zmask => bdytmask 938 CASE( 2 ) ; zmask => bdyumask 939 CASE( 3 ) ; zmask => bdyvmask 940 END SELECT 941 ztmp(:,:) = -999._wp 942 IF( lrim0 ) THEN ! extent of rim 0 943 ibeg = 1 ; iend = idx_bdy(ib_bdy)%nblenrim0(igrd) 944 ELSE ! extent of rim 1 945 ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1 ; iend = idx_bdy(ib_bdy)%nblenrim(igrd) 946 END IF 947 DO ib = ibeg, iend 948 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 949 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 950 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE 951 llnon = zmask(ii ,ij+1) == 1. 952 llson = zmask(ii ,ij-1) == 1. 953 llean = zmask(ii+1,ij ) == 1. 954 llwen = zmask(ii-1,ij ) == 1. 955 inn = COUNT( (/ llnon, llson, llean, llwen /) ) 956 IF( inn == 0 ) THEN ! no neighbours -> interior of a corner or cluster of rim points 957 ! ! ! _____ ! _____ ! __ __ 958 ! 1 | o ! 2 o | ! 3 | x ! 4 x | ! | | -> error 959 ! |_x_ _ ! _ _x_| ! | o ! o | ! |x_x| 960 IF( zmask(ii+1,ij+1) == 1. ) THEN ; ztmp(ii,ij) = 1. 961 ELSEIF( zmask(ii-1,ij+1) == 1. ) THEN ; ztmp(ii,ij) = 2. 962 ELSEIF( zmask(ii+1,ij-1) == 1. ) THEN ; ztmp(ii,ij) = 3. 963 ELSEIF( zmask(ii-1,ij-1) == 1. ) THEN ; ztmp(ii,ij) = 4. 964 ELSE ; ztmp(ii,ij) = -1. 965 WRITE(ctmp1,*) 'Problem with ',cgrid(igrd) ,' grid point', ii, ij, & 966 ' on boundary set ', ib_bdy, ' has no free ocean neighbour' 967 IF( lrim0 ) THEN 968 WRITE(ctmp2,*) ' There seems to be a cluster of rim 0 points.' 969 ELSE 970 WRITE(ctmp2,*) ' There seems to be a cluster of rim 1 points.' 971 END IF 972 CALL ctl_warn( ctmp1, ctmp2 ) 973 END IF 974 END IF 975 IF( inn == 1 ) THEN ! middle of linear bdy or incomplete corner ! ___ o 976 ! | ! | ! o ! ______ ! |x___ 977 ! 5 | x o ! 6 o x | ! 7 __x__ ! 8 x 978 ! | ! | ! ! o 979 IF( llean ) ztmp(ii,ij) = 5. 980 IF( llwen ) ztmp(ii,ij) = 6. 981 IF( llnon ) ztmp(ii,ij) = 7. 982 IF( llson ) ztmp(ii,ij) = 8. 983 END IF 984 IF( inn == 2 ) THEN ! exterior of a corner 985 ! o ! o ! _____| ! |_____ 986 ! 9 ____x o ! 10 o x___ ! 11 x o ! 12 o x 987 ! | ! | ! o ! o 988 IF( llnon .AND. llean ) ztmp(ii,ij) = 9. 989 IF( llnon .AND. llwen ) ztmp(ii,ij) = 10. 990 IF( llson .AND. llean ) ztmp(ii,ij) = 11. 991 IF( llson .AND. llwen ) ztmp(ii,ij) = 12. 992 END IF 993 IF( inn == 3 ) THEN ! 3 neighbours __ __ 994 ! |_ o ! o _| ! |_| ! o 995 ! 13 _| x o ! 14 o x |_ ! 15 o x o ! 16 o x o 996 ! | o ! o | ! o ! __|¨|__ 997 IF( llnon .AND. llean .AND. llson ) ztmp(ii,ij) = 13. 998 IF( llnon .AND. llwen .AND. llson ) ztmp(ii,ij) = 14. 999 IF( llwen .AND. llson .AND. llean ) ztmp(ii,ij) = 15. 1000 IF( llwen .AND. llnon .AND. llean ) ztmp(ii,ij) = 16. 1001 END IF 1002 IF( inn == 4 ) THEN 1003 WRITE(ctmp1,*) 'Problem with ',cgrid(igrd) ,' grid point', ii, ij, & 1004 ' on boundary set ', ib_bdy, ' have 4 neighbours' 1005 CALL ctl_stop( ctmp1 ) 1006 END IF 1007 END DO 1008 SELECT CASE( igrd ) 1009 CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp ) 1010 CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp ) 1011 CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp ) 1012 END SELECT 1013 DO ib = ibeg, iend 1014 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1015 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1016 idx_bdy(ib_bdy)%ntreat(ib,igrd) = NINT(ztmp(ii,ij)) 1017 END DO 1018 END DO 1019 END DO 1020 1021 END SUBROUTINE bdy_rim_treat 1022 1023 1024 SUBROUTINE find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 ) 1025 !!---------------------------------------------------------------------- 1026 !! *** ROUTINE find_neib *** 1027 !! 1028 !! ** Purpose : get ii1, ij1, ii2, ij2, ii3, ij3, the indices of 1029 !! the free ocean neighbours of (ii,ij) for bdy treatment 1030 !! 1031 !! ** Method : use itreat input to select a case 1032 !! N.B. ntreat is defined for all bdy points in routine bdy_rim_treat 1033 !! 1034 !!---------------------------------------------------------------------- 1035 INTEGER, INTENT(in ) :: ii, ij, itreat 1036 INTEGER, INTENT( out) :: ii1, ij1, ii2, ij2, ii3, ij3 1037 !!---------------------------------------------------------------------- 1038 SELECT CASE( itreat ) ! points that will be used by bdy routines, -1 will be discarded 1039 ! ! ! _____ ! _____ 1040 ! 1 | o ! 2 o | ! 3 | x ! 4 x | 1041 ! |_x_ _ ! _ _x_| ! | o ! o | 1042 CASE( 1 ) ; ii1 = ii+1 ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1043 CASE( 2 ) ; ii1 = ii-1 ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1044 CASE( 3 ) ; ii1 = ii+1 ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1045 CASE( 4 ) ; ii1 = ii-1 ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1046 ! | ! | ! o ! ______ ! or incomplete corner 1047 ! 5 | x o ! 6 o x | ! 7 __x__ ! 8 x ! 7 ____ o 1048 ! | ! | ! ! o ! |x___ 1049 CASE( 5 ) ; ii1 = ii+1 ; ij1 = ij ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1050 CASE( 6 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1051 CASE( 7 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1052 CASE( 8 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1053 ! o ! o ! _____| ! |_____ 1054 ! 9 ____x o ! 10 o x___ ! 11 x o ! 12 o x 1055 ! | ! | ! o ! o 1056 CASE( 9 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1 1057 CASE( 10 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1 1058 CASE( 11 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1 1059 CASE( 12 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1 1060 ! |_ o ! o _| ! ¨¨|_|¨¨ ! o 1061 ! 13 _| x o ! 14 o x |_ ! 15 o x o ! 16 o x o 1062 ! | o ! o | ! o ! __|¨|__ 1063 CASE( 13 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = ii ; ij3 = ij-1 1064 CASE( 14 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = ii ; ij3 = ij-1 1065 CASE( 15 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = ii ; ij2 = ij-1 ; ii3 = ii+1 ; ij3 = ij 1066 CASE( 16 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = ii ; ij2 = ij+1 ; ii3 = ii+1 ; ij3 = ij 1067 END SELECT 1068 END SUBROUTINE find_neib 1069 1070 1071 SUBROUTINE bdy_read_seg( kb_bdy, knblendta ) 1072 !!---------------------------------------------------------------------- 1073 !! *** ROUTINE bdy_coords_seg *** 1074 !! 1075 !! ** Purpose : build bdy coordinates with segments defined in namelist 1076 !! 1077 !! ** Method : read namelist nambdy_index blocks 1078 !! 1079 !!---------------------------------------------------------------------- 1080 INTEGER , INTENT (in ) :: kb_bdy ! bdy number 1081 INTEGER, DIMENSION(jpbgrd), INTENT ( out) :: knblendta ! length of index arrays 1082 !! 1083 INTEGER :: ios ! Local integer output status for namelist read 1084 INTEGER :: nbdyind, nbdybeg, nbdyend 1085 INTEGER :: nbdy_count, nbdy_rdstart, nbdy_loc 1086 CHARACTER(LEN=1) :: ctypebdy ! - - 1087 CHARACTER(LEN=50):: cerrmsg ! - - 1088 NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 1089 !!---------------------------------------------------------------------- 1090 ! Need to support possibility of reading more than one nambdy_index from 1091 ! the namelist_cfg internal file. 1092 ! Do this by finding the kb_bdy'th occurence of nambdy_index in the 1093 ! character buffer as the starting point. 1094 nbdy_rdstart = 1 1095 DO nbdy_count = 1, kb_bdy 1096 nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_index' ) 1097 IF( nbdy_loc .GT. 0 ) THEN 1098 nbdy_rdstart = nbdy_rdstart + nbdy_loc 1099 ELSE 1100 WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',kb_bdy,' of nambdy_index not found' 1101 ios = -1 1102 CALL ctl_nam ( ios , cerrmsg ) 1103 ENDIF 1104 END DO 1105 nbdy_rdstart = MAX( 1, nbdy_rdstart - 2 ) 1106 READ ( numnam_cfg( nbdy_rdstart: ), nambdy_index, IOSTAT = ios, ERR = 904) 1107 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_index in configuration namelist' ) 1108 IF(lwm) WRITE ( numond, nambdy_index ) 1109 1110 SELECT CASE ( TRIM(ctypebdy) ) 1111 CASE( 'N' ) 1112 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 1113 nbdyind = jpjglo - 2 ! set boundary to whole side of model domain. 1114 nbdybeg = 2 1115 nbdyend = jpiglo - 1 1116 ENDIF 1117 nbdysegn = nbdysegn + 1 1118 npckgn(nbdysegn) = kb_bdy ! Save bdy package number 1119 jpjnob(nbdysegn) = nbdyind 1120 jpindt(nbdysegn) = nbdybeg 1121 jpinft(nbdysegn) = nbdyend 1122 ! 1123 CASE( 'S' ) 1124 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 1125 nbdyind = 2 ! set boundary to whole side of model domain. 1126 nbdybeg = 2 1127 nbdyend = jpiglo - 1 1128 ENDIF 1129 nbdysegs = nbdysegs + 1 1130 npckgs(nbdysegs) = kb_bdy ! Save bdy package number 1131 jpjsob(nbdysegs) = nbdyind 1132 jpisdt(nbdysegs) = nbdybeg 1133 jpisft(nbdysegs) = nbdyend 1134 ! 1135 CASE( 'E' ) 1136 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 1137 nbdyind = jpiglo - 2 ! set boundary to whole side of model domain. 1138 nbdybeg = 2 1139 nbdyend = jpjglo - 1 1140 ENDIF 1141 nbdysege = nbdysege + 1 1142 npckge(nbdysege) = kb_bdy ! Save bdy package number 1143 jpieob(nbdysege) = nbdyind 1144 jpjedt(nbdysege) = nbdybeg 1145 jpjeft(nbdysege) = nbdyend 1146 ! 1147 CASE( 'W' ) 1148 IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1 1149 nbdyind = 2 ! set boundary to whole side of model domain. 1150 nbdybeg = 2 1151 nbdyend = jpjglo - 1 1152 ENDIF 1153 nbdysegw = nbdysegw + 1 1154 npckgw(nbdysegw) = kb_bdy ! Save bdy package number 1155 jpiwob(nbdysegw) = nbdyind 1156 jpjwdt(nbdysegw) = nbdybeg 1157 jpjwft(nbdysegw) = nbdyend 1158 ! 1159 CASE DEFAULT ; CALL ctl_stop( 'ctypebdy must be N, S, E or W' ) 1160 END SELECT 1161 1162 ! For simplicity we assume that in case of straight bdy, arrays have the same length 1163 ! (even if it is true that last tangential velocity points 1164 ! are useless). This simplifies a little bit boundary data format (and agrees with format 1165 ! used so far in obc package) 1166 1167 knblendta(1:jpbgrd) = (nbdyend - nbdybeg + 1) * nn_rimwidth(kb_bdy) 1168 1169 END SUBROUTINE bdy_read_seg 1170 1171 1252 1172 SUBROUTINE bdy_ctl_seg 1253 1173 !!---------------------------------------------------------------------- … … 1279 1199 &(jpjnob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' ) 1280 1200 IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 1281 IF (jpindt(ib).l e.1 ) CALL ctl_stop( 'Start index out of domain' )1282 IF (jpinft(ib).g e.jpiglo) CALL ctl_stop( 'End index out of domain' )1201 IF (jpindt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' ) 1202 IF (jpinft(ib).gt.jpiglo) CALL ctl_stop( 'End index out of domain' ) 1283 1203 END DO 1284 1204 ! … … 1288 1208 &(jpjsob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' ) 1289 1209 IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 1290 IF (jpisdt(ib).l e.1 ) CALL ctl_stop( 'Start index out of domain' )1291 IF (jpisft(ib).g e.jpiglo) CALL ctl_stop( 'End index out of domain' )1210 IF (jpisdt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' ) 1211 IF (jpisft(ib).gt.jpiglo) CALL ctl_stop( 'End index out of domain' ) 1292 1212 END DO 1293 1213 ! … … 1297 1217 &(jpieob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' ) 1298 1218 IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 1299 IF (jpjedt(ib).l e.1 ) CALL ctl_stop( 'Start index out of domain' )1300 IF (jpjeft(ib).g e.jpjglo) CALL ctl_stop( 'End index out of domain' )1219 IF (jpjedt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' ) 1220 IF (jpjeft(ib).gt.jpjglo) CALL ctl_stop( 'End index out of domain' ) 1301 1221 END DO 1302 1222 ! … … 1306 1226 &(jpiwob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' ) 1307 1227 IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 1308 IF (jpjwdt(ib).l e.1 ) CALL ctl_stop( 'Start index out of domain' )1309 IF (jpjwft(ib).g e.jpjglo) CALL ctl_stop( 'End index out of domain' )1228 IF (jpjwdt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' ) 1229 IF (jpjwft(ib).gt.jpjglo) CALL ctl_stop( 'End index out of domain' ) 1310 1230 ENDDO 1311 1231 ! … … 1336 1256 icorns(ib2,1) = npckgw(ib1) 1337 1257 ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN 1338 WRITE(ctmp1,*) ' E R R O R :Found an acute open boundary corner at point (i,j)= ', &1258 WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 1339 1259 & jpisft(ib2), jpjwft(ib1) 1340 WRITE(ctmp2,*) ' ==========Not allowed yet'1341 WRITE(ctmp3,*) ' 1342 & 1343 CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ')1260 WRITE(ctmp2,*) ' Not allowed yet' 1261 WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), & 1262 & ' and South segment: ',npckgs(ib2) 1263 CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 1344 1264 ELSE 1345 WRITE(ctmp1,*) ' E R R O R :Check South and West Open boundary indices'1346 WRITE(ctmp2,*) ' ==========Crossing problem with West segment: ',npckgw(ib1) , &1347 & 1348 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ')1265 WRITE(ctmp1,*) ' Check South and West Open boundary indices' 1266 WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1) , & 1267 & ' and South segment: ',npckgs(ib2) 1268 CALL ctl_stop( ctmp1, ctmp2 ) 1349 1269 END IF 1350 1270 END IF … … 1368 1288 icorns(ib2,2) = npckge(ib1) 1369 1289 ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN 1370 WRITE(ctmp1,*) ' E R R O R :Found an acute open boundary corner at point (i,j)= ', &1290 WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 1371 1291 & jpisdt(ib2), jpjeft(ib1) 1372 WRITE(ctmp2,*) ' ==========Not allowed yet'1373 WRITE(ctmp3,*) ' 1374 & 1375 CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ')1292 WRITE(ctmp2,*) ' Not allowed yet' 1293 WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), & 1294 & ' and South segment: ',npckgs(ib2) 1295 CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 1376 1296 ELSE 1377 WRITE(ctmp1,*) ' E R R O R :Check South and East Open boundary indices'1378 WRITE(ctmp2,*) ' ==========Crossing problem with East segment: ',npckge(ib1), &1379 & 1380 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ')1297 WRITE(ctmp1,*) ' Check South and East Open boundary indices' 1298 WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), & 1299 & ' and South segment: ',npckgs(ib2) 1300 CALL ctl_stop( ctmp1, ctmp2 ) 1381 1301 END IF 1382 1302 END IF … … 1400 1320 icornn(ib2,1) = npckgw(ib1) 1401 1321 ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN 1402 WRITE(ctmp1,*) ' E R R O R :Found an acute open boundary corner at point (i,j)= ', &1322 WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 1403 1323 & jpinft(ib2), jpjwdt(ib1) 1404 WRITE(ctmp2,*) ' ==========Not allowed yet'1405 WRITE(ctmp3,*) ' 1406 & 1407 CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ')1324 WRITE(ctmp2,*) ' Not allowed yet' 1325 WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), & 1326 & ' and North segment: ',npckgn(ib2) 1327 CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 1408 1328 ELSE 1409 WRITE(ctmp1,*) ' E R R O R :Check North and West Open boundary indices'1410 WRITE(ctmp2,*) ' ==========Crossing problem with West segment: ',npckgw(ib1), &1411 & 1412 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ')1329 WRITE(ctmp1,*) ' Check North and West Open boundary indices' 1330 WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1), & 1331 & ' and North segment: ',npckgn(ib2) 1332 CALL ctl_stop( ctmp1, ctmp2 ) 1413 1333 END IF 1414 1334 END IF … … 1432 1352 icornn(ib2,2) = npckge(ib1) 1433 1353 ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN 1434 WRITE(ctmp1,*) ' E R R O R :Found an acute open boundary corner at point (i,j)= ', &1354 WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 1435 1355 & jpindt(ib2), jpjedt(ib1) 1436 WRITE(ctmp2,*) ' ==========Not allowed yet'1437 WRITE(ctmp3,*) ' 1438 & 1439 CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ')1356 WRITE(ctmp2,*) ' Not allowed yet' 1357 WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), & 1358 & ' and North segment: ',npckgn(ib2) 1359 CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 1440 1360 ELSE 1441 WRITE(ctmp1,*) ' E R R O R :Check North and East Open boundary indices'1442 WRITE(ctmp2,*) ' ==========Crossing problem with East segment: ',npckge(ib1), &1443 & 1444 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ')1361 WRITE(ctmp1,*) ' Check North and East Open boundary indices' 1362 WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), & 1363 & ' and North segment: ',npckgn(ib2) 1364 CALL ctl_stop( ctmp1, ctmp2 ) 1445 1365 END IF 1446 1366 END IF … … 1458 1378 DO ji = 1, jpi 1459 1379 DO jj = 1, jpj 1460 IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 1461 & ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1) 1462 IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 1463 & ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1) 1380 IF( mig(ji) == jpiwob(ib) .AND. mjg(jj) == jpjwdt(ib) ) ztestmask(1) = tmask(ji,jj,1) 1381 IF( mig(ji) == jpiwob(ib) .AND. mjg(jj) == jpjwft(ib) ) ztestmask(2) = tmask(ji,jj,1) 1464 1382 END DO 1465 1383 END DO … … 1468 1386 IF (ztestmask(1)==1) THEN 1469 1387 IF (icornw(ib,1)==0) THEN 1470 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgw(ib) 1471 WRITE(ctmp2,*) ' ========== does not start on land or on a corner' 1472 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1388 WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib) 1389 CALL ctl_stop( ctmp1, ' does not start on land or on a corner' ) 1473 1390 ELSE 1474 1391 ! This is a corner … … 1480 1397 IF (ztestmask(2)==1) THEN 1481 1398 IF (icornw(ib,2)==0) THEN 1482 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgw(ib) 1483 WRITE(ctmp2,*) ' ========== does not end on land or on a corner' 1484 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1399 WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib) 1400 CALL ctl_stop( ' ', ctmp1, ' does not end on land or on a corner' ) 1485 1401 ELSE 1486 1402 ! This is a corner … … 1498 1414 DO ji = 1, jpi 1499 1415 DO jj = 1, jpj 1500 IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 1501 & ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1) 1502 IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 1503 & ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1) 1416 IF( mig(ji) == jpieob(ib)+1 .AND. mjg(jj) == jpjedt(ib) ) ztestmask(1) = tmask(ji,jj,1) 1417 IF( mig(ji) == jpieob(ib)+1 .AND. mjg(jj) == jpjeft(ib) ) ztestmask(2) = tmask(ji,jj,1) 1504 1418 END DO 1505 1419 END DO … … 1508 1422 IF (ztestmask(1)==1) THEN 1509 1423 IF (icorne(ib,1)==0) THEN 1510 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckge(ib) 1511 WRITE(ctmp2,*) ' ========== does not start on land or on a corner' 1512 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1424 WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib) 1425 CALL ctl_stop( ctmp1, ' does not start on land or on a corner' ) 1513 1426 ELSE 1514 1427 ! This is a corner … … 1520 1433 IF (ztestmask(2)==1) THEN 1521 1434 IF (icorne(ib,2)==0) THEN 1522 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckge(ib) 1523 WRITE(ctmp2,*) ' ========== does not end on land or on a corner' 1524 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1435 WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib) 1436 CALL ctl_stop( ctmp1, ' does not end on land or on a corner' ) 1525 1437 ELSE 1526 1438 ! This is a corner … … 1538 1450 DO ji = 1, jpi 1539 1451 DO jj = 1, jpj 1540 IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 1541 & ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1) 1542 IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 1543 & ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1) 1452 IF( mjg(jj) == jpjsob(ib) .AND. mig(ji) == jpisdt(ib) ) ztestmask(1) = tmask(ji,jj,1) 1453 IF( mjg(jj) == jpjsob(ib) .AND. mig(ji) == jpisft(ib) ) ztestmask(2) = tmask(ji,jj,1) 1544 1454 END DO 1545 1455 END DO … … 1547 1457 1548 1458 IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN 1549 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgs(ib) 1550 WRITE(ctmp2,*) ' ========== does not start on land or on a corner' 1551 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1459 WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib) 1460 CALL ctl_stop( ctmp1, ' does not start on land or on a corner' ) 1552 1461 ENDIF 1553 1462 IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN 1554 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgs(ib) 1555 WRITE(ctmp2,*) ' ========== does not end on land or on a corner' 1556 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1463 WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib) 1464 CALL ctl_stop( ctmp1, ' does not end on land or on a corner' ) 1557 1465 ENDIF 1558 1466 END DO … … 1564 1472 DO ji = 1, jpi 1565 1473 DO jj = 1, jpj 1566 IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 1567 & ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1) 1568 IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 1569 & ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1) 1474 IF( mjg(jj) == jpjnob(ib)+1 .AND. mig(ji) == jpindt(ib) ) ztestmask(1) = tmask(ji,jj,1) 1475 IF( mjg(jj) == jpjnob(ib)+1 .AND. mig(ji) == jpinft(ib) ) ztestmask(2) = tmask(ji,jj,1) 1570 1476 END DO 1571 1477 END DO … … 1573 1479 1574 1480 IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN 1575 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgn(ib) 1576 WRITE(ctmp2,*) ' ========== does not start on land' 1577 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1481 WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib) 1482 CALL ctl_stop( ctmp1, ' does not start on land' ) 1578 1483 ENDIF 1579 1484 IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN 1580 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgn(ib) 1581 WRITE(ctmp2,*) ' ========== does not end on land' 1582 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1485 WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib) 1486 CALL ctl_stop( ctmp1, ' does not end on land' ) 1583 1487 ENDIF 1584 1488 END DO … … 1593 1497 END SUBROUTINE bdy_ctl_seg 1594 1498 1595 1499 1500 SUBROUTINE bdy_coords_seg( nbidta, nbjdta, nbrdta ) 1501 !!---------------------------------------------------------------------- 1502 !! *** ROUTINE bdy_coords_seg *** 1503 !! 1504 !! ** Purpose : build nbidta, nbidta, nbrdta for bdy built with segments 1505 !! 1506 !! ** Method : 1507 !! 1508 !!---------------------------------------------------------------------- 1509 INTEGER, DIMENSION(:,:,:), intent( out) :: nbidta, nbjdta, nbrdta ! Index arrays: i and j indices of bdy dta 1510 !! 1511 INTEGER :: ii, ij, ir, iseg 1512 INTEGER :: igrd ! grid type (t=1, u=2, v=3) 1513 INTEGER :: icount ! 1514 INTEGER :: ib_bdy ! bdy number 1515 !!---------------------------------------------------------------------- 1516 1517 ! East 1518 !----- 1519 DO iseg = 1, nbdysege 1520 ib_bdy = npckge(iseg) 1521 ! 1522 ! ------------ T points ------------- 1523 igrd=1 1524 icount=0 1525 DO ir = 1, nn_rimwidth(ib_bdy) 1526 DO ij = jpjedt(iseg), jpjeft(iseg) 1527 icount = icount + 1 1528 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 1529 nbjdta(icount, igrd, ib_bdy) = ij 1530 nbrdta(icount, igrd, ib_bdy) = ir 1531 ENDDO 1532 ENDDO 1533 ! 1534 ! ------------ U points ------------- 1535 igrd=2 1536 icount=0 1537 DO ir = 1, nn_rimwidth(ib_bdy) 1538 DO ij = jpjedt(iseg), jpjeft(iseg) 1539 icount = icount + 1 1540 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir 1541 nbjdta(icount, igrd, ib_bdy) = ij 1542 nbrdta(icount, igrd, ib_bdy) = ir 1543 ENDDO 1544 ENDDO 1545 ! 1546 ! ------------ V points ------------- 1547 igrd=3 1548 icount=0 1549 DO ir = 1, nn_rimwidth(ib_bdy) 1550 ! DO ij = jpjedt(iseg), jpjeft(iseg) - 1 1551 DO ij = jpjedt(iseg), jpjeft(iseg) 1552 icount = icount + 1 1553 nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 1554 nbjdta(icount, igrd, ib_bdy) = ij 1555 nbrdta(icount, igrd, ib_bdy) = ir 1556 ENDDO 1557 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1558 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1559 ENDDO 1560 ENDDO 1561 ! 1562 ! West 1563 !----- 1564 DO iseg = 1, nbdysegw 1565 ib_bdy = npckgw(iseg) 1566 ! 1567 ! ------------ T points ------------- 1568 igrd=1 1569 icount=0 1570 DO ir = 1, nn_rimwidth(ib_bdy) 1571 DO ij = jpjwdt(iseg), jpjwft(iseg) 1572 icount = icount + 1 1573 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 1574 nbjdta(icount, igrd, ib_bdy) = ij 1575 nbrdta(icount, igrd, ib_bdy) = ir 1576 ENDDO 1577 ENDDO 1578 ! 1579 ! ------------ U points ------------- 1580 igrd=2 1581 icount=0 1582 DO ir = 1, nn_rimwidth(ib_bdy) 1583 DO ij = jpjwdt(iseg), jpjwft(iseg) 1584 icount = icount + 1 1585 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 1586 nbjdta(icount, igrd, ib_bdy) = ij 1587 nbrdta(icount, igrd, ib_bdy) = ir 1588 ENDDO 1589 ENDDO 1590 ! 1591 ! ------------ V points ------------- 1592 igrd=3 1593 icount=0 1594 DO ir = 1, nn_rimwidth(ib_bdy) 1595 ! DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 1596 DO ij = jpjwdt(iseg), jpjwft(iseg) 1597 icount = icount + 1 1598 nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 1599 nbjdta(icount, igrd, ib_bdy) = ij 1600 nbrdta(icount, igrd, ib_bdy) = ir 1601 ENDDO 1602 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1603 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1604 ENDDO 1605 ENDDO 1606 ! 1607 ! North 1608 !----- 1609 DO iseg = 1, nbdysegn 1610 ib_bdy = npckgn(iseg) 1611 ! 1612 ! ------------ T points ------------- 1613 igrd=1 1614 icount=0 1615 DO ir = 1, nn_rimwidth(ib_bdy) 1616 DO ii = jpindt(iseg), jpinft(iseg) 1617 icount = icount + 1 1618 nbidta(icount, igrd, ib_bdy) = ii 1619 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 1620 nbrdta(icount, igrd, ib_bdy) = ir 1621 ENDDO 1622 ENDDO 1623 ! 1624 ! ------------ U points ------------- 1625 igrd=2 1626 icount=0 1627 DO ir = 1, nn_rimwidth(ib_bdy) 1628 ! DO ii = jpindt(iseg), jpinft(iseg) - 1 1629 DO ii = jpindt(iseg), jpinft(iseg) 1630 icount = icount + 1 1631 nbidta(icount, igrd, ib_bdy) = ii 1632 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 1633 nbrdta(icount, igrd, ib_bdy) = ir 1634 ENDDO 1635 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1636 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1637 ENDDO 1638 ! 1639 ! ------------ V points ------------- 1640 igrd=3 1641 icount=0 1642 DO ir = 1, nn_rimwidth(ib_bdy) 1643 DO ii = jpindt(iseg), jpinft(iseg) 1644 icount = icount + 1 1645 nbidta(icount, igrd, ib_bdy) = ii 1646 nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir 1647 nbrdta(icount, igrd, ib_bdy) = ir 1648 ENDDO 1649 ENDDO 1650 ENDDO 1651 ! 1652 ! South 1653 !----- 1654 DO iseg = 1, nbdysegs 1655 ib_bdy = npckgs(iseg) 1656 ! 1657 ! ------------ T points ------------- 1658 igrd=1 1659 icount=0 1660 DO ir = 1, nn_rimwidth(ib_bdy) 1661 DO ii = jpisdt(iseg), jpisft(iseg) 1662 icount = icount + 1 1663 nbidta(icount, igrd, ib_bdy) = ii 1664 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 1665 nbrdta(icount, igrd, ib_bdy) = ir 1666 ENDDO 1667 ENDDO 1668 ! 1669 ! ------------ U points ------------- 1670 igrd=2 1671 icount=0 1672 DO ir = 1, nn_rimwidth(ib_bdy) 1673 ! DO ii = jpisdt(iseg), jpisft(iseg) - 1 1674 DO ii = jpisdt(iseg), jpisft(iseg) 1675 icount = icount + 1 1676 nbidta(icount, igrd, ib_bdy) = ii 1677 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 1678 nbrdta(icount, igrd, ib_bdy) = ir 1679 ENDDO 1680 nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1681 nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 1682 ENDDO 1683 ! 1684 ! ------------ V points ------------- 1685 igrd=3 1686 icount=0 1687 DO ir = 1, nn_rimwidth(ib_bdy) 1688 DO ii = jpisdt(iseg), jpisft(iseg) 1689 icount = icount + 1 1690 nbidta(icount, igrd, ib_bdy) = ii 1691 nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 1692 nbrdta(icount, igrd, ib_bdy) = ir 1693 ENDDO 1694 ENDDO 1695 ENDDO 1696 1697 1698 END SUBROUTINE bdy_coords_seg 1699 1700 1596 1701 SUBROUTINE bdy_ctl_corn( ib1, ib2 ) 1597 1702 !!---------------------------------------------------------------------- … … 1619 1724 ! 1620 1725 IF( itest>0 ) THEN 1621 WRITE(ctmp1,*) ' E R R O R : Segments ', ib1, 'and ', ib2 1622 WRITE(ctmp2,*) ' ========== have different open bdy schemes' 1623 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1726 WRITE(ctmp1,*) ' Segments ', ib1, 'and ', ib2 1727 CALL ctl_stop( ctmp1, ' have different open bdy schemes' ) 1624 1728 ENDIF 1625 1729 ! 1626 1730 END SUBROUTINE bdy_ctl_corn 1627 1731 1732 1733 SUBROUTINE bdy_meshwri() 1734 !!---------------------------------------------------------------------- 1735 !! *** ROUTINE bdy_meshwri *** 1736 !! 1737 !! ** Purpose : write netcdf file with nbr, flagu, flagv, ntreat for T, U 1738 !! and V points in 2D arrays for easier visualisation/control 1739 !! 1740 !! ** Method : use iom_rstput as in domwri.F 1741 !!---------------------------------------------------------------------- 1742 INTEGER :: ib_bdy, ii, ij, igrd, ib ! dummy loop indices 1743 INTEGER :: inum ! - - 1744 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! pointer to 2D mask fields 1745 REAL(wp) , DIMENSION(jpi,jpj) :: ztmp 1746 CHARACTER(LEN=1) , DIMENSION(jpbgrd) :: cgrid 1747 !!---------------------------------------------------------------------- 1748 cgrid = (/'t','u','v'/) 1749 CALL iom_open( 'bdy_mesh', inum, ldwrt = .TRUE. ) 1750 DO igrd = 1, jpbgrd 1751 SELECT CASE( igrd ) 1752 CASE( 1 ) ; zmask => tmask(:,:,1) 1753 CASE( 2 ) ; zmask => umask(:,:,1) 1754 CASE( 3 ) ; zmask => vmask(:,:,1) 1755 END SELECT 1756 ztmp(:,:) = zmask(:,:) 1757 DO ib_bdy = 1, nb_bdy 1758 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) ! nbr deined for all rims 1759 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1760 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1761 ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%nbr(ib,igrd), wp) + 10. 1762 IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 1763 END DO 1764 END DO 1765 CALL iom_rstput( 0, 0, inum, 'bdy_nbr_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 1766 ztmp(:,:) = zmask(:,:) 1767 DO ib_bdy = 1, nb_bdy 1768 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! flagu defined only for rims 0 and 1 1769 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1770 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1771 ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagu(ib,igrd), wp) + 10. 1772 IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 1773 END DO 1774 END DO 1775 CALL iom_rstput( 0, 0, inum, 'flagu_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 1776 ztmp(:,:) = zmask(:,:) 1777 DO ib_bdy = 1, nb_bdy 1778 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! flagv defined only for rims 0 and 1 1779 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1780 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1781 ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagv(ib,igrd), wp) + 10. 1782 IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 1783 END DO 1784 END DO 1785 CALL iom_rstput( 0, 0, inum, 'flagv_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 1786 ztmp(:,:) = zmask(:,:) 1787 DO ib_bdy = 1, nb_bdy 1788 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! ntreat defined only for rims 0 and 1 1789 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1790 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1791 ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%ntreat(ib,igrd), wp) + 10. 1792 IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 1793 END DO 1794 END DO 1795 CALL iom_rstput( 0, 0, inum, 'ntreat_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 1796 END DO 1797 CALL iom_close( inum ) 1798 1799 END SUBROUTINE bdy_meshwri 1800 1628 1801 !!================================================================================= 1629 1802 END MODULE bdyini
Note: See TracChangeset
for help on using the changeset viewer.
