subroutine energetics(u,v,w,pd,uc,vc,wc,wnx,wny,wnz, * U0,Lx,Ly,Lz,f,g,rho_0,DGRAD,Re,kappa,nu, * nx,ny,nz,num_dims,bc_flag,force_flag,time, * comm,myid,numprocs,locnx,locnz,ambient_density) c Modified for parallel execution under MPI, KW, 8/24/00 c N.B. All processors recieve the global summation results c when done using MPI_ALL_REDUCE. c c Routine to compute and output physical space, volume averaged c energetics diagnostics. Input fields are dimensionless, outputs c are in dimensional units. Transforms of u,v and w are stored c in uc, vc and wc respectively. c implicit none #include "mpif.h" integer comm,myid,numprocs,locnx,locnz,ierr,ktop integer i,j,k,nx,ny,nz,num_dims,nyplanes,nzplanes real Lx,Ly,Lz,U0,f,g,rho_0,DGRAD,Re,x,y,z,density_scale real dx,dy,dz,dA,dV,Vol,time,kappa,nu,fac,N2 real rho_bar,ddz,junk character*80 bc_flag character*3 force_flag c Physical space storage locations real u(nx+2,ny+1,locnz+1),v(nx+2,ny+1,locnz+1) real w(nx+2,ny+1,locnz+1),pd(nx+2,ny+1,locnz+1) real ambient_density(nx+2,ny+1,locnz+1) c Wavenumber space storage locations complex uc(locnx,ny+1,nz+1),vc(locnx,ny+1,nz+1),wc(locnx,ny+1,nz+1) real wnx(locnx),wny(ny+1),wnz(nz+1) real*8 ke,pe,bf,rho,kappa2,diss,Ep_surf_adv,Ep_surf_dif,rho_top real*8 rho_bottom,zgrad,phi_i,phi_a,global_val,xmag,xnorm real*8 ke_forced,pe_forced real F1,F2,F3,F4,F5 real s1_scale s1_scale=1.0 c THIS NEEDS TO BE UPDATED FROM PD TO T/S LOGIC FOR DENSITY c c preliminaries: c dx = Lx/float(nx) ! [m] dz = Lz/float(nz) ! [m] ke = 0.d0 pe = 0.d0 ke_forced = 0.d0 pe_forced = 0.d0 bf = 0.d0 diss = 0.d0 rho_top = 0.d0 rho_bottom = 0.d0 Ep_surf_dif=0.d0 Ep_surf_adv=0.d0 zgrad=0.d0 density_scale = DGRAD*Lz if( num_dims .eq. 2 ) then nyplanes = 1 dy = 1.0 ! [m] work per unit length in infinite y direction Vol = Lx*1.*Lz ! [m^3] dV = dx*dy*dz ! [m^3] elseif( num_dims .eq. 3 ) then nyplanes = ny dy = Ly/float(ny) ! [m] Vol = Lx*Ly*Lz ! [m^3] dV = dx*dy*dz ! [m^3] endif c N.B. work with dimensionless quantities ONLY within (ijk-xyz) loop c horizontally periodic bcs imply _xy=0 --> perturbation c density can be used rather than rho_total to compute buoy. flux if( bc_flag .eq. 'zperiodic' ) then c volume integrals calculated in physical space, c each processor integrates over its assigned data space c z*pd is not periodic so this form of integration is not quite right for pe c at bottom and at "missing" top", the weighting should be 1/2 c OK to have fac=1 at bottom since z=0 anyway, need to add the top bit c note that pd(top)=pd(1) by periodicity, so add this when k=1, {z}=1 do k=1,locnz z = myid*locnz*dz/Lz + (k-1.)*dz/Lz ! dz is dimensional if( k .eq. 1 ) then fac=0.5 else fac=0.0 endif do j=1,nyplanes do i=1,nx rho = pd(i,j,k) ke = ke + ( u(i,j,k)**2 + v(i,j,k)**2 + w(i,j,k)**2 ) pe = pe + rho*z + fac*rho*1.0 ! 1.0 is dless z at top bf = bf + pd(i,j,k)*w(i,j,k) enddo enddo enddo if(force_flag .eq. 'yes' ) then do k=1,locnz z = myid*locnz*dz/Lz + (k-1.)*dz/Lz ! dz is dimensional do j=1,nyplanes do i=1,nx x=(i-1.)*dx y=(j-1.)*dy call user_defined_forcing(x,y,z,time,F1,F2,F3,F4,F5,Lx,Ly,Lz,f,g,rho_0) ke_forced = ke_forced + ( u(i,j,k)*F1 + v(i,j,k)*F2 + w(i,j,k)*F3 )/(U0**2/Lz) ! make F_i dless pe_forced = pe_forced + ( (z/Lz)*F4/(density_scale/(Lz/U0)) ) enddo enddo enddo endif c Do global sums using MPI_ALLREDUCE call MPI_ALLREDUCE(ke,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) ke=global_val call MPI_ALLREDUCE(pe,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) pe=global_val call MPI_ALLREDUCE(bf,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) bf=global_val if(force_flag .eq. 'yes' ) then call MPI_ALLREDUCE(ke_forced,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) ke_forced=global_val call MPI_ALLREDUCE(pe_forced,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) pe_forced=global_val endif c c surface integrals calculated in physical space c ! for pd, z=Lz, by periodicity do j=1,nyplanes do i=1,nx if( myid .eq. 0 ) then rho_top=rho_top + pd(i,j,1) ! perturbation density at top by periodicity rho_bottom=rho_bottom + ( ambient_density(i,j,1) + pd(i,j,1) ) Ep_surf_adv=Ep_surf_adv + pd(i,j,1)*w(i,j,1) ! z=0 at bottom, --> accumulate ddz(pd) at top only zgrad = zgrad + (pd(i,j,2)-pd(i,j,1))/(dz/Lz) ! =~ ddz (pd) at top by periodicity endif if( myid .eq. numprocs-1 ) then rho_top=rho_top + ambient_density(i,j,locnz+1) ! rho_bar at top rho_bottom=rho_bottom + 0.0 Ep_surf_adv=Ep_surf_adv + 0.0 zgrad = zgrad + (ambient_density(i,j,locnz+1)-ambient_density(i,j,locnz))/(dz/Lz) ! ~ ddz(rho_bar) at top endif enddo enddo c Do global sums using MPI_ALLREDUCE call MPI_ALLREDUCE(rho_top,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) rho_top=global_val call MPI_ALLREDUCE(rho_bottom,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) rho_bottom=global_val call MPI_ALLREDUCE(Ep_surf_adv,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) Ep_surf_adv=global_val call MPI_ALLREDUCE(zgrad,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) zgrad=global_val/(nx*nyplanes) Ep_surf_adv=Ep_surf_adv/(nx*nyplanes) Ep_surf_dif=zgrad phi_i=(rho_top-rho_bottom)/(nx*nyplanes) c volume integrals calculated in spectral space c xnorm = 1./2. ! normalization for fff transforms do i=1,locnx if( i .eq. 1 .and. myid .eq. 0 ) then fac=xnorm*1.0 ! add zero wavenumber component once else fac=xnorm*2.0 ! all other modes have complex conjugate pairs endif do j=1,nyplanes do k=1,nz kappa2 = wnx(i)**2 + wny(j)**2 + wnz(k)**2 xmag = cabs(uc(i,j,k))**2 + cabs(vc(i,j,k))**2 + cabs(wc(i,j,k))**2 diss = diss - fac*(1./Re)*kappa2*xmag enddo enddo enddo c Do global sums using MPI_ALLREDUCE call MPI_ALLREDUCE(diss,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) diss=global_val elseif( bc_flag .eq. 'zslip' ) then c c volume integrals calculated in physical space c do k=1,locnz z = myid*locnz*dz/Lz + (k-1.)*dz/Lz ! dz is dimensional, we need dless z for pe fac=1.0 if( k .eq. 1 .and. myid .eq. 0 ) fac=0.5 ! see below do j=1,nyplanes do i=1,nx rho = pd(i,j,k) ke = ke + fac*( u(i,j,k)**2 + v(i,j,k)**2 + w(i,j,k)**2 ) pe = pe + fac*rho*z bf = bf + fac*pd(i,j,k)*w(i,j,k) enddo enddo enddo if(force_flag .eq. 'yes' ) then do k=1,locnz z = myid*locnz*dz + (k-1.)*dz ! we need dimensional position fac=1.0 if( k .eq. 1 .and. myid .eq. 0 ) fac=0.5 ! see below do j=1,nyplanes do i=1,nx x=(i-1.)*dx y=(j-1.)*dy call user_defined_forcing(x,y,z,time,F1,F2,F3,F4,F5,Lx,Ly,Lz,f,g,rho_0) ke_forced = ke_forced + fac*( u(i,j,k)*F1 + v(i,j,k)*F2 + w(i,j,k)*F3 )/(U0**2/Lz) ! make F_i dless pe_forced = pe_forced + fac*( (z/Lz)*F4/(density_scale/(Lz/U0)) ) enddo enddo enddo endif c Uppermost processor must add in z=Lz contribution c Give only 1/2 weight to end point of closed interval if (myid .eq. numprocs-1) then fac=0.5 k=locnz+1 z = 1.0 ! d'less do j=1,nyplanes do i=1,nx rho = pd(i,j,k) ke = ke + fac*( u(i,j,k)**2 + v(i,j,k)**2 + w(i,j,k)**2 ) pe = pe + fac*rho*z bf = bf + fac*pd(i,j,k)*w(i,j,k) enddo enddo endif if(force_flag .eq. 'yes' ) then c Uppermost processor must add in z=Lz contribution c Give only 1/2 weight to end point of closed interval if (myid .eq. numprocs-1) then fac=0.5 k=locnz+1 z = Lz ! we need dimensional position do j=1,nyplanes do i=1,nx x=(i-1.)*dx y=(j-1.)*dy call user_defined_forcing(x,y,z,time,F1,F2,F3,F4,F5,Lx,Ly,Lz,f,g,rho_0) ke_forced = ke_forced + fac*( u(i,j,k)*F1 + v(i,j,k)*F2 + w(i,j,k)*F3 )/(U0**2/Lz) ! make F_i dless pe_forced = pe_forced + fac*( (z/Lz)*F4/(density_scale/(Lz/U0)) ) enddo enddo endif endif c Do global sums using MPI_ALLREDUCE call MPI_ALLREDUCE(ke,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) ke=global_val call MPI_ALLREDUCE(pe,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) pe=global_val call MPI_ALLREDUCE(bf,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) bf=global_val if( force_flag .eq. 'yes' ) then call MPI_ALLREDUCE(ke_forced,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) ke_forced=global_val call MPI_ALLREDUCE(pe_forced,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) pe_forced=global_val endif c c surface integrals calculated in physical space c do j=1,nyplanes do i=1,nx if(myid .eq. numprocs-1 ) then k=locnz+1 Ep_surf_adv=Ep_surf_adv + pd(i,j,k)*w(i,j,k) ! w=0 by bcs rho_top=rho_top + ( ambient_density(i,j,k) + pd(i,j,k) ) rho_bottom=rho_bottom + 0. zgrad = ( ambient_density(i,j,k)-ambient_density(i,j,k-1) )/(dz/Lz) * + ( pd(i,j,k)-pd(i,j,k-1) )/(dz/Lz) Ep_surf_dif=Ep_surf_dif + zgrad endif if( myid .eq. 0 ) then k=1 Ep_surf_adv=Ep_surf_adv + 0. rho_top=rho_top + 0. Ep_surf_dif=Ep_surf_dif + 0. rho_bottom=rho_bottom + ( ambient_density(i,j,k) + pd(i,j,k) ) endif enddo enddo c Do global sums using MPI_ALLREDUCE call MPI_ALLREDUCE(rho_top,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) rho_top=global_val call MPI_ALLREDUCE(rho_bottom,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) rho_bottom=global_val call MPI_ALLREDUCE(Ep_surf_adv,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) Ep_surf_adv=global_val call MPI_ALLREDUCE(Ep_surf_dif,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) Ep_surf_dif=global_val Ep_surf_adv=Ep_surf_adv/(nx*nyplanes) Ep_surf_dif=Ep_surf_dif/(nx*nyplanes) phi_i=(rho_top-rho_bottom)/(nx*nyplanes) c c volume integrals calculated in spectral space c xnorm = (1./nz)**2 ! normalization for ffc/s transforms do i=1,locnx if(i .eq. 1 .and. myid .eq. 0) then fac=1.0*xnorm ! add zero wavenumber component once else fac=2.0*xnorm ! all other modes have complex conjugate pairs endif do j=1,nyplanes do k=1,nz+1 kappa2 = wnx(i)**2 + wny(j)**2 + wnz(k)**2 xmag = cabs(uc(i,j,k))**2 + cabs(vc(i,j,k))**2 + cabs(wc(i,j,k))**2 diss = diss - fac*(1./Re)*kappa2*xmag enddo enddo enddo c Do global sums using MPI_ALLREDUCE call MPI_ALLREDUCE(diss,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) diss=global_val endif c c dimensionalize the results: energy in [J/kg], transfers in [W/kg] c ke = (0.5)*(dV/Vol)*(U0**2)*ke ! [m2/s2]=[J/kg] ke_forced = (dV/Vol)*(U0**3/Lz)*ke_forced ! [m2/s3]=[W/kg] pe = (dV/Vol)*(g/rho_0)*(density_scale*Lz)*pe ! [J/kg] pe_forced = (dV/Vol)*(g/rho_0)*(density_scale*U0)*pe_forced ! [W/kg] bf = (dV/Vol)*(g/rho_0)*(density_scale*U0)*bf ! [m2/s3]=[W/kg] diss = 2.0*(U0**3/Lz)*diss ! [m2/s3]=[W/kg] Ep_surf_adv=-(g/rho_0)*(density_scale*U0)*Ep_surf_adv ! [m2/s3]=[W/kg] phi_i =-((kappa*g)/(rho_0*Lz))*(density_scale)*phi_i ! [m2/s3]=[W/kg] Ep_surf_dif=(kappa*g/rho_0)*(density_scale/Lz)*Ep_surf_dif ! [W/kg] phi_a = bf c c write the dimensional values to an ascii data file c if( myid .eq. 0 ) then #ifdef F77 open(1,file='output/energetics',status='unknown',access='append') #else open(1,file='output/energetics',status='unknown',position='append') #endif F77 write(1,100) time,ke,pe,phi_a,diss,Ep_surf_adv,Ep_surf_dif,phi_i,ke_forced,pe_forced 100 format(1x,10e16.8) close(1) WRITE(6,*) time,ke,pe,ke+pe endif return end