! THIS VERSION: ! 9: * Removes old del2/del4(ish ... d^4/dx^4 + d^4/dy^4) horizontal diffusion ! that was embedded with advections. Replaces with horizontal_deln subroutine. ! ! THIS VERSION: ! 8: * Adds divergence damping by using a lagged p in p grad on small step ! ! Outputs time series at high time resolution ! ! THIS VERSION: Add sponge layer ! ! 7a:* Just use tavg = horizontal mean t (don't make theta0 the surface temp) ! ... doesn't make any difference, i.e. solution independent of adding a ! constant to theta(x,y,z) ! ! 6: * Read ICs then set tavg(k) (used to treat stratification ! on small step) to area average of theta. This change fixes ! slow instability that appeared in w. (tavg(k) differed from ! avg theta in ICs by factor of 9.81/9.8 ... different g's used.) ! ! 5: * fixed bugs in y BC for diffusion (see !5 lines below) ! ! 4: * changed various bits of periodic conditions; didn't make ! any difference (see !4 lines below) ! ! 2: * Expt 1: original (chnp.new_1.f) ! * Expt 2: in coriolis terms, use f0 instead of fcorv ! (fcorm already replaced by f0 in expt 1) ! * Expt 3: remove following from loop that sets periodicity ! in tendencies ! fu(0,ny,k) = fu(0,1,k) ! ! 1: * Remove special y BC loop for t advections on large step ! * Check for periodicity at beginning of each step ! ! THIS VERSION: ! ---Arrays uinit, tinit removed ! ---New parameters: f0, Nbv, gbyt0, theta0, ! xl, yl, ztop (old) 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=121,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) !!!!! & xl = 1500.e3,! x periodicity (m) !!!!! & yl = 1500.e3,! y periodicity (m) & ztop = 15.e3, ! height of lid (m) !! & ztop = 30.e3, ! height of lid (m) !!!!! & ztop = 7.5e3, ! height of lid (m) & tau = 1.e3, ! 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) !!!! TINIT(NY,NZ1),UINIT(NY,NZ1) 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 real zeta(1:nx-1,1:ny-1), tmean, zetamean ! for calculating vorticity for output real tmp(-1:nx+1, -1:ny+1) integer i_diff_order 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)) 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 write(6,*) 'Enter t_out for time series' read (5,*) t_out2 ! DT=900. ! ns0=4 c-- Viscosities; use +/- depending on 2nd/4th order. gamma1 = 0. ! horizontal filtering coeff for u,v gammaw = 0. ! ... for w gammat = 0. ! ... for t gamma1 = -0.001 ! horizontal filtering coeff for u,v gammaw = -0.001 ! ... for w gammat = -0.001 ! ... for t vsmoth = 0.007 ! vertical filtering coeff for w write(6,*) 'Enter gamma' read (5,*) rjunk ; gamma1= rjunk; gammaw= rjunk; gammat= rjunk; !9: set order of diffusion i_diff_order = 4 write(6,*) ' ... i_diff_order = ', i_diff_order !8: for divergence damping divd_coeff = 0.1 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 write(6,*) 'In main, Nbv=',Nbv 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) time0= time ! horizontal mean theta; set t = t1 = tbar tbar = sum( sum( t, 1), 1) / (nx*ny) ! write(6,*) '***SETTING u,v, and pertn theta to zero !!!!!' ! u = 0; u1 = 0; ! v = 0; v1 = 0; ! do k=1,nz1; t(:,:,k) = tbar(k); t1(:,:,k) = tbar(k); enddo !-- 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) 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 !7a tavg(k) = theta0 -tbar(1) +tbar(k) 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 ! write(15,'('' maximum n2 = '',e12.5)')an2max ! do 26 k=1,nz1 ! write(15,'('' avg n2,t at level '',2e12.5,2x,i4)')an2(k), ! * tavg(k),k !26 continue 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 if (0.eq.1) then !- check for periodicity do k=1,nz1 do i=1,nx if (u(i,1,k).ne.u(i,ny,k)) then write(6,*) 'u not periodic: i,k u=', i,k,u(i,1,k),u(i,ny,k) endif if (v(i,0,k).ne.v(i,ny-1,k)) then write(6,*) 'v0 not periodic: i,k v=',i,k,v(i,0,k),v(i,ny-1,k) endif if (v(i,1,k).ne.v(i,ny,k)) then write(6,*) 'v not periodic: i,k v=', i,k,v(i,1,k),v(i,ny,k) endif if (w(i,1,k).ne.w(i,ny,k)) then write(6,*) 'w not periodic: i,k w=', i,k,w(i,1,k),w(i,ny,k) endif if (t(i,1,k).ne.t(i,ny,k)) then write(6,*) 't not periodic: i,k t=', i,k,t(i,1,k),t(i,ny,k) endif enddo do j=1,ny if (u(0,j,k).ne.u(nx-1,j,k)) then write(6,*) 'u0 not periodic: j,k u=', j,k,u(0,j,k),u(nx-1,j,k) endif if (u(1,j,k).ne.u(nx,j,k)) then write(6,*) 'u1 not periodic: j,k u=', j,k,u(1,j,k),u(nx,j,k) endif if (v(1,j,k).ne.v(nx,j,k)) then write(6,*) 'v not periodic: j,k v=', j,k,v(1,j,k),v(nx,j,k) endif if (w(1,j,k).ne.w(nx,j,k)) then write(6,*) 'w not periodic: j,k w=', j,k,w(1,j,k),w(nx,j,k) endif if (t(1,j,k).ne.t(nx,j,k)) then write(6,*) 't not periodic: j,k t=', j,k,t(1,j,k),t(nx,j,k) endif enddo enddo endif 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 IF( abs( deltat -nint(deltat/t_out2)*t_out2 ) .le. 0.5*dt ) THEN tmean = sum( t(1:nx-1,1:ny-1, 1) ) / ( (nx-1) * (ny-1) ) write(31,*) maxval(t(:,:,1)), minval(t(:,:,1)), & sqrt( sum( (t(:,:,1) - tmean )**2 ) / ( (nx-1) * (ny-1) ) ) write(32,*) maxval(w(:,:,5)), minval(w(:,:,5)), & sqrt( sum( (w(:,:,5) )**2 ) / ( (nx-1) * (ny-1) ) ) write(33,*) maxval(w(:,:,21)), minval(w(:,:,21)), & sqrt( sum( (w(:,:,21) )**2 ) / ( (nx-1) * (ny-1) ) ) write(34,*) maxval(w(:,:,2)), minval(w(:,:,2)), & sqrt( sum( (w(:,:,2) )**2 ) / ( (nx-1) * (ny-1) ) ) zeta = & -( 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 zetamean = sum( zeta ) / ( (nx-1) * (ny-1) ) write(35,*) maxval(zeta), minval(zeta), & sqrt( sum( ( zeta - zetamean )**2 ) / ( (nx-1) * (ny-1) ) ) write(36,*) tmean, zetamean !!---More diagnostics: Would like to calc position and shape of dipole ! at each time: calc center of mass of cycl. and ! anti separately (vals > 0.1 of variation from mean ! temp), then calc their eccentricity and orientation !!---More diagnostics: ??? SLICK WAY TO DEAL WITH WRAP AROUND??? ! t1d = 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 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 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 and 4rth order filter 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: i_diff_order = 2 or 4 ! write(6,*) ' ' ! write(6,*) ' fu(20,25,k) = ' , fu(20,25,k) ! write(6,*) ' ft(20,25,k) = ' , ft(20,25,k) ! write(6,*) ' v1(88, 9,k) = ' , v1(88, 9,k) ! write(6,*) ' w1(12,65,k) = ' , w1(12,65,k) 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 ) 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 !--- 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 !7a ft(:,:,k) = ft(:,:,k) - fac * ( t1(:,:,k) - tbar(k) ) 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 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 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 ; !8: 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 !8 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 !8 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 !8: 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 do 1212 k=1,nz1 !4 do 1212 j=1,ny !4 p1(nx,j,k)=0.5*(p1(nx,j,k)+p1(1,j,k)) !4 p1(1,j,k) = p1(nx,j,k) !4 w1(nx,j,k)=0.5*(w1(nx,j,k)+w1(1,j,k)) !4 w1(1,j,k) = w1(nx,j,k) !4 t1(nx,j,k)=0.5*(t1(nx,j,k)+t1(1,j,k)) !4 t1(1,j,k)=t1(nx,j,k) !41212 continue !4 do 1213 k=1,nz1 !4 do 1213 i=1,nx !4 p1(i,ny,k)=0.5*(p1(i,ny,k)+p1(i,1,k)) !4 p1(i,1,k) = p1(i,ny,k) !4 w1(i,ny,k)=0.5*(w1(i,ny,k)+w1(i,1,k)) !4 w1(i,1,k) = w1(i,ny,k) !4 t1(i,ny,k)=0.5*(t1(i,ny,k)+t1(i,1,k)) !4 t1(i,1,k)= t1(i,ny,k) !41213 continue !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,:) !8: divergence damping p1_old = p1 + divd_coeff * ( p1 - p1_old ) ! p1_old = p1 + divd_coeff * ( p1_old - p1 ) 1099 continue C C********RESET TIME LEVELS AND TIME SMOOTHING C do 222 ijk=0,(nxy+ny)*nz1-1 fu(ijk,1,1) = u(ijk,1,1) + epsts*u1(ijk,1,1) 222 continue do 223 ijk=0,(nxy+ny)*nz1-1 u(ijk,1,1) = u1(ijk,1,1) u1(ijk,1,1) = fu(ijk,1,1) 223 continue do 224 ijk=1,(nxy+nx)*nz1 fv(ijk,0,1) = v(ijk,0,1) + epsts*v1(ijk,0,1) 224 continue do 225 ijk=1,(nxy+nx)*nz1 v(ijk,0,1) = v1(ijk,0,1) v1(ijk,0,1) = fv(ijk,0,1) 225 continue do 226 ijk=1,nxy*nz fw(ijk,1,1) = w(ijk,1,1) + epsts*w1(ijk,1,1) 226 continue do 227 ijk=1,nxy*nz w(ijk,1,1) = w1(ijk,1,1) w1(ijk,1,1) = fw(ijk,1,1) 227 continue do 228 ijk=1,nxy*nz1 fw(ijk,1,1) = p(ijk,1,1) + epsts*p1(ijk,1,1) 228 continue do 229 ijk=1,nxy*nz1 p(ijk,1,1) = p1(ijk,1,1) p1(ijk,1,1) = fw(ijk,1,1) 229 continue do 231 ijk=1,nxy*nz1 c t1(ijk,1,1) = t1(ijk,1,1)+ft(ijk,1,1) fw(ijk,1,1) = t(ijk,1,1) + epsts*t1(ijk,1,1) 231 continue do 232 ijk=1,nxy*nz1 t(ijk,1,1) = t1(ijk,1,1) t1(ijk,1,1) = fw(ijk,1,1) 232 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 221 ijk=1,nx*ny*nz1 pmaxt = amax1(pmaxt,abs(p(ijk,1,1))) wmaxt = amax1(wmaxt,abs(w(ijk,1,1))) 221 continue do 2220 ij=1,nx*(ny+1) vmaxt = amax1(vmaxt,abs(v(ij,0,1))) 2220 continue umaxt = maxval(u) if(mod(nit,iavgg) .eq. 0) then write(6,'('' average growth rate = '',e12.5)') * grwrta grwrta = 0. end if grwrte = 0. if(nit .gt. 1) * grwrte = (86400.*2./dtl)*log(pmaxt/pmax) grwrta = grwrta + grwrte/float(iavgg) pmax = pmaxt wmax = wmaxt write(6,'('' step '', * i7,5(1x,e10.3))')nit,time,pmax,wmax,umaxt ,vmaxt ccc if(vmaxt .ge. 1.70) go to 7001 if(wmax .gt. 1000.) then write(6,'('' blowup!!!, step = ''i5)')nit rewind(9) write(9)nx,ny,nz1,time write(9)u,u1,v,v1,w,w1,p,p1,t,t1 stop end if endif DTL=2.*DT NS=2*NS0 c c take another timestep if necessary c if(nit .lt. nstopc) go to 3000 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) 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) 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; write(6,*) 'k=',k; do j=1,ny; do i=0,nx; read(11,*) u(i,j,k); write(6,*) i,j,k, u(i,j,k) !!!! write(6,*) u(i,j,k); read(5,*) time; enddo; enddo; enddo; 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