MODULE obcini !!====================================================================== !! *** MODULE obcini *** !! Unstructured open boundaries : initialisation !!====================================================================== !! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code !! - ! 2007-01 (D. Storkey) Update to use IOM module !! - ! 2007-01 (D. Storkey) Tidal forcing !! 3.0 ! 2008-04 (NEMO team) add in the reference version !! 3.3 ! 2010-09 (E.O'Dea) updates for Shelf configurations !! 3.3 ! 2010-09 (D.Storkey) add ice boundary conditions !! 3.4 ! 2011 (D. Storkey, J. Chanut) OBC-BDY merge !! ! --- Renamed bdyini.F90 -> obcini.F90 --- !!---------------------------------------------------------------------- #if defined key_obc !!---------------------------------------------------------------------- !! 'key_obc' Unstructured Open Boundary Conditions !!---------------------------------------------------------------------- !! obc_init : Initialization of unstructured open boundaries !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers variables USE dom_oce ! ocean space and time domain USE obc_oce ! unstructured open boundary conditions USE in_out_manager ! I/O units USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE lib_mpp ! for mpp_sum USE iom ! I/O IMPLICIT NONE PRIVATE PUBLIC obc_init ! routine called by opa.F90 !!---------------------------------------------------------------------- !! NEMO/OPA 4.0 , NEMO Consortium (2011) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE obc_init !!---------------------------------------------------------------------- !! *** ROUTINE obc_init *** !! !! ** Purpose : Initialization of the dynamics and tracer fields with !! unstructured open boundaries. !! !! ** Method : Read initialization arrays (mask, indices) to identify !! an unstructured open boundary !! !! ** Input : obc_init.nc, input file for unstructured open boundaries !!---------------------------------------------------------------------- INTEGER :: ib_obc, ii, ij, ik, igrd, ib, ir ! dummy loop indices INTEGER :: icount, icountr, ibr_max, ilen1 ! local integers INTEGER :: iw, ie, is, in, inum, id_dummy ! - - INTEGER :: igrd_start, igrd_end, jpbdta ! - - INTEGER, POINTER :: nbi, nbj, nbr ! short cuts REAL , POINTER :: flagu, flagv ! - - REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars INTEGER, DIMENSION (2) :: kdimsz INTEGER, DIMENSION(jpbgrd,jp_obc) :: nblendta ! Length of index arrays INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of obc dta INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points REAL(wp), DIMENSION(jpidta,jpjdta) :: zmask ! global domain mask CHARACTER(LEN=80),DIMENSION(jpbgrd) :: clfile CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid !! NAMELIST/namobc/ nb_obc, ln_coords_file, cn_coords_file, & & ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn3d, & & nn_tra, & #if defined key_lim2 & nn_ice_lim2, & #endif & nn_tides, ln_vol, ln_clim, nn_dtactl, nn_volctl, & & nn_rimwidth, nn_dmp2d_in, nn_dmp2d_out, & & nn_dmp3d_in, nn_dmp3d_out !!---------------------------------------------------------------------- IF( obc_oce_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'obc_init : unable to allocate oce arrays' ) IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'obc_init : initialization of open boundaries' IF(lwp) WRITE(numout,*) '~~~~~~~~' ! IF( jperio /= 0 ) CALL ctl_stop( 'Cyclic or symmetric,', & & ' and general open boundary condition are not compatible' ) cgrid= (/'T','U','V'/) ! ----------------------------------------- ! Initialise and read namelist parameters ! ----------------------------------------- nb_obc = 0 ln_coords_file(:) = .false. cn_coords_file(:) = '' ln_mask_file = .false. cn_mask_file(:) = '' nn_dyn2d(:) = 0 nn_dyn3d(:) = 0 nn_tra(:) = 0 #if defined key_lim2 nn_ice_lim2(:) = 0 #endif ln_vol = .false. ln_clim(:) = .false. nn_dtactl(:) = -1 ! uninitialised flag nn_tides(:) = 0 ! default to no tidal forcing nn_volctl = -1 ! uninitialised flag nn_rimwidth(:) = -1 ! uninitialised flag nn_dmp2d_in(:) = -1 ! uninitialised flag nn_dmp2d_out(:) = -1 ! uninitialised flag nn_dmp3d_in(:) = -1 ! uninitialised flag nn_dmp3d_out(:) = -1 ! uninitialised flag REWIND( numnam ) READ ( numnam, namobc ) ! ----------------------------------------- ! Check and write out namelist parameters ! ----------------------------------------- ! ! control prints IF(lwp) WRITE(numout,*) ' namobc' IF( nb_obc .eq. 0 ) THEN IF(lwp) WRITE(numout,*) 'nb_obc = 0, NO OPEN BOUNDARIES APPLIED.' ELSE IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_obc ENDIF DO ib_obc = 1,nb_obc IF(lwp) WRITE(numout,*) ' ' IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_obc,'------' ! ! check type of data used (nn_dtactl value) SELECT CASE( nn_dtactl(ib_obc) ) ! CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for obc data' CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file' CASE DEFAULT ; CALL ctl_stop( 'nn_dtactl must be 0 or 1' ) END SELECT IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution: ' SELECT CASE( nn_dyn2d(ib_obc) ) CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Flather radiation condition' CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) END SELECT IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities: ' SELECT CASE( nn_dyn3d(ib_obc) ) CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn3d' ) END SELECT IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity: ' SELECT CASE( nn_tra(ib_obc) ) CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_tra' ) END SELECT IF(lwp) WRITE(numout,*) #if defined key_lim2 IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice: ' SELECT CASE( nn_tra(ib_obc) ) CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_tra' ) END SELECT IF(lwp) WRITE(numout,*) #endif IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nn_rimwidth = ', nn_rimwidth IF(lwp) WRITE(numout,*) SELECT CASE( nn_tides(ib_obc) ) CASE(0) IF(lwp) WRITE(numout,*) 'No tidal harmonic forcing' IF(lwp) WRITE(numout,*) CASE(1) IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ONLY for barotropic solution' IF(lwp) WRITE(numout,*) CASE(2) IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ADDED to other barotropic boundary conditions' IF(lwp) WRITE(numout,*) CASE DEFAULT CALL ctl_stop( 'obc_ini: ERROR: incorrect value for nn_tides ' ) END SELECT ENDDO IF( ln_vol ) THEN ! check volume conservation (nn_volctl value) IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' IF(lwp) WRITE(numout,*) SELECT CASE ( nn_volctl ) CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' The total volume will be constant' CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' The total volume will vary according to the surface E-P flux' CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' ) END SELECT IF(lwp) WRITE(numout,*) ELSE IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' IF(lwp) WRITE(numout,*) ENDIF ! ------------------------------------------------- ! Initialise indices arrays for open boundaries ! ------------------------------------------------- ! Work out global dimensions of boundary data ! --------------------------------------------- DO ib_obc = 1, nb_obc jpbdta = 1 IF( .NOT. ln_coords_file(ib_obc) ) THEN ! Work out size of global arrays from namelist parameters !! 1. Read parameters from namelist !! 2. Work out global size of boundary data arrays nblendta and jpbdta ELSE ! Read size of arrays in boundary coordinates file. CALL iom_open( cn_coords_file(ib_obc), inum ) jpbdta = 1 DO igrd = 1, jpbgrd id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) nblendta(igrd,ib_obc) = kdimsz(1) jpbdta = MAX(jpbdta, kdimsz(1)) ENDDO ENDIF ENDDO ! Allocate arrays !--------------- ALLOCATE( nbidta(jpbdta, jpbgrd, nb_obc), nbjdta(jpbdta, jpbgrd, nb_obc), & & nbrdta(jpbdta, jpbgrd, nb_obc) ) ALLOCATE( dta_global(jpbdta, 1, jpk) ) ! Calculate global boundary index arrays or read in from file !------------------------------------------------------------ DO ib_obc = 1, nb_obc IF( .NOT. ln_coords_file(ib_obc) ) THEN ! Calculate global index arrays from namelist parameters !! Calculate global index arrays from namelist parameters ELSE ! Read global index arrays from boundary coordinates file. DO igrd = 1, jpbgrd CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_obc),:,1) ) DO ii = 1,nblendta(igrd,ib_obc) nbidta(ii,igrd,ib_obc) = INT( dta_global(ii,1,1) ) END DO CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_obc),:,1) ) DO ii = 1,nblendta(igrd,ib_obc) nbjdta(ii,igrd,ib_obc) = INT( dta_global(ii,1,1) ) END DO CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_obc),:,1) ) DO ii = 1,nblendta(igrd,ib_obc) nbrdta(ii,igrd,ib_obc) = INT( dta_global(ii,1,1) ) END DO ibr_max = MAXVAL( nbrdta(:,igrd,ib_obc) ) IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_obc) IF (ibr_max < nn_rimwidth(ib_obc)) & CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_obc) ) END DO CALL iom_close( inum ) ENDIF ENDDO ! Work out dimensions of boundary data on each processor ! ------------------------------------------------------ iw = mig(1) + 1 ! if monotasking and no zoom, iw=2 ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1 is = mjg(1) + 1 ! if monotasking and no zoom, is=2 in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1 DO ib_obc = 1, nb_obc DO igrd = 1, jpbgrd icount = 0 icountr = 0 idx_obc(ib_obc)%nblen(igrd) = 0 idx_obc(ib_obc)%nblenrim(igrd) = 0 DO ib = 1, nblendta(igrd,ib_obc) ! check if point is in local domain IF( nbidta(ib,igrd,ib_obc) >= iw .AND. nbidta(ib,igrd,ib_obc) <= ie .AND. & & nbjdta(ib,igrd,ib_obc) >= is .AND. nbjdta(ib,igrd,ib_obc) <= in ) THEN ! icount = icount + 1 ! IF( nbrdta(ib,igrd,ib_obc) == 1 ) icountr = icountr+1 ENDIF ENDDO idx_obc(ib_obc)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc idx_obc(ib_obc)%nblen (igrd) = icount !: length of boundary data on each proc ENDDO ! igrd ! Allocate index arrays for this boundary set !-------------------------------------------- ilen1 = MAXVAL(idx_obc(ib_obc)%nblen(:)) ALLOCATE( idx_obc(ib_obc)%nbi(ilen1,jpbgrd) ) ALLOCATE( idx_obc(ib_obc)%nbj(ilen1,jpbgrd) ) ALLOCATE( idx_obc(ib_obc)%nbr(ilen1,jpbgrd) ) ALLOCATE( idx_obc(ib_obc)%nbmap(ilen1,jpbgrd) ) ALLOCATE( idx_obc(ib_obc)%nbw(ilen1,jpbgrd) ) ALLOCATE( idx_obc(ib_obc)%flagu(ilen1) ) ALLOCATE( idx_obc(ib_obc)%flagv(ilen1) ) ! Dispatch mapping indices and discrete distances on each processor ! ----------------------------------------------------------------- DO igrd = 1, jpbgrd icount = 0 ! Loop on rimwidth to ensure outermost points come first in the local arrays. DO ir=1, nn_rimwidth(ib_obc) DO ib = 1, nblendta(igrd,ib_obc) ! check if point is in local domain and equals ir IF( nbidta(ib,igrd,ib_obc) >= iw .AND. nbidta(ib,igrd,ib_obc) <= ie .AND. & & nbjdta(ib,igrd,ib_obc) >= is .AND. nbjdta(ib,igrd,ib_obc) <= in .AND. & & nbrdta(ib,igrd,ib_obc) == ir ) THEN ! icount = icount + 1 idx_obc(ib_obc)%nbi(icount,igrd) = nbidta(ib,igrd,ib_obc)- mig(1)+1 idx_obc(ib_obc)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_obc)- mjg(1)+1 idx_obc(ib_obc)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_obc) idx_obc(ib_obc)%nbmap(icount,igrd) = ib ENDIF ENDDO ENDDO ENDDO ! Compute rim weights for FRS scheme ! ---------------------------------- DO igrd = 1, jpbgrd DO ib = 1, idx_obc(ib_obc)%nblen(igrd) nbr => idx_obc(ib_obc)%nbr(ib,igrd) idx_obc(ib_obc)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ! tanh formulation ! idx_obc(ib_obc)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2 ! quadratic ! idx_obc(ib_obc)%nbw(ib,igrd) = FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth) ! linear END DO END DO ENDDO ! ------------------------------------------------------ ! Initialise masks and find normal/tangential directions ! ------------------------------------------------------ ! Read global 2D mask at T-points: obctmask ! ----------------------------------------- ! obctmask = 1 on the computational domain AND on open boundaries ! = 0 elsewhere IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN ! EEL configuration at 5km resolution zmask( : ,:) = 0.e0 zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0 ELSE IF( ln_mask_file ) THEN CALL iom_open( cn_mask_file, inum ) CALL iom_get ( inum, jpdom_data, 'obc_msk', zmask(:,:) ) CALL iom_close( inum ) ELSE zmask(:,:) = 1.e0 ENDIF DO ij = 1, nlcj ! Save mask over local domain DO ii = 1, nlci obctmask(ii,ij) = zmask( mig(ii), mjg(ij) ) END DO END DO ! Derive mask on U and V grid from mask on T grid obcumask(:,:) = 0.e0 obcvmask(:,:) = 0.e0 DO ij=1, jpjm1 DO ii=1, jpim1 obcumask(ii,ij)=obctmask(ii,ij)*obctmask(ii+1, ij ) obcvmask(ii,ij)=obctmask(ii,ij)*obctmask(ii ,ij+1) END DO END DO CALL lbc_lnk( obcumask(:,:), 'U', 1. ) ; CALL lbc_lnk( obcvmask(:,:), 'V', 1. ) ! Lateral boundary cond. ! Mask corrections ! ---------------- DO ik = 1, jpkm1 DO ij = 1, jpj DO ii = 1, jpi tmask(ii,ij,ik) = tmask(ii,ij,ik) * obctmask(ii,ij) umask(ii,ij,ik) = umask(ii,ij,ik) * obcumask(ii,ij) vmask(ii,ij,ik) = vmask(ii,ij,ik) * obcvmask(ii,ij) bmask(ii,ij) = bmask(ii,ij) * obctmask(ii,ij) END DO END DO END DO DO ik = 1, jpkm1 DO ij = 2, jpjm1 DO ii = 2, jpim1 fmask(ii,ij,ik) = fmask(ii,ij,ik) * obctmask(ii,ij ) * obctmask(ii+1,ij ) & & * obctmask(ii,ij+1) * obctmask(ii+1,ij+1) END DO END DO END DO tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:) obctmask(:,:) = tmask(:,:,1) ! obc masks and bmask are now set to zero on boundary points: igrd = 1 ! In the free surface case, bmask is at T-points DO ib_obc = 1, nb_obc DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) bmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0 ENDDO ENDDO ! igrd = 1 DO ib_obc = 1, nb_obc DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) obctmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0 ENDDO ENDDO ! igrd = 2 DO ib_obc = 1, nb_obc DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) obcumask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0 ENDDO ENDDO ! igrd = 3 DO ib_obc = 1, nb_obc DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) obcvmask(idx_obc(ib_obc)%nbi(ib,igrd), idx_obc(ib_obc)%nbj(ib,igrd)) = 0.e0 ENDDO ENDDO ! Lateral boundary conditions CALL lbc_lnk( fmask , 'F', 1. ) ; CALL lbc_lnk( obctmask(:,:), 'T', 1. ) CALL lbc_lnk( obcumask(:,:), 'U', 1. ) ; CALL lbc_lnk( obcvmask(:,:), 'V', 1. ) DO ib_obc = 1, nb_obc ! Indices and directions of rim velocity components idx_obc(ib_obc)%flagu(:) = 0.e0 idx_obc(ib_obc)%flagv(:) = 0.e0 icount = 0 !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward !flagu = 0 : u is tangential !flagu = 1 : u is normal to the boundary and is direction is inward igrd = 2 ! u-component DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) nbi => idx_obc(ib_obc)%nbi(ib,igrd) nbj => idx_obc(ib_obc)%nbj(ib,igrd) zefl = obctmask(nbi ,nbj) zwfl = obctmask(nbi+1,nbj) IF( zefl + zwfl == 2 ) THEN icount = icount + 1 ELSE idx_obc(ib_obc)%flagu(ib)=-zefl+zwfl ENDIF END DO !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward !flagv = 0 : u is tangential !flagv = 1 : u is normal to the boundary and is direction is inward igrd = 3 ! v-component DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) nbi => idx_obc(ib_obc)%nbi(ib,igrd) nbj => idx_obc(ib_obc)%nbj(ib,igrd) znfl = obctmask(nbi,nbj ) zsfl = obctmask(nbi,nbj+1) IF( znfl + zsfl == 2 ) THEN icount = icount + 1 ELSE idx_obc(ib_obc)%flagv(ib) = -znfl + zsfl END IF END DO IF( icount /= 0 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,', & ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_obc IF(lwp) WRITE(numout,*) ' ========== ' IF(lwp) WRITE(numout,*) nstop = nstop + 1 ENDIF ENDDO ! Compute total lateral surface for volume correction: ! ---------------------------------------------------- obcsurftot = 0.e0 IF( ln_vol ) THEN igrd = 2 ! Lateral surface at U-points DO ib_obc = 1, nb_obc DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) nbi => idx_obc(ib_obc)%nbi(ib,igrd) nbj => idx_obc(ib_obc)%nbi(ib,igrd) flagu => idx_obc(ib_obc)%flagu(ib) obcsurftot = obcsurftot + hu (nbi , nbj) & & * e2u (nbi , nbj) * ABS( flagu ) & & * tmask_i(nbi , nbj) & & * tmask_i(nbi+1, nbj) ENDDO ENDDO igrd=3 ! Add lateral surface at V-points DO ib_obc = 1, nb_obc DO ib = 1, idx_obc(ib_obc)%nblenrim(igrd) nbi => idx_obc(ib_obc)%nbi(ib,igrd) nbj => idx_obc(ib_obc)%nbi(ib,igrd) flagv => idx_obc(ib_obc)%flagv(ib) obcsurftot = obcsurftot + hv (nbi, nbj ) & & * e1v (nbi, nbj ) * ABS( flagv ) & & * tmask_i(nbi, nbj ) & & * tmask_i(nbi, nbj+1) ENDDO ENDDO ! IF( lk_mpp ) CALL mpp_sum( obcsurftot ) ! sum over the global domain END IF ! ! Tidy up !-------- DEALLOCATE(nbidta, nbjdta, nbrdta) END SUBROUTINE obc_init #else !!--------------------------------------------------------------------------------- !! Dummy module NO open boundaries !!--------------------------------------------------------------------------------- CONTAINS SUBROUTINE obc_init ! Dummy routine END SUBROUTINE obc_init #endif !!================================================================================= END MODULE obcini