program rdtracep ! ! EMISSIONS ON GRID FOR CARBON MONOXIDE ! get from UCI James Randoson (http://www.ess.uci.edu/~jranders/) ! Satellite fire count for biomass burning (1997-2001 and part of 2002) ! CARBON MONOXIDE ! original data array as: 180W-180E, 90N-90S ! converted data array as: 180W-180E, 90S-90N ! original value's unit: g CO /m2/month ! converted value's unit: kg CO /m2/month ! ! read in 1x1, convert to 2x2.5 integer, parameter :: lon1 = 360 integer, parameter :: lat1 = 180 integer, parameter :: lon2 = 144 integer, parameter :: lat2 = 91 integer, parameter :: mon = 12 real, parameter :: radius = 6.37122d+06 real, parameter :: pi = 3.141592653589793238462643383279 real value11(lon1,lat1,mon),value2(lon2,lat2,mon) real value1(lon1,lat1,mon),sf(lon1,lat1,mon),sf1(lon1,lat1,mon) real qtmp(lon1,lat1) character(len=13) filenm character(len=21) filename character(len=20) filesf character(len=80) dir character(len=3) nmon(mon) data nmon/'jan','feb','mar','apr','may','jun', & 'jul','aug','sep','oct','nov','dec'/ dir='/misc/mcc09/bian/uctm/co/mopitt/Ave-scalingfactor/' do i=1,lon1 do j=1,lat1 do m=1,mon value1(i,j,m)=0.d0 enddo enddo enddo dlambda=2.00*pi/360.00 dphi=pi/180.00 filesf='scalar_cmdl.dat' open(2,file=TRIM(dir)//TRIM(filesf),form='formatted') sf = 0.0 sf1 = 0.0 read(2,*) sf ! 100 format(360(e12.5,1x)) ! shift month ! original scaling factor is from 200004-200103 do nm=1,9 sf1(:,:,nm+3)=sf(:,:,nm) enddo do nm=10,12 sf1(:,:,nm-9)=sf(:,:,nm) enddo do ny = 1,3 filenm='CO_200001.txt' filename = 'co_bb_csf.2x25_200001' write(filenm(7:7),'(i1)') ny-1 totff =0.d0 do nm = 1,12 if(nm.lt.10) write(filenm(9:9),'(i1)') nm if(nm.ge.10) write(filenm(8:9),'(i2)') nm ! print *, 'filenm =',filenm, ny open(1,file=filenm) read(1,*) value1(:,:,nm) ! totff =0.d0 do i=1,360 do j=1,180 if( value1(i,j,nm) .eq. -999 ) value1(i,j,nm)=0.0 if( sf1(i,j,nm) .eq. -999 ) sf1(i,j,nm)=0.0 value11(i,j,nm)=value1(i,j,nm)*sf1(i,j,nm) rlat=89.5-j+1 da=radius*radius*cos(rlat*pi/180.00)*dlambda*dphi totff = totff+value11(i,j,nm)*da enddo enddo ! 100 print *, lon,lat,value1(lon,lat,nm) print *, 'nm =', nm, ' totff =', totff ! reflect in y-direct (lat), and change unit from g to kg do j=1,180 value11(:,j,nm)=value11(:,180-j+1,nm)/1000.d0 enddo !note: value must have a unit of over m2, not over box call map_a2a(lon1,lat1,1,value11(:,:,nm), & lon2,lat2,value2(:,:,nm),0,0,qtmp) ! call shift180(value2,lon2,lat2,1) write(filename(19:19),'(i1)') ny-1 if(nm.lt.10) write(filename(21:21),'(i1)') nm if(nm.ge.10) write(filename(20:21),'(i2)') nm open(2,file=filename) do i=1,lon2 do j=1,lat2 if(value2(i,j,nm)>0.d0) then ! write(2,'(2i4,1x,1p,e13.6)') i,j,value2(i,j) ! gctm write(2,'(2i4,1x,1p,e13.6)') j,i,value2(i,j,nm) ! uctm endif enddo enddo enddo ! nm enddo ! ny stop end !****************************************************************************** subroutine MAP_A2A( im, jm, km, q1, & in, jn, q2, ig, iv, qtmp) ! Horizontal arbitrary grid to arbitrary grid conservative high-order mapping ! S.-J. Lin implicit none integer im, jm, km, kt integer in, jn, jj integer ig ! ig=0: pole to pole; ig=1 j=1 is half-dy ! north of south pole integer iv ! iv=0: scalar; iv=1: vector real pi real q1(im,jm,km) real lon1(im+1) real lon2(in+1) real sin1(jm+1) real sin2(jn+1) real dx1, dx2, dy1, dy2 ! Output real q2(in,jn,km) ! local integer i,j,k,m real qtmp(im,jm) ! km=1 pi=4.*atan(1.) c dx1 = 360./float(im) dx2 = 360./float(in) dy1 = pi/float(jm-1) dy2 = pi/float(jn-1) c do i=1,im+1 lon1(i) = dx1 * (-0.5 + (i-1) ) enddo do i=1,in+1 lon2(i) = dx2 * (-0.5 + (i-1) ) enddo c sin1( 1) = -1. sin1(jm+1) = 1. sin2( 1) = -1. sin2(jn+1) = 1. do j=2,jm sin1(j) = sin( -0.5*pi + dy1*(-0.5+float(j-1)) ) enddo do j=2,jn sin2(j) = sin( -0.5*pi + dy2*(-0.5+float(j-1)) ) enddo do k=1, km kt = k if( im .eq. in ) then do j=1,jm-ig do i=1,im qtmp(i,j+ig) = q1(i,j+ig,kt) enddo enddo else ! Mapping in the E-W direction call xmap(im, jm-ig, lon1, q1(1,1+ig,kt), & in, lon2, qtmp(1,1+ig) ) endif if( jm .eq. jn ) then do j=1,jm-ig do i=1,in q2(i,j+ig,kt) = qtmp(i,j+ig) enddo enddo else ! Mapping in the N-S direction call ymap(in, jm, sin1, qtmp(1,1+ig), & jn, sin2, q2(1,1+ig,kt), ig, iv) endif enddo return end subroutine MAP_A2A !####################################################################### subroutine YMAP(im, jm, sin1, q1, jn, sin2, q2, ig, iv) ! Routine to perform aream preserving mapping in N-S from an arbitrary ! resolution to another. ! ! sin1 (1) = -1 must be south pole; sin1(jm+1)=1 must be N pole. ! ! sin1(1) < sin1(2) < sin1(3) < ... < sin1(jm) < sin1(jm+1) ! sin2(1) < sin2(2) < sin2(3) < ... < sin2(jn) < sin2(jn+1) ! ! Developer: S.-J. Lin ! First version: piece-wise constant mapping ! Apr 1, 2000 ! Last modified: implicit none ! Input integer im ! original E-W dimension integer jm ! original N-S dimension integer jn ! Target N-S dimension integer ig ! ig=0: scalars from S.P. to N. P. ! D-grid v-wind is also ig 0 ! ig=1: D-grid u-wind integer iv ! iv=0 scalar; iv=1: vector real sin1(jm+1-ig) ! original southern edge of the cell ! sin(lat1) ! real q1(im,jm-ig) ! original data at center of the cell real q1(im,jm) ! original data at center of the cell real sin2(jn+1-ig) ! Target cell's southern edge ! sin(lat2) ! Output ! real q2(im,jn-ig) ! Mapped data at the target resolution real q2(im,jn) ! Mapped data at the target resolution ! Local integer i, j0, m, mm integer j ! PPM related arrays real al(im,jm) real ar(im,jm) real a6(im,jm) real dy1(jm) real r3, r23 parameter ( r3 = 1./3., r23 = 2./3. ) real pl, pr, qsum, esl real dy, sum do j=1,jm-ig dy1(j) = sin1(j+1) - sin1(j) enddo ! *********************** ! Area preserving mapping ! *********************** ! Construct subgrid PP distribution call ppm_lat(im, jm, ig, q1, al, ar, a6, 3, iv) do 1000 i=1,im j0 = 1 do 555 j=1,jn-ig do 100 m=j0,jm-ig ! ! locate the southern edge: sin2(i) ! if(sin2(j) .ge. sin1(m) .and. sin2(j) .le. sin1(m+1)) then pl = (sin2(j)-sin1(m)) / dy1(m) if(sin2(j+1) .le. sin1(m+1)) then ! entire new cell is within the original cell pr = (sin2(j+1)-sin1(m)) / dy1(m) q2(i,j) = al(i,m) + 0.5*(a6(i,m)+ar(i,m)-al(i,m)) . *(pr+pl)-a6(i,m)*r3*(pr*(pr+pl)+pl**2) j0 = m goto 555 else ! South most fractional aream qsum = (sin1(m+1)-sin2(j))*(al(i,m)+0.5*(a6(i,m)+ . ar(i,m)-al(i,m))*(1.+pl)-a6(i,m)* . (r3*(1.+pl*(1.+pl)))) do mm=m+1,jm-ig ! locate the eastern edge: sin2(j+1) if(sin2(j+1) .gt. sin1(mm+1) ) then ! Whole layer qsum = qsum + dy1(mm)*q1(i,mm) else ! North most fractional aream dy = sin2(j+1)-sin1(mm) esl = dy / dy1(mm) qsum = qsum + dy*(al(i,mm)+0.5*esl* . (ar(i,mm)-al(i,mm)+a6(i,mm)*(1.-r23*esl))) j0 = mm goto 123 endif enddo goto 123 endif endif 100 continue 123 q2(i,j) = qsum / ( sin2(j+1) - sin2(j) ) 555 continue 1000 continue ! Final processing for poles if ( ig .eq. 0 .and. iv .eq. 0 ) then ! South pole sum = 0. do i=1,im sum = sum + q2(i,1) enddo sum = sum / float(im) do i=1,im q2(i,1) = sum enddo ! North pole: sum = 0. do i=1,im sum = sum + q2(i,jn) enddo sum = sum / float(im) do i=1,im q2(i,jn) = sum enddo endif end subroutine YMAP !######################################################################### subroutine PPM_LAT(im, jm, ig, q, al, ar, a6, jord, iv) implicit none !INPUT ! ig=0: scalar pole to pole ! ig=1: D-grid u-wind; not defined at poles because of staggering integer im, jm ! Dimensions integer ig real q(im,jm-ig) real al(im,jm-ig) real ar(im,jm-ig) real a6(im,jm-ig) integer jord integer iv ! iv=0 scalar ! iv=1 vector ! Local real dm(im,jm-ig) real r3 parameter ( r3 = 1./3. ) integer i, j, im2, iop, jm1 real tmp, qmax, qmin real qop ! Compute dm: linear slope do j=2,jm-1-ig do i=1,im dm(i,j) = 0.25*(q(i,j+1) - q(i,j-1)) qmax = max(q(i,j-1),q(i,j),q(i,j+1)) - q(i,j) qmin = q(i,j) - min(q(i,j-1),q(i,j),q(i,j+1)) dm(i,j) = sign(min(abs(dm(i,j)),qmin,qmax),dm(i,j)) enddo enddo im2 = im/2 jm1 = jm - 1 !Poles: if (iv .eq. 1 ) then ! !********* ! u-wind (ig=1) ! v-wind (ig=0) !********* ! ! SP do i=1,im if( i .le. im2) then qop = -q(i+im2,2-ig) else qop = -q(i-im2,2-ig) endif tmp = 0.25*(q(i,2) - qop) qmax = max(q(i,2),q(i,1), qop) - q(i,1) qmin = q(i,1) - min(q(i,2),q(i,1), qop) dm(i,1) = sign(min(abs(tmp),qmax,qmin),tmp) enddo ! NP do i=1,im if( i .le. im2) then qop = -q(i+im2,jm1) else qop = -q(i-im2,jm1) endif tmp = 0.25*(qop - q(i,jm1-ig)) qmax = max(qop,q(i,jm-ig), q(i,jm1-ig)) - q(i,jm-ig) qmin = q(i,jm-ig) - min(qop,q(i,jm-ig), q(i,jm1-ig)) dm(i,jm-ig) = sign(min(abs(tmp),qmax,qmin),tmp) enddo else ! !********* ! Scalar: !********* ! This code segment currently works only if ig=0 ! SP do i=1,im2 tmp = 0.25*(q(i,2)-q(i+im2,2)) qmax = max(q(i,2),q(i,1), q(i+im2,2)) - q(i,1) qmin = q(i,1) - min(q(i,2),q(i,1), q(i+im2,2)) dm(i,1) = sign(min(abs(tmp),qmax,qmin),tmp) enddo do i=im2+1,im dm(i, 1) = - dm(i-im2, 1) enddo ! NP do i=1,im2 tmp = 0.25*(q(i+im2,jm1)-q(i,jm1)) qmax = max(q(i+im2,jm1),q(i,jm), q(i,jm1)) - q(i,jm) qmin = q(i,jm) - min(q(i+im2,jm1),q(i,jm), q(i,jm1)) dm(i,jm) = sign(min(abs(tmp),qmax,qmin),tmp) enddo do i=im2+1,im dm(i,jm) = - dm(i-im2,jm) enddo endif do j=2,jm-ig do i=1,im al(i,j) = 0.5*(q(i,j-1)+q(i,j)) + r3*(dm(i,j-1) - dm(i,j)) enddo enddo do j=1,jm-1-ig do i=1,im ar(i,j) = al(i,j+1) enddo enddo if ( iv .eq. 1 ) then ! Vector: ! ig=0 if ( ig .eq. 0 ) then do i=1,im2 al(i, 1) = -al(i+im2,2) al(i+im2,1) = -al(i, 2) enddo do i=1,im2 ar(i, jm) = -ar(i+im2,jm1) ar(i+im2,jm) = -ar(i, jm1) enddo else ! ig=1 ! SP do i=1,im if( i .le. im2) then iop = i+im2 else iop = i-im2 endif al(i,1) = 0.5*(q(i,1)-q(iop,1)) - r3*(dm(iop,1) + dm(i,1)) enddo ! NP do i=1,im if( i .le. im2) then iop = i+im2 else iop = i-im2 endif ar(i,jm1) = 0.5*(q(i,jm1)-q(iop,jm1)) - & r3*(dm(iop,jm1) + dm(i,jm1)) enddo endif else ! Scalar (works for ig=0 only): do i=1,im2 al(i, 1) = al(i+im2,2) al(i+im2,1) = al(i, 2) enddo do i=1,im2 ar(i, jm) = ar(i+im2,jm1) ar(i+im2,jm) = ar(i, jm1) enddo endif do j=1,jm-ig do i=1,im a6(i,j) = 3.*(q(i,j)+q(i,j) - (al(i,j)+ar(i,j))) enddo call LMPPM(dm(1,j), a6(1,j), ar(1,j), & al(1,j), q(1,j), im, jord-3) enddo end subroutine PPM_LAT !####################################################################### subroutine XMAP(im, jm, lon1, q1, in, lon2, q2) ! Routine to perform area preserving mapping in E-W from an arbitrary ! resolution to another. ! Periodic domain will be assumed, i.e., the eastern wall bounding cell ! im is lon1(im+1) = lon1(1); Note the equal sign is true geographysically. ! ! lon1(1) < lon1(2) < lon1(3) < ... < lon1(im) < lon1(im+1) ! lon2(1) < lon2(2) < lon2(3) < ... < lon2(in) < lon2(in+1) ! ! Developer: S.-J. Lin ! First version: piece-wise constant mapping ! Apr 1, 2000 ! Last modified: implicit none ! Input integer im ! original E-W dimension integer in ! Target E-W dimension integer jm ! original N-S dimension real lon1(im+1) ! original western edge of the cell real q1(im,jm) ! original data at center of the cell real lon2(in+1) ! Target cell's western edge ! Output real q2(in,jm) ! Mapped data at the target resolution ! Local integer i1, i2 integer i, i0, m, mm integer j ! PPM related arrays real qtmp(-im:im+im) real al(-im:im+im) real ar(-im:im+im) real a6(-im:im+im) real x1(-im:im+im+1) real dx1(-im:im+im) real r3, r23 parameter ( r3 = 1./3., r23 = 2./3. ) real pl, pr, qsum, esl real dx integer iord data iord /3/ logical found do i=1,im+1 x1(i) = lon1(i) enddo do i=1,im dx1(i) = x1(i+1) - x1(i) enddo ! check to see if ghosting is necessary !************** ! Western edge: !************** found = .false. i1 = 1 do while ( .not. found ) if( lon2(1) .ge. x1(i1) ) then found = .true. else i1 = i1 - 1 if (i1 .lt. -im) then write(6,*) 'failed in xmap' stop else x1(i1) = x1(i1+1) - dx1(im+i1) dx1(i1) = dx1(im+i1) endif endif enddo !************** ! Eastern edge: !************** found = .false. i2 = im+1 do while ( .not. found ) if( lon2(in+1) .le. x1(i2) ) then found = .true. else i2 = i2 + 1 if (i2 .gt. 2*im) then write(6,*) 'failed in xmap' stop else dx1(i2-1) = dx1(i2-1-im) x1(i2) = x1(i2-1) + dx1(i2-1) endif endif enddo ! write(6,*) 'i1,i2=',i1,i2 do 1000 j=1,jm ! *********************** ! Area preserving mapping ! *********************** ! Construct subgrid PP distribution call PPM_CYCLE(im, q1(1,j), al(1), ar(1), a6(1), qtmp(0), iord) ! check to see if ghosting is necessary ! Western edge if ( i1 .le. 0 ) then do i=i1,0 qtmp(i) = qtmp(im+i) al(i) = al(im+i) ar(i) = ar(im+i) a6(i) = a6(im+i) enddo endif ! Eastern edge: if ( i2 .gt. im+1 ) then do i=im+1,i2-1 qtmp(i) = qtmp(i-im) al(i) = al(i-im) ar(i) = ar(i-im) a6(i) = a6(i-im) enddo endif i0 = i1 do 555 i=1,in do 100 m=i0,i2-1 ! ! locate the western edge: lon2(i) ! if(lon2(i) .ge. x1(m) .and. lon2(i) .le. x1(m+1)) then pl = (lon2(i)-x1(m)) / dx1(m) if(lon2(i+1) .le. x1(m+1)) then ! entire new grid is within the original grid pr = (lon2(i+1)-x1(m)) / dx1(m) q2(i,j) = al(m) + 0.5*(a6(m)+ar(m)-al(m)) . *(pr+pl)-a6(m)*r3*(pr*(pr+pl)+pl**2) i0 = m goto 555 else ! Left most fractional area qsum = (x1(m+1)-lon2(i))*(al(m)+0.5*(a6(m)+ . ar(m)-al(m))*(1.+pl)-a6(m)* . (r3*(1.+pl*(1.+pl)))) do mm=m+1,i2-1 ! locate the eastern edge: lon2(i+1) if(lon2(i+1) .gt. x1(mm+1) ) then ! Whole layer qsum = qsum + dx1(mm)*qtmp(mm) else ! Right most fractional area dx = lon2(i+1)-x1(mm) esl = dx / dx1(mm) qsum = qsum + dx*(al(mm)+0.5*esl* . (ar(mm)-al(mm)+a6(mm)*(1.-r23*esl))) i0 = mm goto 123 endif enddo goto 123 endif endif 100 continue 123 q2(i,j) = qsum / ( lon2(i+1) - lon2(i) ) 555 continue 1000 continue end subroutine XMAP !####################################################################### subroutine LMPPM(dm, a6, ar, al, p, im, lmt) implicit none real r12 parameter ( r12 = 1./12. ) integer im, lmt integer i real a6(im),ar(im),al(im),p(im),dm(im) real da1, da2, fmin, a6da c LMT = 0: full monotonicity c LMT = 1: semi-monotonic constraint (no undershoot) c LMT = 2: positive-definite constraint if(lmt.eq.0) then c Full constraint do 100 i=1,im if(dm(i) .eq. 0.) then ar(i) = p(i) al(i) = p(i) a6(i) = 0. else da1 = ar(i) - al(i) da2 = da1**2 a6da = a6(i)*da1 if(a6da .lt. -da2) then a6(i) = 3.*(al(i)-p(i)) ar(i) = al(i) - a6(i) elseif(a6da .gt. da2) then a6(i) = 3.*(ar(i)-p(i)) al(i) = ar(i) - a6(i) endif endif 100 continue elseif(lmt.eq.1) then c Semi-monotonic constraint do 150 i=1,im if(abs(ar(i)-al(i)) .ge. -a6(i)) go to 150 if(p(i).lt.ar(i) .and. p(i).lt.al(i)) then ar(i) = p(i) al(i) = p(i) a6(i) = 0. elseif(ar(i) .gt. al(i)) then a6(i) = 3.*(al(i)-p(i)) ar(i) = al(i) - a6(i) else a6(i) = 3.*(ar(i)-p(i)) al(i) = ar(i) - a6(i) endif 150 continue elseif(lmt.eq.2) then c Positive definite constraint do 250 i=1,im if(abs(ar(i)-al(i)) .ge. -a6(i)) go to 250 fmin = p(i) + 0.25*(ar(i)-al(i))**2/a6(i) + a6(i)*r12 if(fmin.ge.0.) go to 250 if(p(i).lt.ar(i) .and. p(i).lt.al(i)) then ar(i) = p(i) al(i) = p(i) a6(i) = 0. elseif(ar(i) .gt. al(i)) then a6(i) = 3.*(al(i)-p(i)) ar(i) = al(i) - a6(i) else a6(i) = 3.*(ar(i)-p(i)) al(i) = ar(i) - a6(i) endif 250 continue endif return end subroutine LMPPM !####################################################################### subroutine PPM_CYCLE(im, q, al, ar, a6, p, iord) implicit none real r3 parameter ( r3 = 1./3. ) ! Input integer im, iord real q(1) ! Output real al(1) real ar(1) real a6(1) real p(0:im+1) ! local real dm(0:im) integer i, lmt real tmp, qmax, qmin p(0) = q(im) do i=1,im p(i) = q(i) enddo p(im+1) = q(1) ! 2nd order slope do i=1,im tmp = 0.25*(p(i+1) - p(i-1)) qmax = max(p(i-1), p(i), p(i+1)) - p(i) qmin = p(i) - min(p(i-1), p(i), p(i+1)) dm(i) = sign(min(abs(tmp),qmax,qmin), tmp) enddo dm(0) = dm(im) do i=1,im al(i) = 0.5*(p(i-1)+p(i)) + (dm(i-1) - dm(i))*r3 enddo do i=1,im-1 ar(i) = al(i+1) enddo ar(im) = al(1) if(iord .le. 6) then do i=1,im a6(i) = 3.*(p(i)+p(i) - (al(i)+ar(i))) enddo lmt = iord - 3 if(lmt.le.2) call lmppm(dm(1),a6(1),ar(1),al(1),p(1),im,lmt) else call HUYNH(im, ar(1), al(1), p(1), a6(1), dm(1)) endif return end subroutine PPM_CYCLE !####################################################################### subroutine HUYNH(im, ar, al, p, d2, d1) ! Enforce Huynh's 2nd constraint in 1D periodic domain implicit none integer im, i real ar(im) real al(im) real p(im) real d2(im) real d1(im) ! Local scalars: real pmp real lac real pmin real pmax ! Compute d1 and d2 d1(1) = p(1) - p(im) do i=2,im d1(i) = p(i) - p(i-1) enddo do i=1,im-1 d2(i) = d1(i+1) - d1(i) enddo d2(im) = d1(1) - d1(im) ! Constraint for AR ! i = 1 pmp = p(1) + 2.0 * d1(1) lac = p(1) + 0.5 * (d1(1)+d2(im)) + d2(im) pmin = min(p(1), pmp, lac) pmax = max(p(1), pmp, lac) ar(1) = min(pmax, max(ar(1), pmin)) do i=2, im pmp = p(i) + 2.0*d1(i) lac = p(i) + 0.5*(d1(i)+d2(i-1)) + d2(i-1) pmin = min(p(i), pmp, lac) pmax = max(p(i), pmp, lac) ar(i) = min(pmax, max(ar(i), pmin)) enddo ! Constraint for AL do i=1, im-1 pmp = p(i) - 2.0*d1(i+1) lac = p(i) + 0.5*(d2(i+1)-d1(i+1)) + d2(i+1) pmin = min(p(i), pmp, lac) pmax = max(p(i), pmp, lac) al(i) = min(pmax, max(al(i), pmin)) enddo ! i=im i = im pmp = p(im) - 2.0*d1(1) lac = p(im) + 0.5*(d2(1)-d1(1)) + d2(1) pmin = min(p(im), pmp, lac) pmax = max(p(im), pmp, lac) al(im) = min(pmax, max(al(im), pmin)) ! compute A6 (d2) do i=1, im d2(i) = 3.*(p(i)+p(i) - (al(i)+ar(i))) enddo return end subroutine HUYNH !####################################################################### subroutine shift180(q,long,lat,levels) real q(long,lat,levels),tmp(long) longhalf=long/2 do k=1,levels do j=1,lat tmp(1:longhalf)=q(longhalf+1:long,j,k) q(longhalf+1:long,j,k)=q(1:longhalf,j,k) q(1:longhalf,j,k)=tmp(1:longhalf) end do end do return end subroutine shift180 !#######################################################################