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 baro_init, added ibar_key c 3/03/93 got rid of the deph_hx and deph_hy arrays c 3/08/93 changed call to baro_init, 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 baro_init 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 baro_init, 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 baro_init c 5/20/93 changed call to baro_solv, including Rayleigh friction c 6/04/93 changed call to baro_init, including nbaro c added subroutine baro_rinit 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 "baro_init" & "baro_rinit" c 7/12/95 added geometric term to time dependent operator, see "mod_rhs" c and "add_friction" 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 and emxy, 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_friction, c mod_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 baro_init, baro_rinit & 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-------------------------------------------------------------------------- subroutine baro_init (i_p, eps, nxp, nyp, nxyc, iox, nbx, lxx, * nby,lyx, alon, blon, alat, blat, x, y, db, glub, mgrid) 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), emxy(1), emx2(1), emy2(1), area(1) pointer (p_f, f), (p_emx, emx), (p_emy, emy), (p_emxy, emxy), * (p_emx2, emx2), (p_emy2, emy2), (p_area, area) common/data_geom/ p_f,p_emx,p_emy,p_emxy,p_emx2,p_emy2, p_area c--------------------------------------------------------------------------- #include "barotropic.h" integer iox(1), lxx(1), lyx(1) real db(1), x(1), y(1) integer itmp(1) pointer (pitmp, itmp) if_per = i_p GLUBINA = glub mod_cosys = mgrid X_MIN = alon X_MAX = blon Y_MIN = alat Y_MAX = blat hx_s = (X_MAX - X_MIN)/(nxp-1) hy_s = (Y_MAX - Y_MIN)/(nyp-1) NX = nxp + 2*if_per do i = 1 , nxp xx(i+if_per) = x(i) enddo NY = nyp do j = 1 , nyp yy(j) = y(j) enddo c add border of land/periodic points if (if_per.eq.1) then iper = NX - 2 xx(1) = X_MIN - hx_s xx(NX) = X_MAX + hx_s endif NXY = NX * NY call mem_alloc (pdeph, NXY, 2, 'baro deph') 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 (pbemxy, NXY, 2, 'baro emxy') call mem_alloc (pbemxx, NXY, 2, 'baro emxx') call mem_alloc (pbemyy, NXY, 2, 'baro emyy') call mem_alloc (pmask, NXY, 0, '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) fcor(i_xy) = f(i) bemx(i_xy) = emx(i) bemy(i_xy) = emy(i) bemxy(i_xy) = emxy(i) bemxx(i_xy) = emx2(i) bemyy(i_xy) = emy2(i) mask(i_xy) = BC_W enddo if (ibar_key .ne. 0) * open (unit = IUNIT_OUT, file = f_bar(1:n_bar)) if (ibar_key .eq. 3) 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. 3) 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 if (ibar_key .eq. 3) 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. 3) write(IUNIT_OUT,*) i, ix, iy, lyx(i) mask(i_xy) = BC_L enddo call mod_open (eps) call mod_init call mod_gra(0) if ( ibar_key .eq. 3 ) then call mem_alloc (pitmp, NXY, 1, '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)*NX),1000), ix=1,NX) enddo call mem_free (pitmp, NXY, 1) endif 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_DSINK) 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_DSINK) 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_DSINK) rely_p(i) = dep/deph(i_xy+NX) if (deph(i_xy-NX).ge.BAR_DSINK) 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') call mem_alloc(pbound_rhs, NXY, 2, 'baro bound rhs') #ifdef ISLAND call mem_alloc(prhs_bc0, NPACK, 2, 'baro rhs_bc0') #endif #ifdef ISLAND1 call mem_alloc(prhs_bc1, NPACK, 2, 'baro rhs_bc1') #endif call mod_mat call mod_lu(0) 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 baro_init 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, (compressed x) c---------------------- #include "barotropic.h" real tb(1), txb(1), tyb(1), psi(1) integer iox(1) call mem_alloc (prhs, NPACK, 2, 'baro rhs') call mem_alloc (ptaux, NPACK, 2, 'baro taux') call mem_alloc (ptauy, NPACK, 2, 'baro tauy') call do_winds(nxp,nxyc,iox,txb,tyb) call mod_rhs (nxp,nxyc,iox,tb) call mod_sol #ifdef ISLAND_SINK do i = 1, NPACK i_xy = list(i) bound_rhs(i_xy) = rhs(i) enddo #endif 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 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_1) dpsi = b_island(1) #ifdef ISLAND if (mask(i_xy).eq.BC_0) dpsi = b_island(0) #endif psi(i) = dpsi bound_rhs(i_xy) = dpsi enddo call mem_free(prhs, NPACK, 2) call mem_free(ptaux, NPACK, 2) call mem_free(ptauy, NPACK, 2) if (ibar_key .ne. 0) call flush(IUNIT_OUT) return end subroutine mod_open (eps) c------------------------------------------- #include "barotropic.h" CNST_EPS = rayl CNST_EPT = eps/nbaro 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 mod_init c------------------------- #include "barotropic.h" R_EARTH = 6378000 CNST_2OMEGA = 2. * 2. * MATH_PI / (24. * 3600.) ! 2 * 2*pi/day phi_0 = to_RAD(Y_MIN + Y_MAX) / 2. CNST_BETA = CNST_2OMEGA * cos(phi_0) CNST_F0 = CNST_2OMEGA * ( sin(phi_0) - phi_0 * cos(phi_0) ) c set various repeated-use vectors call extend if (ibar_key.eq.3) 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 endif #ifdef ISLAND1 call do_island #endif ISLAND1 #ifdef ISLAND_SINK call mod_sink if (ibar_key.eq.3) 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 endif #endif !ISLAND_SINK end subroutine mod_sink c------------------------- c fictitious domain method c------------------------------------- #include "barotropic.h" do i = 1, NXY if (mask(i).eq.BC_S) then deph(i) = BAR_DELTA * GLUBINA mask(i) = BC_W endif enddo end subroutine extend c------------------------- #include "barotropic.h" 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) deph(i_xy) = deph(i_xy+iper) enddo ix = NX do iy = 1, NY i_xy = (iy - 1) * NX + ix mask(i_xy) = mask(i_xy-iper) deph(i_xy) = deph(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 endif if (mask(i_xy).eq.BC_L) then mask(i_xy-NX+1) = BC_L endif enddo endif c make the Southern land mass into an island #ifdef ISLAND do i_xy = 1, NXY ix = mod(i_xy-1, NX) + 1 iy = (i_xy-ix)/NX + 1 dlat = yy(iy) if ((mask(i_xy).eq.BC_L) * .and. dlat.le.-60) mask(i_xy) = BC_0 enddo #endif c mark ICELAND #ifdef ISLAND1 do i_xy = 1, NXY ix = mod(i_xy-1, NX) + 1 iy = (i_xy-ix)/NX + 1 dlon = xx(ix) dlat = yy(iy) if (mask(i_xy).ne.BC_W .and. (dlat .lt. 67 .and. dlat. gt. 62) * .and. (dlon .lt. 347 .and. dlon .gt. 335)) mask(i_xy) = BC_1 enddo #endif #ifdef ISLAND_SINK do i_xy = 1, NXY ix = mod(i_xy-1, NX) + 1 iy = (i_xy-ix)/NX + 1 dlon = xx(ix) dlat = yy(iy) c mark Antilles c if (mask(i_xy).ne.BC_W .and. (dlat .lt. 24 .and. dlat. gt. 15) c * .and. (dlon .lt. 301 .and. dlon .gt. 276)) mask(i_xy) = BC_S c mark New Zealand c if (mask(i_xy).ne.BC_W .and. (dlat .lt. -30 .and. dlat. gt. -50) c * .and. (dlon .lt. 190 .and. dlon .gt. 160)) mask(i_xy) = BC_S c mark Kerguelen islands if (mask(i_xy).ne.BC_W .and. (dlat .lt. -40 .and. dlat. gt. -60) * .and. (dlon .lt. 80 .and. dlon .gt. 60)) mask(i_xy) = BC_S c mark Madagascar if (mask(i_xy).ne.BC_W .and. (dlat .lt. 0 .and. dlat. gt. -30) * .and. (dlon .lt. 55 .and. dlon .gt. 41)) mask(i_xy) = BC_S c mark Australia if (mask(i_xy).ne.BC_W .and. (dlat .lt. -10 .and. dlat. gt. -50) * .and. (dlon .lt. 160 .and. dlon .gt. 110)) mask(i_xy) = BC_S if (mask(i_xy).ne.BC_W .and. (dlat .lt. 5 .and. dlat. gt. -15) * .and. (dlon .lt. 160 .and. dlon .gt. 125)) mask(i_xy) = BC_S enddo #endif end subroutine do_island c------------------------- #include "barotropic.h" i_max = 1 i_min = NX j_max = 1 j_min = NY do i_xy = 1, NXY ix = mod(i_xy-1, NX) + 1 iy = (i_xy-ix)/NX + 1 if (mask(i_xy).eq.BC_1) then i_max = max(i_max,ix) i_min = min(i_min,ix) j_max = max(j_max,iy) j_min = min(j_min,iy) endif enddo i_max1 = i_max + 1 i_min1 = i_min - 1 j_max1 = j_max + 1 j_min1 = j_min - 1 end subroutine do_island_int c------------------------- #include "barotropic.h" #if defined ISLAND || defined ISLAND1 a_island(0,0) = 1. a_island(0,1) = 0. a_island(1,0) = 0. a_island(1,1) = 1. coef0 = 0. coef1 = 0. c line integrals - trapezoid rule on uniform grid #ifdef ISLAND c find the y index for 60 degrees south (a good place to do the line integral) iy = 1 10 iy = iy + 1 if (yy(iy).ge.-60) goto 20 goto 10 20 ib = iy iy = ib do ix = 2, NX-1 i_xy = (iy - 1) * NX + ix fc = fcor(i_xy) ipac = ilst(i_xy) ipyp = ilst(i_xy + NX) if (ix.gt.2) then ipxm = ilst(i_xy - 1) else ipxm = ilst(i_xy - 1 + iper) endif if (ix.lt.NX-1) then ipxp = ilst(i_xy + 1) else ipxp = ilst(i_xy + 1 - iper) endif dep = deph(i_xy) bex = bemx(i_xy) bey = bemy(i_xy) hx = 1./bex d2 = dep**2 eps = CNST_EPS + dep*CNST_EPT coef0 = coef0 + hx*eps*bey*(rhs0(ipyp)-rhs0(ipac))/d2 coef0 = coef0 + fc/dep*(rhs0(ipxp)-rhs0(ipxm))/2. #ifdef ISLAND1 coef1 = coef1 + hx*eps*bey*(rhs1(ipyp)-rhs1(ipac))/d2 coef1 = coef1 + fc/dep*(rhs1(ipxp)-rhs0(ipxm))/2. #endif enddo a_island(0,0) = coef0 a_island(0,1) = coef1 #endif #ifdef ISLAND1 c integrate on rectangular path surrounding island c assuming all points on path are ocean points coef0 = 0. coef1 = 0. i_maxm = i_max1 - 1 j_maxm = j_max1 - 1 iy = j_max1 s = 0.5 do ix = i_min1, i_max1 i_xy = (iy - 1) * NX + ix fc = fcor(i_xy) ipac = ilst(i_xy) 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 d2 = dep**2 eps = CNST_EPS + dep*CNST_EPT coef1 = coef1 + hxs*eps*bey*(rhs1(ipyp)-rhs1(ipac))/d2 coef1 = coef1 + s*fc/dep*(rhs1(ipxp)-rhs1(ipxm))/2. #ifdef ISLAND coef0 = coef0 + hxs*eps*bey*(rhs0(ipyp)-rhs0(ipac))/d2 coef0 = coef0 + s*fc/dep*(rhs0(ipxp)-rhs0(ipxm))/2. #endif s = 1. if (ix.eq.i_maxm) s = 0.5 enddo iy = j_min1 s = 0.5 do ix = i_min1, i_max1 i_xy = (iy - 1) * NX + ix fc = fcor(i_xy) ipac = ilst(i_xy) ipym = 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 d2 = dep**2 eps = CNST_EPS + dep*CNST_EPT coef1 = coef1 - hxs*eps*bey*(rhs1(ipac)-rhs1(ipym))/d2 coef1 = coef1 - s*fc/dep*(rhs1(ipxp)-rhs1(ipxm))/2. #ifdef ISLAND coef0 = coef0 - hxs*eps*bey*(rhs0(ipac)-rhs0(ipym))/d2 coef0 = coef0 - s*fc/dep*(rhs0(ipxp)-rhs0(ipxm))/2. #endif s = 1. if (ix.eq.i_maxm) s = 0.5 enddo ix = i_max1 s = 0.5 do iy = j_min1, j_max1 i_xy = (iy - 1) * NX + ix fc = fcor(i_xy) ipac = ilst(i_xy) ipxp = ilst(i_xy + 1) ipym = ilst(i_xy - NX) ipyp = ilst(i_xy + NX) dep = deph(i_xy) bex = bemx(i_xy) bey = bemy(i_xy) hys = s/bey d2 = dep**2 eps = CNST_EPS + dep*CNST_EPT coef1 = coef1 + hys*eps*bex*(rhs1(ipxp)-rhs1(ipac))/d2 coef1 = coef1 - s*fc/dep*(rhs1(ipyp)-rhs1(ipym))/2. #ifdef ISLAND coef0 = coef0 + hys*eps*bex*(rhs0(ipxp)-rhs0(ipac))/d2 coef0 = coef0 - s*fc/dep*(rhs0(ipyp)-rhs0(ipym))/2. #endif s = 1. if (iy.eq.j_maxm) s = 0.5 enddo ix = i_min1 s = 0.5 do iy = j_min1, j_max1 i_xy = (iy - 1) * NX + ix fc = fcor(i_xy) ipac = ilst(i_xy) ipxm = ilst(i_xy - 1) ipym = ilst(i_xy - NX) ipyp = ilst(i_xy + NX) dep = deph(i_xy) bex = bemx(i_xy) bey = bemy(i_xy) hys = s/bey eps = CNST_EPS + dep*CNST_EPT coef1 = coef1 - hys*eps*bex*(rhs1(ipac)-rhs1(ipxm))/d2 coef1 = coef1 + s*fc/dep*(rhs1(ipyp)-rhs1(ipym))/2. #ifdef ISLAND coef0 = coef0 - hys*eps*bex*(rhs0(ipac)-rhs0(ipxm))/d2 coef0 = coef0 + s*fc/dep*(rhs0(ipyp)-rhs0(ipym))/2. #endif s = 1. if (iy.eq.j_maxm) s = 0.5 enddo a_island(1,0) = coef0 a_island(1,1) = coef1 #endif if ( ibar_key .ge. 1 ) then write(IUNIT_OUT, *) 'a_island(2:2)=' write(IUNIT_OUT, *) a_island(0,0),a_island(0,1) write(IUNIT_OUT, *) a_island(1,0),a_island(1,1) endif #endif end subroutine do_island_int2 c------------------------- #include "barotropic.h" dimension b(0:2) #if defined ISLAND || defined ISLAND1 b(0) = 0. b(1) = 0. coef = 0. wind = 0. c line integrals - trapezoid rule on uniform grid #ifdef ISLAND iy = ib do ix = 2, NX-1 i_xy = (iy - 1) * NX + ix fc = fcor(i_xy) ipac = ilst(i_xy) ipyp = ilst(i_xy + NX) if (ix.gt.2) then ipxm = ilst(i_xy - 1) else ipxm = ilst(i_xy - 1 + iper) endif if (ix.lt.NX-1) then ipxp = ilst(i_xy + 1) else ipxp = ilst(i_xy + 1 - iper) endif dep = deph(i_xy) bex = bemx(i_xy) bey = bemy(i_xy) hx = 1./bex d2 = dep**2 wind = wind + hx * taux(ipac)/ dep eps = CNST_EPS + dep*CNST_EPT coef = coef + hx* eps* bey* (rhs(ipyp)-rhs(ipac))/d2 coef = coef + fc/dep*(rhs(ipxp)-rhs(ipxm))/2. enddo b(0) = - ( wind + coef ) #endif #ifdef ISLAND1 c integrate on rectangular path surrounding island coef = 0. wind = 0. i_maxm = i_max1 - 1 j_maxm = j_max1 - 1 iy = j_max1 s = 0.5 do ix = i_min1, i_max1 i_xy = (iy - 1) * NX + ix fc = fcor(i_xy) ipac = ilst(i_xy) ipxm = ilst(i_xy - 1) ipxp = ilst(i_xy + 1) ipyp = ilst(i_xy + NX) dep = deph(i_xy) bex = bemx(i_xy) bey = bemy(i_xy) hxs = s/bex d2 = dep**2 wind = wind + hxs * taux(ipac)/ dep eps = CNST_EPS + dep*CNST_EPT coef = coef + hxs* eps* bey*(rhs(ipyp)-rhs(ipac))/ d2 coef = coef + s* fc/ dep* (rhs(ipxp)-rhs(ipxm))/2. s = 1. if (ix.eq.i_maxm) s = 0.5 enddo iy = j_min1 s = 0.5 do ix = i_min1, i_max1 i_xy = (iy - 1) * NX + ix fc = fcor(i_xy) ipac = ilst(i_xy) ipxm = ilst(i_xy - 1) ipxp = ilst(i_xy + 1) ipym = ilst(i_xy - NX) dep = deph(i_xy) bex = bemx(i_xy) bey = bemy(i_xy) hxs = s/bex d2 = dep**2 wind = wind - hxs * taux(ipym)/ dep eps = CNST_EPS + dep*CNST_EPT coef = coef - hxs* eps* bey*(rhs(ipac)-rhs(ipym))/ d2 coef = coef - s * fc/dep*(rhs(ipxp)-rhs(ipxm))/2. s = 1. if (ix.eq.i_maxm) s = 0.5 enddo ix = i_max1 s = 0.5 do iy = j_min1, j_max1 i_xy = (iy - 1) * NX + ix fc = fcor(i_xy) ipac = ilst(i_xy) ipxp = ilst(i_xy + 1) ipym = ilst(i_xy - NX) ipyp = ilst(i_xy + NX) dep = deph(i_xy) bex = bemx(i_xy) bey = bemy(i_xy) hys = s/bey d2 = dep**2 wind = wind - hys * tauy(ipac)/ dep eps = CNST_EPS + dep*CNST_EPT coef = coef + hys* eps* bex*(rhs(ipxp)-rhs(ipac))/d2 coef = coef - s * fc/dep*(rhs(ipyp)-rhs(ipym))/2. s = 1. if (iy.eq.j_maxm) s = 0.5 enddo ix = i_min1 s = 0.5 do iy = j_min1, j_max1 i_xy = (iy - 1) * NX + ix fc = fcor(i_xy) ipac = ilst(i_xy) ipxm = ilst(i_xy - 1) ipym = ilst(i_xy - NX) ipyp = ilst(i_xy + NX) dep = deph(i_xy) bex = bemx(i_xy) bey = bemy(i_xy) hys = s/bey d2 = dep**2 wind = wind + hys * tauy(ipxm)/ dep eps = CNST_EPS + dep*CNST_EPT coef = coef - hys* eps* bex*(rhs(ipac)-rhs(ipxm))/d2 coef = coef + s * fc/dep*(rhs(ipyp)-rhs(ipym))/2. s = 1. if (iy.eq.j_maxm) s = 0.5 enddo b(1) = - ( wind + coef ) #endif det = a_island(0,0)*a_island(1,1) - a_island(0,1)*a_island(1,0) b_island(0) = (b(0)*a_island(1,1) - b(1)*a_island(0,1)) / det b_island(1) = (b(1)*a_island(0,0) - b(0)*a_island(1,0)) / det if ( ibar_key .ge. 1 ) then write(IUNIT_OUT, *) 'a_island(2:2) = ' write(IUNIT_OUT, *) a_island(0,0),a_island(0,1),b(0) write(IUNIT_OUT, *) a_island(1,0),a_island(1,1),b(1) write(IUNIT_OUT, *) 'setting CNST_T = ', b_island(0), b_island(1) endif #endif end subroutine mod_mat 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_NORM * = (GLUBINA*CNST_EPT + CNST_EPS) * (1./hx_ave**2 + 1./hy_ave**2) + * CNST_BETA * GLUBINA / hx_ave / R_EARTH do i_xy = 1, NXY sol(i_xy) = 0.0 enddo do i = 1, NPACK rhs_bc(i) = 0. #ifdef ISLAND rhs_bc0(i) = 0. #endif #ifdef ISLAND1 rhs_bc1(i) = 0. #endif 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_friction(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) heps = hinv*CNST_EPS gam = xi/ heps am = 1.0 - (gam-abs(gam))/2. ap = 1.0 + (gam+abs(gam))/2. ac = -(am+ap) c_h = CNST_EPS* hinv* hinv/ CNST_NORM ap = c_h * ap am = c_h * am ac = c_h * ac c ............................. y - direction hinv = bemy(i_xy) heps = hinv*CNST_EPS gam = eta/ heps bm = 1.0 - (gam-abs(gam))/2. bp = 1.0 + (gam+abs(gam))/2. bc = -(bm+bp) c_h = CNST_EPS* 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_friction (i_pac, i_xy) c----------------------------- #include "barotropic.h" common /matr_local/ ce, cw, cn, cs, csum real*8 ce, cw, cn, cs, csum ix = mod (i_xy -1 ,NX) + 1 iy = (i_xy - ix)/NX + 1 dep = deph(i_xy) const = 0.5* dep* CNST_EPT/ CNST_NORM bex = bemx(i_xy) bey = bemy(i_xy) hxe = const*bex**2 hye = const*bey**2 hxy = const*bemxy(i_xy)*bey hxx = const*bemxx(i_xy)*bex hyy = const*bemyy(i_xy)*bey cea = hxe * (1. + relx_p(i_pac)) + hxx cwa = hxe * (1. + relx_m(i_pac)) - hxx cna = hye * (1. + rely_p(i_pac)) + hxy + hyy csa = hye * (1. + rely_m(i_pac)) - hxy - hyy 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) bexy = bemxy(i_xy) bexx = bemxx(i_xy) beyy = bemyy(i_xy) dep = deph(i_xy) if (dep.ge.BAR_DSINK) then deph_x =bex*(relx_p(ip)-relx_m(ip)) * (relx_p(ip)+relx_m(ip)) / 2. deph_y =bey*(rely_p(ip)-rely_m(ip)) * (rely_p(ip)+rely_m(ip)) / 2. if (mask(i_xy-1).eq.BC_P) then f_cmx = fcor(i_xy - 1 + iper) else f_cmx = fcor(i_xy - 1) endif if (mask(i_xy+1).eq.BC_P) then f_cpx = fcor(i_xy + 1 - iper) else f_cpx = fcor(i_xy + 1) endif f_cpy = fcor(i_xy + NX) f_cmy = fcor(i_xy - NX) c calculate H**2 (f/H)x fh_x = dep * bex* (f_cpx*relx_p(ip) - f_cmx*relx_m(ip))/ 2. c calculate H**2 (f/H)y fh_y = dep * bey* (f_cpy*rely_p(ip) - f_cmy*rely_m(ip))/ 2. xi = CNST_EPS * (deph_x + bexx) + fh_y eta= CNST_EPS * (deph_y + beyy + bexy) - fh_x else ! sunken island xi = 0. eta= 0. endif 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 #ifdef ISLAND if (mask(i_xy) .eq. BC_0) then rhs_bc0(i_pac)=rhs_bc0(i_pac) - elem * solu return endif #endif #ifdef ISLAND1 if (mask(i_xy) .eq. BC_1) then rhs_bc1(i_pac)=rhs_bc1(i_pac) - elem * solu return endif #endif end subroutine mod_rhs (nxp,nxyc,iox,tb) c------------------------------------------- #include "barotropic.h" real tb(1) integer iox(1) real tmp1(1) pointer (ptmp1, tmp1) call mem_alloc (ptmp1, NXY, 2, 'baro tmp1') 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 dep = deph(i_xy) tmp1(i_xy) = dep * dep * tb(i) enddo do i_pac = 1, NPACK i_xy = list(i_pac) ix = mod(i_xy-1, NX) + 1 iy = (i_xy-ix)/NX + 1 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) bey = bemy(i_xy) dep = deph(i_xy) const = 0.5* dep* CNST_EPT hxe = const* bex**2 hye = const* bey**2 hxy = const* bemxy(i_xy)*bey hxx = const* bemxx(i_xy)*bex hyy = const* bemyy(i_xy)*bey ce = hxe * (1. + relx_p(i_pac)) + hxx cw = hxe * (1. + relx_m(i_pac)) - hxx cn = hye * (1. + rely_p(i_pac)) + hxy + hyy cs = hye * (1. + rely_m(i_pac)) - hxy - hyy csum = - (ce + cw + cn + cs) btemp= cw*bw + ce*be + cn*bn + cs*bs + csum*bc rhs(i_pac)= btemp enddo do i = 1, NXY bound_rhs(i) = 0. enddo do i_pac = 1, NPACK i = list(i_pac) bound_rhs(i) = (rhs(i_pac) + tmp1(i) ) / CNST_NORM enddo call mem_free (ptmp1, NXY, 2) end subroutine mod_lu(id) 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 real b1(1), sol1(1) pointer (pb1, b1), (psol1,sol1) data M12_HA /13/ call mem_alloc (pb1, NPACK, 2, 'baro b1') call mem_alloc (psol1, NPACK, 2, 'baro sol1') #ifdef ISLAND1 call mem_alloc (prhs1, NPACK, 2, 'baro rhs1') #endif if ( id .eq. 0) then 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') call mem_alloc (prhs0, NPACK, 2, 'baro rhs0') endif c modify rhs to include non-zero dirichlet boundary conditions #ifdef ISLAND do i = 1, NPACK rhs0(i) = rhs_bc0(i) enddo #endif NN12 = NONZ * RM12_NN call mem_realloc (piro, NGRAPH, NN12, NONZ, 1) if (id .eq. 0) then call mem_realloc (pico, NGRAPH, NN12, NONZ, 1) call mem_realloc (paa, NGRAPH, NN12, NONZ, 2) endif itime = icpu_time() if ( ibar_key.ge.1 ) 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 open (unit=12, file = 'matr', form = 'unformatted') write(12) NPACK, NONZ, NN12 write(12) (iflag(i),i=1,8) write(12) (aflag(i),i=1,8) write(12) (aa(i),i=1,NONZ) write(12) (iro(i),i=1,NONZ) 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, rhs0, b1, sol1, pivot, aflag, iflag, ifail) iflag(5) = 3 if (ibar_key.ge.1) 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.1) then write(IUNIT_OUT, *) write(IUNIT_OUT, *) 'CPU time for LU factorization was',itime,' sec' endif #ifdef ISLAND if (ibar_key.ge.1) 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 rhs0(i) = sol1(i) enddo #endif #ifdef ISLAND1 do i = 1, NPACK rhs1(i) = rhs_bc1(i) 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.1) then write(IUNIT_OUT, *) write(IUNIT_OUT, *)"ISLAND1" 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 rhs1(i) = sol1(i) enddo #endif 10 if (ifail .ne. 0) then print*, 'BARO: !!! Error code from Y12M:',ifail stop endif call do_island_int call mem_free (pb1, NPACK, 2) call mem_free (psol1, NPACK, 2) #ifdef ISLAND1 call mem_free (prhs1, NPACK, 2) #endif end subroutine mod_sol c-------------------------- #include "barotropic.h" real b1(1), sol1(1) pointer (pb1, b1), (psol1,sol1) c modify rhs to include non-zero dirichlet boundary conditions do i = 1, NPACK i_xy = list(i) rhs0(i) = bound_rhs(i_xy) + rhs_bc(i) enddo call mem_alloc (pb1, NPACK, 2, 'baro b1') call mem_alloc (psol1, NPACK, 2, 'baro sol1') if (ibar_key.ge.2) 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, rhs0, b1, sol1, pivot, aflag, iflag, ifail) if (ifail .ne. 0) goto 10 if (ibar_key.ge.2) 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) rhs0(i) = b1(i) enddo #if defined ISLAND || defined ISLAND1 c correct for the island influence, using a_island call do_island_int2 c# b_island(0) = 120.e+09 c# write(54,*) b_island(0) #ifdef ISLAND do i = 1, NPACK rhs0(i) = rhs0(i) + rhs_bc0(i)*b_island(0) enddo #endif #ifdef ISLAND1 do i = 1, NPACK rhs0(i) = rhs0(i) + rhs_bc1(i)*b_island(1) enddo #endif call y12mfe (NPACK, aa, ico, NN12, iro, NN12, a1, sn, NONZ, * ha, NPACK, rhs0, 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 call mem_free (pb1, NPACK, 2) call mem_free (psol1, NPACK, 2) 10 if (ifail .ne. 0) then print*, 'BARO: !!! Error code from Y12M:',ifail stop endif end subroutine mod_gra(id) c------------------------- #include "barotropic.h" if (id.eq.0) then call init_pack endif call to_y12m(id) 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_alloc (plist, NPACK, 1, 'baro list') call mem_alloc (pilst, NXY, 1, '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(id) c------------------------ #include "barotropic.h" NGRAPH = NPACK * 5 call mem_alloc (piro, NGRAPH, 1, 'baro iro') if ( id .eq. 0 ) then call mem_alloc (pico, NGRAPH, 1, 'baro ico') call mem_alloc (paa, NGRAPH, 2, 'baro aa') endif 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 baro_rinit (eps) c---------------------- c recompute the coefficients of the discrete barotropic equations and c then perform an approximate lu decomposition c eps = (input) effective friction induced by time step c---------------------- #include "barotropic.h" CNST_EPS = rayl CNST_EPT = eps/nbaro if (ibar_key .ge. 1) then write (IUNIT_OUT, *) write (IUNIT_OUT, *) 'rinit: epsilon =', CNST_EPT write (IUNIT_OUT, *) 'rinit: Rayleigh friction=', CNST_EPS endif call mod_gra(1) call mod_mat call mod_lu(1) call mem_free(piro, NN12, 1) end subroutine do_winds(nxp,nxyc,iox,txb,tyb) c---------------------- #include "barotropic.h" real txb(1), tyb(1) integer iox(1) real tmp1(1) pointer (ptmp1, tmp1) call mem_alloc (ptmp1, NXY, 2, 'baro tmp1') #if defined ISLAND || defined ISLAND1 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 tmp1(i_xy) = txb(i) enddo do i_pac = 1, NPACK i_xy = list(i_pac) bc = tmp1(i_xy) bn = tmp1(i_xy+NX) bnn = bn if (mask(i_xy+NX).eq.BC_W) bnn = tmp1(i_xy+NX+NX) bs = tmp1(i_xy-NX) taux(i_pac) = (7.*(bn + bc) - (bnn + bs) ) / 12. 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 tmp1(i_xy) = tyb(i) enddo do i_pac = 1, NPACK i_xy = list(i_pac) bc = tmp1(i_xy) be = tmp1(i_xy+1) bee = be if (mask(i_xy+1).eq.BC_W) bee = tmp1(i_xy+2) if (mask(i_xy-1).eq.BC_P) bw = tmp1(i_xy-1+iper) if (mask(i_xy+1).eq.BC_P) then be = tmp1(i_xy+1-iper) bee = tmp1(i_xy+2-iper) endif tauy(i_pac) = (7.*(be + bc) - (bee + bw) ) / 12. enddo h2y = 2 * hy do i_pac = 1, NPACK i_xy = list(i_pac) ix = mod(i_xy-1, NX) + 1 iy = (i_xy-ix)/NX + 1 bex = bemx(i_xy) bey = bemy(i_xy) be = bound_rhs(i_xy+1) if (mask(i_xy+1).eq.BC_P) be = bound_rhs(i_xy+1-iper) bn = bound_rhs(i_xy+NX) bc = bound_rhs(i_xy) taux(i_pac) = taux(i_pac) - bey* CNST_EPT*(bn - bc) tauy(i_pac) = tauy(i_pac) + bex* CNST_EPT*(be - bc) enddo #else do i_pac = 1, NPACK taux(i_pac) = 0. tauy(i_pac) = 0. enddo #endif call mem_free (ptmp1, NXY, 2) end