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