MODULE bdyini !!================================================================================= !! *** MODULE bdyini *** !! Initialization of unstructured open boundaries !!================================================================================= #if defined key_bdy || defined key_bdy_tides !!--------------------------------------------------------------------------------- !! 'key_bdy' Unstructured Open Boundary Conditions !!--------------------------------------------------------------------------------- !! bdy_init : Initialization of unstructured open boundaries !!--------------------------------------------------------------------------------- !! * Modules used USE oce ! ocean dynamics and tracers variables USE dom_oce ! ocean space and time domain USE bdy_oce ! unstructured open boundary conditions USE bdytides ! tides at open boundaries initialization (tide_init routine) 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 !! * Routine accessibility PUBLIC bdy_init ! routine called by opa.F90 !! * Substitutions !!--------------------------------------------------------------------------------- !! OPA 9.0 , LODYC-IPSL (2003) !!--------------------------------------------------------------------------------- CONTAINS SUBROUTINE bdy_init !!---------------------------------------------------------------------- !! *** ROUTINE bdy_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 : bdy_init.nc, input file for unstructured open boundaries !! !! History : !! OPA 9.0 ! 05-01 (J. Chanut, A. Sellar) Original code !! ! 07-01 (D. Storkey) Update to use IOM module. !! ! 07-01 (D. Storkey) Tidal forcing. !!---------------------------------------------------------------------- !! * Local declarations INTEGER :: ji, jj, jk, jgrd, & ! dummy loop indices jb, jr, icount, & icountr, nb_rim, nb_len, nbr_max INTEGER :: iw, ie, is, in INTEGER :: inum ! temporary logical unit INTEGER :: & ! temporary integers dummy_id INTEGER, DIMENSION (2) :: kdimsz INTEGER, DIMENSION(jpbdta, jpbgrd) :: & !: Index arrays nbidta, nbjdta, & !: i and j indices of bdy dta nbrdta !: Discrete distance from rim points REAL(wp) :: & efl, wfl, nfl, sfl ! temporary scalars REAL(wp) , DIMENSION(jpidta,jpjdta) :: & tmpmsk ! global domain mask REAL(wp) , DIMENSION(jpbdta,1) :: & ndta ! temporary array CHARACTER(LEN=80),DIMENSION(3) :: bdyfile NAMELIST/nambdy/filbdy_mask, filbdy_data_T, filbdy_data_U, filbdy_data_V, & ln_bdy_clim, ln_bdy_vol, ln_bdy_fla, ln_bdy_mask, & nbdy_dta, nb_rimwidth, volbdy !!---------------------------------------------------------------------- IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'bdy_init : initialization of unstructured open boundaries' IF(lwp) WRITE(numout,*) '~~~~~~~~' IF( jperio /= 0 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' E R R O R : Cyclic or symmetric,', & ' and unstructured open boundary condition are not compatible' IF(lwp) WRITE(numout,*) ' ========== ' IF(lwp) WRITE(numout,*) nstop = nstop + 1 END IF #if defined key_obc IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' E R R O R : Straight open boundaries,', & ' and unstructured open boundaries are not compatible' IF(lwp) WRITE(numout,*) ' ========== ' IF(lwp) WRITE(numout,*) nstop = nstop + 1 #endif # if defined key_dynspg_rl IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' E R R O R : Rigid lid,', & ' and unstructured open boundaries are not compatible' IF(lwp) WRITE(numout,*) ' ========== ' IF(lwp) WRITE(numout,*) nstop = nstop + 1 #endif ! 0. Read namelist parameters ! --------------------------- REWIND( numnam ) READ ( numnam, nambdy ) ! control prints IF(lwp) WRITE(numout,*) ' nambdy' IF ((nbdy_dta/=0).AND.(nbdy_dta/=1)) THEN ! Check nbdy_dta value IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' E R R O R : nbdy_dta =',nbdy_dta, & 'but it should have been 0 or 1' IF(lwp) WRITE(numout,*) ' ========== ' IF(lwp) WRITE(numout,*) nstop = nstop + 1 ELSE IF(lwp) WRITE(numout,*) ' ' IF(lwp) WRITE(numout,*) ' data in file (=1) or nbdy_dta = ', nbdy_dta IF(lwp) WRITE(numout,*) ' initial state used (=0)' END IF IF(lwp) WRITE(numout,*) ' ' IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nb_rimwidth = ', nb_rimwidth IF (ln_bdy_vol) THEN IF (volbdy==1) THEN ! Check volbdy value IF(lwp) WRITE(numout,*) ' ' IF(lwp) WRITE(numout,*) ' volbdy = ', volbdy IF(lwp) WRITE(numout,*) ' The total volume will be constant' ELSEIF (volbdy==0) THEN IF(lwp) WRITE(numout,*) ' ' IF(lwp) WRITE(numout,*) ' volbdy = ', volbdy IF(lwp) WRITE(numout,*) ' The total volume will vary according to & &the surface E-P flux' ELSE IF(lwp) WRITE(numout,*) ' ' IF(lwp) WRITE(numout,*) ' E R R O R : volbdy =',volbdy, & 'but it should have been 0 or 1' IF(lwp) WRITE(numout,*) ' ========== ' IF(lwp) WRITE(numout,*) nstop = nstop + 1 END IF ELSE IF(lwp) WRITE(numout,*) ' ' IF(lwp) WRITE(numout,*) 'No volume correction with unstructured open boundaries' IF(lwp) WRITE(numout,*) ' ' ENDIF IF (ln_bdy_fla) THEN IF(lwp) WRITE(numout,*) ' ' IF(lwp) WRITE(numout,*) 'Flather bc with unstructured open boundaries' IF(lwp) WRITE(numout,*) ' ' ELSE IF(lwp) WRITE(numout,*) ' ' IF(lwp) WRITE(numout,*) 'NO Flather bc with unstructured open boundaries' IF(lwp) WRITE(numout,*) ' ' ENDIF ! 0.5 Read tides namelist ! ------------------------ IF ( lk_bdy_tides ) CALL tide_init ! 1. Read arrays defining unstructured open boundaries ! ---------------------------------------------------- ! 1.1 Read global 2D mask at T-points: bdytmask ! ********************************************* ! bdytmask=1 on the computational domain AND on open boundaries ! =0 elsewhere IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN tmpmsk(: , : ) = 0.e0 tmpmsk(jpizoom+1:jpizoom+jpiglo-2,: ) = 1.e0 ELSE IF ( ln_bdy_mask ) THEN CALL iom_open( filbdy_mask, inum ) CALL iom_get ( inum, jpdom_data, 'bdy_msk', tmpmsk(:,:) ) CALL iom_close( inum ) ELSE tmpmsk(:,:) = 1.0 ENDIF ! Save mask over local domain DO jj = 1, nlcj DO ji = 1, nlci bdytmask(ji,jj) = tmpmsk( mig(ji), mjg(jj)) END DO END DO ! Derive mask on U and V grid from mask on T grid bdyumask(:,:)=0.e0 bdyvmask(:,:)=0.e0 DO jj=1, jpjm1 DO ji=1, jpim1 bdyumask(ji,jj)=bdytmask(ji,jj)*bdytmask(ji+1, jj ) bdyvmask(ji,jj)=bdytmask(ji,jj)*bdytmask(ji ,jj+1) END DO END DO ! Lateral boundary conditions CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) ! 1.2 Read discrete distance and mapping indices ! ********************************************** nbidta(:,:)=0. nbjdta(:,:)=0. nbrdta(:,:)=0. IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN icount = 0 ! Define west boundary (from ji=2 to ji=1+nb_rimwidth): DO jr=1,nb_rimwidth DO jj=3,jpjglo-2 icount=icount+1 nbidta(icount,:) = jr + 1 + (jpizoom-1) nbjdta(icount,:) = jj + (jpjzoom-1) nbrdta(icount,:) = jr END DO END DO ! Define east boundary (from ji=jpiglo-1 to ji=jpiglo-nb_rimwidth): DO jr=1,nb_rimwidth DO jj=3,jpjglo-2 icount=icount+1 nbidta(icount,:) = jpiglo-jr + (jpizoom-1) nbidta(icount,2) = jpiglo-jr-1 + (jpizoom-1) ! special case for u points nbjdta(icount,:) = jj + (jpjzoom-1) nbrdta(icount,:) = jr END DO END DO ELSE ! Read indices and distances in unstructured boundary data files IF ( lk_bdy ) THEN bdyfile(1) = filbdy_data_T bdyfile(2) = filbdy_data_U bdyfile(3) = filbdy_data_V ELSE ! In this case we have tides only at the boundaries ! so read index arrays from tides files for first tidal component bdyfile(1) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_T.nc' bdyfile(2) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_U.nc' bdyfile(3) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_V.nc' ENDIF DO jgrd = 1,3 CALL iom_open( bdyfile(jgrd), inum ) dummy_id = iom_varid( inum, 'nbidta', kdimsz=kdimsz ) WRITE(numout,*) 'kdimsz : ',kdimsz nb_len = kdimsz(1) IF (nb_len > jpbdta) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' E R R O R : jpbdta is too small:' IF(lwp) WRITE(numout,*) ' ========== Boundary array length in file is ', nb_len IF(lwp) WRITE(numout,*) ' But jpbdta is ', jpbdta IF(lwp) WRITE(numout,*) ' File : ', bdyfile(jgrd) IF(lwp) WRITE(numout,*) nstop = nstop + 1 ENDIF CALL iom_get ( inum, jpdom_unknown, 'nbidta', ndta(1:nb_len,:) ) DO ji=1,nb_len nbidta(ji,jgrd) = INT( ndta(ji,1) ) ENDDO CALL iom_get ( inum, jpdom_unknown, 'nbjdta', ndta(1:nb_len,:) ) DO ji=1,nb_len nbjdta(ji,jgrd) = INT( ndta(ji,1) ) ENDDO CALL iom_get ( inum, jpdom_unknown, 'nbrdta', ndta(1:nb_len,:) ) DO ji=1,nb_len nbrdta(ji,jgrd) = INT( ndta(ji,1) ) ENDDO CALL iom_close( inum ) ! Check that rimwidth in file is big enough: nbr_max = MAXVAL(nbrdta(:,jgrd)) IF (nbr_max < nb_rimwidth) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' E R R O R : Maximum rimwidth in file is ', nbr_max IF(lwp) WRITE(numout,*) ' ========== but nb_rimwidth is ', nb_rimwidth IF(lwp) WRITE(numout,*) nstop = nstop + 1 ELSE IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', nbr_max IF(lwp) WRITE(numout,*) END IF ENDDO END IF ! 1.3 Dispatch mapping indices and discrete distances 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 jgrd = 1, jpbgrd icount = 0 icountr = 0 DO nb_rim=1, nb_rimwidth DO jb = 1, jpbdta ! check if point is in local domain and equals nb_rim IF ( (nbidta(jb,jgrd) >= iw ).AND. & (nbidta(jb,jgrd) <= ie ).AND. & (nbjdta(jb,jgrd) >= is ).AND. & (nbjdta(jb,jgrd) <= in ).AND. & (nbrdta(jb,jgrd) == nb_rim ) ) THEN icount = icount + 1 IF (nb_rim==1) icountr = icountr+1 IF (icount > jpbdim) THEN IF(lwp) WRITE(numout,*) 'bdy_ini: jpbdim too small' nstop = nstop + 1 ELSE nbi(icount, jgrd) = nbidta(jb,jgrd)- mig(1)+1 nbj(icount, jgrd) = nbjdta(jb,jgrd)- mjg(1)+1 nbr(icount, jgrd) = nbrdta(jb,jgrd) nbmap(icount,jgrd) = jb END IF END IF END DO END DO nblenrim(jgrd) = icountr !: length of rim boundary data on each proc nblen (jgrd) = icount !: length of boundary data on each proc END DO ! 2. Compute rim weights ! ---------------------- DO jgrd = 1, jpbgrd DO jb = 1, nblen(jgrd) ! tanh formulation nbw(jb,jgrd) = 1.-TANH((FLOAT(nbr(jb,jgrd)-1))/2.) ! quadratic ! nbw(jb,jgrd) = (FLOAT(nb_rimwidth+1-nbr(jb,jgrd))/FLOAT(nb_rimwidth))**2 ! linear ! nbw(jb,jgrd) = FLOAT(nb_rimwidth+1-nbr(jb,jgrd))/FLOAT(nb_rimwidth) END DO END DO ! 3. Mask corrections ! ------------------- DO jk=1, jpkm1 DO jj=1, jpj DO ji=1, jpi tmask(ji,jj,jk)=tmask(ji,jj,jk)*bdytmask(ji,jj) umask(ji,jj,jk)=umask(ji,jj,jk)*bdyumask(ji,jj) vmask(ji,jj,jk)=vmask(ji,jj,jk)*bdyvmask(ji,jj) bmask(ji,jj)=bmask(ji,jj)*bdytmask(ji,jj) END DO END DO END DO ! I am not sure that it is useful: DO jk=1, jpkm1 DO jj=2, jpjm1 DO ji=2, jpim1 fmask(ji,jj,jk) = fmask(ji,jj,jk) & & * bdytmask(ji, jj ) * bdytmask(ji+1, jj ) & & * bdytmask(ji,jj+1) * bdytmask(ji+1,jj+1) END DO END DO END DO tmask_i(:,:) = tmask(:,:,1)*tmask_i(:,:) bdytmask(:,:)=tmask(:,:,1) ! bdy masks and bmask are now set to zero on boundary points: jgrd=1 ! In the free surface case, bmask is at T-points DO jb=1, nblenrim(jgrd) bmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 END DO jgrd=1 DO jb=1, nblenrim(jgrd) bdytmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 END DO jgrd=2 DO jb=1, nblenrim(jgrd) bdyumask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 END DO jgrd=3 DO jb=1, nblenrim(jgrd) bdyvmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0 END DO ! Lateral boundary conditions CALL lbc_lnk( fmask, 'F', 1. ) CALL lbc_lnk( bdytmask(:,:), 'T', 1. ) CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) IF ((ln_bdy_vol).OR.(ln_bdy_fla)) THEN ! 4 Indices and directions of rim velocity components ! --------------------------------------------------- !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 icount = 0 flagu(:)=0.e0 jgrd=2 ! u-component DO jb=1, nblenrim(jgrd) efl=bdytmask(nbi(jb,jgrd) , nbj(jb,jgrd)) wfl=bdytmask(nbi(jb,jgrd)+1, nbj(jb,jgrd)) IF ((efl+wfl)==2) THEN icount = icount +1 ELSE flagu(jb)=-efl+wfl END IF 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 flagv(:)=0.e0 jgrd=3 ! v-component DO jb=1, nblenrim(jgrd) nfl = bdytmask(nbi(jb,jgrd), nbj(jb,jgrd)) sfl = bdytmask(nbi(jb,jgrd), nbj(jb,jgrd)+1) IF ((nfl+sfl)==2) THEN icount = icount +1 ELSE flagv(jb)=-nfl+sfl 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.' IF(lwp) WRITE(numout,*) ' ========== ' IF(lwp) WRITE(numout,*) nstop = nstop + 1 END IF END IF ! 5 Compute total lateral surface for volume correction: ! ------------------------------------------------------ bdysurftot = 0.e0 IF (ln_bdy_vol) THEN jgrd=2 ! Lateral surface at U-points DO jb=1, nblenrim(jgrd) bdysurftot = bdysurftot + & hu(nbi(jb,jgrd), nbj(jb,jgrd)) * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) & * ABS(flagu(jb))*tmask_i(nbi(jb,jgrd) , nbj(jb,jgrd)) & *tmask_i(nbi(jb,jgrd)+1, nbj(jb,jgrd)) END DO jgrd=3 ! Add lateral surface at V-points DO jb=1, nblenrim(jgrd) bdysurftot = bdysurftot + & hv(nbi(jb,jgrd), nbj(jb,jgrd)) * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) & * ABS(flagv(jb))*tmask_i(nbi(jb,jgrd), nbj(jb,jgrd)) & *tmask_i(nbi(jb,jgrd), nbj(jb,jgrd)+1) END DO IF( lk_mpp ) CALL mpp_sum( bdysurftot ) ! sum over the global domain END IF ! 6. Initialise bdy data arrays ! ----------------------------- tbdy(:,:) = 0.e0 sbdy(:,:) = 0.e0 ubdy(:,:) = 0.e0 vbdy(:,:) = 0.e0 sshbdy(:) = 0.e0 ubtbdy(:) = 0.e0 vbtbdy(:) = 0.e0 ! 7. Read in tidal constituents and adjust for model start time ! ------------------------------------------------------------- IF ( lk_bdy_tides ) CALL tide_data END SUBROUTINE bdy_init #else !!--------------------------------------------------------------------------------- !! Dummy module NO unstructured open boundaries !!--------------------------------------------------------------------------------- CONTAINS SUBROUTINE bdy_init ! Dummy routine END SUBROUTINE bdy_init #endif !!================================================================================= END MODULE bdyini