c =========================================================================== c BAROTROPIC.f consists of subroutines which solve the elliptic problem c for the barotropic stream function after a coordinate transform c to a uniform grid involving the metric terms emx,emy,emx2 and emy2. c Islands are included using line integrals (Kamenkovich) to c set constant values of the streamfunction on multiple land masses. c The resulting elliptic equations are solved with Zlatev's Y12M c package, by performing an approximate LU factorization and iterating c to a predetermined tolerance level c Naomi Naik, LDEO, September 17, 1996 c =========================================================================== c c 2/24/93 changed call to baro_solv, added iox c 2/25/93 changed call to baro_solv, added psi c 2/25/93 changed call to init_baro, added ibar_key c 3/03/93 got rid of the deph_hx and deph_hy arrays c 3/08/93 changed call to init_baro, added periodic/non-periodic switch c 3/10/93 major revision: c baro_solv c (input) values of forcing at water and boundary points c (output) values of stream function at water and boundry pnts c init_baro c (input) number and list of water+boundary points c (input) number and list of boundary points c mask set according to the above list, NOT by using the c depth values c 3/24/93 changed call to init_baro, including Rayleigh friction and glubina c deleted the 'key' parameter, now controlled from mod.in c 5/18/93 allocation for bound_rhs moved to init_baro c 5/20/93 changed call to baro_solv, including Rayleigh friction c 6/04/93 changed call to init_baro, including nbaro c 6/29/93 variable depth version c 7/12/93 fixed to stabilize third order scheme c 5/12/94 new variable topography c 3/27/95 fixed 2 memory allocation bugs in "init_baro" c 7/12/95 added geometric term to time dependent operator, see "make_rhs" c and "add_dt" c 7/14/95 all input parameters from "mod.in" and ".y12m" put in common c blocks, for future setting from model_input subroutine c 7/17/95 rewrote so that differencing is on a uniform grid of unit meshsize c with the analytic coordinate transformation factors appearing c in the metric terms emx, emy as in main code c 7/19/95 grid stretching implemented c first order derivatives :the stretch factors appear in emx, emy c second order derivatives:new terms were added in add_dt, c make_rhs and do_xi_eta c 12/19/95 grid stretching invoked; baro-parameters are now read and set c in model_input; changed: calls to init_baro & c baro_solv; debug output has been cleaned /Senya/ c 5/14/96 coriolis function in generalized coordinates from c common/data_geom instead of computed separately here c 6/19/96 fixed a bug in periodic use of coriolis function c 6/19/96 created this stripped down version with reduced memory c (second order, five point stencil only) c 6/26/96 fixed the sunken island which was spoiled when we imported c the geometric terms c 9/17/96 changed all island stuff so it can be specified from input file c10/16/96 changed bottom friction parameterization and cosmetic stuff c12/19/96 interior corners changed to land points c-------------------------------------------------------------------------- subroutine init_baro (i_p, eps, nxp, nyp, nxyc, ncs, iox, nbx, lxx, * nby,lyx, alon, blon, alat, blat, x, y, db, glub, gradH, adelt_io) c---------------------- c compute the coefficients of the discrete barotropic equations and c then perform an approximate lu decomposition c ip = (input) 0: non-periodic in x, 1: periodic in x c eps = (input) effective friction induced by time step c nxp = (input) number of x-coordinate values c nyp = (input) number of y-coordinate values c nxyc = (input) number of compressed (water+boundary) points c iox = (input) compressed -> uncompressed c nbx = (input) number of boundary points c lxx = (input) boundary list position -> compressed position c alon = (input) minimum longitude c blon = (input) maximum longitude c alat = (input) minimum latitude c blat = (input) maximum latitude c x = (input) x coordinate values in degrees c y = (input) y coordinate values in degrees c db = (input) depth values c glub = (input) total depth of ocean c-------------------------------------------------------------- dimension f(1), emx(1), emy(1), emx2(1), emy2(1), area(1) dimension sxm(1) dimension cosl(1),cosli(1) dimension saxmk(1), saxpk(1), saymk(1), saypk(1), sxp(1), syp(1) dimension sponge(1), relax(1) pointer (p_f, f), (p_emx, emx), (p_emy, emy), * (p_emx2, emx2), (p_emy2, emy2), (p_area, area) * , (p_relax, relax), (p_sponge, sponge) * , (p_saxmk, saxmk), (p_saymk, saymk) * , (p_saxpk, saxpk), (p_saypk, saypk) * , (p_sxp, sxp), (p_syp, syp) * , (p_cosl, cosl), (p_cosli, cosli) * , (p_sxm, sxm) common/data_geom/ p_f,p_emx,p_emy,p_emx2,p_emy2,p_area * ,p_relax,p_sponge,p_cosl,p_cosli, * p_saxmk,p_saxpk,p_saymk,p_saypk,p_sxp,p_syp,p_sxm c--------------------------------------------------------------------------- include 'barotropic.h' integer iox(1), lxx(1), lyx(1) real db(1), x(1), y(1), gradH(nxyc,2) integer itmp(1) pointer (pitmp, itmp) save pitmp adelt = adelt_io if_per = i_p GLUBINA = glub X_MIN = alon X_MAX = blon Y_MIN = alat Y_MAX = blat NX = nxp + 2*if_per NY = nyp iper = nxp NXY = NX * NY call mem_alloc (pdeph, NXY, 2, 'baro deph') call mem_alloc (pgradHx, NXY, 2, 'baro gradH_x') call mem_alloc (pgradHy, NXY, 2, 'baro gradH_y') call mem_alloc (pbcosli, NXY, 2, 'baro cosli') call mem_alloc (pfcor, NXY, 2, 'baro fcor') call mem_alloc (pbemx, NXY, 2, 'baro emx') call mem_alloc (pbemy, NXY, 2, 'baro emy') call mem_alloc (pbemxx, NXY, 2, 'baro emxx') call mem_alloc (pbemyy, NXY, 2, 'baro emyy') call mem_alloc (pbound_rhs, NXY, 2, 'baro bound rhs') call mem_alloc0 (pmask, NXY, 'baro mask') do i = 1, NXY mask(i) = BC_L enddo do i = 1, nxyc ixy = iox(i) ix = mod(ixy-1,nxp)+1 iy = (ixy-ix)/nxp + 1 i_xy = (iy - 1)*NX + ix + if_per deph(i_xy) = db(i) gradHx(i_xy) = gradH(i,1) gradHy(i_xy) = gradH(i,2) fcor(i_xy) = f(i) bcosli(i_xy) = cosli(i) bemx(i_xy) = emx(i) bemy(i_xy) = emy(i) bemxx(i_xy) = emx2(i) bemyy(i_xy) = emy2(i) mask(i_xy) = BC_W enddo call mark_islands (nxp,nyp,x,y) if (ibar_key .ne. 0) * open (unit = IUNIT_OUT, file = f_bar(1:n_bar)) if (ibar_key .eq. 4) then write(IUNIT_OUT, *) 'X-boundary points: i, ix,iy, lxx(i)' endif do i = 1, nbx ixy = iox(lxx(i)) ix = mod(ixy-1,nxp)+1 iy = (ixy-ix)/nxp + 1 i_xy = (iy - 1)*NX + ix + if_per if (ibar_key .eq. 4) write(IUNIT_OUT,*) i, ix, iy, lxx(i) if ( (ix.gt.1.and.ix.lt.nxp) .or. if_per.eq.0) mask(i_xy) = BC_L enddo do i = nbx+1, nbx+ncs ixy = iox(lxx(i)) ix = mod(ixy-1,nxp)+1 iy = (ixy-ix)/nxp + 1 i_xy = (iy - 1)*NX + ix + if_per mask(i_xy) = BC_L enddo if (ibar_key .eq. 4) then write(IUNIT_OUT,*)'interior points + X-boundary points 1/2:' do iy = NY, 1, -1 write(IUNIT_OUT,'(200a1)') (mask(ix+(iy-1)*NX), ix=1,NX/2) enddo write(IUNIT_OUT,*)'interior points + X-boundary points 2/2:' do iy = NY, 1, -1 write(IUNIT_OUT,'(200a1)') (mask(ix+(iy-1)*NX), ix=NX/2,NX) enddo write(IUNIT_OUT,*) 'Y-boundary points: i, ix, iy, lyx(i)' endif do i = 1, nby ixy = iox(lyx(i)) ix = mod(ixy-1,nxp)+1 iy = (ixy-ix)/nxp + 1 i_xy = (iy - 1)*NX + ix + if_per if (ibar_key .eq. 4) write(IUNIT_OUT,*) i, ix, iy, lyx(i) mask(i_xy) = BC_L enddo call do_open (eps) call do_init(x,y) call make_gra if ( ibar_key .eq. 4 ) then call mem_alloc1 (pitmp, NXY, 'BARO iox:') do i = 1, nxyc ixy = iox(i) itmp(ixy) = i enddo do iy = NY, 1, -1 write(IUNIT_OUT, '(200i3)') (mod(itmp(ix+(iy-1)*nxp),1000), ix=1,NX) enddo call mem_free (pitmp, NXY, 1) endif do kis=0,n_islands in_rhs(kis) = NPACK*kis enddo call mem_alloc (prelx_p, NPACK, 2, 'baro relx_p') call mem_alloc (prelx_m, NPACK, 2, 'baro relx_m') call mem_alloc (prely_p, NPACK, 2, 'baro rely_p') call mem_alloc (prely_m, NPACK, 2, 'baro rely_m') do i = 1, NPACK i_xy = list(i) ix = mod (i_xy -1 ,NX) + 1 iy = (i_xy - ix)/NX + 1 relx_p(i) = 1. relx_m(i) = 1. rely_p(i) = 1. rely_m(i) = 1. dep = deph(i_xy) if (deph(i_xy+1).ge.BAR_DEPMIN) relx_p(i) = dep/deph(i_xy+1) if (mask(i_xy+1).eq.BC_P) relx_p(i) = dep/deph(i_xy+1-iper) if (deph(i_xy-1).ge.BAR_DEPMIN) relx_m(i) = dep/deph(i_xy-1) if (mask(i_xy-1).eq.BC_P) relx_m(i) = dep/deph(i_xy-1+iper) if (deph(i_xy+NX).ge.BAR_DEPMIN) rely_p(i) = dep/deph(i_xy+NX) if (deph(i_xy-NX).ge.BAR_DEPMIN) rely_m(i) = dep/deph(i_xy-NX) enddo call mem_alloc(psol, NXY, 2, 'baro sol') call mem_alloc(prhs_bc, NPACK, 2, 'baro rhs_bc') nisl = NPACK*(n_islands + 1) call mem_alloc (pb1, NPACK, 2, 'baro b1') call mem_alloc (psol1, NPACK, 2, 'baro sol1') call mem_alloc (ptmp, NXY, 2, 'baro tmp') call mem_alloc (prhs, NPACK, 2, 'baro rhs') call mem_alloc (prhs1, NPACK, 2, 'baro rhs1') call mem_alloc (ptaux, NPACK, 2, 'baro taux') call mem_alloc (ptauy, NPACK, 2, 'baro tauy') call mem_alloc(prhs_bci, nisl, 2, 'baro rhs_bci') call mem_alloc(prhsi, nisl, 2, 'baro rhsi') call make_matr call make_lu call mem_free (prhsi, nisl, 2) call mem_free (piro, NN12, 1) if (ibar_key .ne. 0) call flush(IUNIT_OUT) return end subroutine baro_solv (nxp, nyp, nxyc, iox, tb, txb, tyb, psi) c---------------------- c compute the barotropic transports, ub, vb, given the depth and c the average forcing on a uniform grid, using the lu decomposition c from init_baro c nxp = (input) number of x-coordinate values c nyp = (input) number of y-coordinate values c nxyc = (input) number of compressed (water) points c iox = (input) compressed -> uncompressed c txb = (input) barotropic forcing, x-direction c tyb = (input) barotropic forcing, y-direction c tb = (input) curl of the barotropic forcing c psi =(output) barotropic stream function c---------------------- include 'barotropic.h' real tb(1), txb(1), tyb(1), psi(0:nxp+1,0:nyp+1) integer iox(1) call do_winds(nxp,nyp,nxyc,iox,txb,tyb,psi) call make_rhs (nxp,nxyc,iox,tb) call do_solve do kis = 1, n_islands if (ityp_island(kis).eq.2) then nsk = nsunk(kis) spsi = 0. do i = 1, nsk i_xy = isunk(i,kis) i_pac = ilst(i_xy) spsi = spsi + rhs(i_pac) enddo spsi = spsi/nsk b_island(kis) = spsi do i = 1, nsk i_xy = isunk(i,kis) i_pac = ilst(i_xy) rhs(i_pac) = spsi bound_rhs(i_xy) = spsi enddo endif enddo do iy = 1, ny do ix = 1, nx i_xy = (iy - 1)*NX + ix + if_per if (mask(i_xy).eq.BC_W) then i_pac = ilst(i_xy) dpsi = rhs(i_pac) endif if (mask(i_xy).eq.BC_P) then if (ix.eq.1) then i_pac = ilst(i_xy+iper) dpsi = rhs(i_pac) else i_pac = ilst(i_xy-iper) dpsi = rhs(i_pac) endif endif if (mask(i_xy).eq.BC_L) dpsi = 0. if (mask(i_xy).eq.BC_0) dpsi = b_island(0) do kis = 1, n_islands if (mask(i_xy).eq.BC_I(kis)) dpsi = b_island(kis) enddo psi(ix,iy) = dpsi enddo enddo if (ibar_key .ne. 0) call flush(IUNIT_OUT) return end subroutine do_open (eps) c------------------------------------------- include 'barotropic.h' CNST_EPS = rayl CNST_EPT = eps if ( ibar_key .gt. 0 ) then write (IUNIT_OUT, *) write (IUNIT_OUT, *) ' BAROTROPIC SOLVER INFO' write (IUNIT_OUT, *) ' ----------------------' write (IUNIT_OUT, *) write (IUNIT_OUT, *) 'NX, NY =', NX, NY write (IUNIT_OUT, *) 'X_MIN, X_MAX =', X_MIN,X_MAX write (IUNIT_OUT, *) 'Y_MIN, Y_MAX =', Y_MIN,Y_MAX write (IUNIT_OUT, *) 'epsilon =', CNST_EPT write (IUNIT_OUT, *) 'Rayleigh friction=', CNST_EPS endif return end subroutine do_init(x,y) c------------------------- include 'barotropic.h' dimension x(1),y(1) character*1 ma c set various repeated-use vectors call extend (x,y) if (ibar_key.eq.3) then if (NX.gt.80) then write(IUNIT_OUT, *) 'masks 1/2 :' do iy = NY, 1, -1 write(IUNIT_OUT,'(200a1)') (mask(ix+(iy-1)*NX), ix=1,NX/2) enddo write(IUNIT_OUT, *) 'masks 2/2 :' do iy = NY, 1, -1 write(IUNIT_OUT,'(200a1)') (mask(ix+(iy-1)*NX), ix=NX/2,NX) enddo else write(IUNIT_OUT, *) 'masks:' do iy = NY, 1, -1 write(IUNIT_OUT,'(200a1)') (mask(ix+(iy-1)*NX), ix=1,NX) enddo endif endif if (ibar_key.eq.5) then write(IUNIT_OUT, *) 'masks:' do iy = NY, 1, -1 write(IUNIT_OUT,'(200a1)') (mask(ix+(iy-1)*NX), ix=1,NX) enddo if (use_per_island) then iy = ib do ix = 2, NX-1 mask(ix+(iy-1)*NX) = 'x' enddo endif do kis = 1, n_islands do ipoly = 1,n_poly(kis) iy = j_poly(ipoly,kis) i_min = i_poly(ipoly,kis) i_max = i_poly(ipoly+1,kis) j_min = j_poly(ipoly,kis) j_max = j_poly(ipoly+1,kis) istep = 1 if (i_max.lt.i_min) istep = -1 jstep = 1 if (j_max.lt.j_min) jstep = -1 if (i_min.ne.i_max) then do ix = i_min,i_max,istep ma = mask(ix+(iy-1)*NX) if (ma .eq. 'x'.or.ma .eq. 'o' ) then mask(ix+(iy-1)*NX) = 'o' else mask(ix+(iy-1)*NX) = 'x' endif enddo endif ix = i_min if (j_min.ne.j_max) then do iy = j_min,j_max,jstep ma = mask(ix+(iy-1)*NX) if (ma .eq. 'x'.or.ma .eq. 'o' ) then mask(ix+(iy-1)*NX) = 'o' else mask(ix+(iy-1)*NX) = 'x' endif enddo endif enddo enddo write(IUNIT_OUT, *) 'masks:' do iy = NY, 1, -1 write(IUNIT_OUT,'(200a1)') (mask(ix+(iy-1)*NX), ix=1,NX) enddo stop endif call do_sink if (ibar_key.eq.3) then if (NX.gt.80) then write(IUNIT_OUT, *) 'after sinking, masks 1/2 :' do iy = NY, 1, -1 write(IUNIT_OUT,'(200a1)') (mask(ix+(iy-1)*NX), ix=1,NX/2) enddo write(IUNIT_OUT, *) 'after sinking, masks 2/2 :' do iy = NY, 1, -1 write(IUNIT_OUT,'(200a1)') (mask(ix+(iy-1)*NX), ix=NX/2,NX) enddo else write(IUNIT_OUT, *) 'after sinking, masks:' do iy = NY, 1, -1 write(IUNIT_OUT,'(200a1)') (mask(ix+(iy-1)*NX), ix=1,NX) enddo endif endif end subroutine do_sink c------------------------- c fictitious domain method c------------------------------------- include 'barotropic.h' do kis = 1, n_islands nsunk(kis) = 0 if (ityp_island(kis).eq.2) then bci = BC_I(kis) nsk = 0 do iy = 1, NY do ix = 2, NX i_xy = (iy - 1) * NX + ix if (mask(i_xy).eq.bci) then nsk = nsk+1 isunk(nsk,kis)=i_xy deph(i_xy) = BAR_DSINK mask(i_xy) = BC_W c the rest is a cludge - just to fill in these values fcor(i_xy) = fcor(i_xy - 1) bemx(i_xy) = bemx(i_xy - 1) bemy(i_xy) = bemy(i_xy - 1) bemxx(i_xy) = bemxx(i_xy - 1) bemyy(i_xy) = bemyy(i_xy - 1) endif enddo enddo if (nsk.gt.1000) then print*,'increase the size of isunk in barotropic.h' stop endif nsunk(kis) = nsk endif enddo end subroutine extend (x,y) c------------------------- include 'barotropic.h' dimension x(1),y(1) c extend east/west boundaries if (if_per .eq. 1) then ix = 1 do iy = 1, NY i_xy = (iy - 1) * NX + ix mask(i_xy) = mask(i_xy+iper) enddo ix = NX do iy = 1, NY i_xy = (iy - 1) * NX + ix mask(i_xy) = mask(i_xy-iper) enddo ix = NX do iy = 2, NY-1 i_xy = (iy - 1) * NX + ix if (mask(i_xy).eq.BC_W) then mask(i_xy) = BC_P mask(i_xy-NX+1) = BC_P elseif (mask(i_xy).eq.BC_L) then mask(i_xy-NX+1) = BC_L endif enddo endif if (use_per_island) then ib = 0 do iy = 2, NY if (y(iy-1).lt.per_lat.and.per_lat.le.y(iy)) ib = iy enddo iy = ib do ix = 1, NX i_xy = (iy - 1) * NX + ix if (mask(i_xy).ne.BC_W.and.mask(i_xy).ne.BC_P) then print*,'trouble with line integral for periodic island' print*,'at:',ix,iy stop endif enddo 10 continue do ix = 1, NX do iy = 1, NY i_xy = (iy - 1) * NX + ix dlat = y(iy) if ((mask(i_xy).eq.BC_L) * .and. dlat.lt.per_lat) mask(i_xy) = BC_0 enddo enddo endif if (n_islands.gt.0) then do k_island = 1, n_islands bct = 'T' bci = BC_I(k_island) do ix = 1, NX do iy = 1, NY itp = (iy - 1) * NX + ix ipc = min(itp+1,nx*ny) imc = max(itp-1, 1) icp = min(itp+NX,nx*ny) icm = max(itp-NX, 1) ipp = max(min(itp+1+NX,nx*ny),1) imp = max(min(itp-1+NX,nx*ny),1) ipm = max(min(itp+1-NX,nx*ny),1) imm = max(min(itp-1-NX,nx*ny),1) if (mask(itp).eq.BC_L) then if (mask(ipc).eq.bci.or.mask(imc).eq.bci.or. * mask(icp).eq.bci.or.mask(icm).eq.bci.or. * mask(ipp).eq.bci.or.mask(ipm).eq.bci.or. * mask(imp).eq.bci.or.mask(imm).eq.bci * ) mask(itp) = bct endif enddo enddo do ix = 1, NX do iy = 1, NY itp = (iy - 1) * NX + ix if (mask(itp).eq.bct) then mask(itp) = bci endif enddo enddo enddo endif return end subroutine mark_islands (nxp,nyp,x,y) c------------------------- include 'barotropic.h' dimension x(1),y(1) if (n_islands.gt.0) then do k_island = 1, n_islands npoly = n_poly(k_island) alon_poly(npoly+1,k_island) = alon_poly(1,k_island) alat_poly(npoly+1,k_island) = alat_poly(1,k_island) c find i_poly and j_poly for each segment xis = alon_poly(1,k_island) yis = alat_poly(1,k_island) d_poly = 2*(x(nxp)-x(1)) do ix = 2, nxp dis = abs(x(ix)-xis) if (dis.lt.d_poly) then i_poly(1,k_island) = ix d_poly = dis endif enddo d_poly = 2*(y(nyp)-y(1)) do iy = 2, nyp dis = abs(y(iy)-yis) if (dis.lt.d_poly) then j_poly(1,k_island) = iy d_poly = dis endif enddo do ipoly = 2, npoly yism = yis xis = alon_poly(ipoly,k_island) yis = alat_poly(ipoly,k_island) c find closest i,j value to each segment if (yis.eq.yism) then j_poly(ipoly,k_island) = j_poly(ipoly-1,k_island) i_poly(ipoly,k_island) = 0 do ix = 2, nxp if (x(ix-1).lt.xis.and.xis.le.x(ix)) i_poly(ipoly,k_island) = ix enddo else i_poly(ipoly,k_island) = i_poly(ipoly-1,k_island) j_poly(ipoly,k_island) = 0 do iy = 2, nyp if (y(iy-1).lt.yis.and.yis.le.y(iy)) j_poly(ipoly,k_island) = iy enddo endif enddo i_poly(npoly+1,k_island) = i_poly(1,k_island) j_poly(npoly+1,k_island) = j_poly(1,k_island) c mark all island points xis = alon_island(k_island) yis = alat_island(k_island) iis = 0 do ix = 2, nxp if (x(ix-1).lt.xis.and.xis.le.x(ix)) iis = ix enddo jis = 0 do iy = 2, nyp if (y(iy-1).lt.yis.and.yis.le.y(iy)) jis = iy enddo c print*,'center of island:',k_island,' is:', iis,jis i_xy = (jis - 1) * NX + iis bci = BC_I(k_island) mask(i_xy) = bci i_distance = max(nxp-iis,iis,nyp-jis,jis) do id = 1,i_distance idxm = max(2,iis-id) idxp = min(nxp-1,iis+id) idym = max(2,jis-id) idyp = min(nyp-1,jis+id) do idx = idxm,idxp do idy = idym,idyp itp = (idy - 1) * NX + idx if (mask(itp).eq.BC_L) then if (mask(itp+1).eq.bci.or.mask(itp-1).eq.bci.or. * mask(itp+nx).eq.bci.or.mask(itp-nx).eq.bci) mask(itp) = bci endif enddo enddo enddo enddo endif end subroutine make_matr c------------------------ c computes approximate matrix norm c sets up sparse matrix aa, using stencils computed c in do_elem* c scales all matrix entries in aa by approx. matrix norm c------------------------ include 'barotropic.h' common /band_local/ iband, iaa hx_ave = 0. hy_ave = 0. do i = 1, NPACK i_xy = list(i) hx_ave = hx_ave + 1./bemx(i_xy) hy_ave = hy_ave + 1./bemy(i_xy) enddo hx_ave = hx_ave/NPACK hy_ave = hy_ave/NPACK CNST_2OMEGA = 2. * 2. * PI_MATH / (24. * 3600.) ! 2 * 2*pi/day phi_0 = PI_MATH/ 180.* (Y_MIN + Y_MAX) / 2. CNST_BETA = CNST_2OMEGA * cos(phi_0) CNST_NORM * = (CNST_EPT + CNST_EPS/GLUBINA) * (1./hx_ave**2 + 1./hy_ave**2) + * CNST_BETA / hx_ave / R_EARTH do i_xy = 1, NXY sol(i_xy) = 0.0 enddo do i = 1, NPACK rhs_bc(i) = 0. enddo if (use_per_island) then do i = 1, NPACK rhs_bci(i) = 0. enddo endif do kis = 1, n_islands ind = in_rhs(kis) do i = 1, NPACK rhs_bci(i+ind) = 0. enddo enddo iaa = 0 do i_pac = 1, NPACK i_xy = list(i_pac) call do_xi_eta (i_pac,i_xy,xi,eta) call do_elem(i_pac,i_xy,xi,eta) call add_dt(i_pac, i_xy) iband = 0 call assem (i_pac,i_xy) enddo if (iaa .ne. NONZ) then print*, 'BARO: The matrix has the wrong size !!!' stop endif end subroutine assem (i_pac,i_xy) c------------------------------ include 'barotropic.h' real*8 ce, cw, cn, cs, csum common /matr_local/ ce, cw, cn, cs, csum call coef_matr (i_pac, i_xy, csum) if (mask(i_xy+1).eq.BC_P) then call coef_matr (i_pac, i_xy + 1 - iper, ce) else call coef_matr (i_pac, i_xy + 1, ce) endif if (mask(i_xy-1).eq.BC_P) then call coef_matr (i_pac, i_xy - 1 + iper, cw) else call coef_matr (i_pac, i_xy - 1, cw) endif call coef_matr (i_pac, i_xy + NX, cn) call coef_matr (i_pac, i_xy - NX, cs) end subroutine do_elem(i_pac,i_xy,xi,eta) c----------------------------- include 'barotropic.h' real*8 ce, cw, cn, cs, csum common /matr_local/ ce, cw, cn, cs, csum real*8 hinv, gam, am, ac, ap, heps real*8 bm, bc, bp c ............................. x - direction hinv = bemx(i_xy) gam = xi/ hinv dep = deph(i_xy) edep = CNST_EPS/dep c edep = 0. am = edep - (gam-abs(gam))/2. ap = edep + (gam+abs(gam))/2. ac = -(am+ap) c_h = hinv* hinv/ CNST_NORM ap = c_h * ap am = c_h * am ac = c_h * ac c ............................. y - direction hinv = bemy(i_xy) gam = eta/ hinv bm = edep - (gam-abs(gam))/2. bp = edep + (gam+abs(gam))/2. bc = -(bm+bp) c_h = hinv* hinv/ CNST_NORM bp = c_h * bp bm = c_h * bm bc = c_h * bc c ............................. construct stencil cw = am ce = ap cs = bm cn = bp csum = ac + bc end subroutine add_dt (i_pac, i_xy) c----------------------------- include 'barotropic.h' common /matr_local/ ce, cw, cn, cs, csum real*8 ce, cw, cn, cs, csum dep = deph(i_xy) const = 0.5* CNST_EPT/ CNST_NORM bex = bemx(i_xy) bey = bemy(i_xy) hxe = const*bex**2 hye = const*bey**2 hxx = 2.*const*bemxx(i_xy)*bex hyy = 2.*const*bemyy(i_xy)*bey cea = hxe * (1. + relx_p(i_pac)) cwa = hxe * (1. + relx_m(i_pac)) cna = hye * (1. + rely_p(i_pac)) csa = hye * (1. + rely_m(i_pac)) c if (hxx.gt.0) then c cea = cea + hxx c else c cwa = cwa - hxx c endif c if (hyy.gt.0) then c cna = cna + hyy c else c csa = csa - hyy c endif ce = ce + cea cw = cw + cwa cn = cn + cna cs = cs + csa csum = csum - (cea + cwa + cna + csa) end subroutine do_xi_eta (ip,i_xy,xi,eta) c----------------------------- include 'barotropic.h' bex = bemx(i_xy) bey = bemy(i_xy) dep = deph(i_xy) if (dep.ge.BAR_DEPMIN) then fc = fcor(i_xy) fcm = fcor(i_xy-1) fcp = fcor(i_xy+1) if (mask(i_xy-1).eq.BC_P) fcm = fcor(i_xy-1+iper) if (mask(i_xy+1).eq.BC_P) fcp = fcor(i_xy+1-iper) alfa = bex* (fcp - fcm)/ 2. beta = bey* (fcor(i_xy+NX) - fcor(i_xy-NX))/ 2. xi = beta - fc*gradHy(i_xy) eta = fc*gradHx(i_xy) - alfa else xi = 0. eta= 0. endif return end subroutine coef_matr (i_pac, i_xy, elem) c---------------------------------------- include 'barotropic.h' real*8 elem common /band_local/ iband, iaa solu = 1. if (mask(i_xy) .eq. BC_W) then c if ocean iaa = iaa + 1 aa(iaa) = real(elem) return endif if (mask(i_xy) .eq. BC_L) then c if outer boundary of ocean iband = 1 rhs_bc(i_pac)=rhs_bc(i_pac) - elem * sol(i_xy) return endif if (mask(i_xy) .eq. BC_0) then rhs_bci(i_pac)=rhs_bci(i_pac) - elem * solu return endif do kis = 1, n_islands if (mask(i_xy) .eq. BC_I(kis)) then ind = i_pac + in_rhs(kis) rhs_bci(ind)=rhs_bci(ind) - elem * solu return endif enddo end subroutine make_rhs (nxp,nxyc,iox,tb) c------------------------------------------- include 'barotropic.h' real tb(1) integer iox(1) do i = 1, NXY tmp(i) = 0. enddo do i = 1, nxyc ixy = iox(i) ix = mod(ixy-1,nxp)+1 iy = (ixy-ix)/nxp + 1 i_xy = (iy - 1)*NX + ix + if_per tmp(i_xy) = tb(i) enddo do kis = 1, n_islands do i = 1, nsunk(kis) i_xy = isunk(i,kis) tmp(i_xy) = 0. enddo enddo call add2rhs do i_pac = 1, NPACK i_xy = list(i_pac) bound_rhs(i_xy) = (rhs(i_pac) + tmp(i_xy) ) / CNST_NORM enddo return end subroutine add2rhs c----------------------------- include 'barotropic.h' do ip = 1, NPACK i_xy = list(ip) c construct time derivative term, matches that in add_dt: bw = bound_rhs(i_xy-1) be = bound_rhs(i_xy+1) if (mask(i_xy-1).eq.BC_P) bw = bound_rhs(i_xy-1+iper) if (mask(i_xy+1).eq.BC_P) be = bound_rhs(i_xy+1-iper) bn = bound_rhs(i_xy+NX) bs = bound_rhs(i_xy-NX) bc = bound_rhs(i_xy) bex = bemx(i_xy) bex2 = bex*bex bey = bemy(i_xy) bey2 = bey*bey bexx = bemxx(i_xy) beyy = bemyy(i_xy) dep = deph(i_xy) const = 0.5 hxe = const* bex2 hye = const* bey2 hxx = bexx*bex hyy = beyy*bey ce = hxe * (1. + relx_p(ip)) cw = hxe * (1. + relx_m(ip)) cn = hye * (1. + rely_p(ip)) cs = hye * (1. + rely_m(ip)) c if (hxx.gt.0) then c ce = ce + hxx c else c cw = cw - hxx c endif c if (hyy.gt.0) then c cn = cn + hyy c else c cs = cs - hyy c endif csum = - (ce + cw + cn + cs) btemp= CNST_EPT*(cw*bw + ce*be + cn*bn + cs*bs + csum*bc) c construct bottom drag term, to match implicit term in do_xi_eta,do_elem c if (CNST_EPS.lt.1e-10) then atemp= 0. c else c if (dep.ge.BAR_DEPMIN) then c xi = bexx c eta= beyy + bexy c else ! sunken island c xi = 0. c eta= 0. c endif c c gam = xi/ bex c am = 1. - (gam-abs(gam))/2. c ap = 1. + (gam+abs(gam))/2. c ce = bex2 * ap c cw = bex2 * am c c gam = eta/ bey c am = 1. - (gam-abs(gam))/2. c ap = 1. + (gam+abs(gam))/2. c cn = bey2 * ap c cs = bey2 * am c c csum = - (ce + cw + cn + cs) c atemp= CNST_EPS*(cw*bw + ce*be + cn*bn + cs*bs + csum*bc)/dep c endif c add terms rhs(ip)= btemp + atemp enddo return end subroutine make_lu c-------------------------- c The contents of NPACK, NONZ, aa, ico, NN12, pivot, sn, a1, c columns 1, 3, 4, 6, 7, 8 and 11 of HA, c AFLAG(6), IFLAG(1), IFLAG(4) and IFAIL should c not be altered between calls of y12mfe c-------------------------- include 'barotropic.h' external icpu_time data M12_HA /13/ call mem_alloc (pa1, NONZ, 2, 'baro a1') call mem_alloc (psn, NONZ, 2, 'baro sn') call mem_alloc (pha, NPACK*M12_HA, 2, 'baro ha') call mem_alloc (ppivot, NPACK, 2, 'baro pivot') c modify rhs to include non-zero dirichlet boundary conditions if (use_per_island) then do i = 1, NPACK rhsi(i) = rhs_bci(i) enddo endif NN12 = NONZ * RM12_NN call mem_realloc (piro, NGRAPH, NN12, NONZ, 1) call mem_realloc (pico, NGRAPH, NN12, NONZ, 1) call mem_realloc (paa, NGRAPH, NN12, NONZ, 2) itime = icpu_time() if ( ibar_key.ge.3 ) then write(IUNIT_OUT, *) write(IUNIT_OUT, *) 'LU factorization:' write(IUNIT_OUT, *) 'NPACK, NONZ, NN12 = ', NPACK, NONZ, NN12 endif if (ibar_key .eq. 11) then ! output the matrix for separate analyses c open (unit=12, file = 'matr', form = 'unformatted') open (unit=12, file = 'matr', form = 'formatted') write(12,*) NPACK, NONZ, NN12 write(12,*) (iflag(i),i=1,8) write(12,*) (aflag(i),i=1,8) write(12,*) 'aa:' write(12,*) (aa(i),i=1,NONZ) write(12,*) 'iro:' write(12,*) (iro(i),i=1,NONZ) write(12,*) 'ico:' write(12,*) (ico(i),i=1,NONZ) close(12) stop endif iflag(5) = 2 call y12mfe(NPACK, aa, ico, NN12, iro, NN12, a1, sn, NONZ, * ha, NPACK, rhsi, b1, sol1, pivot, aflag, iflag, ifail) iflag(5) = 3 if (ibar_key.ge.3) then write(IUNIT_OUT, *) 'growth factor =',aflag(5) write(IUNIT_OUT, *) 'minimal pivotal element =',aflag(8) write(IUNIT_OUT, *) 'max number of non-zero elements needed ', * 'in array aa =',iflag(8) write(IUNIT_OUT, *) 'if this last number was much smaller than' * ,NN12,', reduce NN12' endif if (ifail .ne. 0) goto 10 itime = (icpu_time() - itime) / 1000000 if (ibar_key.ge.3) then write(IUNIT_OUT, *) write(IUNIT_OUT, *) 'CPU time for LU factorization was',itime,' sec' endif if (use_per_island) then if (ibar_key.ge.3) then write(IUNIT_OUT, *) write(IUNIT_OUT, *)"ANTARCTICA" write(IUNIT_OUT, *)'number of iterations performed =',iflag(12) write(IUNIT_OUT, *)'max-norm of last correction vector =',aflag(9) write(IUNIT_OUT, *)'max-norm of last residual vector =',aflag(10) write(IUNIT_OUT, *)'max-norm of corrected solution vector =',aflag(11) endif do i = 1, NPACK rhsi(i) = sol1(i) enddo endif do kis = 1, n_islands if (ityp_island(kis).ne.2) then ind = in_rhs(kis) do i = 1, NPACK rhsi(i+ind) = rhs_bci(i+ind) enddo call y12mfe (NPACK, aa, ico, NN12, iro, NN12, a1, sn, NONZ, * ha, NPACK, rhsi(ind+1), b1, sol1, pivot, aflag, iflag, ifail) if (ibar_key.ge.3) then write(IUNIT_OUT, *) write(IUNIT_OUT, *)"ISLAND:",kis write(IUNIT_OUT, *)'number of iterations performed =',iflag(12) write(IUNIT_OUT, *)'max-norm of last correction vector =',aflag(9) write(IUNIT_OUT, *)'max-norm of last residual vector =',aflag(10) write(IUNIT_OUT, *)'max-norm of corrected solution vector =',aflag(11) endif do i = 1, NPACK rhsi(i+ind) = sol1(i) enddo endif enddo 10 if (ifail .ne. 0) then print*, 'BARO (1): !!! Error code from Y12M:',ifail stop endif if (use_per_island.or.use_stan_island) call do_island_init end subroutine do_solve c-------------------------- include 'barotropic.h' c modify rhs to include non-zero dirichlet boundary conditions c on outer boundary do i = 1, NPACK i_xy = list(i) rhs1(i) = bound_rhs(i_xy) + rhs_bc(i) enddo if (ibar_key.ge.3) then write(IUNIT_OUT, *) write(IUNIT_OUT, *) 'Baro_solv:' endif c solve once to get the solution for zero boundary condition, c in order to determine b_island call y12mfe(NPACK, aa, ico, NN12, iro, NN12, a1, sn, NONZ, * ha, NPACK, rhs1, b1, sol1, pivot, aflag, iflag, ifail) if (ifail .ne. 0) goto 10 if (ibar_key.ge.3) then write(IUNIT_OUT, *)'number of iterations performed =', iflag(12) write(IUNIT_OUT, *)'max-norm of last correction vector =', aflag(9) write(IUNIT_OUT, *)'max-norm of last residual vector =', aflag(10) write(IUNIT_OUT, *)'max-norm of corrected solution vector =',aflag(11) endif do i = 1, NPACK rhs(i) = sol1(i) ! the solutions rhs1(i) = b1(i) ! the righthand side enddo c correct for the island influence, using a_island if (use_per_island.or.use_stan_island) then call do_island_integral if (use_per_island) then do i = 1, NPACK rhs1(i) = rhs1(i) + rhs_bci(i)*b_island(0) enddo endif do kis = 1, n_islands if (ityp_island(kis).ne.2) then ind = in_rhs(kis) do i = 1, NPACK rhs1(i) = rhs1(i) + rhs_bci(i+ind)*b_island(kis) enddo endif enddo call y12mfe (NPACK, aa, ico, NN12, iro, NN12, a1, sn, NONZ, * ha, NPACK, rhs1, b1, sol1, pivot, aflag, iflag, ifail) if (ibar_key.ge.2) then write(IUNIT_OUT, *) write(IUNIT_OUT, *)"correcting for island influence" write(IUNIT_OUT, *)'number of iterations performed =', iflag(12) write(IUNIT_OUT, *)'max-norm of last correction vector =', aflag(9) write(IUNIT_OUT, *)'max-norm of last residual vector =', aflag(10) write(IUNIT_OUT, *)'max-norm of corrected solution vector =',aflag(11) endif do i = 1, NPACK rhs(i) = sol1(i) enddo endif if ( ibar_key .ge. 1 ) then write(IUNIT_OUT, *) 'setting CNST_T = ', * (b_island(kis), kis = 0,n_islands) endif 10 if (ifail .ne. 0) then print*, 'BARO (2): !!! Error code from Y12M:',ifail stop endif end subroutine make_gra c------------------------- include 'barotropic.h' call init_pack call to_y12m end subroutine init_pack c--------------------------- include 'barotropic.h' NPACK = 0 do i = 1, NXY if (mask(i) .eq. BC_W) NPACK = NPACK + 1 enddo call mem_alloc1 (plist, NPACK, 'baro list') call mem_alloc1 (pilst, NXY, 'baro ilst') n = 0 do i = 1, NXY if (mask(i) .eq. BC_W) then n = n + 1 list(n) = i ilst(i) = n endif enddo end subroutine to_y12m c------------------------ include 'barotropic.h' NGRAPH = NPACK * 5 call mem_alloc1 (piro, NGRAPH, 'baro iro') call mem_alloc1 (pico, NGRAPH, 'baro ico') call mem_alloc (paa, NGRAPH, 2, 'baro aa') NONZ = 0 do k = 1, NPACK call do_adj1(k) enddo end subroutine do_adj1 (k) c--------------------------------------- include 'barotropic.h' i = list(k) call add_to_graph (k, i) if (mask(i+1) .eq. BC_P) then call add_to_graph (k, i + 1 - iper) else call add_to_graph (k, i + 1) endif if (mask(i-1) .eq. BC_P) then call add_to_graph (k, i - 1 + iper) else call add_to_graph (k, i - 1) endif call add_to_graph (k, i + NX) call add_to_graph (k, i - NX) end subroutine add_to_graph (k, i) c-------------------------------- include 'barotropic.h' ii = ilst(i) if (ii .ne. 0) then NONZ = NONZ + 1 iro(NONZ) = k ico(NONZ) = ii endif end subroutine do_winds (nxp,nyp,nxyc,iox,txb,tyb,psi) c---------------------- include 'barotropic.h' real txb(1), tyb(1), psi(0:nxp+1,0:nyp+1) integer iox(1) do iy = 1, nyp do ix = 1, nxp i_xy = (iy - 1)*NX + ix + if_per bound_rhs(i_xy) = psi(ix,iy) enddo enddo if (use_per_island.or.n_islands.gt.0) then do i = 1, nxyc ixy = iox(i) ix = mod(ixy-1,nxp)+1 iy = (ixy-ix)/nxp + 1 i_xy = (iy - 1)*NX + ix + if_per i_pac = ilst(i_xy) if (i_pac.ne.0) then taux(i_pac) = txb(i) tauy(i_pac) = tyb(i) endif enddo do i_pac = 1, NPACK i_xy = list(i_pac) ix = mod(i_xy-1, NX) + 1 iy = (i_xy-ix)/NX + 1 bex = 0.5*bemx(i_xy) bey = 0.5*bemy(i_xy) be = bound_rhs(i_xy+1) bw = bound_rhs(i_xy-1) if (mask(i_xy+1).eq.BC_P) be = bound_rhs(i_xy+1-iper) if (mask(i_xy-1).eq.BC_P) bw = bound_rhs(i_xy-1+iper) bn = bound_rhs(i_xy+NX) bs = bound_rhs(i_xy-NX) epsd = CNST_EPT taux(i_pac) = taux(i_pac) - bey* epsd*(bn - bs) tauy(i_pac) = tauy(i_pac) + bex* epsd*(be - bw) enddo else do i_pac = 1, NPACK taux(i_pac) = 0. tauy(i_pac) = 0. enddo endif end subroutine do_island_init c------------------------- include 'barotropic.h' character*1 mayp,maym,maxp,maxm,bcij real coef(0:9) common /island_mat/ ipvt(10),z(10) common /test_time/ t_day,t_ramp do jis = 0, 9 do kis = 0, 9 a_island(jis,kis) = 0. enddo enddo do kis = 0, 9 a_island(kis,kis) = 1. coef(kis) = 0. enddo c line integrals - trapezoid rule on uniform grid if (use_per_island) then iy = ib do ix = 2, NX-1 i_xy = (iy - 1) * NX + ix ipac = ilst(i_xy) ipyp = ilst(i_xy + NX) ipym = ilst(i_xy - NX) mayp = mask(i_xy + NX) maym = mask(i_xy - NX) ipxm = ilst(i_xy - 1) if (ix.eq.2) ipxm = ilst(i_xy - 1 + iper) ipxp = ilst(i_xy + 1) if (ix.eq.NX-1) ipxp = ilst(i_xy + 1 - iper) dep = deph(i_xy) cosi = bcosli(i_xy) fc = fcor(i_xy)/dep bex = bemx(i_xy) bey = bemy(i_xy) hx = 1./bex eps = (CNST_EPS/dep + CNST_EPT)/dep do jis = 0, n_islands if (ityp_island(jis).ne.2) then ind = in_rhs(jis) d_0j = jis.eq.0 rhsyp = d_0j if (mayp.eq.BC_W) rhsyp = rhsi(ipyp+ind) rhsym = d_0j if (maym.eq.BC_W) rhsym = rhsi(ipym+ind) coef(jis) = coef(jis) + cosi*hx*eps*bey*(rhsyp-rhsym)*0.5 coef(jis) = coef(jis) + * cosi*fc*(rhsi(ipxp+ind)-rhsi(ipxm+ind))*0.5 endif enddo enddo do jis = 0, n_islands a_island(0,jis) = coef(jis) enddo endif n_per_island = 1 if (use_per_island) n_per_island = 0 do kis = 1, n_islands if (ityp_island(kis).eq.2) goto 100 c integrate on rectangular path surrounding the kth (kis) island ind = in_rhs(kis) do jis = 0, n_islands coef(jis) = 0. enddo do ipoly = 1, n_poly(kis) i_min = i_poly(ipoly,kis) i_max = i_poly(ipoly+1,kis) j_min = j_poly(ipoly,kis) j_max = j_poly(ipoly+1,kis) istep = 1 if (i_max.lt.i_min) istep = -1 jstep = 1 if (j_max.lt.j_min) jstep = -1 i_maxm = i_max - istep j_maxm = j_max - jstep c write(60,*)'ipoly: ',kis,istep,jstep if (i_min.eq.i_max) goto 30 s = 0.5*istep iy = j_min do ix = i_min, i_max, istep i_xy = (iy - 1) * NX + ix if (mask(i_xy).ne.BC_W) then print*,'island integral (E/W) crosses a land or boundary * point at: ',ix,iy stop endif ipac = ilst(i_xy) ipym = ilst(i_xy - NX) ipyp = ilst(i_xy + NX) ipxm = ilst(i_xy - 1) ipxp = ilst(i_xy + 1) dep = deph(i_xy) bex = bemx(i_xy) bey = bemy(i_xy) hxs = s/bex eps = (CNST_EPS/dep + CNST_EPT) / dep maxp = mask(i_xy + 1) maxm = mask(i_xy - 1) mayp = mask(i_xy + NX) maym = mask(i_xy - NX) fc = fcor(i_xy)/dep cosi = bcosli(i_xy) do jis = n_per_island, n_islands if (ityp_island(jis).ne.2) then jnd = in_rhs(jis) bcij = BC_I(jis) rhsyp = 0. rhsym = 0. if (mayp.eq.BC_W) rhsyp = rhsi(ipyp+jnd) if (mayp.eq.bcij) rhsyp = 1. if (maym.eq.BC_W) rhsym = rhsi(ipym+jnd) if (maym.eq.bcij) rhsym = 1. rhsxp = 0. rhsxm = 0. if (maxp.eq.BC_W) rhsxp = rhsi(ipxp+jnd) if (maxp.eq.bcij) rhsxp = 1. if (maxm.eq.BC_W) rhsxm = rhsi(ipxm+jnd) if (maxm.eq.bcij) rhsxm = 1. fric = cosi*hxs*eps*bey*(rhsyp-rhsym)*0.5 cor = cosi*s*fc*(rhsxp-rhsxm)*0.5 coef(jis) = coef(jis) + fric coef(jis) = coef(jis) + cor endif c write(60,*)jis, ix,iy, fric, cor, coef(jis) enddo s = istep if (ix.eq.i_maxm) s = 0.5*istep enddo 30 continue if (j_min.eq.j_max) goto 40 s = 0.5*jstep ix = i_min do iy = j_min, j_max, jstep i_xy = (iy - 1) * NX + ix if (mask(i_xy).ne.BC_W) then print*,'island integral (N/S) crosses a land or boundary * point at: ',ix,iy stop endif ipac = ilst(i_xy) ipym = ilst(i_xy - NX) ipyp = ilst(i_xy + NX) ipxm = ilst(i_xy - 1) ipxp = ilst(i_xy + 1) dep = deph(i_xy) bex = bemx(i_xy) bey = bemy(i_xy) hys = s/bey eps = (CNST_EPS/dep + CNST_EPT) / dep maxp = mask(i_xy + 1) maxm = mask(i_xy - 1) mayp = mask(i_xy + NX) maym = mask(i_xy - NX) fc = fcor(i_xy)/dep do jis = n_per_island, n_islands if (ityp_island(jis).ne.2) then jnd = in_rhs(jis) bcij = BC_I(jis) rhsyp = 0. rhsym = 0. if (mayp.eq.BC_W) rhsyp = rhsi(ipyp+jnd) if (mayp.eq.bcij) rhsyp = 1. if (maym.eq.BC_W) rhsym = rhsi(ipym+jnd) if (maym.eq.bcij) rhsym = 1. rhsxp = 0. rhsxm = 0. if (maxp.eq.BC_W) rhsxp = rhsi(ipxp+jnd) if (maxp.eq.bcij) rhsxp = 1. if (maxm.eq.BC_W) rhsxm = rhsi(ipxm+jnd) if (maxm.eq.bcij) rhsxm = 1. fric = - hys*eps*bex*(rhsxp-rhsxm)*0.5 cor = s*fc*(rhsyp-rhsym)*0.5 coef(jis) = coef(jis) + fric coef(jis) = coef(jis) + cor endif c write(60,*)jis, ix,iy, fric, cor, coef(jis) enddo s = jstep if (iy.eq.j_maxm) s = 0.5*jstep enddo 40 continue enddo do jis = n_per_island, n_islands a_island(kis,jis) = coef(jis) enddo 100 continue enddo do kis = 1, n_islands if (ityp_island(kis).eq.3) then do jis = 1, n_islands if (ityp_island(jis).eq.1) then a_island(jis,kis) = 0. endif enddo endif enddo if ( ibar_key .ge. 2 ) then write(IUNIT_OUT, *) 'a_island=' do kis = 0,n_islands write(IUNIT_OUT, *) (a_island(kis,jis),jis=0,n_islands) enddo endif do jis = 0, 9 do kis = 0, 9 a_fac(jis,kis) = a_island(jis,kis) enddo enddo call sgeco(a_fac(0,0),10,n_islands+1,ipvt,rcond,z) if (rcond.lt.1e-12) then print*,'there is trouble in computing the influence of islands on' print*,' the barotropic transport' stop endif end subroutine do_island_integral c------------------------- include 'barotropic.h' character*1 mayp,maym,maxp,maxm common /island_mat/ ipvt(10),z(10) do kis = 0, n_islands b(kis) = 0. enddo c line integrals - trapezoid rule on uniform grid if (use_per_island) then iy = ib coef0 = 0. coef1 = 0. wind = 0. do ix = 2, NX-1 i_xy = (iy - 1) * NX + ix ipac = ilst(i_xy) ipyp = ilst(i_xy + NX) ipym = ilst(i_xy - NX) mayp = mask(i_xy + NX) maym = mask(i_xy - NX) ipxm = ilst(i_xy - 1) if (ix.eq.2) ipxm = ilst(i_xy - 1 + iper) ipxp = ilst(i_xy + 1) if (ix.eq.NX-1) ipxp = ilst(i_xy + 1 - iper) dep = deph(i_xy) cosi = bcosli(i_xy) fc = fcor(i_xy)/dep bex = bemx(i_xy) bey = bemy(i_xy) hx = 1./bex eps = (CNST_EPS/dep + CNST_EPT)/dep rhsyp = 0. if (mayp.eq.BC_W) rhsyp = rhs(ipyp) rhsym = 0. if (maym.eq.BC_W) rhsym = rhs(ipym) wind = wind + cosi* hx * taux(ipac) / dep coef0 = coef0 + cosi*hx*eps*bey*(rhsyp-rhsym)*0.5 coef1 = coef1 + cosi*fc*(rhs(ipxp)-rhs(ipxm))*0.5 enddo b(0) = - ( wind + coef0 + coef1 ) endif do kis = 1, n_islands c integrate on rectangular path surrounding each island c write(61,*)'island: ',kis if (ityp_island(kis).ne.2) then ind = in_rhs(kis) coef0 = 0. coef1 = 0. wind = 0. do ipoly = 1, n_poly(kis) i_min = i_poly(ipoly,kis) i_max = i_poly(ipoly+1,kis) j_min = j_poly(ipoly,kis) j_max = j_poly(ipoly+1,kis) istep = 0 if (i_max.lt.i_min) istep = -1 if (i_max.gt.i_min) istep = 1 jstep = 0 if (j_max.lt.j_min) jstep = -1 if (j_max.gt.j_min) jstep = 1 i_maxm = i_max - istep j_maxm = j_max - jstep c write(61,*)'ipoly: ',kis,istep,jstep if (istep.eq.0) goto 30 iy = j_min s = 0.5*istep do ix = i_min, i_max, istep i_xy = (iy - 1) * NX + ix ipac = ilst(i_xy) ipym = ilst(i_xy - NX) ipyp = ilst(i_xy + NX) ipxm = ilst(i_xy - 1) ipxp = ilst(i_xy + 1) dep = deph(i_xy) bex = bemx(i_xy) bey = bemy(i_xy) hxs = s/bex eps = (CNST_EPS/dep + CNST_EPT) / dep mayp = mask(i_xy + NX) maym = mask(i_xy - NX) maxp = mask(i_xy + 1) maxm = mask(i_xy - 1) fc = fcor(i_xy)/dep cosi = bcosli(i_xy) rhsxp = 0. rhsxm = 0. if (maxp.eq.BC_W) rhsxp = rhs(ipxp) if (maxm.eq.BC_W) rhsxm = rhs(ipxm) rhsyp = 0. rhsym = 0. if (mayp.eq.BC_W) rhsyp = rhs(ipyp) if (maym.eq.BC_W) rhsym = rhs(ipym) w = cosi*hxs*taux(ipac)/dep fric = cosi*hxs*eps*bey*(rhsyp-rhsym)*0.5 cor = cosi*s* fc* (rhsxp-rhsxm)*0.5 wind = wind + w coef0= coef0 + fric coef1= coef1 + cor s = istep if (ix.eq.i_maxm) s = 0.5*istep c write(61,*)ix,iy,w,fric,cor enddo c write(61,*)wind,coef0,coef1 30 continue if (jstep.eq.0) goto 40 ix = i_min s = 0.5*jstep do iy = j_min, j_max, jstep i_xy = (iy - 1) * NX + ix ipac = ilst(i_xy) ipym = ilst(i_xy - NX) ipyp = ilst(i_xy + NX) ipxm = ilst(i_xy - 1) ipxp = ilst(i_xy + 1) dep = deph(i_xy) bex = bemx(i_xy) bey = bemy(i_xy) hys = s/bey eps = (CNST_EPS/dep + CNST_EPT) / dep mayp = mask(i_xy + NX) maym = mask(i_xy - NX) maxp = mask(i_xy + 1) maxm = mask(i_xy - 1) fc = fcor(i_xy)/dep rhsxp = 0. rhsxm = 0. if (maxp.eq.BC_W) rhsxp = rhs(ipxp) if (maxm.eq.BC_W) rhsxm = rhs(ipxm) rhsyp = 0. rhsym = 0. if (mayp.eq.BC_W) rhsyp = rhs(ipyp) if (maym.eq.BC_W) rhsym = rhs(ipym) w = hys*tauy(ipac)/dep fric = - hys*eps*bex*(rhsxp-rhsxm)*0.5 cor = s*fc*(rhsyp-rhsym)*0.5 wind = wind + w coef0= coef0 + fric coef1= coef1 + cor s = jstep if (iy.eq.j_maxm) s = 0.5*jstep c write(61,*)ix,iy,w,fric,cor enddo c write(61,*)wind,coef0,coef1 40 continue enddo b(kis) = - (wind + coef0 + coef1) c write(61,*)'wind+coef0+coef1: ',b(kis) endif enddo c? do kis = 1, n_islands c? if (ityp_island(kis).eq.3) then c? do jis = 1, n_islands c? if (ityp_island(jis).eq.1) then c? b(jis) = b(jis) - a_island(jis,kis)*b_island(kis) c? endif c? enddo c? endif c? enddo c solve the 10x10 matrix equation for b_island: c a_island * b_island = b n.b. a is not a_island, since sgeco call sgesl(a_fac(0,0),10,n_islands+1,ipvt,b(0),0) c adjust the transports, as required if (itrans.eq.1) then atrans_t = btrans(0)*t_ramp*sin(2.*3.14159*t_day/ftrans) else atrans_t = btrans(0) endif if (ityp_island(0).eq.3) b_island(0) = atrans_t do kis = 1, n_islands if (ityp_island(kis).eq.3) then b_island(kis) = adelt* b_island(kis) + (1.-adelt)* btrans(kis) elseif (ityp_island(kis).eq.1) then b_island(kis) = adelt* b_island(kis) + (1.-adelt)* b(kis) endif enddo c write(61,*)'transports: ',(b(kis),kis = 0,n_islands) c write(61,*) c write(62,*) c write(63,*) c write(64,*) return end subroutine sgesl(a,lda,n,ipvt,b,job) integer lda,n,ipvt(1),job real a(lda,1),b(1) c c sgesl solves the real system c a * x = b or trans(a) * x = b c using the factors computed by sgeco or sgefa. c c on entry c c a real(lda, n) c the output from sgeco or sgefa. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c ipvt integer(n) c the pivot vector from sgeco or sgefa. c c b real(n) c the right hand side vector. c c job integer c = 0 to solve a*x = b , c = nonzero to solve trans(a)*x = b where c trans(a) is the transpose. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains a c zero on the diagonal. technically this indicates singularity c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if sgeco has set rcond .gt. 0.0 c or sgefa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call sgeco(a,lda,n,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call sgesl(a,lda,n,ipvt,c(1,j),0) c 10 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sdot c c internal variables c real sdot,t integer k,kb,l,nm1 c nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call saxpy(n-k,t,a(k+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call saxpy(k-1,t,a(1,k),1,b(1),1) 40 continue go to 100 50 continue c c job = nonzero, solve trans(a) * x = b c first solve trans(u)*y = b c do 60 k = 1, n t = sdot(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/a(k,k) 60 continue c c now solve trans(l)*x = y c if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb b(k) = b(k) + sdot(n-k,a(k+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end subroutine sgeco(a,lda,n,ipvt,rcond,z) integer lda,n,ipvt(1) real a(lda,1),z(1) real rcond c c sgeco factors a real matrix by gaussian elimination c and estimates the condition of the matrix. c c if rcond is not needed, sgefa is slightly faster. c to solve a*x = b , follow sgeco by sgesl. c to compute inverse(a)*c , follow sgeco by sgesl. c to compute determinant(a) , follow sgeco by sgedi. c to compute inverse(a) , follow sgeco by sgedi. c c on entry c c a real(lda, n) c the matrix to be factored. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a an upper triangular matrix and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c rcond real c an estimate of the reciprocal condition of a . c for the system a*x = b , relative perturbations c in a and b of size epsilon may cause c relative perturbations in x of size epsilon/rcond . c if rcond is so small that the logical expression c 1.0 + rcond .eq. 1.0 c is true, then a may be singular to working c precision. in particular, rcond is zero if c exact singularity is detected or the estimate c underflows. c c z real(n) c a work vector whose contents are usually unimportant. c if a is close to a singular matrix, then z is c an approximate null vector in the sense that c norm(a*z) = rcond*norm(a)*norm(z) . c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c linpack sgefa c blas saxpy,sdot,sscal,sasum c fortran abs,amax1,sign c c internal variables c real sdot,ek,t,wk,wkm real anorm,s,sasum,sm,ynorm integer info,j,k,kb,kp1,l c c c compute 1-norm of a c anorm = 0.0e0 do 10 j = 1, n anorm = amax1(anorm,sasum(n,a(1,j),1)) 10 continue c c factor c call sgefa(a,lda,n,ipvt,info) c c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) . c estimate = norm(z)/norm(y) where a*z = y and trans(a)*y = e . c trans(a) is the transpose of a . the components of e are c chosen to cause maximum local growth in the elements of w where c trans(u)*w = e . the vectors are frequently rescaled to avoid c overflow. c c solve trans(u)*w = e c ek = 1.0e0 do 20 j = 1, n z(j) = 0.0e0 20 continue do 100 k = 1, n if (z(k) .ne. 0.0e0) ek = sign(ek,-z(k)) if (abs(ek-z(k)) .le. abs(a(k,k))) go to 30 s = abs(a(k,k))/abs(ek-z(k)) call sscal(n,s,z,1) ek = s*ek 30 continue wk = ek - z(k) wkm = -ek - z(k) s = abs(wk) sm = abs(wkm) if (a(k,k) .eq. 0.0e0) go to 40 wk = wk/a(k,k) wkm = wkm/a(k,k) go to 50 40 continue wk = 1.0e0 wkm = 1.0e0 50 continue kp1 = k + 1 if (kp1 .gt. n) go to 90 do 60 j = kp1, n sm = sm + abs(z(j)+wkm*a(k,j)) z(j) = z(j) + wk*a(k,j) s = s + abs(z(j)) 60 continue if (s .ge. sm) go to 80 t = wkm - wk wk = wkm do 70 j = kp1, n z(j) = z(j) + t*a(k,j) 70 continue 80 continue 90 continue z(k) = wk 100 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c c solve trans(l)*y = w c do 120 kb = 1, n k = n + 1 - kb if (k .lt. n) z(k) = z(k) + sdot(n-k,a(k+1,k),1,z(k+1),1) if (abs(z(k)) .le. 1.0e0) go to 110 s = 1.0e0/abs(z(k)) call sscal(n,s,z,1) 110 continue l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t 120 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) c ynorm = 1.0e0 c c solve l*v = y c do 140 k = 1, n l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t if (k .lt. n) call saxpy(n-k,t,a(k+1,k),1,z(k+1),1) if (abs(z(k)) .le. 1.0e0) go to 130 s = 1.0e0/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 130 continue 140 continue s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c c solve u*z = v c do 160 kb = 1, n k = n + 1 - kb if (abs(z(k)) .le. abs(a(k,k))) go to 150 s = abs(a(k,k))/abs(z(k)) call sscal(n,s,z,1) ynorm = s*ynorm 150 continue if (a(k,k) .ne. 0.0e0) z(k) = z(k)/a(k,k) if (a(k,k) .eq. 0.0e0) z(k) = 1.0e0 t = -z(k) call saxpy(k-1,t,a(1,k),1,z(1),1) 160 continue c make znorm = 1.0 s = 1.0e0/sasum(n,z,1) call sscal(n,s,z,1) ynorm = s*ynorm c if (anorm .ne. 0.0e0) rcond = ynorm/anorm if (anorm .eq. 0.0e0) rcond = 0.0e0 return end integer function isamax(n,sx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c real sx(*),smax integer i,incx,ix,n c isamax = 0 if( n.lt.1 .or. incx.le.0 ) return isamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 smax = abs(sx(1)) ix = ix + incx do 10 i = 2,n if(abs(sx(ix)).le.smax) go to 5 isamax = i smax = abs(sx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 smax = abs(sx(1)) do 30 i = 2,n if(abs(sx(i)).le.smax) go to 30 isamax = i smax = abs(sx(i)) 30 continue return end real function sasum(n,sx,incx) c c takes the sum of the absolute values. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c real sx(*),stemp integer i,incx,m,mp1,n,nincx c sasum = 0.0e0 stemp = 0.0e0 if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx stemp = stemp + abs(sx(i)) 10 continue sasum = stemp return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,6) if( m .eq. 0 ) go to 40 do 30 i = 1,m stemp = stemp + abs(sx(i)) 30 continue if( n .lt. 6 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,6 stemp = stemp + abs(sx(i)) + abs(sx(i + 1)) + abs(sx(i + 2)) * + abs(sx(i + 3)) + abs(sx(i + 4)) + abs(sx(i + 5)) 50 continue 60 sasum = stemp return end subroutine saxpy(n,sa,sx,incx,sy,incy) c c constant times a vector plus a vector. c uses unrolled loop for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c real sx(*),sy(*),sa integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (sa .eq. 0.0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n sy(iy) = sy(iy) + sa*sx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m sy(i) = sy(i) + sa*sx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 sy(i) = sy(i) + sa*sx(i) sy(i + 1) = sy(i + 1) + sa*sx(i + 1) sy(i + 2) = sy(i + 2) + sa*sx(i + 2) sy(i + 3) = sy(i + 3) + sa*sx(i + 3) 50 continue return end real function sdot(n,sx,incx,sy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c real sx(*),sy(*),stemp integer i,incx,incy,ix,iy,m,mp1,n c stemp = 0.0e0 sdot = 0.0e0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n stemp = stemp + sx(ix)*sy(iy) ix = ix + incx iy = iy + incy 10 continue sdot = stemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m stemp = stemp + sx(i)*sy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) + * sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4) 50 continue 60 sdot = stemp return end subroutine sscal(n,sa,sx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to 1. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c real sa,sx(*) integer i,incx,m,mp1,n,nincx c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx sx(i) = sa*sx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m sx(i) = sa*sx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 sx(i) = sa*sx(i) sx(i + 1) = sa*sx(i + 1) sx(i + 2) = sa*sx(i + 2) sx(i + 3) = sa*sx(i + 3) sx(i + 4) = sa*sx(i + 4) 50 continue return end subroutine sgefa(a,lda,n,ipvt,info) integer lda,n,ipvt(1),info real a(lda,1) c c sgefa factors a real matrix by gaussian elimination. c c sgefa is usually called by sgeco, but it can be called c directly with a saving in time if rcond is not needed. c (time for sgeco) = (1 + 9/n)*(time for sgefa) . c c on entry c c a real(lda, n) c the matrix to be factored. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a an upper triangular matrix and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does c indicate that sgesl or sgedi will divide by zero c if called. use rcond in sgeco for a reliable c indication of singularity. c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sscal,isamax c c internal variables c real t integer isamax,j,k,kp1,l,nm1 c c c gaussian elimination with partial pivoting c info = 0 nm1 = n - 1 if (nm1 .lt. 1) go to 70 do 60 k = 1, nm1 kp1 = k + 1 c c find l = pivot index c l = isamax(n-k+1,a(k,k),1) + k - 1 ipvt(k) = l c c zero pivot implies this column already triangularized c if (a(l,k) .eq. 0.0e0) go to 40 c c interchange if necessary c if (l .eq. k) go to 10 t = a(l,k) a(l,k) = a(k,k) a(k,k) = t 10 continue c c compute multipliers c t = -1.0e0/a(k,k) call sscal(n-k,t,a(k+1,k),1) c c row elimination with column indexing c do 30 j = kp1, n t = a(l,j) if (l .eq. k) go to 20 a(l,j) = a(k,j) a(k,j) = t 20 continue call saxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 30 continue go to 50 40 continue info = k 50 continue 60 continue 70 continue ipvt(n) = n if (a(n,n) .eq. 0.0e0) info = n return end