!* 21 Nov 2007: Added divergence damping (!8, !b8, !p8) and real ! del^4 diffusion (via horizontal_deln.f; !9). Not ! related to forcing but keeps perturbation code ! in step with full nonlinear code. ! !* Perturbation eqns forced by "residual" tendencies ! from QG basic state (assumed stationary, but really not). ! ! QG basic state read from file. ! ! Begins from chnp.pert.f (08 Feb 2006). ! ! Begins from chnp.new_3.f (11 Jan 2005). ! ! Computes both nonlinear solution and perturbation solution(s) ! with dynamics linearized about nonlinear solution. ! c c 3-d, boussinesq, elastic model in spherical coordinates c in this code, v is the south-north velocity, u is west-east c rigid lid upper boundary, c periodic in east-west. c advective form of equations, 4rth order differencing of c the horizontal advection terms, second order for the vertical c vertical advection terms c c this version has periodic b.c north-south. c c presently, the spherical terms are off. This version includes c the shift of bouyancy to the small timestep. c last modified 30 march 1990, wcs ! PARAMETER (NZ=121,NX=257,NX1=NX-1,NZ1=NZ-1, ! * NY=257,NY1=NY-1,NXY=NX*NY ) PARAMETER (NZ=61,NX=129,NX1=NX-1,NZ1=NZ-1, * NY=129,NY1=NY-1,NXY=NX*NY ) ! PARAMETER (NZ=31,NX=65,NX1=NX-1,NZ1=NZ-1, ! * NY=65,NY1=NY-1,NXY=NX*NY ) real f0,Nbv, theta0,gbyt0, xl,yl,ztop, tau integer kdepth parameter ( f0 = 1.e-4, ! Coriolis parameter (s^-1) & Nbv = 1.e-2, ! Buoyancy frequency (s^-1) & theta0 = 300., ! Reference theta (K) & gbyt0 = 9.81/theta0, ! g/theta_0 (m s^-1 K^-1) & xl = 3000.e3,! x periodicity (m) & yl = 3000.e3,! y periodicity (m) & ztop = 15.e3, ! height of lid (m) & tau = 1.e4, ! damping time at top of sponge layer (s) & kdepth = 10 ! no. levels in sponge layer & ) c c variables and some work space c DIMENSION U (0:NX,NY,NZ1),W (NX,NY,NZ),P (NX,NY,NZ1), * U1 (0:NX,NY,NZ1),W1 (NX,NY,NZ),P1 (NX,NY,NZ1), * P1_old (NX,NY,NZ1), * FU (0:NX,NY,NZ1),FW (NX,NY,NZ),fpp(nx,ny,nz1), * V (NX,0:NY,NZ1),T (NX,NY,NZ1),ft(nx,ny,nz1), * V1 (NX,0:NY,NZ1),T1 (NX,NY,NZ1),ftt(nx,ny,nz1), * t1t(nx,ny) , * fv (nx,0:ny,nz1), * ZU(NZ1),ZW(NZ),H(NZ),XM(NZ1),X(NX),A(NZ1),B(NZ1),C(NZ1), * ALPHA(NZ1),GAMMA(NZ1),WDZ(NX,NY) dimension tn(nx,2,nz1),tn1(nx,2,nz1),ts(nx,2,nz1),ts1(nx,2,nz1) dimension wduz(nx,ny,2),wdvz(nx,ny,2),wdtz(nx,ny,2),wdwz(nx,ny,2) dimension iph(20),ipv(20),an2(0:nz),tavg(nz1),tplt(ny*nz) real tbar(nz1) ! horizontal average of initial theta real fac ! stuff for sponge calculations ! arrays to save basic-state tendecies real * FU_bst(0:NX,NY,NZ1),FW_bst(NX,NY,NZ), & fv_bst(nx,0:ny,nz1), ft_bst(nx,ny,nz1), & fp_bst(nx,ny,nz1) integer iforce ! perturbation arrays and ancilliaries DIMENSION Up (0:NX,NY,NZ1),Wp (NX,NY,NZ),Pp (NX,NY,NZ1), * U1p(0:NX,NY,NZ1),W1p(NX,NY,NZ),P1p(NX,NY,NZ1), * Vp (NX,0:NY,NZ1),Tp (NX,NY,NZ1), * V1p(NX,0:NY,NZ1),T1p(NX,NY,NZ1) integer ievolve_basic_state, imake_pert_ics, ipert0, jpert0, kpert0, & iseed real cx_phase, pert_ampl, gaussdev logical debug integer idiffuse_basic_state real divd_coeff c c spherical parameters and coriolis parameters c dimension cospm(ny),cospv(0:ny),rcospm(ny),rcospv(0:ny), * sinpm(ny),sinpv(0:ny),rsinpm(ny), * fcorm(ny),fcorv(0:ny) c c data statements to set up the grid c alats, alatn are the latitude of the southern and northern c borders respectively and waveno is the longest periodic wave in the c slice we're simulating c data rearth,omega/6.370e+06,7.292e-05/ data alats,alatn,waveno/15.,85.,6./ data iph/1,19,9,18,16*0/ data ipv/31,19*0/ data nph,npv/2,1/ data irsetp,iavgg/192,24/ logical rstart,rdump c c function for fourth order difference c F4AS(UM1,UC,TM2,TM1,TP1,TP2) = 1 (UC+UM1)*((-TP2+TM2)+8.*(TP1-TM1)) !--Linearized counterpart to F4AS F4ASp(UM1,UC,TM2,TM1,TP1,TP2, UM1p,UCp,TM2p,TM1p,TP1p,TP2p) = & (UCp+UM1p)*((-TP2+TM2)+8.*(TP1-TM1)) & + (UC+UM1)*((-TP2p+TM2p)+8.*(TP1p-TM1p)) c c begin c rstart = .true. rdump = .true. !---Parameters controlling integration: dt = 1. ! time step ns0= 5 ! small time steps per large step t_end = 1000. ! length of integration (s) t_out = 10. ! interval for data output write(6,*) 'Enter dt, ns, t_end, t_out' read (5,*) dt, ns0, t_end, t_out !-- Viscosities; use +/- depending on 2nd/4th order. gamma1 = 0. ! horizontal filtering coeff for u,v gammaw = 0. ! ... for w gammat = 0. ! ... for t vsmoth = 0.007 ! vertical filtering coeff for w write(6,*) 'Enter gamma' read (5,*) rjunk ; gamma1= rjunk; gammaw= rjunk; gammat= rjunk; !-- Inputs for perturbation runs write(6,*) 'Evolve basic state? (Enter 1)' read (5,*) ievolve_basic_state ! if 1, basic state evolves; otherwise, ! basic state is (assumed) steady. write(6,*) 'Initialize perturbations? (Enter 1)' read (5,*) imake_pert_ics ! if 1, perturbations are initialized ! by this code; else read from pert.11 if ( imake_pert_ics .eq. 1 ) then !--- For white-noise ICs: no inputs needed; assume ampl. = 1 iseed = 1 ; pert_ampl = 1. !--- For delta-fn ICs: ! write(6,*) 'Enter i,j,k for initial delta-fn in theta:' ! read (5,*) ipert0, jpert0, kpert0 ! if ( ipert0 .gt. Nx .or. jpert .gt. Ny .or. kpert .gt. Nz ) then ! write(6,*) 'One of i,j,k outside array bounds: check Nx,Ny,Nz' ! stop ! endif endif if (ievolve_basic_state .eq. 0) then !-- Remove constant zonal velocity ... move in frame in which ! dipole is stationary (since we assume it is stationary!). write(6,*) 'Input dipole phase speed' read (5,*) cx_phase endif write(6,*) 'Include forcing? (Enter 1)' read (5,*) iforce ! if 1, linearized equations are forced by residual ! tendencies from QG dipole !9: set order of diffusion i_diff_order = 4 write(6,*) ' ... i_diff_order = ', i_diff_order !8: for divergence damping divd_coeff = 0.0 write(6,*) ' ... divd_coeff = ', divd_coeff !---Derived parameters: dx = xl/real(nx-1) ! x grid length dy = yl/real(ny-1) ! y grid length dz = ztop/real(nz-1) ! z grid length DTS= DT/FLOAT(NS0) ! small time step NS = NS0 nstopc = nint( t_end/dt ) ! number of large time steps rnu = (dx**4)*gamma1/(2.*dt) ! dimensional "viscosity" (not used) rkpa = (dx**4)*gammat/(2.*dt) gamma1y= gamma1 * dx**4 / dy**4 ! filtering coeff. for anisotropic grid gammawy= gammaw * dx**4 / dy**4 ! ... as above gammaty= gammat * dx**4 / dy**4 ! ... as above write(6,*) 'New params OK' c c define dx,dy and dz. dx = rearth*delta(longitude), c dy = rearth*delta(latitude) c twopi = 4.*asin(1.0) c dx = rearth*twopi*(360./waveno)/360./float(nx-1) c dy = rearth*twopi*(alatn-alats)/360./float(ny-1) c c set parameters for grid c deltap = (alatn-alats)/float(ny-1) do 10 j=1,ny c philat = twopi*(alats+float(j-1)*deltap)/360. c sinphl = sin(philat) c cosphl = cos(philat) c cospm(j) = cosphl c rcospm(j) = 1./cosphl c sinpm(j) = sinphl c rsinpm(j) = 1./sinphl c fcorm(j) = 2.*omega*sinphl c philat = 1. sinphl = 1. cosphl = 1. cospm(j) = 1. rcospm(j) = 1. sinpm(j) = 1. rsinpm(j) = 1. fcorm(j) = f0 10 continue do 11 j=0,ny c philat = twopi*(alats+float(j)*deltap-0.5*deltap)/360. c sinphl = sin(philat) c cosphl = cos(philat) c cospv(j) = cosphl c rcospv(j) = 1./cosphl c sinpv(j) = sinphl c fcorv(j) = 0.25*(2.*omega*sinphl) philat = 1. sinphl = 1. cosphl = 1. cospv(j) = 1. rcospv(j) = 1. sinpv(j) = 1. fcorv(j) = f0 11 continue B2= gbyt0 C2=90000. smdiv = .1 NZ2=NZ1-1 RDX=1./DX RDX2=RDX*RDX rdx12=rdx/12. rdx24=rdx/24. rdx48=rdx/48. RDY=1./DY RDY2=RDY*RDY rdy12=rdy/12. rdy24=rdy/24. rdy48=rdy/48. RDZ=1./DZ RDZ2=RDZ*RDZ EPSTS=.1 EPSSM=.2 EPSST=.2 DTL=DT kkk = 0 !--Read initial fields for basic state write(6,*) 'In main, Nbv=',Nbv time0= time call rd_data('fort.11', u,u1,v,v1,w,w1,p,p1,t,t1, time, & nx,ny,nz, xl,yl,ztop, f0,Nbv,gbyt0) if (ievolve_basic_state .eq. 0) then !-- Remove constant zonal velocity u = u - cx_phase ; u1 = u write(6,*) & 'CHANGING basic state u by subtracting ', cx_phase, 'm/s !!!' endif ! horizontal mean theta; set t = t1 = tbar tbar = sum( sum( t, 1), 1) / (nx*ny) !-- Make theta0 the mean surface temperature ! t = t + theta0 - tbar(1); t1 = t1 + theta0 - tbar(1); ! tbar = tbar + theta0 - tbar(1); open(19,file='uvwtp0.dat',form='unformatted') write(19) nx,ny,nz,xl,yl,ztop,dt,ns,gamma1,vsmoth,time write(19) u,v,w,t,p close(19) !-- Initialize perturbations if ( imake_pert_ics .eq. 1 ) then up = 0. ; u1p = 0. ; vp = 0. ; v1p = 0. ; wp = 0. ; w1p = 0. ; tp = 0. ; t1p = 0. ; pp = 0. ; p1p = 0. !--white-noise ICs ! do i=1,nx-1; do j=1,ny-1; do k=1,nz-1; do i=1,nx-1; do j=1,ny-1; do k=1,(nz-1)/4 ; tp(i,j,k) = gaussdev( iseed ) total = total + tp(i,j,k)**2 ! wp(i,j,k) = gaussdev( iseed ) ! total = total + wp(i,j,k)**2 enddo; enddo; enddo write(6,*) 'Initial rms tp =', total ! write(6,*) 'Initial rms wp =', total tp(nx,:,:) = tp(1,:,:) tp(:,ny,:) = tp(:,1,:) ! wp(nx,:,:) = wp(1,:,:) ! wp(:,ny,:) = wp(:,1,:) t1p = tp ! w1p = wp !--delta-fn ICs ! tp ( ipert0, jpert0, kpert0 ) = 1. ! t1p( ipert0, jpert0, kpert0 ) = 1. else !--read ICs from pert.11 call rd_data( & 'pert.11', up,u1p,vp,v1p,wp,w1p,pp,p1p,tp,t1p, time, & nx,ny,nz, xl,yl,ztop, f0,Nbv,gbyt0) endif if((time .gt. 0.1*dt).and. rstart) then ! This version ignores possibility of Euler step at beginning, ! ie., rstart is always .true. dtl = 2*dt ns = 2*ns0 end if write(6,'('' input data '', * i7,5(1x,e10.3))')nit,time, ! maxval(p),maxval(w),grwrte,maxval(v) c c set coefficients for implicit vertical small timesteps c first we'll just have a constant with time coefficient c theta_z= Nbv**2/gbyt0 write(6,*) 'dz, theta_z=', dz, theta_z do 21 k=1,nz1 ! tavg(k) = theta0 +dz*( real(k-1) +0.5 )*theta_z tavg(k) = theta0 -tbar(1) +tbar(k) !7 tavg(k) = tbar(k) write(6,*) tbar(k), tavg(k), & theta0 +dz*( real(k-1) )*theta_z 21 continue do 23 k=2,nz1 an2(k) = gbyt0*(tavg(k)-tavg(k-1))/dz 23 continue an2(1) = an2(2) an2(nz) = an2(nz1) an2(0) = an2(1) an2max = 0. do 24 k=2,nz1 do 24 ij=1,nxy an2max = amax1(an2max,(t(ij,1,k)-t(ij,1,k-1))) 24 continue an2max = an2max*gbyt0/dz c c compute the implicit coefficients for the small timestep c COFW=.25*DTS**2*C2/DZ**2*(1.+EPSSM)**2 COFT=.0625*dts**2*(1.+EPSST)**2 A(1)=0. B(1)=1. C(1)=0. DO 20 K=2,NZ1 A(K)=-COFW+COFT*an2(k-1) B(K)=1.+2.*(COFW+COFT*an2(k)) C(K)=-COFW+COFT*an2(k+1) 20 continue DO 30 K=1,NZ1 ALPHA(K)=1./(B(K)-A(K)*GAMMA(K-1)) GAMMA(K)=C(K)*ALPHA(K) 30 continue NIT=0 NPR=0 SWITCH=1. 3000 continue !-- Begin time loop : return here when nit < nstopc C C********PROCESSING FOR PLOTTING C deltat= time -time0 if (nit.eq.0) then write(21,*) 't', nx,ny, time0, t_out write(22,*) 'w', nx,ny, time0, t_out write(23,*) 'w', nx,ny, time0, t_out write(24,*) 'w', nx,ny, time0, t_out write(25,*) 'vort', nx,ny, time0, t_out endif IF( abs( deltat -nint(deltat/t_out)*t_out ) .le. 0.5*dt ) THEN write(21,*) t(:,:,1 ) write(22,*) w(:,:,5 ) write(23,*) w(:,:,21) write(24,*) w(:,:,2 ) write(25,*) -( u(1:nx-1,2:ny ,1) - u(1:nx-1,1:ny-1,1) ) / dy & +( v(2:nx ,1:ny-1,1) - v(1:nx-1,1:ny-1,1) ) / dx end if rns = 1./float(ns) KKK=KKK+1 NIT=NIT+1 NPR=NPR+1 TIME=TIME+DT if(mod(nit,10) .eq. 0) * write(6,'('' step, time '',i7,1x,e10.3)')nit,time !--------- !-- Basic-state calculations begin here if (ievolve_basic_state .eq. 1 .or. & ( iforce .eq. 1 .and. nit .eq. 1 ) ) ! calc b.st tends if nit = 1 & then ! XXXX Include tendencies of p and w ??? Probably need to test, but seems ! like they would disappear in anelastic, hydrostatic model (though ! there would be surf. p tendency??) ! XXX should check if it matters that diffn is included in b.st. tend XXX ! SEE: idiffuse_basic_state C C********LARGE TIME STEP CALCULATIONS C ku = 1 kd = 2 do 55 ij=1,nxy wduz(ij,1,1) = 0. wdtz(ij,1,1) = 0. wdwz(ij,1,1) = 0. wdvz(ij,1,1) = 0. 55 continue c c loop over z c ! do 60 k=1,nz1 do k=1,nz1 ktmp = kd kd = ku ku = ktmp kp1=min0(nz1,k+1) km1=max0(1,k-1) do 61 ij=1,nxy t1t(ij,1) = t1(ij,1,k) 61 continue c c first, vertical advection c do 70 j=1,ny-1 do 70 i=1,nx-1 wduz(i,j,ku)=.5*(W(i,j,k+1)+W(i+1,j,k+1)) * *RDZ*(U(i,j,kp1)-U(i,j,k)) wdtz(i,j,ku)= 0. * + W(i,j,k+1) * *RDZ*(t(i,j,kp1)-t(i,j,k) * -tavg(kp1)+tavg(k) * ) wdwz(i,j,ku)=.5*(W(i,j,k+1)+W(i,j,k)) * *RDZ*(w(i,j,k+1)-w(i,j,k)) 70 continue do 71 j=1,ny-1 do 71 i=1,nx-1 wdvz(i,j,ku)=.5*(W(i,j,k+1)+W(i,j+1,k+1)) * *RDZ*(v(i,j,kp1)-v(i,j,k)) 71 continue do 72 j=1,ny-1 do 72 i=1,nx-1 fu(i,j,k) = -.5*dts*(wduz(i,j,ku)+wduz(i,j,kd)) fv(i,j,k) = -.5*dts*(wdvz(i,j,ku)+wdvz(i,j,kd)) fw(i,j,k) = -.5*dts*(wdwz(i,j,ku)+wdwz(i,j,kd)) ft(i,j,k) = -.5*dtl*(wdtz(i,j,ku)+wdtz(i,j,kd)) 72 continue c c next, horizontal advection c c first, x-components c do 80 j=1,ny-1 jp1 = min0(j+1,ny) do 80 i=3,nx-2 c fu(i,j,k) = fu(i,j,k) * -rdx12*dts*f4as(0.,u(I,J,K),u(I-2,J,K),u(I-1,J,K), * u(I+1,J,K),u(I+2,J,K)) c fv(i,j,k) = fv(i,j,k) * - rdx48*dts*f4as( U(I-1,jp1,K)+U(I-1,J,K), * U(I,jp1,K)+U(I,J,K), V(I-2,J,K), V(I-1,J,K), * V(I+1,J,K), V(I+2,J,K) ) fw(i,j,k) = fw(i,j,k) * -RDX48*dts*F4AS( U(I-1,J,K)+ u(i-1,j,km1), * U(I,J,K)+ u(i,j,km1), * w(I-2,J,K),w(I-1,J,K), w(I+1,J,K), w(I+2,J,K) ) ft(i,j,k) = ft(i,j,k) * -RDX24*dtl*F4AS( U(I-1,J,K), * U(I,J,K), * t(I-2,J,K),t(I-1,J,K), t(I+1,J,K), t(I+2,J,K) ) 80 continue c c use periodic b.c. where needed c do 81 ii=1,3 i = ii if(ii.eq.3) i = nx-1 im1 = i-1 im2 = i-2 ip1 = i+1 ip2 = i+2 if(im1.lt.1) im1=nx-2+i if(im2.lt.2) im2=nx-3+i if(ip2.gt.nx) ip2 = 2 do 82 j=1,ny-1 c fu(i,j,k) = fu(i,j,k) * -rdx12*dts*f4as(0.,u(I,J,K),u(IM2,J,K),u(IM1,J,K), * u(IP1,J,K),u(IP2,J,K)) c fv(i,j,k) = fv(i,j,k) * - rdx48*dts*f4as( U(IM1,j+1,K)+U(IM1,J,K), * U(I,j+1,K)+U(I,J,K), V(IM2,J,K), V(IM1,J,K), * V(IP1,J,K), V(IP2,J,K) ) c fw(i,j,k) = fw(i,j,k) * -RDX48*dts*F4AS( U(IM1,J,K)+ u(im1,j,km1), * U(I,J,K)+ u(i,j,km1), * w(IM2,J,K),w(IM1,J,K), w(IP1,J,K), w(IP2,J,K) ) c ft(i,j,k) = ft(i,j,k) * -RDX24*dtl*F4AS( U(IM1,J,K), * U(I,J,K), * t(IM2,J,K),t(IM1,J,K), t(IP1,J,K), t(IP2,J,K) ) 82 continue 81 continue c c finished with x terms c next, y terms, first, 4rth order advection and filter in interior c do 90 j=3,ny-2 do 90 i=1,nx-1 fu(i,j,k) = fu(i,j,k) * -rdy48*dts*f4as( v(I,j-1,K)+v(I+1,J-1,K), * v(I,j,K)+v(I+1,J,K), u(I,J-2,K), u(I,J-1,K), * u(I,J+1,K), u(I,J+2,K) ) fv(i,j,k) = fv(i,j,k) * -rdy12*dts*f4as(0.,v(I,J,K),v(I,j-2,K),v(I,j-1,K), * v(I,j+1,K),v(I,j+2,K)) fw(i,j,k) = fw(i,j,k) * -RDy48*dts*F4AS( v(I,J-1,K)+ v(i,j-1,km1), * v(I,J,K)+ v(i,j,km1), * w(I,J-2,K),w(I,J-1,K), w(I,J+1,K), w(I,J+2,K) ) ft(i,j,k) = ft(i,j,k) * -RDy24*dtl*F4AS( v(I,j-1,K), * v(I,J,K), * t(I,j-2,K),t(I,j-1,K), t(I,j+1,K), t(I,j+2,K) ) 90 continue c c apply periodic boundary conditions in y c do 91 jj=1,3 j=jj if(jj .eq. 3) j=ny-1 jm1 = j-1 jm2 = j-2 jp1 = j+1 jp2 = j+2 if (jm1.lt.1) jm1 = ny-2+j if (jm2.lt.1) jm2 = ny-3+j if (jp2.gt.ny) jp2 = 2 do 91 i=1,nx-1 fu(i,j,k) = fu(i,j,k) * -rdy48*dts*f4as( v(I,jm1,K)+v(I+1,Jm1,K), * v(I,j,K)+v(I+1,J,K), u(I,Jm2,K), u(I,Jm1,K), * u(I,Jp1,K), u(I,Jp2,K) ) fv(i,j,k) = fv(i,j,k) * -rdy12*dts*f4as(0.,v(I,J,K),v(I,jm2,K),v(I,jm1,K), * v(I,jp1,K),v(I,jp2,K)) fw(i,j,k) = fw(i,j,k) * -RDy48*dts*F4AS( v(I,Jm1,K)+ v(i,jm1,km1), * v(I,J,K)+ v(i,j,km1), * w(I,Jm2,K),w(I,Jm1,K), w(I,Jp1,K), w(I,Jp2,K) ) ft(i,j,k) = ft(i,j,k) * -RDy24*dtl*F4AS( v(I,jm1,K), * v(I,J,K), * t(I,jm2,K),t(I,jm1,K), t(I,jp1,K), t(I,jp2,K) ) 91 continue !9-- Horizontal diffusion: idiffuse_basic_state = 0 if ( idiffuse_basic_state .eq. 0 ) then call horizontal_deln( i_diff_order, nx,ny,nz, & gamma1, gammaw, gammat, rns, & fu(:,:,k), fv(:,:,k), fw(:,:,k), ft(:,:,k), & u1(:,:,k), v1(:,:,k), w1(:,:,k), t1t ) endif c c the coriolis terms and other terms from the spherical grid. c do 110 j=1,ny-1 do 110 i=1,nx-1 ! fu(i,j,k) = fu(i,j,k) + dts*(fcorm(j)* fu(i,j,k) = fu(i,j,k) + dts*( f0 * * 0.25*(v(i,j,k)+v(i,j-1,k)+v(i+1,j,k)+v(i+1,j-1,k)) ) ! fv(i,j,k) = fv(i,j,k) - dts*fcorv(j)*( fv(i,j,k) = fv(i,j,k) - dts* f0 *( * 0.25*(u(i,j,k)+u(i,j+1,k)+u(i-1,j,k)+u(i-1,j+1,k) ) ) 110 continue ! XXX OMIT SPONGE IN CALCULATING B.ST. TENDENDCIES !!--- Sponge layer (19 Aug 03) ! if ( k .gt. nz1 - kdepth ) then ! fac = dts * ( k - (nz1 - kdepth) ) / kdepth / tau ! ! write(6,*) k, fac, (dtl/dts) * fac ! fu(:,:,k) = fu(:,:,k) - fac * u1(:,:,k) ! fv(:,:,k) = fv(:,:,k) - fac * v1(:,:,k) ! fw(:,:,k) = fw(:,:,k) - fac * w1(:,:,k) ! fac = (dtl/dts) * fac ! ft(:,:,k) = ft(:,:,k) - fac * ( t1(:,:,k) - tbar(k) ) ! !7 ft(:,:,k) = ft(:,:,k) - fac * ( t1(:,:,k) - tavg(k) ) ! endif c c set periodic boundary conditions c do 120 j=1,ny fu(nx,j,k) = fu(1,j,k) fu(0,j,k) = fu(nx-1,j,k) fv(nx,j,k) = fv(1,j,k) fw(nx,j,k) = fw(1,j,k) ft(nx,j,k) = ft(1,j,k) 120 continue do 121 i=1,nx fu(i,ny,k) = fu(i,1,k) fw(i,ny,k) = fw(i,1,k) ft(i,ny,k) = ft(i,1,k) fv(i,ny,k) = fv(i,1,k) fv(i,0,k) = fv(i,ny-1,k) 121 continue ! fu(0,ny,k) = fu(0,1,k) c !60 continue enddo !--- END of loop over z c c done computing largetimestep stuff, reset temp and set up small timestep c ! XXX OMIT TIME FILTER WHEN CALCULATING B.ST. TENDENCIES; ! XXX IF ievolve_basic_state = 1, WILL HAVE TO PUT THIS BACK IN! !c first, do partial time filter for middle level variables u,v,w,p !c ! do 130 ijk=1,nx*ny*nz1 ! w(ijk,1,1) = w(ijk,1,1) + epsts*(w1(ijk,1,1)-2*w(ijk,1,1)) ! p(ijk,1,1) = p(ijk,1,1) + epsts*(p1(ijk,1,1)-2*p(ijk,1,1)) ! t(ijk,1,1) = t(ijk,1,1) + epsts*(t1(ijk,1,1)-2*t(ijk,1,1)) !130 continue ! do 131 ijk=0,(nx+1)*ny*nz1 ! u(ijk,1,1) = u(ijk,1,1) + epsts*(u1(ijk,1,1)-2*u(ijk,1,1)) !131 continue ! do 132 ijk=1,nx*(ny+1)*nz1 ! v(ijk,0,1) = v(ijk,0,1) + epsts*(v1(ijk,0,1)-2*v(ijk,0,1)) !132 continue c c set zero b.c. for w c do 133 ij=1,nx*ny fw(ij,1,1) = 0. fw(ij,1,nz) = 0. w(ij,1,1) = 0. w(ij,1,nz) = 0. w1(ij,1,1) = 0. w1(ij,1,nz) = 0. 133 continue vsm4 = -vsmoth*rns vsm2 = -4.*vsm4 do 134 k=3,nz-2 do 134 ij=1,nxy fw(ij,1,k) = fw(ij,1,k) + vsm4*( 1 w1(ij,1,k+2)+w1(ij,1,k-2)+6*w1(ij,1,k) 2 -4*(w1(ij,1,k-1)+w1(ij,1,k+1)) ) 134 continue k=2 do 135 ij=1,nxy fw(ij,1,k) = fw(ij,1,k) + vsm2*( 1 w1(ij,1,k+1)-2*w1(ij,1,k)+w1(ij,1,k-1) ) 135 continue k=nz-1 do 136 ij=1,nxy fw(ij,1,k) = fw(ij,1,k) + vsm2*( 1 w1(ij,1,k+1)-2*w1(ij,1,k)+w1(ij,1,k-1) ) 136 continue C C********SMALL TIME STEP CALCULATIONS C p1_old = p1 ; !b8: for divergence damping DO 1099 M=1,NS c c advance u c do 1110 k=1,nz1 !4 do 1111 j=1,ny do 1111 j=1,ny-1 do 1111 i=1,nx-1 !b8 U1(i,j,k)=U1(i,j,k)-DTS*RDX*(P1(I+1,j,k)-P1(I,j,k)) U1(i,j,k)=U1(i,j,k)-DTS*RDX*(p1_old(I+1,j,k)-p1_old(I,j,k)) * +FU(I,j,k) 1111 continue do 1112 j=1,ny u1(nx,j,k) = u1(1,j,k) u1(0,j,k) = u1(nx-1,j,k) 1112 continue do 1113 i=0,nx !4 u1(i,ny,k) = 0.5*(u1(i,ny,k)+u1(i,1,k)) !4 u1(i,1,k) = u1(i,ny,k) u1(i,ny,k) = u1(i,1,k) 1113 continue 1110 continue c c advance v c do 1115 k=1,nz1 do 1115 j=1,ny-1 !4 do 1115 i=1,nx do 1115 i=1,nx-1 !b8 v1(i,j,k)=v1(i,j,k)-DTS*RDy*(P1(I,j+1,k)-P1(I,j,k)) v1(i,j,k)=v1(i,j,k)-DTS*RDy*(p1_old(I,j+1,k)-p1_old(I,j,k)) * +Fv(I,j,k) 1115 continue do 1116 k=1,nz1 do 1116 j=0,ny !4 v1(nx,j,k) = 0.5*(v1(1,j,k)+v1(nx,j,k)) !4 v1(1,j,k) = v1(nx,j,k) v1(nx,j,k) = v1(1,j,k) 1116 continue do 1117 k=1,nz1 do 1117 i=1,nx v1(i,0,k) = v1(i,ny-1,k) v1(i,ny,k) = v1(i,1,k) 1117 continue c c advance w and p implicitly c DO 1130 K=1,NZ1 do 1130 ij=1,nxy FPP(ij,1,k)=P1(ij,1,k)-.5*DTS*C2*RDZ*(1.-EPSSM) * *(W1(ij,1,k+1)-W1(ij,1,k)) 1130 continue do 1131 k=1,nz1 km1 = max0(1,k-1) kp1 = min0(nz1,k+1) do 1131 ij=1,nxy c t1(ij,1,k) = t1(ij,1,k)+ft(ij,1,k)*rns ftt(ij,1,k) = t1(ij,1,k)+ft(ij,1,k)*rns * -.25*dts*(1.-epsst)*rdz* * ( w1(ij,1,k+1)*(tavg(kp1)-tavg(k)) * +w1(ij,1,k )*(tavg(k)-tavg(km1)) ) 1131 continue c DO 1140 K=2,NZ1 do 1140 ij=1,nxy W1(ij,1,k)=W1(ij,1,k)+FW(ij,1,k) * -.5*DTS*RDZ*(1.-EPSSM)*(P1(ij,1,k)-P1(ij,1,k-1)) * +.25*dts*b2*(1.-epsst)*(t1(ij,1,k)+t1(ij,1,k-1)) 1140 continue DO 1160 K=1,NZ1 !4 do 1160 j=1,ny !4 do 1160 i=1,nx do 1160 j=1,ny-1 do 1160 i=1,nx-1 FPP(i,j,k)=FPP(i,j,k)-DTS*C2* * (RDX*(U1(i,j,k)-U1(I-1,j,k))+rdy*(v1(i,j,k)-v1(i,j-1,k))) 1160 continue !4-----set periodic points: fpp(nx,:,:) = fpp(1,:,:); fpp(:,ny,:) = fpp(:,1,:) DO 1170 K=2,NZ1 do 1170 ij=1,nxy W1(ij,1,K)=W1(ij,1,K)-.5*DTS*RDZ*(1.+EPSSM) * *(FPP(ij,1,K)-FPP(ij,1,K-1)) * +.25*dts*b2*(1.+epsst) * *(ftt(ij,1,k)+ftt(ij,1,k-1)) 1170 continue DO 1180 K=2,NZ1 do 1180 ij=1,nxy W1(ij,1,K)=(W1(ij,1,K)-A(K)*W1(ij,1,K-1))*ALPHA(K) 1180 continue DO 1210 K=NZ1,2,-1 do 1210 ij=1,nxy W1(ij,1,K)=W1(ij,1,K)-GAMMA(K)*W1(ij,1,K+1) 1210 continue do 1181 k=1,nz1 do 1181 ij=1,nxy t1(ij,1,k) = ftt(ij,1,k) * -.25*dts*(1.+epsst)*rdz* * ( w1(ij,1,k+1)*(tavg(k+1)-tavg(k )) * +w1(ij,1,k )*(tavg(k )-tavg(k-1)) ) 1181 continue !b8: divergence damping p1_old = p1 do 1211 k=nz1,1,-1 do 1211 ij=1,nxy P1(ij,1,K)=FPP(ij,1,K)-.5*DTS*C2*RDZ*(1.+EPSSM) * *(W1(ij,1,K+1)-W1(ij,1,K)) 1211 continue c c reset periodic boundaries c !4----------- p1(nx,:,:) = p1(1,:,:); p1(:,ny,:) = p1(:,1,:) w1(nx,:,:) = w1(1,:,:); w1(:,ny,:) = w1(:,1,:) t1(nx,:,:) = t1(1,:,:); t1(:,ny,:) = t1(:,1,:) !b8: divergence damping p1_old = p1 + divd_coeff * ( p1 - p1_old ) 1099 continue C C********RESET TIME LEVELS AND TIME SMOOTHING C ! XXX at this point, u1 etc are updated values at next time level ! and u etc are values at present time level (IF time smoothing ! is turned off). ! ! Note that vert adv of tavg has been added back in on small step. !-- Compute and save tendencies: fu_bst = ( u1 - u ) / dt ; fv_bst = ( v1 - v ) / dt ; fw_bst = ( w1 - w ) / dt ; ft_bst = ( t1 - t ) / dt ; fp_bst = ( p1 - p ) / dt ; ! XXX SHOULD DO SOME SPATIAL AVERAGING??? open(20,file='basic_state_tendencies.dat',form='unformatted') write(20) nx,ny,nz,xl,yl,ztop,dt,ns,gamma1,vsmoth,time write(20) fu_bst, fv_bst, fw_bst, ft_bst, fp_bst close(20) write(6,*) 'Forcing calculated:' write(6,*) ' fu_bst etc written to basic_state_tendencies.dat' else if ( iforce .ne. 1 .and. nit .eq. 1) then ! iforce not equal to 1; no forcing: write(6,*) 'No forcing: setting fu_bst etc to 0' fu_bst = 0; fv_bst = 0; fw_bst = 0; ft_bst = 0; fp_bst = 0; else if ( nit .gt. 1 ) then ! b.-st. tends already calculated; don't need to do anything !-- Basic-state calculations end here !-------- endif !-------- ! Linearized code begins here: debug = .true. if ( debug ) then write(6,*) '... beginning pertn calculations' write(6,*) ' ... step, time =', nit, time endif !-- write rms tp to fort.41 write(41,*) time, sqrt( sum( sum( sum( tp**2, 1), 1 ), 1 ) ), & sqrt( sum( sum( sum( up**2, 1), 1 ), 1 ) ), & sqrt( sum( sum( sum( vp**2, 1), 1 ), 1 ) ), & sqrt( sum( sum( sum( wp**2, 1), 1 ), 1 ) ) !-- write time series to fort.42: ! tp( Nx/2, Ny/2 + [20, 0, -20], 1 ) ! w ( Nx/4, Ny/2 + [20, 0, -20], 21) write(42,*) time, tp( Nx/2, Ny/2 + 20, 1 ), & tp( Nx/2, Ny/2 , 1 ), & tp( Nx/2, Ny/2 - 20, 1 ), & wp( Nx/4, Ny/2 + 20, 21), & wp( Nx/4, Ny/2 , 21), & wp( Nx/4, Ny/2 - 20, 21) !-- Save selected slices/fields: deltat= time -time0 if (nit.eq.0) then write(31,*) 'tp', nx,ny, time0, t_out write(32,*) 'wp', nx,ny, time0, t_out write(33,*) 'wp', nx,ny, time0, t_out write(34,*) 'wp', nx,ny, time0, t_out write(35,*) 'vortp', nx,ny, time0, t_out endif IF( abs( deltat -nint(deltat/t_out)*t_out ) .le. 0.5*dt ) THEN write(31,*) tp(:,:,1 ) write(32,*) wp(:,:,5 ) write(33,*) wp(:,:,21) write(34,*) wp(:,:,2 ) write(35,*) -( up(1:nx-1,2:ny ,1) - up(1:nx-1,1:ny-1,1) ) / dy & +( vp(2:nx ,1:ny-1,1) - vp(1:nx-1,1:ny-1,1) ) / dx end if !-- Large time step: ku = 1 kd = 2 do 85 ij=1,nxy wduz(ij,1,1) = 0. wdtz(ij,1,1) = 0. wdwz(ij,1,1) = 0. wdvz(ij,1,1) = 0. 85 continue c c loop over z c do 86 k=1,nz1 ktmp = kd kd = ku ku = ktmp kp1=min0(nz1,k+1) km1=max0(1,k-1) do 87 ij=1,nxy t1t(ij,1) = t1p(ij,1,k) 87 continue c c first, vertical advection c do 88 j=1,ny-1 do 88 i=1,nx-1 wduz(i,j,ku)=.5*( (Wp(i,j,k+1)+Wp(i+1,j,k+1)) & *RDZ*(U(i,j,kp1)-U(i,j,k)) & + ( W(i,j,k+1)+ W(i+1,j,k+1)) & *RDZ*(Up(i,j,kp1)-Up(i,j,k)) ) wdtz(i,j,ku)= & Wp(i,j,k+1)*RDZ*( t(i,j,kp1)-t(i,j,k) & -tavg(kp1) +tavg(k) ) & + W(i,j,k+1)*RDZ*( tp(i,j,kp1)-tp(i,j,k) ) wdwz(i,j,ku)= .5 * ( (Wp(i,j,k+1)+Wp(i,j,k)) & *RDZ*(w(i,j,k+1)-w(i,j,k)) & + (W(i,j,k+1)+W(i,j,k)) & *RDZ*(wp(i,j,k+1)-wp(i,j,k)) ) 88 continue do 89 j=1,ny-1 do 89 i=1,nx-1 wdvz(i,j,ku)= .5 *( (Wp(i,j,k+1)+Wp(i,j+1,k+1)) & *RDZ*(v(i,j,kp1)-v(i,j,k)) & + (W(i,j,k+1)+W(i,j+1,k+1)) & *RDZ*(vp(i,j,kp1)-vp(i,j,k)) ) 89 continue do 891 j=1,ny-1 do 891 i=1,nx-1 fu(i,j,k) = -.5*dts*(wduz(i,j,ku)+wduz(i,j,kd)) fv(i,j,k) = -.5*dts*(wdvz(i,j,ku)+wdvz(i,j,kd)) fw(i,j,k) = -.5*dts*(wdwz(i,j,ku)+wdwz(i,j,kd)) ft(i,j,k) = -.5*dtl*(wdtz(i,j,ku)+wdtz(i,j,kd)) 891 continue c c next, horizontal advection and 4rth order filter c c first, x-components c do 880 j=1,ny-1 jp1 = min0(j+1,ny) do 880 i=3,nx-2 c fu(i,j,k) = fu(i,j,k) * -rdx12*dts*f4asp( & 0.,u(I,J,K),u(I-2,J,K),u(I-1,J,K),u(I+1,J,K),u(I+2,J,K), & 0.,up(I,J,K),up(I-2,J,K),up(I-1,J,K),up(I+1,J,K),up(I+2,J,K)) c fv(i,j,k) = fv(i,j,k) * - rdx48*dts*f4asp( & U(I-1,jp1,K)+U(I-1,J,K), U(I,jp1,K)+U(I,J,K), & V(I-2,J,K), V(I-1,J,K), V(I+1,J,K), V(I+2,J,K) , & Up(I-1,jp1,K)+Up(I-1,J,K), Up(I,jp1,K)+Up(I,J,K), & Vp(I-2,J,K), Vp(I-1,J,K), Vp(I+1,J,K), Vp(I+2,J,K) ) fw(i,j,k) = fw(i,j,k) * -RDX48*dts*F4ASp( & U(I-1,J,K)+ u(i-1,j,km1), U(I,J,K)+ u(i,j,km1), & w(I-2,J,K),w(I-1,J,K), w(I+1,J,K), w(I+2,J,K) , & Up(I-1,J,K)+ up(i-1,j,km1), Up(I,J,K)+ up(i,j,km1), & wp(I-2,J,K),wp(I-1,J,K), wp(I+1,J,K), wp(I+2,J,K) ) ft(i,j,k) = ft(i,j,k) * -RDX24*dtl*F4ASp( & U(I-1,J,K), U(I,J,K), & t(I-2,J,K),t(I-1,J,K), t(I+1,J,K), t(I+2,J,K) , & Up(I-1,J,K), Up(I,J,K), & tp(I-2,J,K),tp(I-1,J,K), tp(I+1,J,K), tp(I+2,J,K) ) 880 continue c c use periodic b.c. where needed c do 881 ii=1,3 i = ii if(ii.eq.3) i = nx-1 im1 = i-1 im2 = i-2 ip1 = i+1 ip2 = i+2 if(im1.lt.1) im1=nx-2+i if(im2.lt.2) im2=nx-3+i if(ip2.gt.nx) ip2 = 2 do 882 j=1,ny-1 c fu(i,j,k) = fu(i,j,k) & -rdx12*dts*f4asp( & 0.,u(I,J,K),u(IM2,J,K),u(IM1,J,K), u(IP1,J,K),u(IP2,J,K), & 0.,up(I,J,K),up(IM2,J,K),up(IM1,J,K), up(IP1,J,K),up(IP2,J,K) ) c fv(i,j,k) = fv(i,j,k) * - rdx48*dts*f4asp( & U(IM1,j+1,K)+U(IM1,J,K), U(I,j+1,K)+U(I,J,K), & V(IM2,J,K), V(IM1,J,K), V(IP1,J,K), V(IP2,J,K) , & Up(IM1,j+1,K)+Up(IM1,J,K), Up(I,j+1,K)+Up(I,J,K), & Vp(IM2,J,K), Vp(IM1,J,K), Vp(IP1,J,K), Vp(IP2,J,K) ) c fw(i,j,k) = fw(i,j,k) * -RDX48*dts*F4ASp( & U(IM1,J,K)+ u(im1,j,km1), U(I,J,K)+ u(i,j,km1), & w(IM2,J,K),w(IM1,J,K), w(IP1,J,K), w(IP2,J,K) , & Up(IM1,J,K)+ up(im1,j,km1), Up(I,J,K)+ up(i,j,km1), & wp(IM2,J,K),wp(IM1,J,K), wp(IP1,J,K), wp(IP2,J,K) ) c ft(i,j,k) = ft(i,j,k) * -RDX24*dtl*F4ASp( & U(IM1,J,K), U(I,J,K), & t(IM2,J,K),t(IM1,J,K), t(IP1,J,K), t(IP2,J,K) , & Up(IM1,J,K), Up(I,J,K), & tp(IM2,J,K),tp(IM1,J,K), tp(IP1,J,K), tp(IP2,J,K) ) 882 continue 881 continue c c finished with x terms c next, y terms, first, 4rth order advection and filter in interior c do 890 j=3,ny-2 do 890 i=1,nx-1 fu(i,j,k) = fu(i,j,k) & -rdy48*dts*f4asp( & v(I,j-1,K)+v(I+1,J-1,K), v(I,j,K)+v(I+1,J,K), & u(I,J-2,K), u(I,J-1,K), u(I,J+1,K), u(I,J+2,K) , & vp(I,j-1,K)+vp(I+1,J-1,K), vp(I,j,K)+vp(I+1,J,K), & up(I,J-2,K), up(I,J-1,K), up(I,J+1,K), up(I,J+2,K) ) fv(i,j,k) = fv(i,j,k) * -rdy12*dts*f4asp( & 0.,v(I,J,K),v(I,j-2,K),v(I,j-1,K), v(I,j+1,K),v(I,j+2,K), & 0.,vp(I,J,K),vp(I,j-2,K),vp(I,j-1,K), vp(I,j+1,K),vp(I,j+2,K)) fw(i,j,k) = fw(i,j,k) * -RDy48*dts*F4ASp( & v(I,J-1,K)+ v(i,j-1,km1), v(I,J,K)+ v(i,j,km1), & w(I,J-2,K),w(I,J-1,K), w(I,J+1,K), w(I,J+2,K) , & vp(I,J-1,K)+ vp(i,j-1,km1), vp(I,J,K)+ vp(i,j,km1), & wp(I,J-2,K),wp(I,J-1,K), wp(I,J+1,K), wp(I,J+2,K) ) ft(i,j,k) = ft(i,j,k) * -RDy24*dtl*F4ASp( & v(I,j-1,K), v(I,J,K), & t(I,j-2,K),t(I,j-1,K), t(I,j+1,K), t(I,j+2,K) , & vp(I,j-1,K), vp(I,J,K), & tp(I,j-2,K),tp(I,j-1,K), tp(I,j+1,K), tp(I,j+2,K) ) 890 continue c c apply periodic boundary conditions in y c do 8911 jj=1,3 j=jj if(jj .eq. 3) j=ny-1 jm1 = j-1 jm2 = j-2 jp1 = j+1 jp2 = j+2 if (jm1.lt.1) jm1 = ny-2+j if (jm2.lt.1) jm2 = ny-3+j if (jp2.gt.ny) jp2 = 2 do 8911 i=1,nx-1 fu(i,j,k) = fu(i,j,k) * -rdy48*dts*f4asp( & v(I,jm1,K)+v(I+1,Jm1,K), v(I,j,K)+v(I+1,J,K), & u(I,Jm2,K), u(I,Jm1,K), u(I,Jp1,K), u(I,Jp2,K) , & vp(I,jm1,K)+vp(I+1,Jm1,K), vp(I,j,K)+vp(I+1,J,K), & up(I,Jm2,K), up(I,Jm1,K), up(I,Jp1,K), up(I,Jp2,K) ) fv(i,j,k) = fv(i,j,k) * -rdy12*dts*f4asp( & 0.,v(I,J,K),v(I,jm2,K),v(I,jm1,K), v(I,jp1,K),v(I,jp2,K), & 0.,vp(I,J,K),vp(I,jm2,K),vp(I,jm1,K), vp(I,jp1,K),vp(I,jp2,K)) fw(i,j,k) = fw(i,j,k) * -RDy48*dts*F4ASp( & v(I,Jm1,K)+ v(i,jm1,km1), v(I,J,K)+ v(i,j,km1), & w(I,Jm2,K),w(I,Jm1,K), w(I,Jp1,K), w(I,Jp2,K) , & vp(I,Jm1,K)+ vp(i,jm1,km1), vp(I,J,K)+ vp(i,j,km1), & wp(I,Jm2,K),wp(I,Jm1,K), wp(I,Jp1,K), wp(I,Jp2,K) ) ft(i,j,k) = ft(i,j,k) * -RDy24*dtl*F4ASp( & v(I,jm1,K), v(I,J,K), & t(I,jm2,K),t(I,jm1,K), t(I,jp1,K), t(I,jp2,K) , & vp(I,jm1,K), vp(I,J,K), & tp(I,jm2,K),tp(I,jm1,K), tp(I,jp1,K), tp(I,jp2,K) ) 8911 continue !9-- Horizontal diffusion: call horizontal_deln( i_diff_order, nx,ny,nz, & gamma1, gammaw, gammat, rns, & fu(:,:,k), fv(:,:,k), fw(:,:,k), ft(:,:,k), & u1p(:,:,k), v1p(:,:,k), w1p(:,:,k), t1t ) c c the coriolis terms and other terms from the spherical grid. c do 8110 j=1,ny-1 do 8110 i=1,nx-1 ! fu(i,j,k) = fu(i,j,k) + dts*(fcorm(j)* fu(i,j,k) = fu(i,j,k) + dts*( f0 * * 0.25*(vp(i,j,k)+vp(i,j-1,k)+vp(i+1,j,k)+vp(i+1,j-1,k)) ) ! fv(i,j,k) = fv(i,j,k) - dts*fcorv(j)*( fv(i,j,k) = fv(i,j,k) - dts* f0 *( * 0.25*(up(i,j,k)+up(i,j+1,k)+up(i-1,j,k)+up(i-1,j+1,k)) ) 8110 continue !--- Sponge layer (19 Aug 03) if ( k .gt. nz1 - kdepth ) then fac = dts * ( k - (nz1 - kdepth) ) / kdepth / tau ! write(6,*) k, fac, (dtl/dts) * fac fu(:,:,k) = fu(:,:,k) - fac * u1p(:,:,k) fv(:,:,k) = fv(:,:,k) - fac * v1p(:,:,k) fw(:,:,k) = fw(:,:,k) - fac * w1p(:,:,k) fac = (dtl/dts) * fac ft(:,:,k) = ft(:,:,k) - fac * ( t1p(:,:,k) ) endif !--- Forcing by "residual" basic-state tendencies ! (Note that fu_bst etc are really tendencies, and so get multiplied by dtl or dts) fu = fu - dts * fu_bst fv = fv - dts * fv_bst fw = fw - dts * fw_bst ft = ft - dtl * ft_bst ! "fp": pressure forcing fp_bst added on small step c c set periodic boundary conditions c do 8120 j=1,ny fu(nx,j,k) = fu(1,j,k) fu(0,j,k) = fu(nx-1,j,k) fv(nx,j,k) = fv(1,j,k) fw(nx,j,k) = fw(1,j,k) ft(nx,j,k) = ft(1,j,k) 8120 continue do 8121 i=1,nx fu(i,ny,k) = fu(i,1,k) fw(i,ny,k) = fw(i,1,k) ft(i,ny,k) = ft(i,1,k) fv(i,ny,k) = fv(i,1,k) fv(i,0,k) = fv(i,ny-1,k) 8121 continue ! fu(0,ny,k) = fu(0,1,k) c 86 continue c c done computing largetimestep stuff, reset temp and set up small timestep c c first, do partial time filter for middle level variables u,v,w,p c do 8130 ijk=1,nx*ny*nz1 wp(ijk,1,1) = wp(ijk,1,1) + epsts*(w1p(ijk,1,1)-2*wp(ijk,1,1)) pp(ijk,1,1) = pp(ijk,1,1) + epsts*(p1p(ijk,1,1)-2*pp(ijk,1,1)) tp(ijk,1,1) = tp(ijk,1,1) + epsts*(t1p(ijk,1,1)-2*tp(ijk,1,1)) 8130 continue do 8131 ijk=0,(nx+1)*ny*nz1 up(ijk,1,1) = up(ijk,1,1) + epsts*(u1p(ijk,1,1)-2*up(ijk,1,1)) 8131 continue do 8132 ijk=1,nx*(ny+1)*nz1 vp(ijk,0,1) = vp(ijk,0,1) + epsts*(v1p(ijk,0,1)-2*vp(ijk,0,1)) 8132 continue c c set zero b.c. for w c do 8133 ij=1,nx*ny fw(ij,1,1) = 0. fw(ij,1,nz) = 0. wp(ij,1,1) = 0. wp(ij,1,nz) = 0. w1p(ij,1,1) = 0. w1p(ij,1,nz) = 0. 8133 continue vsm4 = -vsmoth*rns vsm2 = -4.*vsm4 do 8134 k=3,nz-2 do 8134 ij=1,nxy fw(ij,1,k) = fw(ij,1,k) + vsm4*( 1 w1p(ij,1,k+2)+w1p(ij,1,k-2)+6*w1p(ij,1,k) 2 -4*(w1p(ij,1,k-1)+w1p(ij,1,k+1)) ) 8134 continue k=2 do 8135 ij=1,nxy fw(ij,1,k) = fw(ij,1,k) + vsm2*( 1 w1p(ij,1,k+1)-2*w1p(ij,1,k)+w1p(ij,1,k-1) ) 8135 continue k=nz-1 do 8136 ij=1,nxy fw(ij,1,k) = fw(ij,1,k) + vsm2*( 1 w1p(ij,1,k+1)-2*w1p(ij,1,k)+w1p(ij,1,k-1) ) 8136 continue C C********SMALL TIME STEP CALCULATIONS C p1_old = p1p ; !p8: for divergence damping DO 8099 M=1,NS c c advance u c do 8210 k=1,nz1 do 8211 j=1,ny-1 do 8211 i=1,nx-1 !p8 U1p(i,j,k)=U1p(i,j,k)-DTS*RDX*(P1p(I+1,j,k)-P1p(I,j,k)) U1p(i,j,k)=U1p(i,j,k)-DTS*RDX*(p1_old(I+1,j,k)-p1_old(I,j,k)) * +FU(I,j,k) 8211 continue do 8212 j=1,ny u1p(nx,j,k) = u1p(1,j,k) u1p(0,j,k) = u1p(nx-1,j,k) 8212 continue do 8213 i=0,nx u1p(i,ny,k) = u1p(i,1,k) 8213 continue 8210 continue c c advance v c do 8215 k=1,nz1 do 8215 j=1,ny-1 do 8215 i=1,nx-1 !p8 v1p(i,j,k)=v1p(i,j,k)-DTS*RDy*(P1p(I,j+1,k)-P1p(I,j,k)) v1p(i,j,k)=v1p(i,j,k)-DTS*RDy*(p1_old(I,j+1,k)-p1_old(I,j,k)) * +Fv(I,j,k) 8215 continue do 8216 k=1,nz1 do 8216 j=0,ny v1p(nx,j,k) = v1p(1,j,k) 8216 continue do 8217 k=1,nz1 do 8217 i=1,nx v1p(i,0,k) = v1p(i,ny-1,k) v1p(i,ny,k) = v1p(i,1,k) 8217 continue c c advance w and p implicitly c DO 8230 K=1,NZ1 do 8230 ij=1,nxy FPP(ij,1,k)=P1p(ij,1,k)-.5*DTS*C2*RDZ*(1.-EPSSM) * *(W1p(ij,1,k+1)-W1p(ij,1,k)) & - dts * fp_bst(ij,1,k) 8230 continue do 8231 k=1,nz1 km1 = max0(1,k-1) kp1 = min0(nz1,k+1) do 8231 ij=1,nxy ftt(ij,1,k) = t1p(ij,1,k)+ft(ij,1,k)*rns * -.25*dts*(1.-epsst)*rdz* * ( w1p(ij,1,k+1)*(tavg(kp1)-tavg(k)) * +w1p(ij,1,k )*(tavg(k)-tavg(km1)) ) 8231 continue c DO 8240 K=2,NZ1 do 8240 ij=1,nxy W1p(ij,1,k)=W1(ij,1,k)+FW(ij,1,k) * -.5*DTS*RDZ*(1.-EPSSM)*(P1p(ij,1,k)-P1p(ij,1,k-1)) * +.25*dts*b2*(1.-epsst)*(t1p(ij,1,k)+t1p(ij,1,k-1)) 8240 continue DO 8260 K=1,NZ1 do 8260 j=1,ny-1 do 8260 i=1,nx-1 FPP(i,j,k)=FPP(i,j,k)-DTS*C2* * (RDX*(U1p(i,j,k)-U1p(I-1,j,k))+rdy*(v1p(i,j,k)-v1p(i,j-1,k))) 8260 continue !4-----set periodic points: fpp(nx,:,:) = fpp(1,:,:); fpp(:,ny,:) = fpp(:,1,:) DO 8270 K=2,NZ1 do 8270 ij=1,nxy W1p(ij,1,K)=W1p(ij,1,K)-.5*DTS*RDZ*(1.+EPSSM) * *(FPP(ij,1,K)-FPP(ij,1,K-1)) * +.25*dts*b2*(1.+epsst) * *(ftt(ij,1,k)+ftt(ij,1,k-1)) 8270 continue DO 8280 K=2,NZ1 do 8280 ij=1,nxy W1p(ij,1,K)=(W1p(ij,1,K)-A(K)*W1p(ij,1,K-1))*ALPHA(K) 8280 continue DO 8281 K=NZ1,2,-1 do 8281 ij=1,nxy W1p(ij,1,K)=W1p(ij,1,K)-GAMMA(K)*W1p(ij,1,K+1) 8281 continue do 8282 k=1,nz1 do 8282 ij=1,nxy t1p(ij,1,k) = ftt(ij,1,k) * -.25*dts*(1.+epsst)*rdz* * ( w1p(ij,1,k+1)*(tavg(k+1)-tavg(k )) * +w1p(ij,1,k )*(tavg(k )-tavg(k-1)) ) 8282 continue !p8: divergence damping p1_old = p1p do 8311 k=nz1,1,-1 do 8311 ij=1,nxy P1p(ij,1,K)=FPP(ij,1,K)-.5*DTS*C2*RDZ*(1.+EPSSM) * *(W1p(ij,1,K+1)-W1p(ij,1,K)) 8311 continue c c reset periodic boundaries c p1p(nx,:,:) = p1p(1,:,:); p1p(:,ny,:) = p1p(:,1,:) w1p(nx,:,:) = w1p(1,:,:); w1p(:,ny,:) = w1p(:,1,:) t1p(nx,:,:) = t1p(1,:,:); t1p(:,ny,:) = t1p(:,1,:) !p8: divergence damping p1_old = p1p + divd_coeff * ( p1p - p1_old ) 8099 continue C C********RESET TIME LEVELS AND TIME SMOOTHING C do 822 ijk=0,(nxy+ny)*nz1-1 fu(ijk,1,1) = up(ijk,1,1) + epsts*u1p(ijk,1,1) 822 continue do 823 ijk=0,(nxy+ny)*nz1-1 up(ijk,1,1) = u1p(ijk,1,1) u1p(ijk,1,1) = fu(ijk,1,1) 823 continue do 824 ijk=1,(nxy+nx)*nz1 fv(ijk,0,1) = vp(ijk,0,1) + epsts*v1p(ijk,0,1) 824 continue do 825 ijk=1,(nxy+nx)*nz1 vp(ijk,0,1) = v1p(ijk,0,1) v1p(ijk,0,1) = fv(ijk,0,1) 825 continue do 826 ijk=1,nxy*nz fw(ijk,1,1) = wp(ijk,1,1) + epsts*w1p(ijk,1,1) 826 continue do 827 ijk=1,nxy*nz wp(ijk,1,1) = w1p(ijk,1,1) w1p(ijk,1,1) = fw(ijk,1,1) 827 continue do 828 ijk=1,nxy*nz1 fw(ijk,1,1) = pp(ijk,1,1) + epsts*p1p(ijk,1,1) 828 continue do 829 ijk=1,nxy*nz1 pp(ijk,1,1) = p1p(ijk,1,1) p1p(ijk,1,1) = fw(ijk,1,1) 829 continue do 831 ijk=1,nxy*nz1 fw(ijk,1,1) = tp(ijk,1,1) + epsts*t1p(ijk,1,1) 831 continue do 832 ijk=1,nxy*nz1 tp(ijk,1,1) = t1p(ijk,1,1) t1p(ijk,1,1) = fw(ijk,1,1) 832 continue c c check pmax and wmax c ! if(mod(nit,10) .eq. 1) then if(mod(nit,1) .eq. 0) then pmaxt = -1.e+20 wmaxt = -1.e+20 vmaxt = -1.e+20 umaxt = -1.e+20 do 821 ijk=1,nx*ny*nz1 pmaxt = amax1(pmaxt,abs(pp(ijk,1,1))) wmaxt = amax1(wmaxt,abs(wp(ijk,1,1))) 821 continue do 8220 ij=1,nx*(ny+1) vmaxt = amax1(vmaxt,abs(vp(ij,0,1))) 8220 continue umaxt = maxval(up) pmax = pmaxt wmax = wmaxt write(6,'('' pert: '', * 5(1x,e10.3))') time,pmax,wmax,umaxt ,vmaxt if(wmax .gt. 1000.) then write(6,*) ' !! perturbation blow up !! step = ', nit rewind(9) write(9)nx,ny,nz1,time write(9)up,u1p,vp,v1p,wp,w1p,pp,p1p,tp,t1p stop end if endif DTL=2.*DT NS=2*NS0 ! Linearized code end here !-------- c c take another timestep if necessary c if(nit .lt. nstopc) go to 3000 ! End time loop 7001 continue call wrt_data('fort.9', u,u1,v,v1,w,w1,p,p1,t,t1, time, & nx,ny,nz, xl,yl,ztop, f0,Nbv,gbyt0) call wrt_data('pert.9', up,u1p,vp,v1p,wp,w1p,pp,p1p,tp,t1p, time, & nx,ny,nz, xl,yl,ztop, f0,Nbv,gbyt0) open(20,file='uvwtp.dat',form='unformatted') write(20) nx,ny,nz,xl,yl,ztop,dt,ns,gamma1,vsmoth,time write(20) u,v,w,t,p close(20) open(20,file='pert_uvwtp.dat',form='unformatted') write(20) nx,ny,nz,xl,yl,ztop,dt,ns,gamma1,vsmoth,time write(20) up,vp,wp,tp,pp close(20) END !--------------------------------------------------------------- subroutine rd_data(file, u,u1,v,v1,w,w1,p,p1,t,t1, time, & nx,ny,nz, xl,yl,ztop, f0,Nbv,gbyt0) ! ! Reads state, time, parameters from specified file. ! implicit none character(*), intent(in) :: file integer, intent(in) :: nx,ny,nz real, intent(in) :: xl,yl,ztop, f0,Nbv,gbyt0 real, intent(out) :: time, ! u(0:nx,ny,nz-1),u1(0:nx,ny,nz-1), * v(nx,0:ny,nz-1),v1(nx,0:ny,nz-1), * w(nx, ny,nz ),w1(nx, ny,nz ), * p(nx, ny,nz-1),p1(nx, ny,nz-1), * t(nx, ny,nz-1),t1(nx, ny,nz-1) !! real, intent(out) :: time !* Local vars integer nxt,nyt,nzt, i,j,k real xlt,ylt,ztopt, f0t,Nbvt,gbyt0t write(6,*) 'In rd_data, Nbv=',Nbv open(11,file=file, status='old', err=100) write(6,*) 'Opened file ', file read(11,*) nxt,nyt,nzt write(6,*) nxt,nyt,nzt if( (nxt.ne.nx) .or. (nyt.ne.ny) .or. (nzt.ne.nz) ) then write(6,*) 'CHECK DIMENSIONS IN ', file write(6,*) ' nxt,nyt,nzt=', nxt,nyt,nzt write(6,*) ' nx ,ny ,nz =', nx ,ny ,nz stop end if read(11,*) xlt,ylt,ztopt, f0t,Nbvt,gbyt0t if( (xlt.ne.xl) .or. (ylt.ne.yl) .or. (ztopt.ne.ztop) & .or. (f0t.ne.f0) .or. (Nbvt.ne.Nbv) ) & then !!! & .or. (f0t.ne.f0) .or. (Nbvt.ne.Nbv) .or. (gbyt0t.ne.gbyt0) ) write(6,*) 'CHECK PARAMETERS IN ', file write(6,*) ' xlt,ylt,ztopt, f0t,Nbvt,gbyt0t=', & xlt,ylt,ztopt, f0t,Nbvt,gbyt0t write(6,*) ' xl ,yl ,ztop , f0 ,Nbv ,gbyt0 =', & xl ,yl ,ztop , f0 ,Nbv ,gbyt0 !! stop end if read(11,*) time write(6,*) '...reading data for time ', time do k=1,nz-1; do j=1,ny; do i=0,nx; read(11,*) u(i,j,k); enddo; enddo; enddo; !!! read(11,*) v(j,i,k); enddo; enddo; enddo; !! read(11,*) u(i,j,k); write(6,*) i,j,k,u(i,j,k) !! enddo; enddo; enddo; write(6,*) ' u read' do k=1,nz-1; do j=0,ny; do i=1,nx; read(11,*) v(i,j,k); enddo; enddo; enddo; !!! read(11,*) u(j,i,k); enddo; enddo; enddo; write(6,*) ' v read' do k=1,nz; do j=1,ny; do i=1,nx; read(11,*) w(i,j,k); enddo; enddo; enddo; !!! read(11,*) w(j,i,k); enddo; enddo; enddo; write(6,*) ' w read' do k=1,nz-1; do j=1,ny; do i=1,nx; read(11,*) t(i,j,k); enddo; enddo; enddo; !!! read(11,*) t(j,i,k); enddo; enddo; enddo; write(6,*) ' t read' do k=1,nz-1; do j=1,ny; do i=1,nx; read(11,*) p(i,j,k); enddo; enddo; enddo; !!! read(11,*) p(j,i,k); enddo; enddo; enddo; write(6,*) ' p read' do k=1,nz-1; do j=1,ny; do i=0,nx; read(11,*,end=99) u1(i,j,k); enddo; enddo; enddo; !* If EOF found, then missing 2nd time level; go to 99 do k=1,nz-1; do j=0,ny; do i=1,nx; read(11,*) v1(i,j,k); enddo; enddo; enddo; do k=1,nz; do j=1,ny; do i=1,nx; read(11,*) w1(i,j,k); enddo; enddo; enddo; do k=1,nz-1; do j=1,ny; do i=1,nx; read(11,*) t1(i,j,k); enddo; enddo; enddo; do k=1,nz-1; do j=1,ny; do i=1,nx; read(11,*) p1(i,j,k); enddo; enddo; enddo; close(11) return 99 write(6,*) & 'Only single time level available; initializing u1=u, etc' u1=u; v1=v; w1=w; p1=p; t1=t; return 100 write(6,*) file, ' NOT FOUND' stop end !--------------------------------------------------------------- subroutine wrt_data(file, u,u1,v,v1,w,w1,p,p1,t,t1, time, & nx,ny,nz, xl,yl,ztop, f0,Nbv,gbyt0) ! ! Writes state, time, parameters to specified file. ! implicit none character(*), intent(in) :: file integer, intent(in) :: nx,ny,nz real, intent(in) :: xl,yl,ztop, f0,Nbv,gbyt0, time, & u(0:nx,ny,nz-1),u1(0:nx,ny,nz-1), * v(nx,0:ny,nz-1),v1(nx,0:ny,nz-1), * w(nx, ny,nz ),w1(nx, ny,nz ), * p(nx, ny,nz-1),p1(nx, ny,nz-1), * t(nx, ny,nz-1),t1(nx, ny,nz-1) !* Local vars integer i,j,k open(11,file=file, status='replace') write(11,*) nx,ny,nz write(11,*) xl,yl,ztop, f0,Nbv,gbyt0 write(11,*) time write(6,*) '...writing data for time ', time do k=1,nz-1; do j=1,ny; do i=0,nx; write(11,*) u(i,j,k); enddo; enddo; enddo; do k=1,nz-1; do j=0,ny; do i=1,nx; write(11,*) v(i,j,k); enddo; enddo; enddo; do k=1,nz; do j=1,ny; do i=1,nx; write(11,*) w(i,j,k); enddo; enddo; enddo; do k=1,nz-1; do j=1,ny; do i=1,nx; write(11,*) t(i,j,k); enddo; enddo; enddo; do k=1,nz-1; do j=1,ny; do i=1,nx; write(11,*) p(i,j,k); enddo; enddo; enddo; do k=1,nz-1; do j=1,ny; do i=0,nx; write(11,*) u1(i,j,k); enddo; enddo; enddo; do k=1,nz-1; do j=0,ny; do i=1,nx; write(11,*) v1(i,j,k); enddo; enddo; enddo; do k=1,nz; do j=1,ny; do i=1,nx; write(11,*) w1(i,j,k); enddo; enddo; enddo; do k=1,nz-1; do j=1,ny; do i=1,nx; write(11,*) t1(i,j,k); enddo; enddo; enddo; do k=1,nz-1; do j=1,ny; do i=1,nx; write(11,*) p1(i,j,k); enddo; enddo; enddo; close(11) return end !------------------------------------------------------------------------ function gaussdev(idum) ! Returns a normally distributed deviate with 0 mean and unit variance, ! using ran3(idum) as the source of uniform deviates; see Num'l Recipes, ! ch 7.2. integer idum real gaussdev !- local variables integer iset real fac,gset,rsq,v1,v2,ran1 save iset, gset data iset/0/ if (iset.eq.0) then 10 v1 = 2.*ran3(idum)-1. v2 = 2.*ran3(idum)-1. rsq = v1**2 + v2**2 if (rsq.ge.1. .or. rsq.eq.0.) goto 10 fac = sqrt( -2.*log(rsq)/rsq ) gset = v1*fac gaussdev = v2*fac iset= 1 else gaussdev = gset iset=0 end if return end !------------------------------------------------------------------------ FUNCTION ran3(idum) INTEGER idum INTEGER MBIG,MSEED,MZ ! REAL MBIG,MSEED,MZ REAL ran3,FAC PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG) ! PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG) INTEGER i,iff,ii,inext,inextp,k INTEGER mj,mk,ma(55) ! REAL mj,mk,ma(55) SAVE iff,inext,inextp,ma DATA iff /0/ if(idum.lt.0.or.iff.eq.0)then iff=1 mj=MSEED-iabs(idum) mj=mod(mj,MBIG) ma(55)=mj mk=1 do 11 i=1,54 ii=mod(21*i,55) ma(ii)=mk mk=mj-mk if(mk.lt.MZ)mk=mk+MBIG mj=ma(ii) 11 continue do 13 k=1,4 do 12 i=1,55 ma(i)=ma(i)-ma(1+mod(i+30,55)) if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG 12 continue 13 continue inext=0 inextp=31 idum=1 endif inext=inext+1 if(inext.eq.56)inext=1 inextp=inextp+1 if(inextp.eq.56)inextp=1 mj=ma(inext)-ma(inextp) if(mj.lt.MZ)mj=mj+MBIG ma(inext)=mj ran3=mj*FAC return END