- Timestamp:
- 2011-07-11T12:53:56+02:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90
r2715 r2797 1 1 MODULE obcini 2 2 !!====================================================================== 3 3 !! *** MODULE obcini *** 4 !! OBC initial state : Open boundary initial state4 !! Unstructured open boundaries : initialisation 5 5 !!====================================================================== 6 !! History : 8.0 ! 97-07 (J.M. Molines, G. Madec) Original code 7 !! NEMO 1.0 ! 02-11 (C. Talandier, A-M. Treguier) Free surface, F90 8 !! 2.0 ! 05-11 (V. Garnier) Surface pressure gradient organization 6 !! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code 7 !! - ! 2007-01 (D. Storkey) Update to use IOM module 8 !! - ! 2007-01 (D. Storkey) Tidal forcing 9 !! 3.0 ! 2008-04 (NEMO team) add in the reference version 10 !! 3.3 ! 2010-09 (E.O'Dea) updates for Shelf configurations 11 !! 3.3 ! 2010-09 (D.Storkey) add ice boundary conditions 12 !! 3.4 ! 2011 (D. Storkey, J. Chanut) OBC-BDY merge 13 !! ! --- Renamed bdyini.F90 -> obcini.F90 --- 9 14 !!---------------------------------------------------------------------- 10 15 #if defined key_obc 11 16 !!---------------------------------------------------------------------- 12 !! 'key_obc' 17 !! 'key_obc' Unstructured Open Boundary Conditions 13 18 !!---------------------------------------------------------------------- 14 !! obc_init : initialization for the open boundary condition19 !! obc_init : Initialization of unstructured open boundaries 15 20 !!---------------------------------------------------------------------- 16 21 USE oce ! ocean dynamics and tracers variables 17 USE dom_oce ! ocean space and time domain variables 22 USE dom_oce ! ocean space and time domain 23 USE obc_oce ! unstructured open boundary conditions 24 USE obctides ! tides at open boundaries initialization (tide_init routine) 25 USE in_out_manager ! I/O units 18 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 19 USE phycst ! physical constants 20 USE obc_oce ! open boundary condition: ocean 21 USE obcdta ! open boundary condition: data 22 USE in_out_manager ! I/O units 23 USE lib_mpp ! MPP library 24 USE dynspg_oce ! flag lk_dynspg_flt 27 USE lib_mpp ! for mpp_sum 28 USE iom ! I/O 25 29 26 30 IMPLICIT NONE … … 29 33 PUBLIC obc_init ! routine called by opa.F90 30 34 31 !! * Substitutions32 # include "obc_vectopt_loop_substitute.h90"33 35 !!---------------------------------------------------------------------- 34 !! NEMO/OPA 3.3 , NEMO Consortium (2010)36 !! NEMO/OPA 4.0 , NEMO Consortium (2011) 35 37 !! $Id$ 36 38 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 42 44 !! *** ROUTINE obc_init *** 43 45 !! 44 !! ** Purpose : Initialization of the dynamics and tracer fields at45 !! theopen boundaries.46 !! ** Purpose : Initialization of the dynamics and tracer fields with 47 !! unstructured open boundaries. 46 48 !! 47 !! ** Method : initialization of open boundary variables 48 !! (u, v) over 3 time step and 3 rows 49 !! (t, s) over 2 time step and 2 rows 50 !! if ln_rstart = .FALSE. : no restart, fields set to zero 51 !! if ln_rstart = .TRUE. : restart, fields are read in a file 52 !! if rdpxxx = 0 then lfbc is set true for this boundary. 49 !! ** Method : Read initialization arrays (mask, indices) to identify 50 !! an unstructured open boundary 53 51 !! 54 !! ** Input : restart.obc file, restart file for open boundaries 52 !! ** Input : obc_init.nc, input file for unstructured open boundaries 53 !!---------------------------------------------------------------------- 54 INTEGER :: ib_obc, ii, ij, ik, igrd, ib, ir ! dummy loop indices 55 INTEGER :: icount, icountr, ibr_max, ilen1 ! local integers 56 INTEGER :: iw, ie, is, in, inum, id_dummy ! - - 57 INTEGER :: igrd_start, igrd_end, jpbdta ! - - 58 INTEGER, POINTER :: nbi, nbj, nbr ! short cuts 59 REAL , POINTER :: flagu, flagv ! - - 60 REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars 61 INTEGER, DIMENSION (2) :: kdimsz 62 INTEGER, DIMENSION(jpbgrd,jp_obc) :: nblendta ! Length of index arrays 63 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of obc dta 64 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points 65 REAL(wp), DIMENSION(jpidta,jpjdta) :: zmask ! global domain mask 66 CHARACTER(LEN=80),DIMENSION(jpbgrd) :: clfile 67 CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid 68 !! 69 NAMELIST/namobc/ nb_obc, ln_coords_file, cn_coords_file, & 70 & ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn3d, & 71 & nn_tra, & 72 #if defined key_lim2 73 & nn_ice_lim2, & 74 #endif 75 & ln_tides, ln_vol, ln_clim, nn_dtactl, nn_volctl, & 76 & nn_rimwidth, nn_dmp2d_in, nn_dmp2d_out, & 77 & nn_dmp3d_in, nn_dmp3d_out 55 78 !!---------------------------------------------------------------------- 56 USE obcrst, ONLY : obc_rst_read ! Make obc_rst_read routine available 57 !! 58 INTEGER :: ji, jj, istop , inumfbc 59 INTEGER, DIMENSION(4) :: icorner 60 REAL(wp), DIMENSION(2) :: ztestmask 61 !! 62 NAMELIST/namobc/ rn_dpein, rn_dpwin, rn_dpnin, rn_dpsin, & 63 & rn_dpeob, rn_dpwob, rn_dpnob, rn_dpsob, & 64 & rn_volemp, nn_obcdta, cn_obcdta, & 65 & ln_obc_clim, ln_vol_cst, ln_obc_fla 66 !!---------------------------------------------------------------------- 67 68 REWIND( numnam ) ! Namelist namobc : open boundaries 69 READ ( numnam, namobc ) 70 71 ! convert DOCTOR namelist name into the OLD names 72 nobc_dta = nn_obcdta 73 cffile = cn_obcdta 74 rdpein = rn_dpein 75 rdpwin = rn_dpwin 76 rdpsin = rn_dpsin 77 rdpnin = rn_dpnin 78 rdpeob = rn_dpeob 79 rdpwob = rn_dpwob 80 rdpsob = rn_dpsob 81 rdpnob = rn_dpnob 82 volemp = rn_volemp 83 84 ! ! allocate obc arrays 85 IF( obc_oce_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'obc_init : unable to allocate obc_oce arrays' ) 86 IF( obc_dta_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'obc_init : unable to allocate obc_dta arrays' ) 87 88 ! By security we set rdpxin and rdpxob respectively to 1. and 15. if the corresponding OBC is not activated 89 IF( .NOT.lp_obc_east ) THEN ; rdpein = 1. ; rdpeob = 15. ; END IF 90 IF( .NOT.lp_obc_west ) THEN ; rdpwin = 1. ; rdpwob = 15. ; END IF 91 IF( .NOT.lp_obc_north ) THEN ; rdpnin = 1. ; rdpnob = 15. ; END IF 92 IF( .NOT.lp_obc_south ) THEN ; rdpsin = 1. ; rdpsob = 15. ; END IF 93 94 ! number of open boudaries and open boundary indicators 95 nbobc = 0 96 IF( lp_obc_east ) nbobc = nbobc + 1 97 IF( lp_obc_west ) nbobc = nbobc + 1 98 IF( lp_obc_north ) nbobc = nbobc + 1 99 IF( lp_obc_south ) nbobc = nbobc + 1 79 80 IF( obc_oce_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'obc_init : unable to allocate oce arrays' ) 100 81 101 82 IF(lwp) WRITE(numout,*) 102 83 IF(lwp) WRITE(numout,*) 'obc_init : initialization of open boundaries' 103 84 IF(lwp) WRITE(numout,*) '~~~~~~~~' 104 IF(lwp) WRITE(numout,*) ' Number of open boundaries nbobc = ', nbobc 105 IF(lwp) WRITE(numout,*) 106 107 ! control prints 108 IF(lwp) WRITE(numout,*) ' Namelist namobc' 109 IF(lwp) WRITE(numout,*) ' data in file (=1) or initial state used (=0) nn_obcdta = ', nn_obcdta 110 IF(lwp) WRITE(numout,*) ' climatology (true) or not ln_obc_clim = ', ln_obc_clim 111 IF(lwp) WRITE(numout,*) ' vol_cst (true) or not: ln_vol_cst = ', ln_vol_cst 112 IF(lwp) WRITE(numout,*) ' ' 113 IF(lwp) WRITE(numout,*) ' WARNING ' 114 IF(lwp) WRITE(numout,*) ' Flather"s algorithm is applied with explicit free surface scheme ' 115 IF(lwp) WRITE(numout,*) ' or with free surface time-splitting scheme ' 116 IF(lwp) WRITE(numout,*) ' Nor radiation neither relaxation is allowed with explicit free surface scheme: ' 117 IF(lwp) WRITE(numout,*) ' Radiation and/or relaxation is allowed with free surface time-splitting scheme ' 118 IF(lwp) WRITE(numout,*) ' depending of the choice of rdpXin = rdpXob = 0. for open boundaries ' 119 IF(lwp) WRITE(numout,*) 120 IF(lwp) WRITE(numout,*) ' For the filtered free surface case, ' 121 IF(lwp) WRITE(numout,*) ' radiation, relaxation or presciption of data can be applied ' 122 IF(lwp) WRITE(numout,*) 123 124 IF( lwp.AND.lp_obc_east ) THEN 125 WRITE(numout,*) ' East open boundary :' 126 WRITE(numout,*) ' i index jpieob = ', jpieob 127 WRITE(numout,*) ' damping time scale (days) rn_dpeob = ', rn_dpeob 128 WRITE(numout,*) ' damping time scale (days) rn_dpein = ', rn_dpein 85 ! 86 87 IF( jperio /= 0 ) CALL ctl_stop( 'Cyclic or symmetric,', & 88 & ' and general open boundary condition are not compatible' ) 89 90 cgrid= (/'T','U','V'/) 91 92 ! ----------------------------------------- 93 ! Initialise and read namelist parameters 94 ! ----------------------------------------- 95 96 nb_obc = 0 97 ln_coords_file(:) = .false. 98 cn_coords_file(:) = '' 99 ln_mask_file = .false. 100 cn_mask_file(:) = '' 101 nn_dyn2d(:) = 0 102 nn_dyn3d(:) = 0 103 nn_tra(:) = 0 104 #if defined key_lim2 105 nn_ice_lim2(:) = 0 106 #endif 107 ln_tides(:) = .false. 108 ln_vol = .false. 109 ln_clim(:) = .false. 110 nn_dtactl(:) = -1 ! uninitialised flag 111 nn_volctl = -1 ! uninitialised flag 112 nn_rimwidth(:) = -1 ! uninitialised flag 113 nn_dmp2d_in(:) = -1 ! uninitialised flag 114 nn_dmp2d_out(:) = -1 ! uninitialised flag 115 nn_dmp3d_in(:) = -1 ! uninitialised flag 116 nn_dmp3d_out(:) = -1 ! uninitialised flag 117 118 REWIND( numnam ) 119 READ ( numnam, namobc ) 120 121 ! ----------------------------------------- 122 ! Check and write out namelist parameters 123 ! ----------------------------------------- 124 125 ! ! control prints 126 IF(lwp) WRITE(numout,*) ' namobc' 127 128 IF( nb_obc .eq. 0 ) THEN 129 IF(lwp) WRITE(numout,*) 'nb_obc = 0, NO OPEN BOUNDARIES APPLIED.' 130 ELSE 131 IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_obc 129 132 ENDIF 130 133 131 IF( lwp.AND.lp_obc_west ) THEN 132 WRITE(numout,*) ' West open boundary :' 133 WRITE(numout,*) ' i index jpiwob = ', jpiwob 134 WRITE(numout,*) ' damping time scale (days) rn_dpwob = ', rn_dpwob 135 WRITE(numout,*) ' damping time scale (days) rn_dpwin = ', rn_dpwin 134 DO ib_obc = 1,nb_obc 135 IF(lwp) WRITE(numout,*) ' ' 136 IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_obc,'------' 137 138 ! ! check type of data used (nn_dtactl value) 139 SELECT CASE( nn_dtactl(ib_obc) ) ! 140 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for obc data' 141 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file' 142 CASE DEFAULT ; CALL ctl_stop( 'nn_dtactl must be 0 or 1' ) 143 END SELECT 144 IF(lwp) WRITE(numout,*) 145 146 IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution: ' 147 SELECT CASE( nn_dyn2d(ib_obc) ) 148 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 149 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 150 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Flather radiation condition' 151 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) 152 END SELECT 153 IF(lwp) WRITE(numout,*) 154 155 IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities: ' 156 SELECT CASE( nn_dyn3d(ib_obc) ) 157 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 158 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 159 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn3d' ) 160 END SELECT 161 IF(lwp) WRITE(numout,*) 162 163 IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity: ' 164 SELECT CASE( nn_tra(ib_obc) ) 165 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 166 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 167 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_tra' ) 168 END SELECT 169 IF(lwp) WRITE(numout,*) 170 171 #if defined key_lim2 172 IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice: ' 173 SELECT CASE( nn_tra(ib_obc) ) 174 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 175 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 176 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_tra' ) 177 END SELECT 178 IF(lwp) WRITE(numout,*) 179 #endif 180 181 IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nn_rimwidth = ', nn_rimwidth 182 IF(lwp) WRITE(numout,*) 183 184 IF( ln_tides(ib_obc) ) THEN 185 IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing at unstructured open boundaries' 186 IF(lwp) WRITE(numout,*) 187 ENDIF 188 189 !!$ ! Presumably will need to read in a separate namelist for each boundary that includes tides??? 190 !!$ IF( ln_tides ) CALL tide_init( ib_obc ) ! Read tides namelist 191 192 193 ENDDO 194 195 IF( ln_vol ) THEN ! check volume conservation (nn_volctl value) 196 IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 197 IF(lwp) WRITE(numout,*) 198 SELECT CASE ( nn_volctl ) 199 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' The total volume will be constant' 200 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' The total volume will vary according to the surface E-P flux' 201 CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 202 END SELECT 203 IF(lwp) WRITE(numout,*) 204 ELSE 205 IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 206 IF(lwp) WRITE(numout,*) 207 ENDIF 208 209 ! ------------------------------------------------- 210 ! Initialise indices arrays for open boundaries 211 ! ------------------------------------------------- 212 213 ! Work out global dimensions of boundary data 214 ! --------------------------------------------- 215 DO ib_obc = 1, nb_obc 216 217 jpbdta = 1 218 IF( .NOT. ln_coords_file(ib_obc) ) THEN ! Work out size of global arrays from namelist parameters 219 220 221 !! 1. Read parameters from namelist 222 !! 2. Work out global size of boundary data arrays nblendta and jpbdta 223 224 225 ELSE ! Read size of arrays in boundary coordinates file. 226 227 228 CALL iom_open( cn_coords_file(ib_obc), inum ) 229 jpbdta = 1 230 DO igrd = 1, jpbgrd 231 id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 232 nblendta(igrd,ib_obc) = kdimsz(1) 233 jpbdta = MAX(jpbdta, kdimsz(1)) 234 ENDDO 235 236 ENDIF 237 238 ENDDO 239 240 ! Allocate arrays 241 !--------------- 242 ALLOCATE( nbidta(jpbdta, jpbgrd, nb_obc), nbjdta(jpbdta, jpbgrd, nb_obc), & 243 & nbrdta(jpbdta, jpbgrd, nb_obc) ) 244 245 ALLOCATE( dta_global(jpbdta, 1, jpk) ) 246 247 ! Calculate global boundary index arrays or read in from file 248 !------------------------------------------------------------ 249 DO ib_obc = 1, nb_obc 250 251 IF( .NOT. ln_coords_file(ib_obc) ) THEN ! Calculate global index arrays from namelist parameters 252 253 !! Calculate global index arrays from namelist parameters 254 255 ELSE ! Read global index arrays from boundary coordinates file. 256 257 DO igrd = 1, jpbgrd 258 CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_obc),:,1) ) 259 DO ii = 1,nblendta(igrd,ib_obc) 260 nbidta(ii,igrd,ib_obc) = INT( dta_global(ii,1,1) ) 261 END DO 262 CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_obc),:,1) ) 263 DO ii = 1,nblendta(igrd,ib_obc) 264 nbjdta(ii,igrd,ib_obc) = INT( dta_global(ii,1,1) ) 265 END DO 266 CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_obc),:,1) ) 267 DO ii = 1,nblendta(igrd,ib_obc) 268 nbrdta(ii,igrd,ib_obc) = INT( dta_global(ii,1,1) ) 269 END DO 270 271 ibr_max = MAXVAL( nbrdta(:,igrd,ib_obc) ) 272 IF(lwp) WRITE(numout,*) 273 IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max 274 IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_obc) 275 IF (ibr_max < nn_rimwidth(ib_obc)) & 276 CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_obc) ) 277 278 END DO 279 CALL iom_close( inum ) 280 281 ENDIF 282 283 ENDDO 284 285 ! Work out dimensions of boundary data on each processor 286 ! ------------------------------------------------------ 287 288 iw = mig(1) + 1 ! if monotasking and no zoom, iw=2 289 ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1 290 is = mjg(1) + 1 ! if monotasking and no zoom, is=2 291 in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1 292 293 DO ib_obc = 1, nb_obc 294 DO igrd = 1, jpbgrd 295 icount = 0 296 icountr = 0 297 idx_obc(ib_obc)%nblen(igrd) = 0 298 idx_obc(ib_obc)%nblenrim(igrd) = 0 299 DO ib = 1, nblendta(igrd,ib_obc) 300 ! check if point is in local domain 301 IF( nbidta(ib,igrd,ib_obc) >= iw .AND. nbidta(ib,igrd,ib_obc) <= ie .AND. & 302 & nbjdta(ib,igrd,ib_obc) >= is .AND. nbjdta(ib,igrd,ib_obc) <= in ) THEN 303 ! 304 icount = icount + 1 305 ! 306 IF( nbrdta(ib,igrd,ib_obc) == 1 ) icountr = icountr+1 307 ENDIF 308 ENDDO 309 idx_obc(ib_obc)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc 310 idx_obc(ib_obc)%nblen (igrd) = icount !: length of boundary data on each proc 311 ENDDO ! igrd 312 313 ! Allocate index arrays for this boundary set 314 !-------------------------------------------- 315 ilen1 = MAXVAL(idx_obc(ib_obc)%nblen(:)) 316 ALLOCATE( idx_obc(ib_obc)%nbi(ilen1,jpbgrd) ) 317 ALLOCATE( idx_obc(ib_obc)%nbj(ilen1,jpbgrd) ) 318 ALLOCATE( idx_obc(ib_obc)%nbr(ilen1,jpbgrd) ) 319 ALLOCATE( idx_obc(ib_obc)%nbmap(ilen1,jpbgrd) ) 320 ALLOCATE( idx_obc(ib_obc)%nbw(ilen1,jpbgrd) ) 321 ALLOCATE( idx_obc(ib_obc)%flagu(ilen1) ) 322 ALLOCATE( idx_obc(ib_obc)%flagv(ilen1) ) 323 324 ! Dispatch mapping indices and discrete distances on each processor 325 ! ----------------------------------------------------------------- 326 327 DO igrd = 1, jpbgrd 328 icount = 0 329 ! Loop on rimwidth to ensure outermost points come first in the local arrays. 330 DO ir=1, nn_rimwidth(ib_obc) 331 DO ib = 1, nblendta(igrd,ib_obc) 332 ! check if point is in local domain and equals ir 333 IF( nbidta(ib,igrd,ib_obc) >= iw .AND. nbidta(ib,igrd,ib_obc) <= ie .AND. & 334 & nbjdta(ib,igrd,ib_obc) >= is .AND. nbjdta(ib,igrd,ib_obc) <= in .AND. & 335 & nbrdta(ib,igrd,ib_obc) == ir ) THEN 336 ! 337 icount = icount + 1 338 idx_obc(ib_obc)%nbi(icount,igrd) = nbidta(ib,igrd,ib_obc)- mig(1)+1 339 idx_obc(ib_obc)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_obc)- mjg(1)+1 340 idx_obc(ib_obc)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_obc) 341 idx_obc(ib_obc)%nbmap(icount,igrd) = ib 342 ENDIF 343 ENDDO 344 ENDDO 345 ENDDO 346 347 ! Compute rim weights 348 ! ------------------- 349 DO igrd = 1, jpbgrd 350 DO ib = 1, idx_obc(ib_obc)%nblen(igrd) 351 nbr => idx_obc(ib_obc)%nbr(ib,igrd) 352 idx_obc(ib_obc)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ! tanh formulation 353 ! idx_obc(ib_obc)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2 ! quadratic 354 ! idx_obc(ib_obc)%nbw(ib,igrd) = FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth) ! linear 355 END DO 356 END DO 357 358 ENDDO 359 360 ! ------------------------------------------------------ 361 ! Initialise masks and find normal/tangential directions 362 ! ------------------------------------------------------ 363 364 ! Read global 2D mask at T-points: obctmask 365 ! ----------------------------------------- 366 ! obctmask = 1 on the computational domain AND on open boundaries 367 ! = 0 elsewhere 368 369 IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN ! EEL configuration at 5km resolution 370 zmask( : ,:) = 0.e0 371 zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0 372 ELSE IF( ln_mask_file ) THEN 373 CALL iom_open( cn_mask_file, inum ) 374 CALL iom_get ( inum, jpdom_data, 'obc_msk', zmask(:,:) ) 375 CALL iom_close( inum ) 376 ELSE 377 zmask(:,:) = 1.e0 136 378 ENDIF 137 379 138 IF( lwp.AND.lp_obc_north ) THEN 139 WRITE(numout,*) ' North open boundary :' 140 WRITE(numout,*) ' j index jpjnob = ', jpjnob 141 WRITE(numout,*) ' damping time scale (days) rn_dpnob = ', rn_dpnob 142 WRITE(numout,*) ' damping time scale (days) rn_dpnin = ', rn_dpnin 143 ENDIF 144 145 IF( lwp.AND.lp_obc_south ) THEN 146 WRITE(numout,*) ' South open boundary :' 147 WRITE(numout,*) ' j index jpjsob = ', jpjsob 148 WRITE(numout,*) ' damping time scale (days) rn_dpsob = ', rn_dpsob 149 WRITE(numout,*) ' damping time scale (days) rn_dpsin = ', rn_dpsin 150 WRITE(numout,*) 151 ENDIF 152 153 IF( nbobc >= 2 .AND. jperio /= 0 ) & 154 & CALL ctl_stop( ' Cyclic or symmetric, and open boundary condition are not compatible' ) 155 156 ! 1. Initialisation of constants 157 ! ------------------------------ 158 ! ... convert rdp$ob in seconds 159 ! Fixed Bdy flag inbound outbound 160 lfbceast = .FALSE. ; rdpein = rdpein * rday ; rdpeob = rdpeob * rday 161 lfbcwest = .FALSE. ; rdpwin = rdpwin * rday ; rdpwob = rdpwob * rday 162 lfbcnorth = .FALSE. ; rdpnin = rdpnin * rday ; rdpnob = rdpnob * rday 163 lfbcsouth = .FALSE. ; rdpsin = rdpsin * rday ; rdpsob = rdpsob * rday 164 inumfbc = 0 165 ! ... look for Fixed Boundaries (rdp = 0 ) 166 ! ... When specified, lbcxxx flags are set to TRUE and rdpxxx are set to 167 ! ... a small arbitrary value, (to avoid division by zero further on). 168 ! ... rdpxxx is not used anymore. 169 IF( lp_obc_east ) THEN 170 IF( (rdpein+rdpeob) == 0 ) THEN 171 lfbceast = .TRUE. ; rdpein = 1e-3 ; rdpeob = 1e-3 172 inumfbc = inumfbc+1 173 ELSEIF ( (rdpein*rdpeob) == 0 ) THEN 174 CALL ctl_stop( 'obc_init : rn_dpein & rn_dpeob must be both zero or non zero' ) 175 END IF 176 END IF 177 178 IF( lp_obc_west ) THEN 179 IF( (rdpwin + rdpwob) == 0 ) THEN 180 lfbcwest = .TRUE. ; rdpwin = 1e-3 ; rdpwob = 1e-3 181 inumfbc = inumfbc+1 182 ELSEIF ( (rdpwin*rdpwob) == 0 ) THEN 183 CALL ctl_stop( 'obc_init : rn_dpwin & rn_dpwob must be both zero or non zero' ) 184 END IF 185 END IF 186 IF( lp_obc_north ) THEN 187 IF( (rdpnin + rdpnob) == 0 ) THEN 188 lfbcnorth = .TRUE. ; rdpnin = 1e-3 ; rdpnob = 1e-3 189 inumfbc = inumfbc+1 190 ELSEIF ( (rdpnin*rdpnob) == 0 ) THEN 191 CALL ctl_stop( 'obc_init : rn_dpnin & rn_dpnob must be both zero or non zero' ) 192 END IF 193 END IF 194 IF( lp_obc_south ) THEN 195 IF( (rdpsin + rdpsob) == 0 ) THEN 196 lfbcsouth = .TRUE. ; rdpsin = 1e-3 ; rdpsob = 1e-3 197 inumfbc = inumfbc+1 198 ELSEIF ( (rdpsin*rdpsob) == 0 ) THEN 199 CALL ctl_stop( 'obc_init : rn_dpsin & rn_dpsob must be both zero or non zero' ) 200 END IF 201 END IF 202 203 ! 2. Clever mpp indices for loops on the open boundaries. 204 ! The loops will be performed only on the processors 205 ! that contain a given open boundary. 206 ! -------------------------------------------------------- 207 208 IF( lp_obc_east ) THEN 209 ! ... mpp initialization 210 nie0 = max( 1, min(jpieob - nimpp+1, jpi ) ) 211 nie1 = max( 0, min(jpieob - nimpp+1, jpi - 1 ) ) 212 nie0p1 = max( 1, min(jpieob+1 - nimpp+1, jpi ) ) 213 nie1p1 = max( 0, min(jpieob+1 - nimpp+1, jpi - 1 ) ) 214 nie0m1 = max( 1, min(jpieob-1 - nimpp+1, jpi ) ) 215 nie1m1 = max( 0, min(jpieob-1 - nimpp+1, jpi - 1 ) ) 216 nje0 = max( 2, min(jpjed - njmpp+1, jpj ) ) 217 nje1 = max( 0, min(jpjef - njmpp+1, jpj - 1 ) ) 218 nje0p1 = max( 1, min(jpjedp1 - njmpp+1, jpj ) ) 219 nje0m1 = max( 1, min(jpjed - njmpp+1, jpj ) ) 220 nje1m1 = max( 0, min(jpjefm1 - njmpp+1, jpj - 1 ) ) 221 nje1m2 = max( 0, min(jpjefm1-1- njmpp+1, jpj - 1 ) ) 222 IF(lwp) THEN 223 IF( lfbceast ) THEN 224 WRITE(numout,*)' ' 225 WRITE(numout,*)' Specified East Open Boundary' 380 DO ij = 1, nlcj ! Save mask over local domain 381 DO ii = 1, nlci 382 obctmask(ii,ij) = zmask( mig(ii), mjg(ij) ) 383 END DO 384 END DO 385 386 ! Derive mask on U and V grid from mask on T grid 387 obcumask(:,:) = 0.e0 388 obcvmask(:,:) = 0.e0 389 DO ij=1, jpjm1 390 DO ii=1, jpim1 391 obcumask(ii,ij)=obctmask(ii,ij)*obctmask(ii+1, ij ) 392 obcvmask(ii,ij)=obctmask(ii,ij)*obctmask(ii ,ij+1) 393 END DO 394 END DO 395 CALL lbc_lnk( obcumask(:,:), 'U', 1. ) ; CALL lbc_lnk( obcvmask(:,:), 'V', 1. ) ! Lateral boundary cond. 396 397 398 ! Mask corrections 399 ! ---------------- 400 DO ik = 1, jpkm1 401 DO ij = 1, jpj 402 DO ii = 1, jpi 403 tmask(ii,ij,ik) = tmask(ii,ij,ik) * obctmask(ii,ij) 404 umask(ii,ij,ik) = umask(ii,ij,ik) * obcumask(ii,ij) 405 vmask(ii,ij,ik) = vmask(ii,ij,ik) * obcvmask(ii,ij) 406 bmask(ii,ij) = bmask(ii,ij) * obctmask(ii,ij) 407 END DO 408 END DO 409 END DO 410 411 DO ik = 1, jpkm1 412 DO ij = 2, jpjm1 413 DO ii = 2, jpim1 414 fmask(ii,ij,ik) = fmask(ii,ij,ik) * obctmask(ii,ij ) * obctmask(ii+1,ij ) & 415 & * obctmask(ii,ij+1) * obctmask(ii+1,ij+1) 416 END DO 417 END DO 418 END DO 419 420 tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:) 421 obctmask(:,:) = tmask(:,:,1) 422 423 ! obc masks and bmask are now set to zero on boundary points: 424 igrd = 1 ! In the free surface case, bmask is at T-points 425 DO ib_obc = 1, nb_obc 426 DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) 427 bmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0 428 ENDDO 429 ENDDO 430 ! 431 igrd = 1 432 DO ib_obc = 1, nb_obc 433 DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) 434 obctmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0 435 ENDDO 436 ENDDO 437 ! 438 igrd = 2 439 DO ib_obc = 1, nb_obc 440 DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) 441 obcumask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0 442 ENDDO 443 ENDDO 444 ! 445 igrd = 3 446 DO ib_obc = 1, nb_obc 447 DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) 448 obcvmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0 449 ENDDO 450 ENDDO 451 452 ! Lateral boundary conditions 453 CALL lbc_lnk( fmask , 'F', 1. ) ; CALL lbc_lnk( obctmask(:,:), 'T', 1. ) 454 CALL lbc_lnk( obcumask(:,:), 'U', 1. ) ; CALL lbc_lnk( obcvmask(:,:), 'V', 1. ) 455 456 DO ib_obc = 1, nb_obc ! Indices and directions of rim velocity components 457 458 idx_obc(ib_obc)%flagu(:) = 0.e0 459 idx_obc(ib_obc)%flagv(:) = 0.e0 460 icount = 0 461 462 !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward 463 !flagu = 0 : u is tangential 464 !flagu = 1 : u is normal to the boundary and is direction is inward 465 466 igrd = 2 ! u-component 467 DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) 468 nbi => idx_obc(ib_obc)%nbi(ib,igrd) 469 nbj => idx_obc(ib_obc)%nbj(ib,igrd) 470 zefl = obctmask(nbi ,nbj) 471 zwfl = obctmask(nbi+1,nbj) 472 IF( zefl + zwfl == 2 ) THEN 473 icount = icount + 1 226 474 ELSE 227 WRITE(numout,*)' ' 228 WRITE(numout,*)' Radiative East Open Boundary' 475 idx_obc(ib_obc)%flagu(ib)=-zefl+zwfl 476 ENDIF 477 END DO 478 479 !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward 480 !flagv = 0 : u is tangential 481 !flagv = 1 : u is normal to the boundary and is direction is inward 482 483 igrd = 3 ! v-component 484 DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) 485 nbi => idx_obc(ib_obc)%nbi(ib,igrd) 486 nbj => idx_obc(ib_obc)%nbj(ib,igrd) 487 znfl = obctmask(nbi,nbj ) 488 zsfl = obctmask(nbi,nbj+1) 489 IF( znfl + zsfl == 2 ) THEN 490 icount = icount + 1 491 ELSE 492 idx_obc(ib_obc)%flagv(ib) = -znfl + zsfl 229 493 END IF 230 END IF 231 END IF 232 233 IF( lp_obc_west ) THEN 234 ! ... mpp initialization 235 niw0 = max( 1, min(jpiwob - nimpp+1, jpi ) ) 236 niw1 = max( 0, min(jpiwob - nimpp+1, jpi - 1 ) ) 237 niw0p1 = max( 1, min(jpiwob+1 - nimpp+1, jpi ) ) 238 niw1p1 = max( 0, min(jpiwob+1 - nimpp+1, jpi - 1 ) ) 239 njw0 = max( 2, min(jpjwd - njmpp+1, jpj ) ) 240 njw1 = max( 0, min(jpjwf - njmpp+1, jpj - 1 ) ) 241 njw0p1 = max( 1, min(jpjwdp1 - njmpp+1, jpj ) ) 242 njw0m1 = max( 1, min(jpjwd - njmpp+1, jpj ) ) 243 njw1m1 = max( 0, min(jpjwfm1 - njmpp+1, jpj - 1 ) ) 244 njw1m2 = max( 0, min(jpjwfm1-1- njmpp+1, jpj - 1 ) ) 245 IF(lwp) THEN 246 IF( lfbcwest ) THEN 247 WRITE(numout,*)' ' 248 WRITE(numout,*)' Specified West Open Boundary' 249 ELSE 250 WRITE(numout,*)' ' 251 WRITE(numout,*)' Radiative West Open Boundary' 252 END IF 253 END IF 254 END IF 494 END DO 255 495 256 IF( lp_obc_north ) THEN 257 ! ... mpp initialization 258 nin0 = max( 2, min(jpind - nimpp+1, jpi ) ) 259 nin1 = max( 0, min(jpinf - nimpp+1, jpi - 1 ) ) 260 nin0p1 = max( 1, min(jpindp1 - nimpp+1, jpi ) ) 261 nin0m1 = max( 1, min(jpind - nimpp+1, jpi ) ) 262 nin1m1 = max( 0, min(jpinfm1 - nimpp+1, jpi - 1 ) ) 263 nin1m2 = max( 0, min(jpinfm1-1- nimpp+1, jpi - 1 ) ) 264 njn0 = max( 1, min(jpjnob - njmpp+1, jpj ) ) 265 njn1 = max( 0, min(jpjnob - njmpp+1, jpj - 1 ) ) 266 njn0p1 = max( 1, min(jpjnob+1 - njmpp+1, jpj ) ) 267 njn1p1 = max( 0, min(jpjnob+1 - njmpp+1, jpj - 1 ) ) 268 njn0m1 = max( 1, min(jpjnob-1 - njmpp+1, jpj ) ) 269 njn1m1 = max( 0, min(jpjnob-1 - njmpp+1, jpj - 1 ) ) 270 IF(lwp) THEN 271 IF( lfbcnorth ) THEN 272 WRITE(numout,*)' ' 273 WRITE(numout,*)' Specified North Open Boundary' 274 ELSE 275 WRITE(numout,*)' ' 276 WRITE(numout,*)' Radiative North Open Boundary' 277 END IF 278 END IF 279 END IF 280 281 IF( lp_obc_south ) THEN 282 ! ... mpp initialization 283 nis0 = max( 2, min(jpisd - nimpp+1, jpi ) ) 284 nis1 = max( 0, min(jpisf - nimpp+1, jpi - 1 ) ) 285 nis0p1 = max( 1, min(jpisdp1 - nimpp+1, jpi ) ) 286 nis0m1 = max( 1, min(jpisd - nimpp+1, jpi ) ) 287 nis1m1 = max( 0, min(jpisfm1 - nimpp+1, jpi - 1 ) ) 288 nis1m2 = max( 0, min(jpisfm1-1- nimpp+1, jpi - 1 ) ) 289 njs0 = max( 1, min(jpjsob - njmpp+1, jpj ) ) 290 njs1 = max( 0, min(jpjsob - njmpp+1, jpj - 1 ) ) 291 njs0p1 = max( 1, min(jpjsob+1 - njmpp+1, jpj ) ) 292 njs1p1 = max( 0, min(jpjsob+1 - njmpp+1, jpj - 1 ) ) 293 IF(lwp) THEN 294 IF( lfbcsouth ) THEN 295 WRITE(numout,*)' ' 296 WRITE(numout,*)' Specified South Open Boundary' 297 ELSE 298 WRITE(numout,*)' ' 299 WRITE(numout,*)' Radiative South Open Boundary' 300 END IF 301 END IF 302 END IF 303 304 ! 3. mask correction for OBCs 305 ! --------------------------- 306 307 IF( lp_obc_east ) THEN 308 !... (jpjed,jpjefm1),jpieob 309 bmask(nie0p1:nie1p1,nje0:nje1m1) = 0.e0 310 311 ! ... initilization to zero 312 uemsk(:,:) = 0.e0 ; vemsk(:,:) = 0.e0 ; temsk(:,:) = 0.e0 313 314 ! ... set 2D mask on East OBC, Vopt 315 DO ji = fs_nie0, fs_nie1 316 DO jj = nje0, nje1 317 uemsk(jj,:) = umask(ji, jj,:) * tmask_i(ji,jj) * tmask_i(ji+1,jj) 318 vemsk(jj,:) = vmask(ji+1,jj,:) * tmask_i(ji+1,jj) 319 temsk(jj,:) = tmask(ji+1,jj,:) * tmask_i(ji+1,jj) 320 END DO 321 END DO 322 323 END IF 324 325 IF( lp_obc_west ) THEN 326 ! ... (jpjwd,jpjwfm1),jpiwob 327 bmask(niw0:niw1,njw0:njw1m1) = 0.e0 328 329 ! ... initilization to zero 330 uwmsk(:,:) = 0.e0 ; vwmsk(:,:) = 0.e0 ; twmsk(:,:) = 0.e0 331 332 ! ... set 2D mask on West OBC, Vopt 333 DO ji = fs_niw0, fs_niw1 334 DO jj = njw0, njw1 335 uwmsk(jj,:) = umask(ji,jj,:) * tmask_i(ji,jj) * tmask_i(ji+1,jj) 336 vwmsk(jj,:) = vmask(ji,jj,:) * tmask_i(ji,jj) 337 twmsk(jj,:) = tmask(ji,jj,:) * tmask_i(ji,jj) 338 END DO 339 END DO 340 341 END IF 342 343 IF( lp_obc_north ) THEN 344 ! ... jpjnob,(jpind,jpisfm1) 345 bmask(nin0:nin1m1,njn0p1:njn1p1) = 0.e0 346 347 ! ... initilization to zero 348 unmsk(:,:) = 0.e0 ; vnmsk(:,:) = 0.e0 ; tnmsk(:,:) = 0.e0 349 350 ! ... set 2D mask on North OBC, Vopt 351 DO jj = fs_njn0, fs_njn1 352 DO ji = nin0, nin1 353 unmsk(ji,:) = umask(ji,jj+1,:) * tmask_i(ji,jj+1) 354 vnmsk(ji,:) = vmask(ji,jj ,:) * tmask_i(ji,jj) * tmask_i(ji,jj+1) 355 tnmsk(ji,:) = tmask(ji,jj+1,:) * tmask_i(ji,jj+1) 356 END DO 357 END DO 358 359 END IF 360 361 IF( lp_obc_south ) THEN 362 ! ... jpjsob,(jpisd,jpisfm1) 363 bmask(nis0:nis1m1,njs0:njs1) = 0.e0 364 365 ! ... initilization to zero 366 usmsk(:,:) = 0.e0 ; vsmsk(:,:) = 0.e0 ; tsmsk(:,:) = 0.e0 367 368 ! ... set 2D mask on South OBC, Vopt 369 DO jj = fs_njs0, fs_njs1 370 DO ji = nis0, nis1 371 usmsk(ji,:) = umask(ji,jj,:) * tmask_i(ji,jj) 372 vsmsk(ji,:) = vmask(ji,jj,:) * tmask_i(ji,jj) * tmask_i(ji,jj+1) 373 tsmsk(ji,:) = tmask(ji,jj,:) * tmask_i(ji,jj) 374 END DO 375 END DO 376 377 END IF 378 379 ! ... Initialize obcumask and obcvmask for the Force filtering 380 ! boundary condition in dynspg_flt 381 obcumask(:,:) = umask(:,:,1) 382 obcvmask(:,:) = vmask(:,:,1) 383 384 ! ... Initialize obctmsk on overlap region and obcs. This mask 385 ! is used in obcvol.F90 to calculate cumulate flux E-P. 386 ! obc Tracer point are outside the domain ( U/V obc points) ==> masked by obctmsk 387 ! - no flux E-P on obcs and overlap region (jpreci = jprecj = 1) 388 obctmsk(:,:) = tmask_i(:,:) 389 390 IF( lp_obc_east ) THEN 391 ! ... East obc Force filtering mask for the grad D 392 obcumask(nie0 :nie1 ,nje0p1:nje1m1) = 0.e0 393 obcvmask(nie0p1:nie1p1,nje0p1:nje1m1) = 0.e0 394 ! ... set to 0 on East OBC 395 obctmsk(nie0p1:nie1p1,nje0:nje1) = 0.e0 396 END IF 397 398 IF( lp_obc_west ) THEN 399 ! ... West obc Force filtering mask for the grad D 400 obcumask(niw0:niw1,njw0:njw1) = 0.e0 401 obcvmask(niw0:niw1,njw0:njw1) = 0.e0 402 ! ... set to 0 on West OBC 403 obctmsk(niw0:niw1,njw0:njw1) = 0.e0 404 END IF 405 406 IF( lp_obc_north ) THEN 407 ! ... North obc Force filtering mask for the grad D 408 obcumask(nin0p1:nin1m1,njn0p1:njn1p1) = 0.e0 409 obcvmask(nin0p1:nin1m1,njn0 :njn1 ) = 0.e0 410 ! ... set to 0 on North OBC 411 obctmsk(nin0:nin1,njn0p1:njn1p1) = 0.e0 412 END IF 413 414 IF( lp_obc_south ) THEN 415 ! ... South obc Force filtering mask for the grad D 416 obcumask(nis0p1:nis1m1,njs0:njs1) = 0.e0 417 obcvmask(nis0p1:nis1m1,njs0:njs1) = 0.e0 418 ! ... set to 0 on South OBC 419 obctmsk(nis0:nis1,njs0:njs1) = 0.e0 420 END IF 421 422 ! 3.1 Total lateral surface 423 ! ------------------------- 424 obcsurftot = 0.e0 425 426 IF( lp_obc_east ) THEN ! ... East open boundary lateral surface 427 DO ji = nie0, nie1 428 DO jj = 1, jpj 429 obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uemsk(jj,1) * MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 430 END DO 431 END DO 432 END IF 433 434 IF( lp_obc_west ) THEN ! ... West open boundary lateral surface 435 DO ji = niw0, niw1 436 DO jj = 1, jpj 437 obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uwmsk(jj,1) * MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 438 END DO 439 END DO 440 END IF 441 442 IF( lp_obc_north ) THEN ! ... North open boundary lateral surface 443 DO jj = njn0, njn1 444 DO ji = 1, jpi 445 obcsurftot = obcsurftot+hv(ji,jj)*e1v(ji,jj)*vnmsk(ji,1) * MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) ) 446 END DO 447 END DO 448 END IF 449 450 IF( lp_obc_south ) THEN ! ... South open boundary lateral surface 451 DO jj = njs0, njs1 452 DO ji = 1, jpi 453 obcsurftot = obcsurftot+hv(ji,jj)*e1v(ji,jj)*vsmsk(ji,1) * MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) ) 454 END DO 455 END DO 456 END IF 457 458 IF( lk_mpp ) CALL mpp_sum( obcsurftot ) ! sum over the global domain 459 460 ! 5. Control print on mask 461 ! The extremities of the open boundaries must be in land 462 ! or else correspond to an "ocean corner" between two open boundaries. 463 ! corner 1 is southwest, 2 is south east, 3 is northeast, 4 is northwest. 464 ! -------------------------------------------------------------------------- 465 466 icorner(:)=0 467 468 ! ... control of the west boundary 469 IF( lp_obc_west ) THEN 470 IF( jpiwob < 2 .OR. jpiwob >= jpiglo-2 ) THEN 471 WRITE(ctmp1,*) ' jpiwob exceed ', jpiglo-2, 'or less than 2' 472 CALL ctl_stop( ctmp1 ) 473 END IF 474 ztestmask(:)=0. 475 DO ji=niw0,niw1 476 IF( (njw0 + njmpp - 1) == jpjwd ) ztestmask(1)=ztestmask(1)+ tmask(ji,njw0,1) 477 IF( (njw1 + njmpp - 1) == jpjwf ) ztestmask(2)=ztestmask(2)+ tmask(ji,njw1,1) 478 END DO 479 IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain 480 481 IF( ztestmask(1) /= 0. ) icorner(1)=icorner(1)+1 482 IF( ztestmask(2) /= 0. ) icorner(4)=icorner(4)+1 483 END IF 484 485 ! ... control of the east boundary 486 IF( lp_obc_east ) THEN 487 IF( jpieob < 4 .OR. jpieob >= jpiglo ) THEN 488 WRITE(ctmp1,*) ' jpieob exceed ', jpiglo, ' or less than 4' 489 CALL ctl_stop( ctmp1 ) 490 END IF 491 ztestmask(:)=0. 492 DO ji=nie0p1,nie1p1 493 IF( (nje0 + njmpp - 1) == jpjed ) ztestmask(1)=ztestmask(1)+ tmask(ji,nje0,1) 494 IF( (nje1 + njmpp - 1) == jpjef ) ztestmask(2)=ztestmask(2)+ tmask(ji,nje1,1) 495 END DO 496 IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain 497 498 IF( ztestmask(1) /= 0. ) icorner(2)=icorner(2)+1 499 IF( ztestmask(2) /= 0. ) icorner(3)=icorner(3)+1 500 END IF 501 502 ! ... control of the north boundary 503 IF( lp_obc_north ) THEN 504 IF( jpjnob < 4 .OR. jpjnob >= jpjglo ) THEN 505 WRITE(ctmp1,*) 'jpjnob exceed ', jpjglo, ' or less than 4' 506 CALL ctl_stop( ctmp1 ) 507 END IF 508 ztestmask(:)=0. 509 DO jj=njn0p1,njn1p1 510 IF( (nin0 + nimpp - 1) == jpind ) ztestmask(1)=ztestmask(1)+ tmask(nin0,jj,1) 511 IF( (nin1 + nimpp - 1) == jpinf ) ztestmask(2)=ztestmask(2)+ tmask(nin1,jj,1) 512 END DO 513 IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain 514 515 IF( ztestmask(1) /= 0. ) icorner(4)=icorner(4)+1 516 IF( ztestmask(2) /= 0. ) icorner(3)=icorner(3)+1 517 END IF 518 519 ! ... control of the south boundary 520 IF( lp_obc_south ) THEN 521 IF( jpjsob < 2 .OR. jpjsob >= jpjglo-2 ) THEN 522 WRITE(ctmp1,*) ' jpjsob exceed ', jpjglo-2, ' or less than 2' 523 CALL ctl_stop( ctmp1 ) 524 END IF 525 ztestmask(:)=0. 526 DO jj=njs0,njs1 527 IF( (nis0 + nimpp - 1) == jpisd ) ztestmask(1)=ztestmask(1)+ tmask(nis0,jj,1) 528 IF( (nis1 + nimpp - 1) == jpisf ) ztestmask(2)=ztestmask(2)+ tmask(nis1,jj,1) 529 END DO 530 IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain 531 532 IF( ztestmask(1) /= 0. ) icorner(1)=icorner(1)+1 533 IF( ztestmask(2) /= 0. ) icorner(2)=icorner(2)+1 534 END IF 535 536 IF( icorner(1) == 2 ) THEN 537 IF(lwp) WRITE(numout,*) 538 IF(lwp) WRITE(numout,*) ' South West ocean corner, two open boudaries' 539 IF(lwp) WRITE(numout,*) ' ========== ' 540 IF(lwp) WRITE(numout,*) 541 IF( jpisd /= jpiwob.OR.jpjsob /= jpjwd ) & 542 & CALL ctl_stop( ' Open boundaries do not fit, we stop' ) 543 544 ELSE IF( icorner(1) == 1 ) THEN 545 CALL ctl_stop( ' Open boundaries do not fit at SW corner, we stop' ) 546 END IF 547 548 IF( icorner(2) == 2 ) THEN 549 IF(lwp) WRITE(numout,*) 550 IF(lwp) WRITE(numout,*) ' South East ocean corner, two open boudaries' 551 IF(lwp) WRITE(numout,*) ' ========== ' 552 IF(lwp) WRITE(numout,*) 553 IF( jpisf /= jpieob+1.OR.jpjsob /= jpjed ) & 554 & CALL ctl_stop( ' Open boundaries do not fit, we stop' ) 555 ELSE IF( icorner(2) == 1 ) THEN 556 CALL ctl_stop( ' Open boundaries do not fit at SE corner, we stop' ) 557 END IF 558 559 IF( icorner(3) == 2 ) THEN 560 IF(lwp) WRITE(numout,*) 561 IF(lwp) WRITE(numout,*) ' North East ocean corner, two open boudaries' 562 IF(lwp) WRITE(numout,*) ' ========== ' 563 IF(lwp) WRITE(numout,*) 564 IF( jpinf /= jpieob+1 .OR. jpjnob+1 /= jpjef ) & 565 & CALL ctl_stop( ' Open boundaries do not fit, we stop' ) 566 ELSE IF( icorner(3) == 1 ) THEN 567 CALL ctl_stop( ' Open boundaries do not fit at NE corner, we stop' ) 568 END IF 569 570 IF( icorner(4) == 2 ) THEN 571 IF(lwp) WRITE(numout,*) 572 IF(lwp) WRITE(numout,*) ' North West ocean corner, two open boudaries' 573 IF(lwp) WRITE(numout,*) ' ========== ' 574 IF(lwp) WRITE(numout,*) 575 IF( jpind /= jpiwob.OR.jpjnob+1 /= jpjwf ) & 576 & CALL ctl_stop( ' Open boundaries do not fit, we stop' ) 577 ELSE IF( icorner(4) == 1 ) THEN 578 CALL ctl_stop( ' Open boundaries do not fit at NW corner, we stop' ) 579 END IF 580 581 ! 6. Initialization of open boundary variables (u, v, t, s) 582 ! -------------------------------------------------------------- 583 ! only if at least one boundary is radiative 584 IF ( inumfbc < nbobc .AND. ln_rstart ) THEN 585 ! Restart from restart.obc 586 CALL obc_rst_read 587 ELSE 588 589 ! ! ... Initialization to zero of radiation arrays. 590 ! ! Those have dimensions of local subdomains 591 592 uebnd(:,:,:,:) = 0.e0 ; unbnd(:,:,:,:) = 0.e0 593 vebnd(:,:,:,:) = 0.e0 ; vnbnd(:,:,:,:) = 0.e0 594 tebnd(:,:,:,:) = 0.e0 ; tnbnd(:,:,:,:) = 0.e0 595 sebnd(:,:,:,:) = 0.e0 ; snbnd(:,:,:,:) = 0.e0 596 597 uwbnd(:,:,:,:) = 0.e0 ; usbnd(:,:,:,:) = 0.e0 598 vwbnd(:,:,:,:) = 0.e0 ; vsbnd(:,:,:,:) = 0.e0 599 twbnd(:,:,:,:) = 0.e0 ; tsbnd(:,:,:,:) = 0.e0 600 swbnd(:,:,:,:) = 0.e0 ; ssbnd(:,:,:,:) = 0.e0 601 602 END IF 603 604 ! 7. Control print 605 ! ----------------------------------------------------------------- 606 607 ! ... control of the east boundary 608 IF( lp_obc_east ) THEN 609 istop = 0 610 IF( jpieob < 4 .OR. jpieob >= jpiglo ) THEN 611 IF(lwp) WRITE(numout,cform_err) 612 IF(lwp) WRITE(numout,*) ' jpieob exceed ', jpim1, ' or less than 4' 613 istop = istop + 1 614 END IF 615 616 IF( lk_mpp ) THEN 617 ! ... 618 IF( nimpp > jpieob-5) THEN 619 IF(lwp) WRITE(numout,cform_err) 620 IF(lwp) WRITE(numout,*) ' A sub-domain is too close to the East OBC' 621 IF(lwp) WRITE(numout,*) ' nimpp must be < jpieob-5' 622 istop = istop + 1 623 ENDIF 624 ELSE 625 626 ! ... stop if e r r o r (s) detected 627 IF( istop /= 0 ) THEN 628 WRITE(ctmp1,*) istop,' obcini : E R R O R (S) detected : stop' 629 CALL ctl_stop( ctmp1 ) 630 ENDIF 631 ENDIF 632 ENDIF 633 634 ! ... control of the west boundary 635 IF( lp_obc_west ) THEN 636 istop = 0 637 IF( jpiwob < 2 .OR. jpiwob >= jpiglo ) THEN 638 IF(lwp) WRITE(numout,cform_err) 639 IF(lwp) WRITE(numout,*) ' jpiwob exceed ', jpim1, ' or less than 2' 640 istop = istop + 1 641 END IF 642 643 IF( lk_mpp ) THEN 644 IF( (nimpp < jpiwob+5) .AND. (nimpp > 1) ) THEN 645 IF(lwp) WRITE(numout,cform_err) 646 IF(lwp) WRITE(numout,*) ' A sub-domain is too close to the West OBC' 647 IF(lwp) WRITE(numout,*) ' nimpp must be > jpiwob-5 or =1' 648 istop = istop + 1 649 ENDIF 650 ELSE 651 652 ! ... stop if e r r o r (s) detected 653 IF( istop /= 0 ) THEN 654 WRITE(ctmp1,*) istop,' obcini : E R R O R (S) detected : stop' 655 CALL ctl_stop( ctmp1 ) 656 ENDIF 657 ENDIF 658 ENDIF 659 660 ! control of the north boundary 661 IF( lp_obc_north ) THEN 662 istop = 0 663 IF( jpjnob < 4 .OR. jpjnob >= jpjglo ) THEN 664 IF(lwp) WRITE(numout,cform_err) 665 IF(lwp) WRITE(numout,*) ' jpjnob exceed ', jpjm1,' or less than 4' 666 istop = istop + 1 667 END IF 668 669 IF( lk_mpp ) THEN 670 IF( njmpp > jpjnob-5) THEN 671 IF(lwp) WRITE(numout,cform_err) 672 IF(lwp) WRITE(numout,*) ' A sub-domain is too close to the North OBC' 673 IF(lwp) WRITE(numout,*) ' njmpp must be < jpjnob-5' 674 istop = istop + 1 675 ENDIF 676 ELSE 677 678 ! ... stop if e r r o r (s) detected 679 IF( istop /= 0 ) THEN 680 WRITE(ctmp1,*) istop,' obcini : E R R O R (S) detected : stop' 681 CALL ctl_stop( ctmp1 ) 682 ENDIF 683 ENDIF 684 ENDIF 685 686 ! control of the south boundary 687 IF( lp_obc_south ) THEN 688 istop = 0 689 IF( jpjsob < 2 .OR. jpjsob >= jpjglo ) THEN 690 IF(lwp) WRITE(numout,cform_err) 691 IF(lwp) WRITE(numout,*) ' jpjsob exceed ', jpjm1,' or less than 2' 692 istop = istop + 1 693 END IF 694 695 IF( lk_mpp ) THEN 696 IF( (njmpp < jpjsob+5) .AND. (njmpp > 1) ) THEN 697 IF(lwp) WRITE(numout,cform_err) 698 IF(lwp) WRITE(numout,*) ' A sub-domain is too close to the South OBC' 699 IF(lwp) WRITE(numout,*) ' njmpp must be > jpjsob+5 or =1' 700 istop = istop + 1 701 ENDIF 702 ELSE 703 704 ! ... stop if e r r o r (s) detected 705 IF( istop /= 0 ) THEN 706 WRITE(ctmp1,*) istop,' obcini : E R R O R (S) detected : stop' 707 CALL ctl_stop( ctmp1 ) 708 ENDIF 709 ENDIF 710 ENDIF 496 IF( icount /= 0 ) THEN 497 IF(lwp) WRITE(numout,*) 498 IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,', & 499 ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_obc 500 IF(lwp) WRITE(numout,*) ' ========== ' 501 IF(lwp) WRITE(numout,*) 502 nstop = nstop + 1 503 ENDIF 504 505 ENDDO 506 507 ! Compute total lateral surface for volume correction: 508 ! ---------------------------------------------------- 509 obcsurftot = 0.e0 510 IF( ln_vol ) THEN 511 igrd = 2 ! Lateral surface at U-points 512 DO ib_obc = 1, nb_obc 513 DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) 514 nbi => idx_obc(ib_obc)%nbi(ib,igrd) 515 nbj => idx_obc(ib_obc)%nbi(ib,igrd) 516 flagu => idx_obc(ib_obc)%flagu(ib) 517 obcsurftot = obcsurftot + hu (nbi , nbj) & 518 & * e2u (nbi , nbj) * ABS( flagu ) & 519 & * tmask_i(nbi , nbj) & 520 & * tmask_i(nbi+1, nbj) 521 ENDDO 522 ENDDO 523 524 igrd=3 ! Add lateral surface at V-points 525 DO ib_obc = 1, nb_obc 526 DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) 527 nbi => idx_obc(ib_obc)%nbi(ib,igrd) 528 nbj => idx_obc(ib_obc)%nbi(ib,igrd) 529 flagv => idx_obc(ib_obc)%flagv(ib) 530 obcsurftot = obcsurftot + hv (nbi, nbj ) & 531 & * e1v (nbi, nbj ) * ABS( flagv ) & 532 & * tmask_i(nbi, nbj ) & 533 & * tmask_i(nbi, nbj+1) 534 ENDDO 535 ENDDO 536 ! 537 IF( lk_mpp ) CALL mpp_sum( obcsurftot ) ! sum over the global domain 538 END IF 539 540 ! Read in tidal constituents and adjust for model start time 541 ! ---------------------------------------------------------- 542 !!$ IF( ln_tides ) CALL tide_data 543 ! 544 ! Tidy up 545 !-------- 546 DEALLOCATE(nbidta, nbjdta, nbrdta) 711 547 712 548 END SUBROUTINE obc_init … … 714 550 #else 715 551 !!--------------------------------------------------------------------------------- 716 !! Dummy module NOopen boundaries552 !! Dummy module NO unstructured open boundaries 717 553 !!--------------------------------------------------------------------------------- 718 554 CONTAINS
Note: See TracChangeset
for help on using the changeset viewer.
