source: trunk/SpectralModelF90/PROJECTS/trivial/energetics.f @ 4

Last change on this file since 4 was 4, checked in by xlvlod, 17 years ago

import initial des vraies codes.
File size: 13.2 KB
Line 
1      subroutine energetics(u,v,w,pd,uc,vc,wc,wnx,wny,wnz,
2     *                      U0,Lx,Ly,Lz,f,g,rho_0,DGRAD,Re,kappa,nu,
3     *                      nx,ny,nz,num_dims,bc_flag,force_flag,time,
4     *                      comm,myid,numprocs,locnx,locnz,ambient_density)
5
6c     Modified for parallel execution under MPI, KW, 8/24/00
7c     N.B. All processors recieve the global summation results
8c     when done using MPI_ALL_REDUCE.
9c
10c     Routine to compute and output physical space, volume averaged
11c     energetics diagnostics. Input fields are dimensionless, outputs
12c     are in dimensional units. Transforms of u,v and w are stored
13c     in uc, vc and wc respectively.
14c
15      implicit none
16#include "mpif.h"
17      integer comm,myid,numprocs,locnx,locnz,ierr,ktop
18      integer i,j,k,nx,ny,nz,num_dims,nyplanes,nzplanes
19      real Lx,Ly,Lz,U0,f,g,rho_0,DGRAD,Re,x,y,z,density_scale
20      real dx,dy,dz,dA,dV,Vol,time,kappa,nu,fac,N2
21      real rho_bar,ddz,junk
22      character*80 bc_flag
23      character*3 force_flag
24c     Physical space storage locations
25      real u(nx+2,ny+1,locnz+1),v(nx+2,ny+1,locnz+1)
26      real w(nx+2,ny+1,locnz+1),pd(nx+2,ny+1,locnz+1)
27      real ambient_density(nx+2,ny+1,locnz+1)
28c     Wavenumber space storage locations
29      complex uc(locnx,ny+1,nz+1),vc(locnx,ny+1,nz+1),wc(locnx,ny+1,nz+1)
30      real wnx(locnx),wny(ny+1),wnz(nz+1)
31      real*8 ke,pe,bf,rho,kappa2,diss,Ep_surf_adv,Ep_surf_dif,rho_top
32      real*8 rho_bottom,zgrad,phi_i,phi_a,global_val,xmag,xnorm
33      real*8 ke_forced,pe_forced
34      real F1,F2,F3,F4,F5
35
36      real s1_scale
37      s1_scale=1.0
38
39c     THIS NEEDS TO BE UPDATED FROM PD TO T/S LOGIC FOR DENSITY
40
41c
42c     preliminaries:
43c
44
45      dx = Lx/float(nx)   ! [m]
46      dz = Lz/float(nz)   ! [m]
47      ke = 0.d0
48      pe = 0.d0
49      ke_forced = 0.d0
50      pe_forced = 0.d0
51      bf = 0.d0
52      diss = 0.d0
53      rho_top = 0.d0
54      rho_bottom = 0.d0
55      Ep_surf_dif=0.d0
56      Ep_surf_adv=0.d0
57      zgrad=0.d0
58      density_scale = DGRAD*Lz
59
60      if( num_dims .eq. 2 ) then
61       nyplanes = 1
62       dy = 1.0            ! [m] work per unit length in infinite y direction
63       Vol = Lx*1.*Lz      ! [m^3]
64       dV = dx*dy*dz       ! [m^3]
65      elseif( num_dims .eq. 3 ) then
66       nyplanes = ny
67       dy = Ly/float(ny)   ! [m]
68       Vol = Lx*Ly*Lz      ! [m^3]
69       dV = dx*dy*dz       ! [m^3]
70      endif
71     
72
73c   N.B. work with dimensionless quantities ONLY within (ijk-xyz) loop
74c        horizontally periodic bcs imply <w>_xy=0 --> perturbation
75c        density can be used rather than rho_total to compute buoy. flux
76
77      if( bc_flag .eq. 'zperiodic' ) then
78       
79c     volume integrals calculated in physical space,
80c     each processor integrates over its assigned data space
81 
82c     z*pd is not periodic so this form of integration is not quite right for pe
83c     at bottom and at "missing" top", the weighting should be 1/2
84c     OK to have fac=1 at bottom since z=0 anyway, need to add the top bit
85c     note that pd(top)=pd(1) by periodicity, so add this when k=1, {z}=1
86
87       do k=1,locnz
88        z = myid*locnz*dz/Lz + (k-1.)*dz/Lz   ! dz is dimensional
89        if( k .eq. 1 ) then
90         fac=0.5
91        else
92         fac=0.0
93        endif
94        do j=1,nyplanes
95         do i=1,nx
96          rho = pd(i,j,k) 
97          ke = ke + ( u(i,j,k)**2 + v(i,j,k)**2 + w(i,j,k)**2 )
98          pe = pe + rho*z + fac*rho*1.0  ! 1.0 is dless z at top
99          bf = bf + pd(i,j,k)*w(i,j,k)
100         enddo
101        enddo
102       enddo
103
104      if(force_flag .eq. 'yes' ) then
105       do k=1,locnz
106        z = myid*locnz*dz/Lz + (k-1.)*dz/Lz   ! dz is dimensional
107        do j=1,nyplanes
108         do i=1,nx
109          x=(i-1.)*dx
110          y=(j-1.)*dy
111          call user_defined_forcing(x,y,z,time,F1,F2,F3,F4,F5,Lx,Ly,Lz,f,g,rho_0)
112          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
113          pe_forced = pe_forced + ( (z/Lz)*F4/(density_scale/(Lz/U0)) )
114         enddo
115        enddo
116       enddo
117      endif
118
119c     Do global sums using MPI_ALLREDUCE
120      call MPI_ALLREDUCE(ke,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
121      ke=global_val
122      call MPI_ALLREDUCE(pe,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
123      pe=global_val
124      call MPI_ALLREDUCE(bf,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
125      bf=global_val
126
127      if(force_flag .eq. 'yes' ) then
128       call MPI_ALLREDUCE(ke_forced,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
129       ke_forced=global_val
130       call MPI_ALLREDUCE(pe_forced,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
131       pe_forced=global_val
132      endif
133
134c     
135c     surface integrals calculated in physical space
136c
137      ! for pd, z=Lz, by periodicity
138      do j=1,nyplanes
139       do i=1,nx
140        if( myid .eq. 0 ) then
141         rho_top=rho_top + pd(i,j,1)  ! perturbation density at top by periodicity
142         rho_bottom=rho_bottom + ( ambient_density(i,j,1) + pd(i,j,1) )
143         Ep_surf_adv=Ep_surf_adv + pd(i,j,1)*w(i,j,1)   ! z=0 at bottom, --> accumulate ddz(pd) at top only
144         zgrad = zgrad + (pd(i,j,2)-pd(i,j,1))/(dz/Lz)  ! =~ ddz (pd) at top by periodicity
145        endif
146        if( myid .eq. numprocs-1 ) then
147         rho_top=rho_top + ambient_density(i,j,locnz+1)   ! rho_bar at top
148         rho_bottom=rho_bottom + 0.0
149         Ep_surf_adv=Ep_surf_adv + 0.0
150         zgrad = zgrad + (ambient_density(i,j,locnz+1)-ambient_density(i,j,locnz))/(dz/Lz) ! ~ ddz(rho_bar) at top
151        endif
152       enddo
153      enddo
154
155c     Do global sums using MPI_ALLREDUCE
156      call MPI_ALLREDUCE(rho_top,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
157      rho_top=global_val
158      call MPI_ALLREDUCE(rho_bottom,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
159      rho_bottom=global_val
160      call MPI_ALLREDUCE(Ep_surf_adv,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
161      Ep_surf_adv=global_val
162
163      call MPI_ALLREDUCE(zgrad,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
164      zgrad=global_val/(nx*nyplanes)
165
166      Ep_surf_adv=Ep_surf_adv/(nx*nyplanes)
167      Ep_surf_dif=zgrad
168      phi_i=(rho_top-rho_bottom)/(nx*nyplanes)
169
170 
171     
172c     volume integrals calculated in spectral space
173c
174       xnorm = 1./2.   ! normalization for fff transforms
175       do i=1,locnx
176        if( i .eq. 1 .and. myid .eq. 0 ) then
177         fac=xnorm*1.0    ! add zero wavenumber component once
178        else
179         fac=xnorm*2.0    ! all other modes have complex conjugate pairs
180        endif
181        do j=1,nyplanes
182         do k=1,nz
183          kappa2 = wnx(i)**2 + wny(j)**2 + wnz(k)**2
184          xmag = cabs(uc(i,j,k))**2 + cabs(vc(i,j,k))**2 + cabs(wc(i,j,k))**2
185          diss = diss - fac*(1./Re)*kappa2*xmag
186         enddo
187        enddo
188       enddo
189                 
190c     Do global sums using MPI_ALLREDUCE
191      call MPI_ALLREDUCE(diss,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
192      diss=global_val
193
194      elseif( bc_flag .eq. 'zslip' ) then
195c
196c     volume integrals calculated in physical space
197c
198       do k=1,locnz
199        z = myid*locnz*dz/Lz + (k-1.)*dz/Lz   ! dz is dimensional, we need dless z for pe
200        fac=1.0
201        if( k .eq. 1 .and. myid .eq. 0 ) fac=0.5   ! see below
202        do j=1,nyplanes
203         do i=1,nx
204          rho = pd(i,j,k)   
205          ke = ke + fac*( u(i,j,k)**2 + v(i,j,k)**2 + w(i,j,k)**2 )
206          pe = pe + fac*rho*z
207          bf = bf + fac*pd(i,j,k)*w(i,j,k)
208         enddo
209        enddo
210       enddo
211
212      if(force_flag .eq. 'yes' ) then
213       do k=1,locnz
214        z = myid*locnz*dz + (k-1.)*dz   ! we need dimensional position
215        fac=1.0
216        if( k .eq. 1 .and. myid .eq. 0 ) fac=0.5   ! see below
217        do j=1,nyplanes
218         do i=1,nx
219          x=(i-1.)*dx
220          y=(j-1.)*dy
221          call user_defined_forcing(x,y,z,time,F1,F2,F3,F4,F5,Lx,Ly,Lz,f,g,rho_0)
222          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
223          pe_forced = pe_forced + fac*( (z/Lz)*F4/(density_scale/(Lz/U0)) )
224         enddo
225        enddo
226       enddo
227      endif
228
229c      Uppermost processor must add in z=Lz contribution
230c      Give only 1/2 weight to end point of closed interval
231       if (myid .eq. numprocs-1) then
232        fac=0.5
233        k=locnz+1 
234         z = 1.0    ! d'less
235         do j=1,nyplanes
236          do i=1,nx
237           rho = pd(i,j,k)
238           ke = ke + fac*( u(i,j,k)**2 + v(i,j,k)**2 + w(i,j,k)**2 )
239           pe = pe + fac*rho*z
240           bf = bf + fac*pd(i,j,k)*w(i,j,k) 
241          enddo
242         enddo
243       endif
244
245      if(force_flag .eq. 'yes' ) then
246c      Uppermost processor must add in z=Lz contribution
247c      Give only 1/2 weight to end point of closed interval
248       if (myid .eq. numprocs-1) then
249        fac=0.5
250        k=locnz+1 
251         z = Lz   ! we need dimensional position
252         do j=1,nyplanes
253          do i=1,nx
254           x=(i-1.)*dx
255           y=(j-1.)*dy
256           call user_defined_forcing(x,y,z,time,F1,F2,F3,F4,F5,Lx,Ly,Lz,f,g,rho_0)
257           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
258           pe_forced = pe_forced + fac*( (z/Lz)*F4/(density_scale/(Lz/U0)) )
259          enddo
260         enddo
261       endif
262      endif
263
264c     Do global sums using MPI_ALLREDUCE
265      call MPI_ALLREDUCE(ke,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
266      ke=global_val
267      call MPI_ALLREDUCE(pe,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
268      pe=global_val
269      call MPI_ALLREDUCE(bf,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
270      bf=global_val
271
272      if( force_flag .eq. 'yes' ) then
273       call MPI_ALLREDUCE(ke_forced,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
274       ke_forced=global_val
275       call MPI_ALLREDUCE(pe_forced,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
276       pe_forced=global_val
277      endif
278
279c     
280c     surface integrals calculated in physical space
281c
282
283       do j=1,nyplanes
284        do i=1,nx
285         if(myid .eq. numprocs-1 ) then
286          k=locnz+1
287          Ep_surf_adv=Ep_surf_adv + pd(i,j,k)*w(i,j,k) ! w=0 by bcs
288          rho_top=rho_top + ( ambient_density(i,j,k) + pd(i,j,k) )
289          rho_bottom=rho_bottom + 0.
290          zgrad = ( ambient_density(i,j,k)-ambient_density(i,j,k-1) )/(dz/Lz) 
291     *                + ( pd(i,j,k)-pd(i,j,k-1) )/(dz/Lz) 
292          Ep_surf_dif=Ep_surf_dif + zgrad
293         endif
294         if( myid .eq. 0 ) then
295          k=1
296          Ep_surf_adv=Ep_surf_adv + 0.
297          rho_top=rho_top + 0.
298          Ep_surf_dif=Ep_surf_dif + 0.
299          rho_bottom=rho_bottom + ( ambient_density(i,j,k) + pd(i,j,k) )
300         endif
301        enddo
302       enddo
303
304     
305c     Do global sums using MPI_ALLREDUCE
306      call MPI_ALLREDUCE(rho_top,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
307      rho_top=global_val
308      call MPI_ALLREDUCE(rho_bottom,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
309      rho_bottom=global_val
310      call MPI_ALLREDUCE(Ep_surf_adv,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
311      Ep_surf_adv=global_val
312      call MPI_ALLREDUCE(Ep_surf_dif,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
313      Ep_surf_dif=global_val
314
315
316      Ep_surf_adv=Ep_surf_adv/(nx*nyplanes)
317      Ep_surf_dif=Ep_surf_dif/(nx*nyplanes)
318      phi_i=(rho_top-rho_bottom)/(nx*nyplanes)
319c     
320c     volume integrals calculated in spectral space
321c
322       xnorm = (1./nz)**2  ! normalization for ffc/s transforms   
323       do i=1,locnx
324        if(i .eq. 1 .and. myid .eq. 0) then
325         fac=1.0*xnorm    ! add zero wavenumber component once
326        else
327         fac=2.0*xnorm    ! all other modes have complex conjugate pairs
328        endif
329        do j=1,nyplanes
330         do k=1,nz+1
331          kappa2 = wnx(i)**2 + wny(j)**2 + wnz(k)**2
332          xmag = cabs(uc(i,j,k))**2 + cabs(vc(i,j,k))**2 + cabs(wc(i,j,k))**2
333          diss = diss - fac*(1./Re)*kappa2*xmag
334         enddo
335        enddo
336       enddo
337                 
338c     Do global sums using MPI_ALLREDUCE
339      call MPI_ALLREDUCE(diss,global_val,1,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
340      diss=global_val
341
342      endif
343
344c
345c     dimensionalize the results: energy in [J/kg], transfers in [W/kg]
346c
347      ke = (0.5)*(dV/Vol)*(U0**2)*ke                                  ! [m2/s2]=[J/kg]
348      ke_forced = (dV/Vol)*(U0**3/Lz)*ke_forced                       ! [m2/s3]=[W/kg]
349      pe = (dV/Vol)*(g/rho_0)*(density_scale*Lz)*pe                   ! [J/kg]
350      pe_forced = (dV/Vol)*(g/rho_0)*(density_scale*U0)*pe_forced     ! [W/kg]
351      bf = (dV/Vol)*(g/rho_0)*(density_scale*U0)*bf                   ! [m2/s3]=[W/kg]
352      diss = 2.0*(U0**3/Lz)*diss                                      ! [m2/s3]=[W/kg]
353      Ep_surf_adv=-(g/rho_0)*(density_scale*U0)*Ep_surf_adv           ! [m2/s3]=[W/kg]
354      phi_i =-((kappa*g)/(rho_0*Lz))*(density_scale)*phi_i            ! [m2/s3]=[W/kg]
355      Ep_surf_dif=(kappa*g/rho_0)*(density_scale/Lz)*Ep_surf_dif      ! [W/kg]
356      phi_a = bf
357
358c
359c     write the dimensional values to an ascii data file
360c
361      if( myid .eq. 0 ) then
362#ifdef F77
363       open(1,file='output/energetics',status='unknown',access='append')
364#else
365       open(1,file='output/energetics',status='unknown',position='append')
366#endif F77
367        write(1,100) time,ke,pe,phi_a,diss,Ep_surf_adv,Ep_surf_dif,phi_i,ke_forced,pe_forced
368100     format(1x,10e16.8)
369       close(1)
370       WRITE(6,*) time,ke,pe,ke+pe
371      endif
372
373      return
374      end
Note: See TracBrowser for help on using the repository browser.